From efdf0f71ed0a6a406d296a2c1f166e8a8476b256 Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 11 May 2018 12:18:30 +0100 Subject: [PATCH 01/69] First attempt at grey mars. Need to add timings before it can really be tested. --- exp/grey_mars/grey_mars.py | 221 +++++++++++++++++++++++++++++++++++++ 1 file changed, 221 insertions(+) create mode 100644 exp/grey_mars/grey_mars.py diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py new file mode 100644 index 000000000..26d612f23 --- /dev/null +++ b/exp/grey_mars/grey_mars.py @@ -0,0 +1,221 @@ +import numpy as np + +from isca import IscaCodeBase, GreyCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE +from isca.util import exp_progress +from ntfy import notify + +NCORES = 16 + +# a CodeBase can be a directory on the computer, +# useful for iterative development +cb = GreyCodeBase.from_directory(GFDL_BASE) +# cb = GreyCodeBase.from_repo(repo='git@github.com:sit23/Isca.git', commit='7145be9') + +# or it can point to a specific git repo and commit id. +# This method should ensure future, independent, reproducibility of results. +# cb = DryCodeBase.from_repo(repo='https://github.com/isca/isca', commit='isca1.1') + +# compilation depends on computer specific settings. The $GFDL_ENV +# environment variable is used to determine which `$GFDL_BASE/src/extra/env` file +# is used to load the correct compilers. The env file is always loaded from +# $GFDL_BASE and not the checked out git repo. + +cb.compile() # compile the source code to working directory $GFDL_WORK/codebase + +# create a diagnostics output file for daily snapshots +diag = DiagTable() +diag.add_file('atmos_monthly', 30, 'days', time_units='days') +diag.add_file('atmos_daily', 1, 'days', time_units='days') + +# Define diag table entries +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk', time_avg=True) +diag.add_field('dynamics', 'pk', time_avg=True) +diag.add_field('dynamics', 'vor', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'div', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'vcomp', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'ucomp_vcomp', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'temp', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'slp', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'omega', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'height', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'sphum', time_avg=True, files=['atmos_monthly']) +diag.add_field('atmosphere', 'precipitation', time_avg=True, files=['atmos_monthly']) +diag.add_field('atmosphere', 'rh', time_avg=True, files=['atmos_monthly']) + +diag.add_field('mixed_layer', 't_surf', time_avg=True, files=['atmos_monthly']) +diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) +diag.add_field('mixed_layer', 'flux_t', time_avg=True) + + +# define namelist values as python dictionary +namelist = Namelist({ + 'main_nml': { + 'dt_atmos': 150, + 'days': 30., + 'seconds': 0, + 'calendar': 'no_calendar' + }, + + 'idealized_moist_phys_nml': { + 'do_damping': True, + 'turb':True, + 'mixed_layer_bc':True, + 'do_virtual' :False, + 'do_simple': True, + 'roughness_mom':3.21e-05, + 'roughness_heat':3.21e-05, + 'roughness_moist':0., + 'two_stream_gray': True, #Use grey radiation + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, # default: True + 'do_diffusivity': True, # default: False + 'do_simple': True, # default: False + 'constant_gust': 0.0, # default: 1.0 + 'use_tau': False + }, + + 'diffusivity_nml': { + 'do_entrain':False, + 'do_simple': True, + }, + + 'surface_flux_nml': { + 'use_virtual_temp': False, + 'do_simple': True, + 'old_dtaudv': True, + 'rh_target': 50., + 'delta_t_relax': 7200., + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + 'mixed_layer_nml': { + 'tconst' : 285., + 'prescribe_initial_dist':True, + 'evaporation':False, + 'albedo': 0.25, + }, + + 'qe_moist_convection_nml': { + 'rhbm':0.0, + 'tau_bm':3600., + }, + + 'dry_convection_nml': { + 'tau':7200, + }, + + 'betts_miller_nml': { + 'rhbm': .7 , + 'do_simp': False, + 'do_shallower': True, + }, + + 'lscale_cond_nml': { + 'do_simple':True, + 'do_evap':True + }, + + 'sat_vapor_pres_nml': { + 'do_simple':True + }, + + 'damping_driver_nml': { + 'do_rayleigh': True, + 'trayfric': -0.25, # neg. value: time in *days* + 'sponge_pbottom': 0.5, #Bottom of the model's sponge down to 50hPa (units are Pa) + 'do_conserve_energy': True, + }, + + 'spectral_dynamics_nml': { + 'num_levels': 30, + 'exponent': 2.5, + 'scale_heights': 4, + 'surf_res': 0.1, + 'robert_coeff': 4e-2, + 'do_water_correction': False, + 'vert_coord_option': 'even_sigma', + 'initial_sphum': 0., + 'valid_range_T': [0, 700] + }, + + 'two_stream_gray_rad_nml': { + 'rad_scheme': 'frierson', #Select radiation scheme to use, which in this case is Frierson + 'do_seasonal': True, + 'atm_abs': 0.2, + 'sw_diff':0.0, + 'ir_tau_eq':0.2, + 'ir_tau_pole':0.2, + 'linear_tau': 1.0, + }, + + +# configure the relaxation profile +# 'hs_forcing_nml': { +# 'equilibrium_t_option' : 'top_down', +# 'ml_depth': 10., +# 'spinup_time': 108000, +# 'ka': -2., +# 'ks': -2., +# 'tau_s': 0.2, +# 'calculate_insolation_from_orbit' : True, +# 'do_rayleigh_damping':False, +# 'albedo':0.25, +# 'pure_rad_equil':True, +# 'stratosphere_t_option':'pure_rad_equil', +# 'peri_time': 0., +# 'h_a': 10.8, +# 'use_olr_from_t_surf':True, +# 'frac_of_year_ae':0.3032894472101727, +# }, + + 'astronomy_nml': { + 'ecc':0.0935, + 'obliq':25.19, + 'per':336.0, + }, + + 'constants_nml': { + 'orbital_period': 686.980, + 'solar_const':589.0, + 'radius':3396.0e3, + 'omega':7.088e-5, + 'rdgas':192.0, + 'kappa':0.22727, + }, + +}) + +if __name__=="__main__": + + conv_schemes = ['dry'] + + scales = [ 1.] + + for conv in conv_schemes: + for scale in scales: + exp = Experiment('grey_mars_mk1', codebase=cb) + exp.clear_rundir() + + exp.diag_table = diag + exp.namelist = namelist.copy() + exp.namelist['constants_nml']['grav'] = scale * 3.71 + exp.namelist['constants_nml']['pstd'] = scale * 6100.0 + exp.namelist['constants_nml']['pstd_mks'] = scale * 610.0 + exp.namelist['spectral_dynamics_nml']['reference_sea_level_press'] = scale * 610.0 + exp.namelist['idealized_moist_phys_nml']['convection_scheme'] = conv + + notify('top down with conv scheme = '+conv+' has started', 'isca') + + with exp_progress(exp, description='o%.0f d{day}' % scale): + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2, 241): + with exp_progress(exp, description='o%.0f d{day}' % scale): + exp.run(i, num_cores=NCORES) + notify('top down with conv scheme = '+conv+' has completed', 'isca') From a58dafda66c59b3eda3976ac8ee3652b9f3d055e Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 11 May 2018 13:01:16 +0100 Subject: [PATCH 02/69] Added Alex Paterson's calculation of eccentric anomaly etc to astronomy module in order that we can correctly calculate rrsun using the true anomaly rather than the mean anomaly. Have also added multiplication by rrsun to two-stream-grey-rad, which makes no difference to results in initial tests when ecc = 0.0 --- exp/grey_mars/grey_mars.py | 20 +++---- .../two_stream_gray_rad.F90 | 2 +- src/shared/astronomy/astronomy.f90 | 57 +++++++++++++++++-- 3 files changed, 64 insertions(+), 15 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index 26d612f23..8fe9d0549 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -48,12 +48,14 @@ diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) diag.add_field('mixed_layer', 'flux_t', time_avg=True) +# diag.add_field('two_stream', 'coszen', time_avg=True) + # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { 'dt_atmos': 150, - 'days': 30., + 'days': 1., 'seconds': 0, 'calendar': 'no_calendar' }, @@ -99,7 +101,7 @@ 'tconst' : 285., 'prescribe_initial_dist':True, 'evaporation':False, - 'albedo': 0.25, + 'albedo_value': 0.25, }, 'qe_moist_convection_nml': { @@ -176,7 +178,7 @@ # }, 'astronomy_nml': { - 'ecc':0.0935, +# 'ecc':0.0935, 'obliq':25.19, 'per':336.0, }, @@ -200,7 +202,7 @@ for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk1', codebase=cb) + exp = Experiment('grey_mars_mk4_without_rrsun_multiplying', codebase=cb) exp.clear_rundir() exp.diag_table = diag @@ -211,11 +213,9 @@ exp.namelist['spectral_dynamics_nml']['reference_sea_level_press'] = scale * 610.0 exp.namelist['idealized_moist_phys_nml']['convection_scheme'] = conv - notify('top down with conv scheme = '+conv+' has started', 'isca') - with exp_progress(exp, description='o%.0f d{day}' % scale): exp.run(1, use_restart=False, num_cores=NCORES) - for i in range(2, 241): - with exp_progress(exp, description='o%.0f d{day}' % scale): - exp.run(i, num_cores=NCORES) - notify('top down with conv scheme = '+conv+' has completed', 'isca') +# for i in range(2, 241): +# with exp_progress(exp, description='o%.0f d{day}' % scale): +# exp.run(i, num_cores=NCORES) +# notify('top down with conv scheme = '+conv+' has completed', 'isca') diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 84215b6e2..3fbe16690 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -444,7 +444,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun) end if - insolation = solar_constant * coszen + insolation = solar_constant * coszen * rrsun else if (sw_scheme==B_SCHNEIDER_LIU) then insolation = (solar_constant/pi)*cos(lat) diff --git a/src/shared/astronomy/astronomy.f90 b/src/shared/astronomy/astronomy.f90 index b5181b56d..a2dc2adc2 100644 --- a/src/shared/astronomy/astronomy.f90 +++ b/src/shared/astronomy/astronomy.f90 @@ -160,11 +160,12 @@ module astronomy_mod integer :: num_angles = 3600 ! number of intervals into which the year ! is divided to compute orbital positions +logical :: use_mean_anom_in_rrsun_calc = .TRUE. !It appears the standard astronomy module uses the mean anomaly for calculating the orbital distances on eccentric orbits, rather than the true anomaly. Keeping this behaviour default for legacy, but .FALSE. seems more correct. namelist /astronomy_nml/ ecc, obliq, per, period, & year_ae, month_ae, day_ae, & hour_ae, minute_ae, second_ae, & - num_angles + num_angles, use_mean_anom_in_rrsun_calc !-------------------------------------------------------------------- !------ public data ---------- @@ -3229,7 +3230,7 @@ function r_inv_squared (ang) !--------------------------------------------------------------------- ! local variables - real :: r_inv_squared, r, rad_per + real :: r_inv_squared, r, rad_per, mean_anomaly, true_anomaly, ecc_anomaly !--------------------------------------------------------------------- ! local variables: @@ -3249,13 +3250,61 @@ function r_inv_squared (ang) ! its square (r_inv_squared) to the calling routine. !-------------------------------------------------------------------- rad_per = per*deg_to_rad - r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per)) - r_inv_squared = r**(-2) + + if (ecc.eq.0.) then + r = 1. + r_inv_squared = 1. + else + + if (use_mean_anom_in_rrsun_calc) then + r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per)) + else + mean_anomaly = ang - rad_per + call calc_ecc_anomaly(mean_anomaly, ecc, ecc_anomaly) + true_anomaly = 2*atan(((1 + ecc)/(1 - ecc))**0.5 * tan(ecc_anomaly/2)) + r = (1. - ecc**2)/(1. + ecc*cos(true_anomaly)) + endif + + r_inv_squared = r**(-2) + + endif + end function r_inv_squared +!####################################################################### +!Written by Alex Paterson in hs_forcing.f90 and copied to astronomy by SIT 11th May '18 + + subroutine calc_ecc_anomaly( mean_anomaly, ecc, ecc_anomaly) + + real, intent(in) :: mean_anomaly, ecc + real, intent(out) :: ecc_anomaly + real :: dE, d + integer, parameter :: maxiter = 30 + real, parameter :: tol = 1.d-10 + integer :: k + + ecc_anomaly = mean_anomaly + d = ecc_anomaly - ecc*sin(ecc_anomaly) - mean_anomaly + do k=1,maxiter + dE = d/(1 - ecc*cos(ecc_anomaly)) + ecc_anomaly = ecc_anomaly - dE + d = ecc_anomaly - ecc*sin(ecc_anomaly) - mean_anomaly + if (abs(d) < tol) then + exit + endif + enddo + + if (k > maxiter) then + if (abs(d) > tol) then + print *, '*** Warning: eccentric anomaly has not converged' + endif + endif + + end subroutine calc_ecc_anomaly +!################################################################### !#################################################################### From 07bb7be3c7d489cc80a5a5a33408adb4582d992e Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 11 May 2018 14:49:18 +0100 Subject: [PATCH 03/69] Made alternative r_inv_squared subroutine so that we can pass back true anomaly, and therefore be able to do mars time-telling. --- exp/grey_mars/grey_mars.py | 14 ++- .../two_stream_gray_rad.F90 | 19 ++-- src/shared/astronomy/astronomy.f90 | 98 +++++++++++++++++-- 3 files changed, 111 insertions(+), 20 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index 8fe9d0549..daacb73e7 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -48,14 +48,16 @@ diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) diag.add_field('mixed_layer', 'flux_t', time_avg=True) -# diag.add_field('two_stream', 'coszen', time_avg=True) +# diag.add_field('two_stream', 'mars_solar_long', time_avg=True) +diag.add_field('two_stream', 'coszen', time_avg=True) +diag.add_field('two_stream', 'rrsun', time_avg=True) # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { 'dt_atmos': 150, - 'days': 1., + 'days': 2., 'seconds': 0, 'calendar': 'no_calendar' }, @@ -178,9 +180,11 @@ # }, 'astronomy_nml': { -# 'ecc':0.0935, + 'ecc':0.0935, 'obliq':25.19, - 'per':336.0, + 'per':250.815799, #to be equivalent to peri_time = 0. in hs_forcing version of mars + 'use_mean_anom_in_rrsun_calc':False, + 'use_old_r_inv_squared':False }, 'constants_nml': { @@ -202,7 +206,7 @@ for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk4_without_rrsun_multiplying', codebase=cb) + exp = Experiment('grey_mars_mk13', codebase=cb) exp.clear_rundir() exp.diag_table = diag diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 3fbe16690..60cfc6fee 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -154,7 +154,7 @@ module two_stream_gray_rad_mod integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, & - id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2 + id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_mars_solar_long, id_rrsun character(len=10), parameter :: mod_name = 'two_stream' @@ -377,6 +377,12 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb register_diag_field ( mod_name, 'lw_dtrans', axes(1:3), Time, & 'LW transmission (non window)', & 'none', missing_value=missing_value ) + + id_mars_solar_long = register_diag_field ( mod_name, 'mars_solar_long', & + Time, 'Martian solar longitude', 'deg') + + id_rrsun = register_diag_field ( mod_name, 'rrsun', & + Time, 'inverse planet sun distance', 'none') return end subroutine two_stream_gray_rad_init @@ -397,7 +403,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, integer :: i, j, k, n, dyofyr integer :: seconds, year_in_s, days -real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians +real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians, true_anomaly logical :: used @@ -405,10 +411,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, n = size(t,3) - - -! albedo(:,:) = albedo_value !s albedo now set in mixed_layer_init. - ! ================================================================================= ! SHORTWAVE RADIATION @@ -446,6 +448,11 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, insolation = solar_constant * coszen * rrsun + if (id_mars_solar_long > 0) used = send_data ( id_mars_solar_long, modulo((180./pi)*(true_anomaly-1.905637),360.), Time_diag) + + if (id_rrsun > 0) used = send_data ( id_rrsun, rrsun, Time_diag) + + else if (sw_scheme==B_SCHNEIDER_LIU) then insolation = (solar_constant/pi)*cos(lat) else diff --git a/src/shared/astronomy/astronomy.f90 b/src/shared/astronomy/astronomy.f90 index a2dc2adc2..1982d571d 100644 --- a/src/shared/astronomy/astronomy.f90 +++ b/src/shared/astronomy/astronomy.f90 @@ -161,11 +161,13 @@ module astronomy_mod ! is divided to compute orbital positions logical :: use_mean_anom_in_rrsun_calc = .TRUE. !It appears the standard astronomy module uses the mean anomaly for calculating the orbital distances on eccentric orbits, rather than the true anomaly. Keeping this behaviour default for legacy, but .FALSE. seems more correct. +logical :: use_old_r_inv_squared namelist /astronomy_nml/ ecc, obliq, per, period, & year_ae, month_ae, day_ae, & hour_ae, minute_ae, second_ae, & - num_angles, use_mean_anom_in_rrsun_calc + num_angles, use_mean_anom_in_rrsun_calc, & + use_old_r_inv_squared !-------------------------------------------------------------------- !------ public data ---------- @@ -1122,7 +1124,7 @@ end subroutine get_ref_date_of_ae ! subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, dt, allow_negative_cosz, & - half_day_out) + half_day_out, true_anom) !--------------------------------------------------------------------- ! diurnal_solar_2d returns 2d fields of cosine of zenith angle, @@ -1138,7 +1140,7 @@ subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & real, intent(in), optional :: dt logical, intent(in), optional :: allow_negative_cosz real, dimension(:,:), intent(out), optional :: half_day_out - +real, intent(out), optional :: true_anom !--------------------------------------------------------------------- ! local variables @@ -1187,7 +1189,16 @@ subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & !--------------------------------------------------------------------- ang = angle(time_since_ae) dec = declination(ang) - rrsun = r_inv_squared(ang) + if (use_old_r_inv_squared) then + rrsun = r_inv_squared(ang) + else + if (present(true_anom)) then + call r_inv_squared_alt(ang,rrsun, true_anom) + else + call r_inv_squared_alt(ang,rrsun) + endif + + endif !--------------------------------------------------------------------- ! define terms needed in the cosine zenith angle equation. @@ -2080,7 +2091,7 @@ subroutine daily_mean_solar_2d (lat, time_since_ae, cosz, h_out, rr_out) ! local variables real, dimension(size(lat,1),size(lat,2)) :: h - real :: ang, dec, rr + real :: ang, dec, rr, ta !-------------------------------------------------------------------- ! be sure the time in the annual cycle is legitimate. @@ -2097,7 +2108,11 @@ subroutine daily_mean_solar_2d (lat, time_since_ae, cosz, h_out, rr_out) ang = angle (time_since_ae) dec = declination(ang) h = half_day (lat, dec) - rr = r_inv_squared (ang) + if (use_old_r_inv_squared) then + rr = r_inv_squared(ang) + else + call r_inv_squared_alt(ang,rr, ta) + endif !--------------------------------------------------------------------- ! where the entire day is dark, define cosz to be zero. otherwise @@ -3230,7 +3245,7 @@ function r_inv_squared (ang) !--------------------------------------------------------------------- ! local variables - real :: r_inv_squared, r, rad_per, mean_anomaly, true_anomaly, ecc_anomaly + real :: r_inv_squared, r, rad_per !--------------------------------------------------------------------- ! local variables: @@ -3255,6 +3270,67 @@ function r_inv_squared (ang) r = 1. r_inv_squared = 1. else + r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per)) + r_inv_squared = r**(-2) + endif + + + +end function r_inv_squared + +subroutine r_inv_squared_alt (ang, r_inv_squared_out, true_anomaly_out) + +!-------------------------------------------------------------------- +! r_inv_squared returns the inverse of the square of the earth-sun +! distance relative to the mean distance at angle ang in the earth's +! orbit. +!-------------------------------------------------------------------- + +!-------------------------------------------------------------------- +real, intent(in) :: ang +real, intent(out) :: r_inv_squared_out +real, intent(out), optional:: true_anomaly_out +!-------------------------------------------------------------------- + +!--------------------------------------------------------------------- +! +! intent(in) variables: +! +! ang angular position of earth in its orbit, relative to a +! value of 0.0 at the NH autumnal equinox, value between +! 0.0 and 2 * pi +! [ radians ] +! +!--------------------------------------------------------------------- + +!--------------------------------------------------------------------- +! local variables + + real :: r, rad_per, mean_anomaly, ecc_anomaly, true_anomaly + +!--------------------------------------------------------------------- +! local variables: +! +! r_inv_squared the inverse of the square of the earth-sun +! distance relative to the mean distance +! [ dimensionless ] +! r earth-sun distance relative to mean distance +! [ dimensionless ] +! rad_per angular position of perihelion +! [ radians ] +! +!-------------------------------------------------------------------- + +!-------------------------------------------------------------------- +! define the earth-sun distance (r) and then return the inverse of +! its square (r_inv_squared) to the calling routine. +!-------------------------------------------------------------------- + rad_per = per*deg_to_rad + + if (ecc.eq.0.) then + r = 1. + r_inv_squared_out = 1. + else if (use_mean_anom_in_rrsun_calc) then r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per)) @@ -3265,13 +3341,17 @@ function r_inv_squared (ang) r = (1. - ecc**2)/(1. + ecc*cos(true_anomaly)) endif - r_inv_squared = r**(-2) + r_inv_squared_out = r**(-2) endif + if (present(true_anomaly_out)) then + true_anomaly_out = true_anomaly + + endif +end subroutine -end function r_inv_squared !####################################################################### !Written by Alex Paterson in hs_forcing.f90 and copied to astronomy by SIT 11th May '18 From 2b166e43834b8df0f21df798f4af9c97e56e1e2f Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 11 May 2018 15:02:39 +0100 Subject: [PATCH 04/69] Small changes to make mars-solar-long calc more consistent across different settings. --- exp/grey_mars/grey_mars.py | 8 +++++--- .../two_stream_gray_rad/two_stream_gray_rad.F90 | 4 ++-- src/shared/astronomy/astronomy.f90 | 1 + 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index daacb73e7..ff15228ca 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -48,16 +48,17 @@ diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) diag.add_field('mixed_layer', 'flux_t', time_avg=True) -# diag.add_field('two_stream', 'mars_solar_long', time_avg=True) +diag.add_field('two_stream', 'mars_solar_long', time_avg=True) diag.add_field('two_stream', 'coszen', time_avg=True) diag.add_field('two_stream', 'rrsun', time_avg=True) +diag.add_field('two_stream', 'swdn_toa', time_avg=True) # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { 'dt_atmos': 150, - 'days': 2., + 'days': 30., 'seconds': 0, 'calendar': 'no_calendar' }, @@ -104,6 +105,7 @@ 'prescribe_initial_dist':True, 'evaporation':False, 'albedo_value': 0.25, + 'depth':10., }, 'qe_moist_convection_nml': { @@ -206,7 +208,7 @@ for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk13', codebase=cb) + exp = Experiment('grey_mars_mk14_first_long_run', codebase=cb) exp.clear_rundir() exp.diag_table = diag diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 60cfc6fee..406e29db1 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -441,9 +441,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, r_dt_rad_avg=real(dt_rad_avg) dt_rad_radians = (r_dt_rad_avg/day_in_s)*2.0*pi - call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, dt_rad_radians) + call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, dt_rad_radians, true_anom=true_anomaly) else - call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun) + call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, true_anom=true_anomaly) end if insolation = solar_constant * coszen * rrsun diff --git a/src/shared/astronomy/astronomy.f90 b/src/shared/astronomy/astronomy.f90 index 1982d571d..417c92c16 100644 --- a/src/shared/astronomy/astronomy.f90 +++ b/src/shared/astronomy/astronomy.f90 @@ -3334,6 +3334,7 @@ subroutine r_inv_squared_alt (ang, r_inv_squared_out, true_anomaly_out) if (use_mean_anom_in_rrsun_calc) then r = (1. - ecc**2)/(1. + ecc*cos(ang - rad_per)) + true_anomaly = ang-rad_per !This obviously isn't true, but is the approximation made by the old astronomy code. else mean_anomaly = ang - rad_per call calc_ecc_anomaly(mean_anomaly, ecc, ecc_anomaly) From 387addcc435548cdc2ecd4d8f21566cc6b753e83 Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 11 May 2018 15:13:54 +0100 Subject: [PATCH 05/69] A few extra changes before full long run. --- exp/grey_mars/grey_mars.py | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index ff15228ca..6dc49928c 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -31,20 +31,16 @@ diag.add_field('dynamics', 'ps', time_avg=True) diag.add_field('dynamics', 'bk', time_avg=True) diag.add_field('dynamics', 'pk', time_avg=True) -diag.add_field('dynamics', 'vor', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'div', time_avg=True, files=['atmos_monthly']) diag.add_field('dynamics', 'ucomp', time_avg=True) -diag.add_field('dynamics', 'vcomp', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'ucomp_vcomp', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'temp', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'slp', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'omega', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'height', time_avg=True, files=['atmos_monthly']) -diag.add_field('dynamics', 'sphum', time_avg=True, files=['atmos_monthly']) -diag.add_field('atmosphere', 'precipitation', time_avg=True, files=['atmos_monthly']) -diag.add_field('atmosphere', 'rh', time_avg=True, files=['atmos_monthly']) - -diag.add_field('mixed_layer', 't_surf', time_avg=True, files=['atmos_monthly']) +diag.add_field('dynamics', 'vcomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) + +diag.add_field('dynamics', 'omega', time_avg=True) +diag.add_field('dynamics', 'height', time_avg=True) +diag.add_field('dynamics', 'sphum', time_avg=True) +diag.add_field('atmosphere', 'precipitation', time_avg=True) + +diag.add_field('mixed_layer', 't_surf', time_avg=True) diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) diag.add_field('mixed_layer', 'flux_t', time_avg=True) @@ -221,7 +217,7 @@ with exp_progress(exp, description='o%.0f d{day}' % scale): exp.run(1, use_restart=False, num_cores=NCORES) -# for i in range(2, 241): -# with exp_progress(exp, description='o%.0f d{day}' % scale): -# exp.run(i, num_cores=NCORES) -# notify('top down with conv scheme = '+conv+' has completed', 'isca') + for i in range(2, 241): + with exp_progress(exp, description='o%.0f d{day}' % scale): + exp.run(i, num_cores=NCORES) + notify('top down with conv scheme = '+conv+' has completed', 'isca') From 6708bfc18e47323976ed39d0f81b7a77af4ef494 Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 11 May 2018 16:43:39 +0100 Subject: [PATCH 06/69] Added extra diagnostics to help with debugging. Problem is that model expects orbital_period in seconds, but I had supplied it in days. This is because Alex's code expects it in days. Will fix Alex code to make it consistent. --- exp/grey_mars/grey_mars.py | 7 +++++-- .../two_stream_gray_rad.F90 | 20 ++++++++++++++++--- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index 6dc49928c..d636af2d0 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -48,6 +48,8 @@ diag.add_field('two_stream', 'coszen', time_avg=True) diag.add_field('two_stream', 'rrsun', time_avg=True) diag.add_field('two_stream', 'swdn_toa', time_avg=True) +diag.add_field('two_stream', 'time_since_ae', time_avg=True) +diag.add_field('two_stream', 'true_anomaly', time_avg=True) # define namelist values as python dictionary @@ -155,6 +157,7 @@ 'ir_tau_eq':0.2, 'ir_tau_pole':0.2, 'linear_tau': 1.0, + 'equinox_day':0.3032894472101727, }, @@ -186,7 +189,7 @@ }, 'constants_nml': { - 'orbital_period': 686.980, + 'orbital_period': 686.980*86400., 'solar_const':589.0, 'radius':3396.0e3, 'omega':7.088e-5, @@ -204,7 +207,7 @@ for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk14_first_long_run', codebase=cb) + exp = Experiment('grey_mars_mk15_fixed_year_length', codebase=cb) exp.clear_rundir() exp.diag_table = diag diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 406e29db1..d2221fc81 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -154,7 +154,8 @@ module two_stream_gray_rad_mod integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, & - id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_mars_solar_long, id_rrsun + id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_mars_solar_long, id_rrsun, id_true_anom, & + id_time_since_ae character(len=10), parameter :: mod_name = 'two_stream' @@ -380,9 +381,16 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb id_mars_solar_long = register_diag_field ( mod_name, 'mars_solar_long', & Time, 'Martian solar longitude', 'deg') + + id_true_anom = register_diag_field ( mod_name, 'true_anomaly', & + Time, 'True anomaly', 'deg') id_rrsun = register_diag_field ( mod_name, 'rrsun', & - Time, 'inverse planet sun distance', 'none') + Time, 'inverse planet sun distance', 'none') + + id_time_since_ae = register_diag_field ( mod_name, 'time_since_ae', & + Time, 'time since ae', 'none') + return end subroutine two_stream_gray_rad_init @@ -433,7 +441,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, endif gmt = abs(mod(frac_of_day, 1.0)) * 2.0 * pi - + time_since_ae = modulo(frac_of_year-equinox_day, 1.0) * 2.0 * pi if(use_time_average_coszen) then @@ -450,6 +458,12 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, if (id_mars_solar_long > 0) used = send_data ( id_mars_solar_long, modulo((180./pi)*(true_anomaly-1.905637),360.), Time_diag) + + if (id_time_since_ae > 0) used = send_data ( id_time_since_ae, time_since_ae, Time_diag) + + + if (id_true_anom > 0) used = send_data ( id_true_anom, true_anomaly, Time_diag) + if (id_rrsun > 0) used = send_data ( id_rrsun, rrsun, Time_diag) From 21a7e88b7e6880a537b1e592c44973b18fbe6f01 Mon Sep 17 00:00:00 2001 From: sit23 Date: Tue, 15 May 2018 17:00:24 +0100 Subject: [PATCH 07/69] Script for calculating the equilibrium surface temperature for a range of surface optical depths and albedos, assuming transparancy in the visible. Useful for working out rough values for optical depths and albedos for a simple mars model. --- .../python/scripts/grey_rad_parm_check.py | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 src/extra/python/scripts/grey_rad_parm_check.py diff --git a/src/extra/python/scripts/grey_rad_parm_check.py b/src/extra/python/scripts/grey_rad_parm_check.py new file mode 100644 index 000000000..ee12d247a --- /dev/null +++ b/src/extra/python/scripts/grey_rad_parm_check.py @@ -0,0 +1,40 @@ +import numpy as np +import matplotlib.pyplot as plt +import pdb + +sigma=5.67e-8 +solar_constant = 589.0 + + +tau_s_range = np.arange(0.01, 2.0, 0.01) +albedo_range = np.arange(0.1, 0.8, 0.01) + +n_tau_s = np.shape(tau_s_range)[0] +n_albedo = np.shape(albedo_range)[0] + +tau_s_arr = np.zeros((n_tau_s, n_albedo)) +albedo_arr = np.zeros((n_tau_s, n_albedo)) + +for al in range(n_albedo): + for ts in range(n_tau_s): + tau_s_arr[ts, al] = tau_s_range[ts] + albedo_arr[ts,al] = albedo_range[al] + +def calc_t_surf(albedo_in, tau_s_in): + olr = (1.-albedo_in)* solar_constant + + t_surf = (olr*(tau_s_in+1.)/(2.*sigma))**0.25 + + return t_surf + +t_surf_arr = calc_t_surf(albedo_arr, tau_s_arr) + +cs = plt.contourf(albedo_range, tau_s_range, t_surf_arr) +cn = plt.contour(albedo_range, tau_s_range, t_surf_arr,levels=[210., 260.], colors=['k']) +#plt.clabel(cn, inline=1, fontsize=10.) +fig = plt.gcf() + +cb = fig.colorbar(cs) + +t_surf_value = calc_t_surf(0.25, 0.2) + From 9df8af6e1078c11cafeb0ffd055c1c6e30d26274 Mon Sep 17 00:00:00 2001 From: sit23 Date: Tue, 15 May 2018 17:00:41 +0100 Subject: [PATCH 08/69] Test mars script using lower optical depths. --- exp/grey_mars/grey_mars.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index d636af2d0..50db3b827 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -24,7 +24,7 @@ # create a diagnostics output file for daily snapshots diag = DiagTable() -diag.add_file('atmos_monthly', 30, 'days', time_units='days') +#diag.add_file('atmos_monthly', 30, 'days', time_units='days') diag.add_file('atmos_daily', 1, 'days', time_units='days') # Define diag table entries @@ -40,6 +40,8 @@ diag.add_field('dynamics', 'sphum', time_avg=True) diag.add_field('atmosphere', 'precipitation', time_avg=True) +diag.add_field('atmosphere', 'dt_tg_convection', time_avg=True) + diag.add_field('mixed_layer', 't_surf', time_avg=True) diag.add_field('mixed_layer', 'flux_lhe', time_avg=True) diag.add_field('mixed_layer', 'flux_t', time_avg=True) @@ -102,7 +104,7 @@ 'tconst' : 285., 'prescribe_initial_dist':True, 'evaporation':False, - 'albedo_value': 0.25, + 'albedo_value': 0.7, 'depth':10., }, @@ -157,7 +159,7 @@ 'ir_tau_eq':0.2, 'ir_tau_pole':0.2, 'linear_tau': 1.0, - 'equinox_day':0.3032894472101727, + 'equinox_day':0.8032894472101727, }, @@ -183,7 +185,7 @@ 'astronomy_nml': { 'ecc':0.0935, 'obliq':25.19, - 'per':250.815799, #to be equivalent to peri_time = 0. in hs_forcing version of mars + 'per':109.18420099566217, #to be equivalent to peri_time = 0. in hs_forcing version of mars 'use_mean_anom_in_rrsun_calc':False, 'use_old_r_inv_squared':False }, @@ -201,13 +203,13 @@ if __name__=="__main__": - conv_schemes = ['dry'] + conv_schemes = ['dry', 'none'] scales = [ 1.] for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk15_fixed_year_length', codebase=cb) + exp = Experiment('grey_mars_mk16_fixed_year_length_shifted_thin_10m_conv_outputs_shiny_'+conv, codebase=cb) exp.clear_rundir() exp.diag_table = diag @@ -220,7 +222,7 @@ with exp_progress(exp, description='o%.0f d{day}' % scale): exp.run(1, use_restart=False, num_cores=NCORES) - for i in range(2, 241): + for i in range(2, 121): with exp_progress(exp, description='o%.0f d{day}' % scale): exp.run(i, num_cores=NCORES) notify('top down with conv scheme = '+conv+' has completed', 'isca') From 0593a8a3b7d9c085054edc4ebdb5d4c405f2d348 Mon Sep 17 00:00:00 2001 From: sit23 Date: Wed, 16 May 2018 16:47:13 +0100 Subject: [PATCH 09/69] Was having problems with output from daily averaging as obviously 1 day is not the same as one mars day. Have altered rotation rate and orbital rate such that each can be set as an integer, and we end up with a length of sol which is also an integer. This can then be used as the averaging period. --- exp/grey_mars/grey_mars.py | 32 +++++++++++-------- .../python/scripts/grey_rad_parm_check.py | 3 +- .../python/scripts/len_of_day_year_check.py | 30 +++++++++++++++++ src/shared/constants/constants.F90 | 11 ++++++- 4 files changed, 60 insertions(+), 16 deletions(-) create mode 100644 src/extra/python/scripts/len_of_day_year_check.py diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index 50db3b827..eec290a2d 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -25,7 +25,9 @@ # create a diagnostics output file for daily snapshots diag = DiagTable() #diag.add_file('atmos_monthly', 30, 'days', time_units='days') -diag.add_file('atmos_daily', 1, 'days', time_units='days') +diag.add_file('atmos_daily',88440 , 'seconds', time_units='days') +diag.add_file('atmos_e_daily',1 , 'days', time_units='days') +diag.add_file('atmos_timestep', 240, 'seconds', time_units='days') # Define diag table entries diag.add_field('dynamics', 'ps', time_avg=True) @@ -57,9 +59,9 @@ # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { - 'dt_atmos': 150, - 'days': 30., - 'seconds': 0, + 'dt_atmos': 220, + 'days': 0., + 'seconds': 88440., 'calendar': 'no_calendar' }, @@ -104,7 +106,7 @@ 'tconst' : 285., 'prescribe_initial_dist':True, 'evaporation':False, - 'albedo_value': 0.7, + 'albedo_value': 0.3, 'depth':10., }, @@ -160,6 +162,8 @@ 'ir_tau_pole':0.2, 'linear_tau': 1.0, 'equinox_day':0.8032894472101727, + 'use_time_average_coszen':True, + 'solar_constant':589.0, }, @@ -191,25 +195,25 @@ }, 'constants_nml': { - 'orbital_period': 686.980*86400., + 'orbital_period': 59166360, 'solar_const':589.0, 'radius':3396.0e3, - 'omega':7.088e-5, 'rdgas':192.0, 'kappa':0.22727, + 'rotation_period':88308, }, }) if __name__=="__main__": - conv_schemes = ['dry', 'none'] + conv_schemes = ['dry'] scales = [ 1.] for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk16_fixed_year_length_shifted_thin_10m_conv_outputs_shiny_'+conv, codebase=cb) + exp = Experiment('grey_mars_mk22_no_solar_abs_alb_0_3_av_test_'+conv, codebase=cb) exp.clear_rundir() exp.diag_table = diag @@ -220,9 +224,9 @@ exp.namelist['spectral_dynamics_nml']['reference_sea_level_press'] = scale * 610.0 exp.namelist['idealized_moist_phys_nml']['convection_scheme'] = conv - with exp_progress(exp, description='o%.0f d{day}' % scale): - exp.run(1, use_restart=False, num_cores=NCORES) - for i in range(2, 121): - with exp_progress(exp, description='o%.0f d{day}' % scale): - exp.run(i, num_cores=NCORES) +# with exp_progress(exp, description='o%.0f d{day}' % scale): + exp.run(1, use_restart=False, num_cores=NCORES) +# for i in range(2, 121): +# with exp_progress(exp, description='o%.0f d{day}' % scale): +# exp.run(i, num_cores=NCORES) notify('top down with conv scheme = '+conv+' has completed', 'isca') diff --git a/src/extra/python/scripts/grey_rad_parm_check.py b/src/extra/python/scripts/grey_rad_parm_check.py index ee12d247a..197af812b 100644 --- a/src/extra/python/scripts/grey_rad_parm_check.py +++ b/src/extra/python/scripts/grey_rad_parm_check.py @@ -31,10 +31,11 @@ def calc_t_surf(albedo_in, tau_s_in): cs = plt.contourf(albedo_range, tau_s_range, t_surf_arr) cn = plt.contour(albedo_range, tau_s_range, t_surf_arr,levels=[210., 260.], colors=['k']) -#plt.clabel(cn, inline=1, fontsize=10.) +plt.clabel(cn, inline=1, fontsize=10.) fig = plt.gcf() cb = fig.colorbar(cs) t_surf_value = calc_t_surf(0.25, 0.2) +t_surf_shiny = calc_t_surf(0.7, 0.2) diff --git a/src/extra/python/scripts/len_of_day_year_check.py b/src/extra/python/scripts/len_of_day_year_check.py new file mode 100644 index 000000000..2469a666c --- /dev/null +++ b/src/extra/python/scripts/len_of_day_year_check.py @@ -0,0 +1,30 @@ +import numpy as np +import matplotlib.pyplot as plt +import pdb + +n_day_range = np.arange(660, 680, 1) +len_day_range = np.arange(88000, 89000, 1) + +n_days = np.shape(n_day_range)[0] +n_len_day = np.shape(len_day_range)[0] + +n_day_arr = np.zeros((n_days, n_len_day)) +n_len_day_arr = np.zeros((n_days, n_len_day)) + +for al in range(n_len_day): + for ts in range(n_days): + n_day_arr[ts, al] = n_day_range[ts] + n_len_day_arr[ts,al] = len_day_range[al] + +integer_values = n_day_arr * n_len_day_arr / (n_day_arr -1) +res_arr = np.mod(integer_values,1) +where_integer = np.where(res_arr==0.0) + +cs = plt.contourf(n_len_day_arr, n_day_arr, integer_values) +#cn = plt.contour(albedo_range, tau_s_range, t_surf_arr,levels=[210., 260.], colors=['k']) +#plt.clabel(cn, inline=1, fontsize=10.) +fig = plt.gcf() + +cb = fig.colorbar(cs) + + diff --git a/src/shared/constants/constants.F90 b/src/shared/constants/constants.F90 index e3978c1f3..d63ebdf93 100644 --- a/src/shared/constants/constants.F90 +++ b/src/shared/constants/constants.F90 @@ -250,6 +250,7 @@ module constants_mod real, public :: RADIUS = 6376.0e3 real, public :: GRAV = EARTH_GRAV real, public :: OMEGA = EARTH_OMEGA ! planetary rotation rate in s^-1 +real, public :: ROTATION_PERIOD = -1.0 real, public :: ORBITAL_PERIOD = EARTH_ORBITAL_PERIOD ! orbital period (periapse -> periapse) in s real, public :: SECONDS_PER_SOL = SECONDS_PER_DAY real, public :: orbital_rate ! this is calculated from 2pi/orbital_period @@ -262,7 +263,7 @@ module constants_mod real, public :: CP_AIR = EARTH_CP_AIR logical :: earthday_multiple = .false. -namelist/constants_nml/ radius, grav, omega, orbital_period, pstd, pstd_mks, rdgas, kappa, solar_const, earthday_multiple +namelist/constants_nml/ radius, grav, omega, orbital_period, pstd, pstd_mks, rdgas, kappa, solar_const, earthday_multiple, rotation_period !----------------------------------------------------------------------- ! version and tagname published @@ -294,6 +295,14 @@ subroutine constants_init if (mpp_pe() == mpp_root_pe()) write (stdlog(),nml=constants_nml) + if (rotation_period.gt.0) then + omega = 2.*pi/rotation_period + + call error_mesg( 'constants_init','rotation_period has been specified, so omega calculated from rotation_period', NOTE) + + endif + + !> SECONDS_PER_SOL is the exoplanet equivalent of seconds_per_day. !! It is the number of seconds between sucessive solar zeniths at longitude 0. !! (The concept of seconds_per_day is kept as seconds per earth day From af6df4ad3bbbdcd9299750b59008fc7a5f7ef34f Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 18 May 2018 12:46:17 +0100 Subject: [PATCH 10/69] Making option for large scale condensation to be turned on and off independently. Then also only calculating rh when rh is asked for. Model now runs and looks vaguely mars-like, but dates of equinox etc still not right. --- exp/grey_mars/grey_mars.py | 30 ++++++++++++------- .../driver/solo/idealized_moist_phys.F90 | 17 +++++++---- 2 files changed, 31 insertions(+), 16 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index eec290a2d..2dbda39e3 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -26,8 +26,8 @@ diag = DiagTable() #diag.add_file('atmos_monthly', 30, 'days', time_units='days') diag.add_file('atmos_daily',88440 , 'seconds', time_units='days') -diag.add_file('atmos_e_daily',1 , 'days', time_units='days') -diag.add_file('atmos_timestep', 240, 'seconds', time_units='days') +#diag.add_file('atmos_e_daily',1 , 'days', time_units='days') +#diag.add_file('atmos_timestep', 240, 'seconds', time_units='days') # Define diag table entries diag.add_field('dynamics', 'ps', time_avg=True) @@ -59,9 +59,9 @@ # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { - 'dt_atmos': 220, + 'dt_atmos': 55, 'days': 0., - 'seconds': 88440., + 'seconds': 30.*88440., 'calendar': 'no_calendar' }, @@ -75,6 +75,7 @@ 'roughness_heat':3.21e-05, 'roughness_moist':0., 'two_stream_gray': True, #Use grey radiation + 'do_lscale_cond': False, }, 'vert_turb_driver_nml': { @@ -131,7 +132,9 @@ }, 'sat_vapor_pres_nml': { - 'do_simple':True + 'do_simple':True, + 'tcmin': -223, #Make sure low temperature limit of saturation vapour pressure is low enough that it doesn't cause an error (note that this giant planet has no moisture anyway, so doesn't directly affect calculation. + 'tcmax': 350., }, 'damping_driver_nml': { @@ -142,17 +145,22 @@ }, 'spectral_dynamics_nml': { - 'num_levels': 30, + 'num_levels': 25, 'exponent': 2.5, 'scale_heights': 4, 'surf_res': 0.1, 'robert_coeff': 4e-2, 'do_water_correction': False, - 'vert_coord_option': 'even_sigma', + 'vert_coord_option': 'input', 'initial_sphum': 0., 'valid_range_T': [0, 700] }, + 'vert_coordinate_nml': { + 'bk': [ 0.00000000e+00, 1.53008955e-04, 4.63790800e-04, 1.10977640e-03, 2.37044150e-03, 4.74479200e-03, 9.12245300e-03, 1.70677050e-02, 3.12516100e-02, 5.59939500e-02, 9.76165000e-02, 1.63754900e-01, 2.60315150e-01, 3.85974250e-01, 5.28084800e-01, 6.65956600e-01, 7.81088000e-01, 8.65400050e-01, 9.21109250e-01, 9.55343500e-01, 9.75416950e-01, 9.86856800e-01, 9.93269300e-01, 9.96830200e-01, 9.98799150e-01, 1.00000000e+00], + 'pk': [0.]*26, + }, + 'two_stream_gray_rad_nml': { 'rad_scheme': 'frierson', #Select radiation scheme to use, which in this case is Frierson 'do_seasonal': True, @@ -207,13 +215,13 @@ if __name__=="__main__": - conv_schemes = ['dry'] + conv_schemes = ['none'] scales = [ 1.] for conv in conv_schemes: for scale in scales: - exp = Experiment('grey_mars_mk22_no_solar_abs_alb_0_3_av_test_'+conv, codebase=cb) + exp = Experiment('grey_mars_mk24_macda_levels_'+conv, codebase=cb) exp.clear_rundir() exp.diag_table = diag @@ -226,7 +234,7 @@ # with exp_progress(exp, description='o%.0f d{day}' % scale): exp.run(1, use_restart=False, num_cores=NCORES) -# for i in range(2, 121): + for i in range(2, 121): # with exp_progress(exp, description='o%.0f d{day}' % scale): -# exp.run(i, num_cores=NCORES) + exp.run(i, num_cores=NCORES) notify('top down with conv scheme = '+conv+' has completed', 'isca') diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 273e6f061..9fe0c0f23 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -6,7 +6,7 @@ module idealized_moist_phys_mod use fms_mod, only: open_namelist_file, close_file #endif -use fms_mod, only: write_version_number, file_exist, close_file, stdlog, error_mesg, NOTE, FATAL, read_data, field_size, uppercase, mpp_pe +use fms_mod, only: write_version_number, file_exist, close_file, stdlog, error_mesg, NOTE, FATAL, WARNING, read_data, field_size, uppercase, mpp_pe use constants_mod, only: grav, rdgas, rvgas, cp_air, PSTD_MKS, dens_h2o !mj cp_air needed for rrtmg !s pstd_mks needed for pref calculation @@ -101,6 +101,8 @@ module idealized_moist_phys_mod logical :: do_bm = .false. logical :: do_ras = .false. +logical :: do_lscale_cond = .true. + !s Radiation options logical :: two_stream_gray = .true. logical :: do_rrtm_radiation = .false. @@ -142,7 +144,7 @@ module idealized_moist_phys_mod land_roughness_prefactor, & gp_surface, convection_scheme, & bucket, init_bucket_depth, init_bucket_depth_land, & !RG Add bucket - max_bucket_depth_land, robert_bucket, raw_bucket + max_bucket_depth_land, robert_bucket, raw_bucket, do_lscale_cond integer, parameter :: num_time_levels = 2 !RG Add bucket - number of time levels added to allow timestepping in this module @@ -674,6 +676,9 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l axes(1:2), Time, 'Rain from convection','kg/m/m/s') !endif +if (r_conv_scheme .eq. DRY_CONV .and. do_lscale_cond .eq. .true.) then + call error_mesg('idealized_moist_phys','do_lscale_cond is .true. but r_conv_scheme is dry. These options may not be consistent.', WARNING) +endif if(two_stream_gray) call two_stream_gray_rad_init(is, ie, js, je, num_levels, get_axis_id(), Time, rad_lonb_2d, rad_latb_2d, dt_real) @@ -863,7 +868,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg ! Perform large scale convection -if (r_conv_scheme .ne. DRY_CONV) then +if ( do_lscale_cond .eq. .true.) then ! Large scale convection is a function of humidity only. This is ! inconsistent with the dry convection scheme, don't run it! rain = 0.0; snow = 0.0 @@ -1131,8 +1136,10 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg endif ! if(turb) then !s Adding relative humidity calculation so as to allow comparison with Frierson's thesis. - call rh_calc (p_full(:,:,:,previous),tg_tmp,qg_tmp,RH) - if(id_rh >0) used = send_data(id_rh, RH*100., Time) + if (id_rh > 0) then + call rh_calc (p_full(:,:,:,previous),tg_tmp,qg_tmp,RH) + used = send_data(id_rh, RH*100., Time) + endif ! RG Add bucket From 1ca7b41f85e13970bf190a93c5155c11a006e9f0 Mon Sep 17 00:00:00 2001 From: sit23 Date: Tue, 29 May 2018 10:56:34 +0100 Subject: [PATCH 11/69] Various updates to mars_dev, includng error checking of surface flux nml, and ability to use specified temperatures for evaporation calculations in surface flux. This means that, when using a dry model where temperatures are outside standard range for sat specific humidity calculations, then false temps can be used to stop it failing. --- exp/grey_mars/grey_mars.py | 53 ++++++++++--------- .../two_stream_gray_rad.F90 | 16 ++++-- .../driver/solo/idealized_moist_phys.F90 | 5 +- src/coupler/surface_flux.F90 | 20 ++++++- src/shared/astronomy/astronomy.f90 | 9 ++-- 5 files changed, 64 insertions(+), 39 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index 2dbda39e3..fa92fee5d 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -54,12 +54,13 @@ diag.add_field('two_stream', 'swdn_toa', time_avg=True) diag.add_field('two_stream', 'time_since_ae', time_avg=True) diag.add_field('two_stream', 'true_anomaly', time_avg=True) - +diag.add_field('two_stream', 'dec', time_avg=True) +diag.add_field('two_stream', 'ang', time_avg=True) # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { - 'dt_atmos': 55, + 'dt_atmos': 220, 'days': 0., 'seconds': 30.*88440., 'calendar': 'no_calendar' @@ -95,8 +96,7 @@ 'use_virtual_temp': False, 'do_simple': True, 'old_dtaudv': True, - 'rh_target': 50., - 'delta_t_relax': 7200., + 'use_actual_surface_temperatures':False, }, 'atmosphere_nml': { @@ -108,7 +108,6 @@ 'prescribe_initial_dist':True, 'evaporation':False, 'albedo_value': 0.3, - 'depth':10., }, 'qe_moist_convection_nml': { @@ -169,7 +168,7 @@ 'ir_tau_eq':0.2, 'ir_tau_pole':0.2, 'linear_tau': 1.0, - 'equinox_day':0.8032894472101727, + 'equinox_day':0.0, 'use_time_average_coszen':True, 'solar_constant':589.0, }, @@ -197,8 +196,7 @@ 'astronomy_nml': { 'ecc':0.0935, 'obliq':25.19, - 'per':109.18420099566217, #to be equivalent to peri_time = 0. in hs_forcing version of mars - 'use_mean_anom_in_rrsun_calc':False, + 'use_mean_anom_in_rrsun_calc':True, 'use_old_r_inv_squared':False }, @@ -217,24 +215,31 @@ conv_schemes = ['none'] - scales = [ 1.] + depths = [10.,2., 0.5, 0.1] + + pers = [70.85] + + scale = 1. for conv in conv_schemes: - for scale in scales: - exp = Experiment('grey_mars_mk24_macda_levels_'+conv, codebase=cb) - exp.clear_rundir() - - exp.diag_table = diag - exp.namelist = namelist.copy() - exp.namelist['constants_nml']['grav'] = scale * 3.71 - exp.namelist['constants_nml']['pstd'] = scale * 6100.0 - exp.namelist['constants_nml']['pstd_mks'] = scale * 610.0 - exp.namelist['spectral_dynamics_nml']['reference_sea_level_press'] = scale * 610.0 - exp.namelist['idealized_moist_phys_nml']['convection_scheme'] = conv + for depth_val in depths: + for per_value in pers: + exp = Experiment('grey_mars_mk34_per_value'+str((per_value))+'_'+conv+'_mld_'+str(depth_val), codebase=cb) + exp.clear_rundir() + + exp.diag_table = diag + exp.namelist = namelist.copy() + exp.namelist['constants_nml']['grav'] = scale * 3.71 + exp.namelist['constants_nml']['pstd'] = scale * 6100.0 + exp.namelist['constants_nml']['pstd_mks'] = scale * 610.0 + exp.namelist['spectral_dynamics_nml']['reference_sea_level_press'] = scale * 610.0 + exp.namelist['idealized_moist_phys_nml']['convection_scheme'] = conv + exp.namelist['mixed_layer_nml']['depth'] = depth_val + exp.namelist['astronomy_nml']['per'] = per_value # with exp_progress(exp, description='o%.0f d{day}' % scale): - exp.run(1, use_restart=False, num_cores=NCORES) - for i in range(2, 121): + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2, 241): # with exp_progress(exp, description='o%.0f d{day}' % scale): - exp.run(i, num_cores=NCORES) - notify('top down with conv scheme = '+conv+' has completed', 'isca') + exp.run(i, num_cores=NCORES) + notify('top down with conv scheme = '+conv+' has completed', 'isca') diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index d2221fc81..62ddbce84 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -155,7 +155,7 @@ module two_stream_gray_rad_mod integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_net_lw_surf, id_lwdn_sfc, id_lwup_sfc, & id_tdt_rad, id_tdt_solar, id_flux_rad, id_flux_lw, id_flux_sw, id_coszen, id_fracsun, & id_lw_dtrans, id_lw_dtrans_win, id_sw_dtrans, id_co2, id_mars_solar_long, id_rrsun, id_true_anom, & - id_time_since_ae + id_time_since_ae, id_dec, id_ang character(len=10), parameter :: mod_name = 'two_stream' @@ -391,6 +391,12 @@ subroutine two_stream_gray_rad_init(is, ie, js, je, num_levels, axes, Time, lonb id_time_since_ae = register_diag_field ( mod_name, 'time_since_ae', & Time, 'time since ae', 'none') + id_dec = register_diag_field ( mod_name, 'dec', & + Time, 'dec', 'none') + + id_ang = register_diag_field ( mod_name, 'ang', & + Time, 'ang', 'none') + return end subroutine two_stream_gray_rad_init @@ -411,7 +417,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, integer :: i, j, k, n, dyofyr integer :: seconds, year_in_s, days -real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians, true_anomaly +real :: r_seconds, frac_of_day, frac_of_year, gmt, time_since_ae, rrsun, day_in_s, r_solday, r_total_seconds, r_days, r_dt_rad_avg, dt_rad_radians, true_anomaly, dec, ang_out logical :: used @@ -449,9 +455,9 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, r_dt_rad_avg=real(dt_rad_avg) dt_rad_radians = (r_dt_rad_avg/day_in_s)*2.0*pi - call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, dt_rad_radians, true_anom=true_anomaly) + call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, dt_rad_radians, true_anom=true_anomaly, dec=dec, ang=ang_out) else - call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, true_anom=true_anomaly) + call diurnal_solar(lat, lon, gmt, time_since_ae, coszen, fracsun, rrsun, true_anom=true_anomaly, dec=dec, ang=ang_out) end if insolation = solar_constant * coszen * rrsun @@ -460,6 +466,8 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, if (id_time_since_ae > 0) used = send_data ( id_time_since_ae, time_since_ae, Time_diag) + if (id_dec > 0) used = send_data ( id_dec, dec, Time_diag) + if (id_ang > 0) used = send_data ( id_ang, ang_out, Time_diag) if (id_true_anom > 0) used = send_data ( id_true_anom, true_anomaly, Time_diag) diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 9fe0c0f23..0d75a8514 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -1016,7 +1016,6 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg - !---------------------------------------------------------------------- ! Copied from MiMA physics_driver.f90 ! call damping_driver to calculate the various model dampings that @@ -1035,9 +1034,6 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg z_pbl) !s have taken the names of arrays etc from vert_turb_driver below. Watch ntp from 2006 call to this routine? endif - - - if(turb) then call vert_turb_driver( 1, 1, & @@ -1106,6 +1102,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg ! ! update surface temperature ! + if(mixed_layer_bc) then call mixed_layer( & Time, Time+Time_step, & diff --git a/src/coupler/surface_flux.F90 b/src/coupler/surface_flux.F90 index 048404e84..09c9edb83 100644 --- a/src/coupler/surface_flux.F90 +++ b/src/coupler/surface_flux.F90 @@ -264,6 +264,8 @@ module surface_flux_mod real :: land_humidity_prefactor = 1.0 !s Default is that land makes no difference to evaporative fluxes real :: land_evap_prefactor = 1.0 !s Default is that land makes no difference to evaporative fluxes +logical :: use_actual_surface_temperatures = .true. !Always true, apart from when running a dry model, where you can set this to false so that escomp is called with temperatures of 200k always, preventing bad temperature errors. + real :: flux_heat_gp = 5.7 !s Default value for Jupiter of 5.7 Wm^-2 real :: diabatic_acce = 1.0 !s Diabatic acceleration?? @@ -282,7 +284,7 @@ module surface_flux_mod land_humidity_prefactor, & !s Added to make land 'dry', i.e. to decrease the evaporative heat flux in areas of land. land_evap_prefactor, & !s Added to make land 'dry', i.e. to decrease the evaporative heat flux in areas of land. flux_heat_gp, & !s prescribed lower boundary heat flux on a giant planet - diabatic_acce + diabatic_acce, use_actual_surface_temperatures @@ -386,12 +388,13 @@ subroutine surface_flux_1d ( & integer :: i, nbad - if (do_init) call surface_flux_init !---- use local value of surf temp ---- t_surf0 = 200. ! avoids out-of-bounds in es lookup + +if (use_actual_surface_temperatures) then where (avail) where (land) t_surf0 = t_ca @@ -399,12 +402,24 @@ subroutine surface_flux_1d ( & t_surf0 = t_surf endwhere endwhere +endif t_surf1 = t_surf0 + del_temp call escomp ( t_surf0, e_sat ) ! saturation vapor pressure call escomp ( t_surf1, e_sat1 ) ! perturbed vapor pressure +if (.not. use_actual_surface_temperatures) then + where (avail) + where (land) + t_surf0 = t_ca + elsewhere + t_surf0 = t_surf + endwhere + endwhere +endif + + if(use_mixing_ratio) then ! surface mixing ratio at saturation q_sat = d622*e_sat /(p_surf-e_sat ) @@ -814,6 +829,7 @@ subroutine surface_flux_init ! read namelist #ifdef INTERNAL_FILE_NML read (input_nml_file, surface_flux_nml, iostat=io) + ierr = check_nml_error(io, 'surface_flux_nml') #else if ( file_exist('input.nml')) then unit = open_namelist_file () diff --git a/src/shared/astronomy/astronomy.f90 b/src/shared/astronomy/astronomy.f90 index 417c92c16..f281992a0 100644 --- a/src/shared/astronomy/astronomy.f90 +++ b/src/shared/astronomy/astronomy.f90 @@ -1124,7 +1124,7 @@ end subroutine get_ref_date_of_ae ! subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & fracday, rrsun, dt, allow_negative_cosz, & - half_day_out, true_anom) + half_day_out, true_anom, dec, ang) !--------------------------------------------------------------------- ! diurnal_solar_2d returns 2d fields of cosine of zenith angle, @@ -1140,14 +1140,13 @@ subroutine diurnal_solar_2d (lat, lon, gmt, time_since_ae, cosz, & real, intent(in), optional :: dt logical, intent(in), optional :: allow_negative_cosz real, dimension(:,:), intent(out), optional :: half_day_out -real, intent(out), optional :: true_anom +real, intent(out), optional :: true_anom, dec, ang !--------------------------------------------------------------------- ! local variables real, dimension(size(lat,1),size(lat,2)) :: t, tt, h, aa, bb, & st, stt, sh - real :: ang, dec logical :: Lallow_negative !--------------------------------------------------------------------- @@ -3347,7 +3346,7 @@ subroutine r_inv_squared_alt (ang, r_inv_squared_out, true_anomaly_out) endif if (present(true_anomaly_out)) then - true_anomaly_out = true_anomaly + true_anomaly_out = modulo(true_anomaly, twopi) endif @@ -3849,4 +3848,4 @@ end subroutine diurnal_exoplanet - end module astronomy_mod \ No newline at end of file + end module astronomy_mod From 30d31d930fcce26cf16dfbf990de2dd0392c19f9 Mon Sep 17 00:00:00 2001 From: sit23 Date: Wed, 30 May 2018 17:52:04 +0100 Subject: [PATCH 12/69] Simple script to interpolate mola mars altitude data into topography input file. --- src/extra/python/scripts/mars_topo_maker.py | 37 +++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 src/extra/python/scripts/mars_topo_maker.py diff --git a/src/extra/python/scripts/mars_topo_maker.py b/src/extra/python/scripts/mars_topo_maker.py new file mode 100644 index 000000000..830bb0574 --- /dev/null +++ b/src/extra/python/scripts/mars_topo_maker.py @@ -0,0 +1,37 @@ +import xarray as xar +import gauss_grid as gg +import numpy as np +import mpl_toolkits.basemap as basemap +import matplotlib.pyplot as plt +from netCDF4 import Dataset + +num_lat_out = 64 +num_lon_out = 128 + + +mola_original = xar.open_dataset('/scratch/sit204/planet_topo_files/mars/mola32.nc') + +longitudes_out = np.arange(0., 360., (360./num_lon_out)) +latitudes_out = gg.gaussian_latitudes(int(num_lat_out/2.))[0] + +lon_array, lat_array = np.meshgrid(longitudes_out, latitudes_out) + +output_array = basemap.interp(mola_original.alt.values[::-1,:], mola_original.longitude.values, mola_original.latitude.values[::-1], lon_array, lat_array, order=1) + + +nlat = latitudes_out.shape[0] +nlon = longitudes_out.shape[0] + +#Write land and topography arrays to file +topo_filename = './t42_mola_mars.nc' +topo_file = Dataset(topo_filename, 'w', format='NETCDF3_CLASSIC') +lat = topo_file.createDimension('lat', nlat) +lon = topo_file.createDimension('lon', nlon) +latitudes = topo_file.createVariable('lat','f4',('lat',)) +longitudes = topo_file.createVariable('lon','f4',('lon',)) +topo_array_netcdf = topo_file.createVariable('zsurf','f4',('lat','lon',)) +latitudes[:] = latitudes_out +longitudes[:] = longitudes_out +topo_array_netcdf[:] = output_array +topo_file.close() + From 32914527a1e6c1509ea3248e6ef37520154d0a40 Mon Sep 17 00:00:00 2001 From: sit23 Date: Wed, 30 May 2018 18:10:26 +0100 Subject: [PATCH 13/69] Model didn't like not having a land mask. Added this now. --- src/extra/python/scripts/mars_topo_maker.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/extra/python/scripts/mars_topo_maker.py b/src/extra/python/scripts/mars_topo_maker.py index 830bb0574..2f7e1c623 100644 --- a/src/extra/python/scripts/mars_topo_maker.py +++ b/src/extra/python/scripts/mars_topo_maker.py @@ -18,6 +18,7 @@ output_array = basemap.interp(mola_original.alt.values[::-1,:], mola_original.longitude.values, mola_original.latitude.values[::-1], lon_array, lat_array, order=1) +land_array = np.zeros_like(output_array) nlat = latitudes_out.shape[0] nlon = longitudes_out.shape[0] @@ -30,8 +31,10 @@ latitudes = topo_file.createVariable('lat','f4',('lat',)) longitudes = topo_file.createVariable('lon','f4',('lon',)) topo_array_netcdf = topo_file.createVariable('zsurf','f4',('lat','lon',)) +land_array_netcdf = topo_file.createVariable('land_mask','f4',('lat','lon',)) latitudes[:] = latitudes_out longitudes[:] = longitudes_out topo_array_netcdf[:] = output_array +land_array_netcdf[:] = land_array topo_file.close() From 4459f5ca261779c3bd7e2d5e1f2f8b9d0d04dada Mon Sep 17 00:00:00 2001 From: sit23 Date: Mon, 18 Jun 2018 07:55:28 +0100 Subject: [PATCH 14/69] Modernising run-plevel.py --- .../scripts/run_plevel.py | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/postprocessing/plevel_interpolation/scripts/run_plevel.py b/postprocessing/plevel_interpolation/scripts/run_plevel.py index 294ecb53a..addba9747 100644 --- a/postprocessing/plevel_interpolation/scripts/run_plevel.py +++ b/postprocessing/plevel_interpolation/scripts/run_plevel.py @@ -7,11 +7,11 @@ import subprocess start_time=time.time() -base_dir='/scratch/sit204/Data_2013/' -exp_name_list = ['no_ice_flux_lhe_exps_q_flux_hadgem_anoms_3'] +base_dir='/scratch/sit204/data_isca/' +exp_name_list = ['project_3_omega_normal'] avg_or_daily_list=['monthly'] -start_file=287 -end_file=288 +start_file=13 +end_file=24 nfiles=(end_file-start_file)+1 do_extra_averaging=False #If true, then 6hourly data is averaged into daily data using cdo @@ -64,15 +64,8 @@ for avg_or_daily in avg_or_daily_list: print(n+start_file) - number_prefix='' - - if n+start_file < 100: - number_prefix='0' - if n+start_file < 10: - number_prefix='00' - - nc_file_in = base_dir+'/'+exp_name+'/run'+number_prefix+str(n+start_file)+'/atmos_'+avg_or_daily+'.nc' - nc_file_out = out_dir+'/'+exp_name+'/run'+number_prefix+str(n+start_file)+'/atmos_'+avg_or_daily+file_suffix+'.nc' + nc_file_in = base_dir+'/'+exp_name+'/run%04d'%(n+start_file)+'/atmos_'+avg_or_daily+'.nc' + nc_file_out = out_dir+'/'+exp_name+'/run%04d'%(n+start_file)+'/atmos_'+avg_or_daily+file_suffix+'.nc' if not os.path.isfile(nc_file_out): plevel_call(nc_file_in,nc_file_out, var_names = var_names[avg_or_daily], p_levels = plevs[avg_or_daily], mask_below_surface_option=mask_below_surface_set) From a217c65e64e70b04c4b23d6be02fd34fde913d06 Mon Sep 17 00:00:00 2001 From: James Penn Date: Tue, 17 Jul 2018 19:30:53 +0100 Subject: [PATCH 15/69] Compilation on Mac Cores on a Mac cannot be tied in the same way that they can on linux, so core affinity is not possible. Here I am just stubbing out the core affinity functions to mimic the linux case where it cannot get the affinity info. More info in links inline in `affinity.c`. --- src/extra/env/docker | 1 + src/extra/env/isca | 9 +++++---- src/extra/env/mac_brew | 7 +++++++ src/extra/python/isca/templates/compile.sh | 6 +++--- .../python/isca/templates/mkmf.template.gfort | 3 +-- src/shared/mpp/affinity.c | 14 +++++++++++++- 6 files changed, 30 insertions(+), 10 deletions(-) create mode 100644 src/extra/env/mac_brew diff --git a/src/extra/env/docker b/src/extra/env/docker index eefeef6f0..0aa43a992 100755 --- a/src/extra/env/docker +++ b/src/extra/env/docker @@ -4,3 +4,4 @@ export GFDL_COMPILER=gfort export F90=mpifort export CC=mpicc +export NETCDF_LIBS = `nf-config --fflags --flibs` \ No newline at end of file diff --git a/src/extra/env/isca b/src/extra/env/isca index 9ceb3dbf8..a4cf1f14c 100644 --- a/src/extra/env/isca +++ b/src/extra/env/isca @@ -1,9 +1,10 @@ echo loadmodules for Isca -export F90=mpiifort -export CC=mpiicc - module purge module load intel/2016b module load HDF5/1.8.16-intel-2016b -module load netCDF-Fortran/4.4.2-intel-2016b \ No newline at end of file +module load netCDF-Fortran/4.4.2-intel-2016b + +export F90=mpiifort +export CC=mpiicc +export NETCDF_LIBS = `nf-config --fflags --flibs` \ No newline at end of file diff --git a/src/extra/env/mac_brew b/src/extra/env/mac_brew new file mode 100644 index 000000000..4beff9c56 --- /dev/null +++ b/src/extra/env/mac_brew @@ -0,0 +1,7 @@ +echo loadmodules for Isca on mac using homebrew + +export GFDL_COMPILER=gfort +export F90=mpifort +export CC=gcc-8 +export CDEFS=" -DMACOS " +export NETCDF_LIBS=" -I/usr/local/include " \ No newline at end of file diff --git a/src/extra/python/isca/templates/compile.sh b/src/extra/python/isca/templates/compile.sh index 11c67c1d4..66cd750bb 100755 --- a/src/extra/python/isca/templates/compile.sh +++ b/src/extra/python/isca/templates/compile.sh @@ -19,9 +19,9 @@ template_debug={{ template_dir }}/mkmf.template.debug execdir={{ execdir }} # where code is compiled and executable is created executable={{ executable_name }} -netcdf_flags=`nf-config --fflags --flibs` +#netcdf_flags=`nf-config --fflags --flibs` -ulimit -s unlimited # Set stack size to unlimited +#ulimit -s unlimited # Set stack size to unlimited export MALLOC_CHECK_=0 # 3. compile the mppncombine tool if it hasn't yet been done. @@ -49,7 +49,7 @@ cd $execdir echo $pathnames # execute mkmf to create makefile -cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 {{compile_flags}}" +cppDefs="-Duse_libMPI -Duse_netCDF -Duse_LARGEFILE -DINTERNAL_FILE_NML -DOVERLOAD_C8 ${CDEFS} {{compile_flags}}" $mkmf -a $sourcedir -t $template -p $executable -c "$cppDefs" $pathnames $sourcedir/shared/include $sourcedir/shared/mpp/include make diff --git a/src/extra/python/isca/templates/mkmf.template.gfort b/src/extra/python/isca/templates/mkmf.template.gfort index b35336dbc..f74bcaf93 100755 --- a/src/extra/python/isca/templates/mkmf.template.gfort +++ b/src/extra/python/isca/templates/mkmf.template.gfort @@ -2,7 +2,6 @@ # typical use with mkmf # mkmf -t template.ifc -c"-Duse_libMPI -Duse_netCDF" path_names /usr/local/include CPPFLAGS = -I/usr/local/include -NETCDF_LIBS = `nf-config --fflags --flibs` # FFLAGS: # -cpp: Use the fortran preprocessor @@ -23,5 +22,5 @@ FFLAGS = $(CPPFLAGS) $(NETCDF_LIBS) -cpp -fcray-pointer \ FC = $(F90) LD = $(F90) $(NETCDF_LIBS) -LDFLAGS = -lnetcdff -lnetcdf -lmpi +LDFLAGS = -lnetcdff -lnetcdf -lmpi CFLAGS = -D__IFC diff --git a/src/shared/mpp/affinity.c b/src/shared/mpp/affinity.c index f0eb78149..9ea3ba3a4 100644 --- a/src/shared/mpp/affinity.c +++ b/src/shared/mpp/affinity.c @@ -21,6 +21,16 @@ #define _GNU_SOURCE +#ifdef MACOS +// Mac OS doesn't permit thread pinning. So we just ignore it... +// https://developer.apple.com/library/archive/releasenotes/Performance/RN-AffinityAPI/index.html +// see also Neil's question on Stack Overflow +// https://stackoverflow.com/questions/45227009/alternative-to-shed-getaffinity-cpu-set-t-etc +int get_cpu_affinity(void) { return -1; }; +void set_cpu_affinity( int cpu ) {}; + +#else + #include #include #include @@ -64,7 +74,6 @@ int get_cpu_affinity(void) return (last_cpu == -1) ? first_cpu : -1; } -int get_cpu_affinity_(void) { return get_cpu_affinity(); } /* Fortran interface */ /* @@ -81,4 +90,7 @@ void set_cpu_affinity( int cpu ) } } +#endif + +int get_cpu_affinity_(void) { return get_cpu_affinity(); } /* Fortran interface */ void set_cpu_affinity_(int *cpu) { set_cpu_affinity(*cpu); } /* Fortran interface */ From 9caa54f72c2be88846752326747792d9e0a05d05 Mon Sep 17 00:00:00 2001 From: sit23 Date: Tue, 31 Jul 2018 12:31:31 +0100 Subject: [PATCH 16/69] Changes necessary to run on stephen macbook pro. --- postprocessing/compile_mppn.sh | 4 +++- .../two_stream_gray_rad/two_stream_gray_rad.F90 | 2 +- src/extra/python/isca/experiment.py | 7 +++++++ src/extra/python/isca/templates/mkmf.template.gfort | 2 +- 4 files changed, 12 insertions(+), 3 deletions(-) diff --git a/postprocessing/compile_mppn.sh b/postprocessing/compile_mppn.sh index 981ee9950..58e6b7eef 100755 --- a/postprocessing/compile_mppn.sh +++ b/postprocessing/compile_mppn.sh @@ -11,7 +11,9 @@ ppdir=./ # path to directory containing the tool for combining d # 2. Load the necessary tools into the environment source $GFDL_BASE/src/extra/env/$GFDL_ENV -netcdf_flags=`nc-config --cflags --libs` +netcdf_flags=`nc-config --cflags` +netcdf_flags+=' -L/usr/local/Cellar/netcdf/4.6.1_2/lib -lnetcdf' +echo ${netcdf_flags} #-------------------------------------------------------------------------------------------------------- # compile combine tool #cd $ppdir diff --git a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 index 84215b6e2..bff89c1ba 100644 --- a/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 +++ b/src/atmos_param/two_stream_gray_rad/two_stream_gray_rad.F90 @@ -514,7 +514,6 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, ! longwave source function b = stefan*t**4 -lw_dtrans_win = 1. if(do_read_co2)then call interpolator( co2_interp, Time_diag, p_half, co2f, trim(co2_variable_name)) @@ -526,6 +525,7 @@ subroutine two_stream_gray_rad_down (is, js, Time_diag, lat, lon, p_half, t, case(B_GEEN) ! split LW in 2 bands: water-vapour window and remaining = non-window ! ref: Ruth Geen etal, GRL 2016 (supp. information). + lw_dtrans_win = 1. do k = 1, n lw_del_tau = ( ir_tau_co2 + 0.2023 * log(carbon_conc/360.) & + ir_tau_wv1*log(ir_tau_wv2*q(:,:,k) + 1) ) & diff --git a/src/extra/python/isca/experiment.py b/src/extra/python/isca/experiment.py index 5176bc98f..ad941517b 100755 --- a/src/extra/python/isca/experiment.py +++ b/src/extra/python/isca/experiment.py @@ -320,6 +320,13 @@ def _outhandler(line): self.log.debug("Restart file %s combined" % restartfile) self.emit('run:combined', self, i) + else: + for file in self.diag_table.files: + netcdf_file = '%s.nc' % file + filebase = P(self.rundir, netcdf_file) + sh.cp(filebase, P(outdir, netcdf_file)) + sh.rm(glob.glob(filebase+'*')) + self.log.debug('%s copied to data directory' % netcdf_file) # make the restart archive and delete the restart files self.make_restart_archive(self.get_restart_file(i), resdir) diff --git a/src/extra/python/isca/templates/mkmf.template.gfort b/src/extra/python/isca/templates/mkmf.template.gfort index f74bcaf93..63500db3b 100755 --- a/src/extra/python/isca/templates/mkmf.template.gfort +++ b/src/extra/python/isca/templates/mkmf.template.gfort @@ -22,5 +22,5 @@ FFLAGS = $(CPPFLAGS) $(NETCDF_LIBS) -cpp -fcray-pointer \ FC = $(F90) LD = $(F90) $(NETCDF_LIBS) -LDFLAGS = -lnetcdff -lnetcdf -lmpi +LDFLAGS = -lnetcdff -lnetcdf -lmpi -L/usr/local/Cellar/netcdf/4.6.1_2/lib CFLAGS = -D__IFC From 0be5bed8774b505f4cd9ccc6b747d346aa29fc11 Mon Sep 17 00:00:00 2001 From: sit23 Date: Mon, 6 Aug 2018 09:51:49 +0100 Subject: [PATCH 17/69] Found stable configuration of Mars model with stronger sponge. --- exp/grey_mars/grey_mars.py | 8 ++++---- src/extra/python/scripts/mars_topo_maker.py | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/exp/grey_mars/grey_mars.py b/exp/grey_mars/grey_mars.py index fa92fee5d..0b45ada0b 100644 --- a/exp/grey_mars/grey_mars.py +++ b/exp/grey_mars/grey_mars.py @@ -60,7 +60,7 @@ # define namelist values as python dictionary namelist = Namelist({ 'main_nml': { - 'dt_atmos': 220, + 'dt_atmos': 110, 'days': 0., 'seconds': 30.*88440., 'calendar': 'no_calendar' @@ -138,7 +138,7 @@ 'damping_driver_nml': { 'do_rayleigh': True, - 'trayfric': -0.25, # neg. value: time in *days* + 'trayfric': -0.125, # neg. value: time in *days* 'sponge_pbottom': 0.5, #Bottom of the model's sponge down to 50hPa (units are Pa) 'do_conserve_energy': True, }, @@ -215,7 +215,7 @@ conv_schemes = ['none'] - depths = [10.,2., 0.5, 0.1] + depths = [2.] pers = [70.85] @@ -224,7 +224,7 @@ for conv in conv_schemes: for depth_val in depths: for per_value in pers: - exp = Experiment('grey_mars_mk34_per_value'+str((per_value))+'_'+conv+'_mld_'+str(depth_val), codebase=cb) + exp = Experiment('grey_mars_mk36_per_value'+str((per_value))+'_'+conv+'_mld_'+str(depth_val), codebase=cb) exp.clear_rundir() exp.diag_table = diag diff --git a/src/extra/python/scripts/mars_topo_maker.py b/src/extra/python/scripts/mars_topo_maker.py index 2f7e1c623..2f14eee27 100644 --- a/src/extra/python/scripts/mars_topo_maker.py +++ b/src/extra/python/scripts/mars_topo_maker.py @@ -5,8 +5,8 @@ import matplotlib.pyplot as plt from netCDF4 import Dataset -num_lat_out = 64 -num_lon_out = 128 +num_lat_out = 128 +num_lon_out = 256 mola_original = xar.open_dataset('/scratch/sit204/planet_topo_files/mars/mola32.nc') @@ -24,7 +24,7 @@ nlon = longitudes_out.shape[0] #Write land and topography arrays to file -topo_filename = './t42_mola_mars.nc' +topo_filename = './t85_mola_mars.nc' topo_file = Dataset(topo_filename, 'w', format='NETCDF3_CLASSIC') lat = topo_file.createDimension('lat', nlat) lon = topo_file.createDimension('lon', nlon) From 19bc997b4c18fec8852ef5420cd25d93ed9a419d Mon Sep 17 00:00:00 2001 From: sit23 Date: Tue, 7 Aug 2018 14:35:44 +0100 Subject: [PATCH 18/69] Added optin to plevel_fn to add back variables that have the scalar_axis. For some reason, the interpolator will get rid of such variables when used on its own, so I modified the python to add them back. --- .../plevel_interpolation/scripts/plevel_fn.py | 35 ++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/postprocessing/plevel_interpolation/scripts/plevel_fn.py b/postprocessing/plevel_interpolation/scripts/plevel_fn.py index d6266cc3d..31f5d03fd 100644 --- a/postprocessing/plevel_interpolation/scripts/plevel_fn.py +++ b/postprocessing/plevel_interpolation/scripts/plevel_fn.py @@ -5,8 +5,10 @@ import sys import os import sh +import xarray as xar +import pdb -def plevel_call(nc_file_in,nc_file_out, var_names = '-a', p_levels='default', mask_below_surface_option=' '): +def plevel_call(nc_file_in,nc_file_out, var_names = '-a', p_levels='default', mask_below_surface_option=' ', add_back_scalar_axis_vars=False): check_gfdl_directories_set() @@ -23,6 +25,37 @@ def plevel_call(nc_file_in,nc_file_out, var_names = '-a', p_levels='default', ma command = interper + nc_file + out_file + plev +' '+mask_below_surface_option+ var_names print(command) subprocess.call([command], shell=True) + + if add_back_scalar_axis_vars: + add_back_scalar_axis_vars_fn(nc_file_in, nc_file_out) + +def add_back_scalar_axis_vars_fn(file_in, file_out): + ''' For some reason the plevel interpolator will not put variables + with `scalar_axis` as a dimension in the interpolated output file. + This piece of python adds them back after the interpolation. Particularly + important for Mars work. + ''' + + ds_in = xar.open_dataset(file_in, decode_times=False) + ds_out = xar.open_dataset(file_out, decode_times=False) + + try: + ds_in.dims['scalar_axis'] + except KeyError: + return + + list_of_vars_to_copy=[] + + for name in ds_in.var().keys(): + if 'scalar_axis' in ds_in[name].dims and name not in ds_out.var().keys(): + list_of_vars_to_copy.append(name) + + ds_out.coords['scalar_axis'] = ('scalar_axis', ds_in['scalar_axis'].values) + + for out_name in list_of_vars_to_copy: + ds_out[out_name] = (ds_in[out_name].dims, ds_in[out_name].values) + + ds_out.to_netcdf(path=file_out) def daily_average(nc_file_in, nc_file_out): subprocess.call('cdo daymean '+nc_file_in+' '+nc_file_out, shell=True) From c829ed8fce3dc8800e26c07f1991c6d7db53867d Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 9 Aug 2018 16:08:20 +0100 Subject: [PATCH 19/69] Managed to resurrect the local local_heating implementation I put into idealized moist phys back in teh old GFDLMoistmodel repo on the local_heating_dev branch. This seems to work now, and should hopefully be useful. Only tried the Isidoro option, and not the input file option, but this should be fine. --- src/atmos_param/hs_forcing/hs_forcing.F90 | 10 ++++---- .../driver/solo/idealized_moist_phys.F90 | 23 +++++++++++++++++-- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/atmos_param/hs_forcing/hs_forcing.F90 b/src/atmos_param/hs_forcing/hs_forcing.F90 index 04f2bdf43..226dafe3c 100644 --- a/src/atmos_param/hs_forcing/hs_forcing.F90 +++ b/src/atmos_param/hs_forcing/hs_forcing.F90 @@ -58,7 +58,7 @@ module hs_forcing_mod !----------------------------------------------------------------------- !---------- interfaces ------------ - public :: hs_forcing, hs_forcing_init, hs_forcing_end + public :: hs_forcing, hs_forcing_init, hs_forcing_end, local_heating type(interpolate_type),save :: heating_source_interp type(interpolate_type),save :: u_interp @@ -229,8 +229,6 @@ subroutine hs_forcing ( is, ie, js, je, dt, Time, lon, lat, p_half, p_full, & if(trim(local_heating_option) /= '') then call local_heating ( Time, is, js, lon, lat, ps, p_full, p_half, ttnd ) tdt = tdt + ttnd -! if (id_local_heating > 0) used = send_data ( id_local_heating, ttnd, Time, is, js) - if (id_local_heating > 0) used = send_data ( id_local_heating, ttnd, Time) endif ! if (id_tdt > 0) used = send_data ( id_tdt, tdt, Time, is, js) @@ -735,6 +733,8 @@ subroutine local_heating ( Time, is, js, lon, lat, ps, p_full, p_half, tdt ) real, dimension(size(lon,1),size(lon,2)) :: lon_factor real, dimension(size(lat,1),size(lat,2)) :: lat_factor real, dimension(size(p_half,1),size(p_half,2),size(p_half,3)) :: p_half2 +logical :: used + do i=1,size(p_half,3) p_half2(:,:,i)=p_half(:,:,size(p_half,3)-i+1) enddo @@ -762,6 +762,8 @@ subroutine local_heating ( Time, is, js, lon, lat, ps, p_full, p_half, tdt ) call error_mesg ('hs_forcing_nml','"'//trim(local_heating_option)//'" is not a valid value for local_heating_option',FATAL) endif +if (id_local_heating > 0) used = send_data ( id_local_heating, tdt, Time) + end subroutine local_heating !####################################################################### @@ -1021,4 +1023,4 @@ subroutine top_down_newtonian_damping ( Time, lat, ps, p_full, p_half, t, tdt, t end subroutine top_down_newtonian_damping -end module hs_forcing_mod \ No newline at end of file +end module hs_forcing_mod diff --git a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 index 273e6f061..20082758e 100644 --- a/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 +++ b/src/atmos_spectral/driver/solo/idealized_moist_phys.F90 @@ -54,6 +54,8 @@ module idealized_moist_phys_mod use rayleigh_bottom_drag_mod, only: rayleigh_bottom_drag_init, compute_rayleigh_bottom_drag +use hs_forcing_mod, only: hs_forcing_init, local_heating + #ifdef RRTM_NO_COMPILE ! RRTM_NO_COMPILE not included #else @@ -134,6 +136,10 @@ module idealized_moist_phys_mod real :: raw_bucket = 0.53 ! default raw coefficient for bucket depth LJJ ! end RG Add bucket +!s Adding localised heating option from Held-Suarez +logical :: do_local_heating = .false. +!s end Adding localised heating option from Held-Suarez + namelist / idealized_moist_phys_nml / turb, lwet_convection, do_bm, do_ras, roughness_heat, & two_stream_gray, do_rrtm_radiation, do_damping,& mixed_layer_bc, do_simple, & @@ -142,7 +148,8 @@ module idealized_moist_phys_mod land_roughness_prefactor, & gp_surface, convection_scheme, & bucket, init_bucket_depth, init_bucket_depth_land, & !RG Add bucket - max_bucket_depth_land, robert_bucket, raw_bucket + max_bucket_depth_land, robert_bucket, raw_bucket, & + do_local_heating integer, parameter :: num_time_levels = 2 !RG Add bucket - number of time levels added to allow timestepping in this module @@ -712,6 +719,11 @@ subroutine idealized_moist_phys_init(Time, Time_step_in, nhum, rad_lon_2d, rad_l id_rh = register_diag_field ( mod_name, 'rh', & axes(1:3), Time, 'relative humidity', 'percent') +if (do_local_heating) then + call hs_forcing_init(get_axis_id(), Time, rad_lonb_2d, rad_latb_2d, rad_lat_2d) +endif + + end subroutine idealized_moist_phys_init !================================================================================================================================= subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg, grid_tracers, & @@ -725,7 +737,7 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg real, dimension(:,:,:,:), intent(inout) :: dt_tracers real :: delta_t -real, dimension(size(ug,1), size(ug,2), size(ug,3)) :: tg_tmp, qg_tmp, RH,tg_interp, mc, dt_ug_conv, dt_vg_conv +real, dimension(size(ug,1), size(ug,2), size(ug,3)) :: tg_tmp, qg_tmp, RH,tg_interp, mc, dt_ug_conv, dt_vg_conv, tdt_local_heating real, intent(in) , dimension(:,:,:), optional :: mask @@ -991,6 +1003,13 @@ subroutine idealized_moist_phys(Time, p_half, p_full, z_half, z_full, ug, vg, tg endif #endif +if (do_local_heating) then + call local_heating ( Time, is, js, rad_lon, rad_lat, & + p_half(:,:,num_levels+1,current), p_full(:,:,:,current), & + p_half(:,:,:,current), tdt_local_heating ) + dt_tg = dt_tg + tdt_local_heating +endif + if(gp_surface) then From 4901b0b05b1b9efccc4cf698fd2cbf2ecd0cee7b Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 23 Aug 2018 17:03:44 +0100 Subject: [PATCH 20/69] Modified ozone cmip timeseries script to allow for the creation of local prescribed heating input files. Seems to work alright, but problem is that local heating code always reads zeros no matter what I seem to do. Reading the created files into RRTM as an ozone does work, and reading ozone-1990 into local heating does not. Very odd. Tried everything I can think of, but we'll have to carry on testing. --- .../create_local_heating_input_file.py | 99 ++++ src/extra/python/scripts/finite_difference.py | 422 ++++++++++++++++++ 2 files changed, 521 insertions(+) create mode 100644 src/extra/python/scripts/create_local_heating_input_file.py create mode 100644 src/extra/python/scripts/finite_difference.py diff --git a/src/extra/python/scripts/create_local_heating_input_file.py b/src/extra/python/scripts/create_local_heating_input_file.py new file mode 100644 index 000000000..3f0360671 --- /dev/null +++ b/src/extra/python/scripts/create_local_heating_input_file.py @@ -0,0 +1,99 @@ +# -*- coding: utf-8 -*-s +import numpy as np +from calendar_calc import day_number_to_date +from netCDF4 import date2num +import xarray as xar +import pdb +import create_timeseries as cts +from finite_difference import diffz + +grav = 9.81 +c_p = 287.04 / (2./7.) + +resolution_file = xar.open_dataset('~/Desktop/full_continents_global_monthly_control_test.nc', decode_times=False) + +lons = resolution_file['lon'] +lats = resolution_file['lat'] + +try: + resolution_file['htr_lw'] + heating_rates_present=True +except KeyError: + heating_rates_present=False + +if not heating_rates_present: + total_flux = resolution_file['lw_recalculated']+resolution_file['sw_recalculated'] + dflux_dp = diffz(total_flux.values, resolution_file['phalf'].values*100., axis=1) + heating_rate = (grav/c_p)*dflux_dp +else: + heating_rate = resolution_file['htr_lw']+resolution_file['htr_sw'] + +# ozone_in = resolution_file.variables['O3'][:] + +latb_temp=np.zeros(lats.shape[0]+1) + +for tick in np.arange(1,lats.shape[0]): + latb_temp[tick]=(lats[tick-1]+lats[tick])/2. + +lonb_temp=np.zeros(lons.shape[0]+1) + +for tick in np.arange(1,lons.shape[0]): + lonb_temp[tick]=(lons[tick-1]+lons[tick])/2. + +latb_temp[0]=-90. +latb_temp[-1]=90. + +lonb_temp[0]=0. +lonb_temp[-1]=360. + +latbs=latb_temp +lonbs=lonb_temp + +nlon=lons.shape[0] +nlat=lats.shape[0] + +nlonb=len(lonbs) +nlatb=latbs.shape[0] + +p_full = resolution_file['pfull'][::-1] + +p_half = resolution_file['phalf'][::-1] + +time_arr = resolution_file['time'].values + +heating_rate_new = np.zeros((12, heating_rate.shape[1], nlat, nlon)) +time_arr_new = np.arange(0,12,1)*30. +time_arr = time_arr_new + +for t in range(len(time_arr)): + heating_rate_new[t,...] = heating_rate + +heating_rate = heating_rate_new[:,1::,...] + +date_arr=np.mod((time_arr-time_arr[0]),12) +date_arr_new=np.zeros(12) + +#Find grid and time numbers + +npfull=p_full.shape[0] +nphalf=p_half.shape[0] + +ntime=len(time_arr) + +#Output it to a netcdf file. +file_name='full_continents_heating_input_file.nc' +variable_name='heating_rate' + +number_dict={} +number_dict['nlat']=nlat +number_dict['nlon']=nlon +number_dict['nlatb']=nlatb +number_dict['nlonb']=nlonb +number_dict['npfull']=npfull +number_dict['nphalf']=nphalf +number_dict['ntime']=ntime + +time_units='days since 0000-01-01 00:00:00.0' + +cts.output_to_file(heating_rate[:,::-1,...],lats,lons,latbs,lonbs,p_full,p_half,time_arr,time_units,file_name,variable_name,number_dict) + diff --git a/src/extra/python/scripts/finite_difference.py b/src/extra/python/scripts/finite_difference.py new file mode 100644 index 000000000..fa97f3504 --- /dev/null +++ b/src/extra/python/scripts/finite_difference.py @@ -0,0 +1,422 @@ +""" +Routines associated with finite differencing on the sphere. +""" +import numpy as np + +# Static parameters +RAD = np.pi / 180.0 +EARTH_R = 6.371e6 + + +class NDSlicer(object): + """N-Dimensional slice class for numpy arrays.""" + + def __init__(self, axis, ndim, start=None, stop=None, step=None): + """ + Create an n-dimensional slice list. + + Parameters + ---------- + axis : integer + Axis on which to apply the slice + ndim : integer + Total number of dimensions of array to be sliced + start, stop, step : integer, optional + Index of beginning, stop and step width of the slice [start:stop:step] + default for each is None. + + """ + self.axis = axis + self.ndim = ndim + self.start = start + self.stop = stop + self.step = step + self.slicer = None + self.__getitem__(slice(start, stop, step)) + + def __getitem__(self, key): + """ + Create an n-dimensional slice list. + + Parameters + ---------- + axis : integer + Axis on which to apply the slice + ndim : integer + Total number of dimensions of array to be sliced + start, stop, step : integer, optional + Index of beginning, stop and step width of the slice [start:stop:step] + default for each is None. + + Returns + ------- + slicer : list + list of slices such that all data at other axes are kept, one axis is sliced + + Examples + -------- + Create random array, slice it:: + + x = np.random.randn(5, 3) + + # Create slicer equivalent to [1:-1, :] + slc = NDSlicer(0, x.ndim) + print(x) + [[ 0.68470539 0.87880216 -0.45086367] + [ 1.06804045 0.63094676 -0.76633033] + [-1.69841915 0.35207064 -0.4582049 ] + [-0.56431067 0.62833728 -0.04101542] + [-0.02760744 2.02814338 0.13195714]] + print(x[slc[1:-1]]) + [[ 1.06804045 0.63094676 -0.76633033] + [-1.69841915 0.35207064 -0.4582049 ] + [-0.56431067 0.62833728 -0.04101542]] + + """ + if isinstance(key, slice): + self.start = key.start + self.stop = key.stop + self.step = key.step + elif isinstance(key, int): + self.start = key + self.stop = key + 1 + self.step = None + + self.slicer = [slice(None)] * self.ndim + self.slicer[self.axis] = slice(self.start, self.stop, self.step) + return self.slicer + + def slice(self, start=None, stop=None, step=None): + """Legacy compatibility method, calls `__getitem__`.""" + self.__getitem__(slice(start, stop, step)) + + +def cfd(data, lons, lats, cyclic=False): + """ + Vectorized central finite difference for N-D array data by lon / lat. + + Parameters + ---------- + data : array_like + N-dimensional array to be differentiated + lons, lats : array_like + 1D coordinate arrays for longitude and latitude dimensions + cyclic : Boolean, optional + Data is cyclic on if true + + Returns + ------- + diff_lon, diff_lat : array_like + ND arrays of longitudinal, latitudinal centered finite differences of `data` + respectively, with same dimensionality as `data` + + """ + # Find the axis where `data` shape matches longitudes + axis_x = np.where(np.array(data.shape) == lons.shape[0])[0][0] + + # Find the axis where `data` shape matches latitudes + axis_y = np.where(np.array(data.shape) == lats.shape[0])[0][0] + + dlong, dlatg = dlon_dlat(lons, lats, cyclic=cyclic) + diff_x = diff_cfd(data, axis_x, cyclic=cyclic) / dlong + diff_y = diff_cfd(data, axis_y, cyclic=False) / dlatg + + return diff_x, diff_y + + +def convert_radians_latlon(lat, lon): + """ + Convert input lat/lon array to radians if input is degrees, do nothing if radians. + + Parameters + ---------- + lat : array_like + ND array of latitude + lon : array_like + ND array of longitude + + Returns + ---------- + lat : array_like + ND array of latitude in radians + lon : array_like + ND array of longitude in radians + + """ + if (np.max(np.abs(lat)) - np.pi / 2.0) > 1.0: + lat_out = lat * RAD + else: + lat_out = lat + + if(np.min(lon) < 0 and np.max(lon) > 0 and + np.abs(np.max(np.abs(lon)) - np.pi) > np.pi): + lon_out = lon * RAD + elif np.abs(np.max(np.abs(lon)) - np.pi * 2) > np.pi: + lon_out = lon * RAD + else: + lon_out = lon + + return lat_out, lon_out + + +def dlon_dlat(lon, lat, cyclic=True): + """ + Compute horizontal center finite differences of latitude/longitude on Earth [m]. + + Parameters + ---------- + lon, lat : array_like + 1D array of longitudes and latitudes respectively + cyclic : boolean + If True, longitudes are cyclic + + Returns + ------- + dlong, dlatg : array_like + Longitudinal, latitudinal centered finite differences in m + Shapes are (lat.shape[0], lon.shape[0]) + + """ + # Check that lat/lon are in radians + lat, lon = convert_radians_latlon(lat, lon) + + # Calculate centre finite difference of lon / lat + dlon = lon[2:] - lon[:-2] + dlat = lat[2:] - lat[:-2] + + # If we want cyclic data, repeat dlon[0] and dlon[-1] at edges + if cyclic: + dlon = np.append(dlon[0], dlon) # cyclic boundary in East + dlon = np.append(dlon, dlon[-1]) # cyclic boundary in West + _, lat2d = np.meshgrid(lon, lat) + else: + _, lat2d = np.meshgrid(lon[1:-1], lat) + + dlat = np.append(lat[1] - lat[0], dlat) # boundary in South + dlat = np.append(dlat, lat[-1] - lat[-2]) # boundary in North + dlong, dlatg = np.meshgrid(dlon, dlat) + + # Lon/Lat differences in spherical coords + dlong *= EARTH_R * np.cos(lat2d) + dlatg *= EARTH_R + + return dlong, dlatg + + +def diff_cfd(data, axis=-1, cyclic=False, ndiff=1): + """ + Calculate centered finite difference of a field along an axis with *even* spacing. + + Parameters + ---------- + data : array_like + ND array of data of which to calculate the differences + axis : integer + Axis of `data` on which differences are calculated + cyclic : bool + Flag to indicate whether `data` is cyclic on `axis` + + Returns + ------- + diff : array_like + ND array of central finite differences of `data` along `axis` + + Notes + ----- + The output of this must be divided by 2x the grid spacing. For example, say `x` is + the coordinate, and y is the `data` + + >>> x = np.array([-3, -2, -1, 0, 1, 2]) + >>> y = np.array([0.0, -0.9, -0.9, 0.0, 0.9, 0.9]) + >>> diff_y = diff_cfd(y, cyclic=True) + >>> dydx = diff_y / (x[2] - x[0]) + + Since the coordinate itself isn't cyclic, and the spacing is constant using + the difference across a point works (shown is the difference across index == 1) + + Also, note that if data is *not* cyclic, then differences on the boundaries are + forward / backward differences, thus the coordinate must be differenced in the + same fashion (e.g. using this function on the coordinate variable) + + """ + # Calculate centred differences along longitude direction + # Equivalent to: diff = data[..., 2:] - data[..., :-2] for axis == -1 + slc = NDSlicer(axis, data.ndim) + if ndiff == 1: + diff = data[slc[2:]] - data[slc[:-2]] + elif ndiff == 2: + diff = data[slc[2:]] - 2 * data[slc[1:-1]] + data[slc[:-2]] + + if cyclic: + # Cyclic boundary in "East" + if ndiff == 1: + # Equiv to diff[..., 0] = data[..., 1:2] - data[..., -1:] + d_1 = data[slc[1:2]] - data[slc[-1:]] + # Cyclic boundary in "West" + # Equiv to diff[..., -1] = data[..., 0:1] - data[..., -2:-1] + d_2 = data[slc[0:1]] - data[slc[-2:-1]] + + elif ndiff == 2: + d_1 = data[slc[1:2]] - 2 * data[slc[0:1]] + data[slc[-1:]] + d_2 = data[slc[0:1]] - 2 * data[slc[-1:]] + data[slc[-2:-1]] + + else: + if ndiff == 1: + # Otherwise edges are forward/backward differences + # Boundary in "South", (data[..., 1:2] - data[..., 0:1]) + d_1 = data[slc[1:2]] - data[slc[0:1]] + + # Boundary in "North" (data[..., -1:] - data[..., -2:-1]) + d_2 = data[slc[-1:]] - data[slc[-2:-1]] + + elif ndiff == 2: + d_1 = data[slc[2:3]] - 2 * data[slc[1:2]] + data[slc[0:1]] + d_2 = data[slc[-3:-2]] - 2 * data[slc[-2:-1]] + data[slc[-1:None]] + + diff = np.append(d_1, diff, axis=axis) + diff = np.append(diff, d_2, axis=axis) + + return diff + + +def diffz(data, vcoord, axis=None): + """ + Calculate vertical derivative for data on uneven vertical levels. + + Parameters + ---------- + data : array_like + N-D array of input data to be differentiated, where + data.shape[axis] == vcoord.shape[0] + vcoord : array_like + Vertical coordinate, 1D + axis : integer + Axis where data.shape[axis] == vcoord.shape[0] + + Returns + ------- + dxdz : array_like + N-D array of d(data)/d(vcoord), same shape as input `data` + + """ + if axis is None: + # Find matching axis between data and vcoord + axis = np.where(np.array(data.shape) == vcoord.shape[0])[0][0] + + # Create array to hold vertical derivative + dxdz = np.ones(data.shape) + + slc = NDSlicer(axis, data.ndim) + # Create an n-dimensional broadcast along matching axis, same as [None, :, None, None] + # for axis=1, ndim=4 + bcast = [np.newaxis] * data.ndim + bcast[axis] = slice(None) + + dz1 = vcoord[1:-1] - vcoord[:-2] # z[i] - z[i - 1] + dz2 = vcoord[2:] - vcoord[1:-1] # z[i + 1] - z[i] + dz1 = dz1[bcast] + dz2 = dz2[bcast] + + dxdz[slc[1:-1]] = ((dz1**2 * data[slc[2:]] + (dz2**2 - dz1**2) * data[slc[1:-1]] - + dz2**2 * data[slc[:-2]]) / (dz1 * dz2 * (dz2 + dz1))) + + # Do forward difference at 0th level [:, 1, :, :] - [:, 0, :, :] + i = 0 + dz1 = vcoord[i + 1] - vcoord[i] + dz2 = vcoord[i + 2] - vcoord[i + 1] + dxdz[slc[i]] = (-dz1**2 * data[slc[i + 2]] + (dz1 + dz2)**2 * data[slc[i + 1]] + - (dz2**2 + 2 * dz1 * dz2) * data[slc[i]]) / (dz1 * dz2 * (dz1 + dz2)) + + # Do backward difference at Nth level [:, -1, :, :] - [:, -2, :, :] + i = data.shape[axis] - 1 + dz1 = vcoord[i - 1] - vcoord[i - 2] + dz2 = vcoord[i] - vcoord[i - 1] + dxdz[slc[i]] = (((dz1**2 + 2 * dz1 * dz2) * data[slc[i]] - + (dz1 + dz2)**2 * data[slc[i - 1]] + dz2**2 * data[slc[i - 2]]) / + (dz1 * dz2 * (dz1 + dz2))) + + return dxdz + + +def diff2z(data, vcoord, axis=None): + """ + Calculate 2nd order vertical derivative for data on uneven vertical levels. + + Parameters + ---------- + data : array_like + N-D array of input data to be differentiated, where + data.shape[axis] == vcoord.shape[0] + vcoord : array_like + Vertical coordinate, 1D + axis : integer + Axis where data.shape[axis] == vcoord.shape[0] + + Returns + ------- + d2xdz2 : array_like + N-D array of d**2(data)/d(vcoord)**2, same shape as input `data` + + """ + if axis is None: + # Find matching axis between data and vcoord + axis = np.where(np.array(data.shape) == vcoord.shape[0])[0][0] + + # Create array to hold vertical derivative + d2xdz2 = np.ones(data.shape) + + slc = NDSlicer(axis, data.ndim) + # Create an n-dimensional broadcast along matching axis, + # this is the same as [None, :, None, None] for axis=1, ndim=4 + bcast = [np.newaxis] * data.ndim + bcast[axis] = slice(None) + + dz1 = vcoord[1:-1] - vcoord[:-2] # z[i] - z[i - 1] + dz2 = vcoord[2:] - vcoord[1:-1] # z[i + 1] - z[i] + dz1 = dz1[bcast] + dz2 = dz2[bcast] + + d2xdz2[slc[1:-1]] = (2 * (dz1 * data[slc[2:]] - (dz1 + dz2) * data[slc[1:-1]] + + dz2 * data[slc[:-2]]) / (dz1 * dz2 * (dz1 + dz2))) + + i = 0 + dz1 = vcoord[i + 1] - vcoord[i] + dz2 = vcoord[i + 2] - vcoord[i + 1] + d2xdz2[slc[i]] = 2 * (dz1 * data[slc[i + 2]] - (dz1 + dz2) * data[slc[i + 1]] + + dz2 * data[slc[i]]) / (dz1 * dz2 * (dz1 + dz2)) + + i = data.shape[axis] - 1 + dz1 = vcoord[i - 1] - vcoord[i - 2] + dz2 = vcoord[i] - vcoord[i - 1] + d2xdz2[slc[i]] = 2 * (dz1 * data[slc[i]] - (dz1 + dz2) * data[slc[i - 1]] + + dz2 * data[slc[i - 2]]) / (dz1 * dz2 * (dz1 + dz2)) + + return d2xdz2 + + +def dp(pres): + """ + Calculate centered finite difference on pressure field with evenly spaced levels. + + Note + ---- + Difference will be negative if input pressure is monotonic decreasing + (surface is 0th element) + + Parameters + ---------- + pres : array_like + 1D array of pressure level data + + Returns + ------- + dp : array_like + 1D array of same shape as `pres`, that is the centered finite difference for + elements [1:-1] and forward/backward difference for elements [0] and [-1] resp. + + """ + d_p = pres[2:] - pres[:-2] + d_p = np.append(pres[1] - pres[0], dp) # boundary at surface + d_p = np.append(dp, pres[-1] - pres[-2]) # boundary aloft + + return dp \ No newline at end of file From 5be137a3bd49c4fdbc73ba61bc7ac7fb6a4db777 Mon Sep 17 00:00:00 2001 From: sit23 Date: Thu, 23 Aug 2018 17:06:50 +0100 Subject: [PATCH 21/69] Adding test case for looking at local heating rate. --- .../frierson_test_case.py | 187 ++++++++++++++++++ 1 file changed, 187 insertions(+) create mode 100644 exp/frierson_dry_heating/frierson_test_case.py diff --git a/exp/frierson_dry_heating/frierson_test_case.py b/exp/frierson_dry_heating/frierson_test_case.py new file mode 100644 index 000000000..5f88ef289 --- /dev/null +++ b/exp/frierson_dry_heating/frierson_test_case.py @@ -0,0 +1,187 @@ +import os + +import numpy as np + +from isca import IscaCodeBase, DiagTable, Experiment, Namelist, GFDL_BASE + +NCORES = 1 +base_dir = os.path.dirname(os.path.realpath(__file__)) +# a CodeBase can be a directory on the computer, +# useful for iterative development +cb = IscaCodeBase.from_directory(GFDL_BASE) + +# or it can point to a specific git repo and commit id. +# This method should ensure future, independent, reproducibility of results. +# cb = DryCodeBase.from_repo(repo='https://github.com/isca/isca', commit='isca1.1') + +# compilation depends on computer specific settings. The $GFDL_ENV +# environment variable is used to determine which `$GFDL_BASE/src/extra/env` file +# is used to load the correct compilers. The env file is always loaded from +# $GFDL_BASE and not the checked out git repo. + +cb.compile() # compile the source code to working directory $GFDL_WORK/codebase + +# create an Experiment object to handle the configuration of model parameters +# and output diagnostics +exp = Experiment('frierson_test_experiment_dry_heating_mk5', codebase=cb) + +exp.inputfiles = [ os.path.join(base_dir,'input/heating_rate.nc'),os.path.join(GFDL_BASE,'input/rrtm_input_files/ozone_1990.nc')] + +#Tell model how to write diagnostics +diag = DiagTable() +diag.add_file('atmos_daily', 1, 'days', time_units='days') + +#Tell model which diagnostics to write +diag.add_field('dynamics', 'ps', time_avg=True) +diag.add_field('dynamics', 'bk') +diag.add_field('dynamics', 'pk') +diag.add_field('atmosphere', 'precipitation', time_avg=True) +diag.add_field('mixed_layer', 't_surf', time_avg=True) +diag.add_field('dynamics', 'sphum', time_avg=True) +diag.add_field('dynamics', 'ucomp', time_avg=True) +diag.add_field('dynamics', 'vcomp', time_avg=True) +diag.add_field('dynamics', 'temp', time_avg=True) +diag.add_field('dynamics', 'vor', time_avg=True) +diag.add_field('dynamics', 'div', time_avg=True) +diag.add_field('hs_forcing', 'local_heating', time_avg=True) + +exp.diag_table = diag + +#Empty the run directory ready to run +exp.clear_rundir() + +#Define values for the 'core' namelist +exp.namelist = namelist = Namelist({ + 'main_nml':{ + 'days' : 2, + 'hours' : 0, + 'minutes': 0, + 'seconds': 0, + 'dt_atmos':720, + 'current_date' : [1,1,1,0,0,0], + 'calendar' : 'thirty_day' + }, + + 'idealized_moist_phys_nml': { + 'do_damping': True, + 'turb':True, + 'mixed_layer_bc':True, + 'do_virtual' :False, + 'do_simple': True, + 'roughness_mom':3.21e-05, + 'roughness_heat':3.21e-05, + 'roughness_moist':3.21e-05, + 'two_stream_gray': True, #Use grey radiation + 'convection_scheme': 'SIMPLE_BETTS_MILLER', #Use the simple Betts Miller convection scheme from Frierson + 'do_local_heating':True, + }, + + 'vert_turb_driver_nml': { + 'do_mellor_yamada': False, # default: True + 'do_diffusivity': True, # default: False + 'do_simple': True, # default: False + 'constant_gust': 0.0, # default: 1.0 + 'use_tau': False + }, + + 'diffusivity_nml': { + 'do_entrain':False, + 'do_simple': True, + }, + + 'surface_flux_nml': { + 'use_virtual_temp': False, + 'do_simple': True, + 'old_dtaudv': True + }, + + 'atmosphere_nml': { + 'idealized_moist_model': True + }, + + #Use a large mixed-layer depth, and the Albedo of the CTRL case in Jucker & Gerber, 2017 + 'mixed_layer_nml': { + 'tconst' : 285., + 'prescribe_initial_dist':True, + 'evaporation':True, + 'depth': 2.5, #Depth of mixed layer used + 'albedo_value': 0.31, #Albedo value used + }, + + 'qe_moist_convection_nml': { + 'rhbm':0.7, + 'Tmin':160., + 'Tmax':350. + }, + + 'betts_miller_nml': { + 'rhbm': .7 , + 'do_simp': False, + 'do_shallower': True, + }, + + 'lscale_cond_nml': { + 'do_simple':True, + 'do_evap':True + }, + + 'sat_vapor_pres_nml': { + 'do_simple':True + }, + + 'damping_driver_nml': { + 'do_rayleigh': True, + 'trayfric': -0.25, # neg. value: time in *days* + 'sponge_pbottom': 5000., #Bottom of the model's sponge down to 50hPa (units are Pa) + 'do_conserve_energy': True, + }, + + 'two_stream_gray_rad_nml': { + 'rad_scheme': 'frierson', #Select radiation scheme to use, which in this case is Frierson + 'do_seasonal': False, #do_seasonal=false uses the p2 insolation profile from Frierson 2006. do_seasonal=True uses the GFDL astronomy module to calculate seasonally-varying insolation. + 'atm_abs': 0.2, # default: 0.0 + }, + + # FMS Framework configuration + 'diag_manager_nml': { + 'mix_snapshot_average_fields': False # time avg fields are labelled with time in middle of window + }, + + 'fms_nml': { + 'domains_stack_size': 600000 # default: 0 + }, + + 'fms_io_nml': { + 'threading_write': 'single', # default: multi + 'fileset_write': 'single', # default: multi + }, + + 'spectral_dynamics_nml': { + 'damping_order': 4, + 'water_correction_limit': 200.e2, + 'reference_sea_level_press':1.0e5, + 'num_levels':25, #How many model pressure levels to use + 'valid_range_t':[100.,800.], + 'initial_sphum':[2.e-6], + 'vert_coord_option':'input', #Use the vertical levels from Frierson 2006 + 'surf_res':0.5, + 'scale_heights' : 11.0, + 'exponent':7.0, + 'robert_coeff':0.03 + }, + 'vert_coordinate_nml': { + 'bk': [0.000000, 0.0117665, 0.0196679, 0.0315244, 0.0485411, 0.0719344, 0.1027829, 0.1418581, 0.1894648, 0.2453219, 0.3085103, 0.3775033, 0.4502789, 0.5244989, 0.5977253, 0.6676441, 0.7322627, 0.7900587, 0.8400683, 0.8819111, 0.9157609, 0.9422770, 0.9625127, 0.9778177, 0.9897489, 1.0000000], + 'pk': [0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000], + }, + + 'hs_forcing_nml' :{ + 'local_heating_option':'from_file', + 'local_heating_file':'ozone_1990' + }, +}) + +#Lets do a run! +if __name__=="__main__": + exp.run(1, use_restart=False, num_cores=NCORES) + for i in range(2,121): + exp.run(i, num_cores=NCORES) From 2dffa4cbd3faaa083c64344f0608de1d168db05d Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 24 Aug 2018 12:02:03 +0100 Subject: [PATCH 22/69] Finally managed to fix the interpolator within local heating. It was not being passed time, and therefore was not being fed to the correct interpolator within the interface type structure. --- src/atmos_param/hs_forcing/hs_forcing.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/atmos_param/hs_forcing/hs_forcing.F90 b/src/atmos_param/hs_forcing/hs_forcing.F90 index 226dafe3c..211dcf173 100644 --- a/src/atmos_param/hs_forcing/hs_forcing.F90 +++ b/src/atmos_param/hs_forcing/hs_forcing.F90 @@ -46,7 +46,7 @@ module hs_forcing_mod use tracer_manager_mod, only: query_method, get_number_tracers use interpolator_mod, only: interpolate_type, interpolator_init, & interpolator, interpolator_end, & - CONSTANT, INTERP_WEIGHTED_P + CONSTANT, INTERP_WEIGHTED_P, ZERO use astronomy_mod, only: diurnal_exoplanet, astronomy_init, obliq, ecc use transforms_mod, only: grid_domain, get_grid_domain @@ -279,7 +279,7 @@ subroutine hs_forcing_init ( axes, Time, lonb, latb, lat ) integer, intent(in) :: axes(4) type(time_type), intent(in) :: Time real, intent(in), dimension(:,:) :: lat - real, intent(in), optional, dimension(:,:) :: lonb, latb + real, intent(in), dimension(:,:) :: lonb, latb !----------------------------------------------------------------------- @@ -447,7 +447,7 @@ subroutine hs_forcing_init ( axes, Time, lonb, latb, lat ) endif if(trim(local_heating_option) == 'from_file') then - call interpolator_init(heating_source_interp, trim(local_heating_file)//'.nc', lonb, latb, data_out_of_bounds=(/CONSTANT/)) + call interpolator_init(heating_source_interp, trim(local_heating_file)//'.nc', lonb, latb, data_out_of_bounds=(/ZERO/)) endif if(trim(equilibrium_t_option) == 'from_file') then call interpolator_init (temp_interp, trim(equilibrium_t_file)//'.nc', lonb, latb, data_out_of_bounds=(/CONSTANT/)) @@ -742,7 +742,7 @@ subroutine local_heating ( Time, is, js, lon, lat, ps, p_full, p_half, tdt ) tdt(:,:,:)=0. if(trim(local_heating_option) == 'from_file') then - call interpolator( heating_source_interp, p_half, tdt, trim(local_heating_file)) + call interpolator( heating_source_interp, Time, p_half, tdt, trim(local_heating_file)) else if(trim(local_heating_option) == 'Isidoro') then do j=1,size(lon,2) do i=1,size(lon,1) From 9644406e0e67a87432830bc5c6c2d3b422573962 Mon Sep 17 00:00:00 2001 From: sit23 Date: Fri, 24 Aug 2018 12:10:05 +0100 Subject: [PATCH 23/69] Other changes to make local heating work. --- .../frierson_test_case.py | 4 ++-- .../input/heating_rate.nc | Bin 0 -> 15733784 bytes .../driver/solo/idealized_moist_phys.F90 | 4 +++- src/extra/python/scripts/create_timeseries.py | 11 +++++++++-- 4 files changed, 14 insertions(+), 5 deletions(-) create mode 100644 exp/frierson_dry_heating/input/heating_rate.nc diff --git a/exp/frierson_dry_heating/frierson_test_case.py b/exp/frierson_dry_heating/frierson_test_case.py index 5f88ef289..77d3887f0 100644 --- a/exp/frierson_dry_heating/frierson_test_case.py +++ b/exp/frierson_dry_heating/frierson_test_case.py @@ -23,7 +23,7 @@ # create an Experiment object to handle the configuration of model parameters # and output diagnostics -exp = Experiment('frierson_test_experiment_dry_heating_mk5', codebase=cb) +exp = Experiment('frierson_test_experiment_dry_heating_mk7', codebase=cb) exp.inputfiles = [ os.path.join(base_dir,'input/heating_rate.nc'),os.path.join(GFDL_BASE,'input/rrtm_input_files/ozone_1990.nc')] @@ -176,7 +176,7 @@ 'hs_forcing_nml' :{ 'local_heating_option':'from_file', - 'local_heating_file':'ozone_1990' + 'local_heating_file':'heating_rate' }, }) diff --git a/exp/frierson_dry_heating/input/heating_rate.nc b/exp/frierson_dry_heating/input/heating_rate.nc new file mode 100644 index 0000000000000000000000000000000000000000..fffc0d0af45c134a197498dc32e4d0e762735343 GIT binary patch literal 15733784 zcmce;by(D0+cs)TO*aeO*uw}0D<7v;1!5&4(QTCrdO8~mw9?w<;{l{_6TnYM`yL`n=pISaMXZyI<6rg#=k-rw_ z^Ks2|%A8noV))<2(u+#``LD5zbrEgK6lM> zxZ`y3ocR8)`!OfJoQeMyW8OJKr*Znf*KmB!Ecfb?#mnX{pSv<=ALkH}Z{}jZ91{Mg z+Lq1vIMiC=>$7ax3SXaE3w&1kETC#G#+X}{D_3cEjQ{I1k`rTo;{VqebKL)z82jeR zRb~m}`OoJzCqBmhU*oGabK0t9wU*6YG{dJ>GyXJg)|~%iv%|mtYm$?1p&8Q__$-<^ z&6n*Q{z{MR+i&>DaZ`GB8#iUpkijGSj^H4FivLG%%JNl9a>(Y-|C9cI{p^i&ug!^@ zpQ*B-EH=l72LO)+-)YNpOy}&2ocSD?^OKY3e}3}kDwoqNYc!7D`0to;l;-O14lk9) zk(pI$KlsPrII=^Gck-3v#u3|QE!*t($2dGD?aKSXNyg!xlgD*!`r0^biRsJi!wrl> zZ_4J?{cy)PwCI_el}i*e4vuW;-9d5HIA~BIZDN<{#sTkMt$c7g$=E;OTp@W>DPzCt z8bv5?Euo>ziSBej+099z z8L!S8yF`w;>s&M5*zw(`xM>%58QWLiz02>O!PsU*gOBl*w;5Y*yI-+ct1`xBH>c&> z8Fa)mtXj|35r*{j%}jrmR$wTUM*( zmY(O`vi?xFZ1UbMTkzw=A6e^4ZrL{7Ej#3FNsf2k=9aQXZrSaCTXz5DmOW>>Wv{1h z*{79T_EouMfA{v44fywr$t4?TbB_=1`GSvOn5ek}utIOo&^K?JmhLZaGfpmg6V8<%BYlyziD1SG(n;CX#&XmXi;; z<&^%C{Nt8W-8Z^ys@~n7R$P*^B^l+G)2m8yg(Rchaz@S*gz?YVEXf47oY_{AdnM^| zOP^knJSxfWZaK?+Ldj;GarbBEJe6p7_B2UebIUp9CAnCVkKA&u`=cbAyH4uIx#hg( zlH4iDRJWYpRgy;}`N=I843gvtw_NBxwqy$@Nc~H0xu~Qh=SuR9TQ06HNk2)N-Ev9! zjX4e3lJ1gREJ>9lA4$^v0gx^AboZAImE<}}>LeNGmcIEU*<6y7CAm|QkYuV`E-NF+ zu993J$s>}yFUe1CxxAJn2T5|ZBu_}v=9VjROR|Y1CrEOeBri$wty``vDalTf{8!{M zj{4enb!MwfI}fV!&hPmkV~2pK2Llc~Hnzp2Vd6`eu~p``10^~w zG`8^EuNdCYU~CdOF(Itkeq;UW`&x(BP#8V8Js7mH&_-j`ci-l&Jo?yJ>gKdh_oM6H z`ZA*N=me47RgyAkTqfzsq;)dMj!d#Ell;gezcMKfnG}~yic==VEtB$)NqNbnJY`be zGKmM7#EVSgNha|olX#R#yvih=WfJc)sSYx!E;6Z3GO2Ddsg5$Kt}?04GO6ylq#s!U z_kNI7kYuf-e`jRgQom6xcfVPZTeg($hc+AC{r1v5(W#cZ-$lAFWJ&IRcj?~fA+75r z-58Cc6bKS;V)21|B_O83jqB=`I<>7E%b`57+VHzOo}Bc*$1 zl;n4mbpMRr=w3fYx`)O{ag3Gjqj6F^3*76%e{W0bWcr^;+!nqSCgf9 zr%3nKR4MMMTiolWCb{J_=^mRVJGO_%bTA>D5?r2J+|_uNb=&zaJF z=Og9oKhBl-nJe9s^CX_; zN%!SEiLd!(-0SB{yv@Jw?k}h%N$DP4An~|Rx=$BMd@g+J)>|a;x=6ZT7fJjsmhRcb z63>gJ`*w-M_Y&#eT_W+mM7n>MB)RP_EiOr^4ohcCeW@NxrF(g)R2N_Ae)g5><15|M zzEYihrTf}fs@F2<-d-lvZJBg`FO%xGEWzz>xm3sHZKb|c&*giizEs!c(*3?ds_zQv zo?ju=d4+V}uaN4!Lb~@?N_Ai9zW=)_qlUg|I;nrxyNi4NC~>Pw*QBl?V}?Di+x1I< zL%tnz1(P=%7)FHXPOr~ST_2L?Dg&KuVs@;&u*MF`h;xSowDs8FIg$`DL1P5pfB5HbLQPv zbUSF2&A*bka(nL|vPD0ZT`ykyh-_)OuO9v@cguXkuR7yr%Vo=(Ec@I_xm~v6;+GA% z%RZ3#byR$hZg^SdfA7;v@585LtH&Ivd+ysa*_zJ<1`dCfOSVq4rCoWmr)AG0k}KBALgxH)eCPA^ zZ9^A$*Yz8t>KgfKPuabXdvrB@4Q{mR8dRLv#yJ5WcgFYgN|R9 zjr;O<+?C*avZ-GC^VM&)TsC9T{_?XNS7o#AO?}bj+GN?h$}t58?Vl@KIK5T+;{qwN zB?BGLsc6K!g{Q z6GcRL(al6Y5nf~$QALCoK1kaY%$WBf_C09)Y>GKD{d>cm*Wf{*n#}v0tBl>gQt2^k26Eb0F$4 zKCfaT+>^Da{hhgHi13<8%sZC(KNI2AXus-xp}k*6Xz%z5tua$*Zw3kNb(zp!ttGVZ z5LyG<7gq{x9tX`$5^i1|b_BA?Kz(}|Zv0`Z<>#fWZ1DD|5XPTF`y zn+0i`>yDMAKiWH0n{gauOrsg!490qcaZh0m-I)*N-Ah9IsH4z6OA%VjF`Y>D7y-*K)BGmmeg}Pr2q3*Ll zsC!Yj=X{~=QA?=1{ULG-wXB{{cU>mbU6P2^L?^Zd054y=1@-UWTLOtp&<5LUuq&Gr6JyNLWcnI~PR?O*z zP_JPOn`v+7KB3;vzG|>ghuehu%s`E+ul>jM2$>--QWvN@Jl;YADop#v8|&EIow!*(TchDAad37F9y1Z{`!~ zYmD>iP@%pOO-Y?=%;&}{p*E&c?}Si4%FlI)LT#-e)GxR`;S>FL6Y6v#gL%HcE!3Y2 z3H6sHLjConP=9MfxlXA6VXPnD2z55+oE;cXh)};SFVwMFLLJ>(s2{YXfAVVR!d%Y@ z^>K2p$rfrk>$8{f?$|HXfpdg<{ZXM_)m11S-4TjAXNBVCc%isr6beIWp@{q~6zAs% z#o3NRaeACk=yXDHs)|sY2p5XuErcTcvQUU_LJ?{e3hjEKI7%Ijk5H%;LZSLVOc4sD zndl%CieRBQa!)7@XA|XwLf%6t4(%l>3&lYfWd$OPd;;Z+RjHGj79tq zievxKzmsv${^{C`ZI4h~m?RVzPYMP03B`>K%ww=nJT5F0=JP`FvXf9Gb3AjWP<;I% zlzGkwW$^_%a6loQ7b<)n#3S)rV; zM<~Z;3gxJQLOJxFP!8BBl)at{W!Jl0uN2Cb+k~>QUMTC06G~5=P}Y1bl+_Z1vdSN# zBXNZN%0gM=yij@!5K3>8P&TSAlr4%0WxI8>Gg2sfPZ7#NON4SHZH#vc<+OrAIV(UY z=T{TTMQ4R_X$_(DsI;&6)pmHibhL>qTxfLqEIwo?do?Iih8V_H|^E=LJTFU6N@=UKXo&yL(UtpW(_}beIv@P z9LwUIuTXf65sF%@VRfxgRJkS;m8e@$CKMG?C}*?JJyeB0s(A}V4aQjWuTa$HzVsZ; zwb$wA8|`>8hPkx)mA1J*nw)0bZ5jU)q39GO6g~WeVqgoQ7+qf|rg;j*f@GmsRYoX+ z1_;Go#-m!my+nVJJA~reGwvJiznpy&U0x_GuZ1FpdpfQJ!ST2WLJ`ZoZX-vQvefA# z6i*uq#lv4haqkn?l1t-LetwwK4f1!LcCSaWU6hCyid(d4j1Y?3-1B#D2*m@&^XMvV ztajV76r$}H^mB=ORxT6?^H?9|_m*o?`1wlX=jYum)~z!o^GH*2e1uTE$tM)A&U4)` z){%91WTB7T{0tx9XZJAktIN+Lbq}b8VrN&O*tAC|{21RtVg~JvC%?nW301B4LRI;c zP?ZfAsuF*Ms%Vf<6>ch21*3`LLRDa=Q00#ms(gEdD$ik|{F^P5zb*>p_s&B3rI1j5 zIw+JMxobZ;$}E>qI`Rl* zN?V~!;^!;TLnvQ|QF^oAPblA<6H1=P$`pPMQ~CKwt3w!tlIOGX-3y^i*K++6p>%Db ztz1I+VFv9c{cB&9=Lz#@D^%r}Tdg@l)#$lUwHHFwd$v#w*9+BT<~K(rR7;uH>I9(* zxFS^B#|u@+XQA59ww$>uM+lW>xKL^1g(@sXsKWh(>ev~fI<6L~6O)AMq?b^gswq^u zAwm^#Q>f0&6RLBah3Y~zp^8iss*Axwb?KK-T`opk6)NPWq@By;!f=3dQ-n&-br-;O zK0V|-3C4GfzC}x+3Z?&}_Zi0)p;G$_m5SqvpF(w{lTaOEtOr=XJ-dWz7xUUm zegjyewaG%Yk~Wu6cOK)L!TcxI6{?ZUZ(v!W>OubTX?q{t!=R#AmqtKKpCN!m|2~F9FLQ|ot&{XEy8YAid zmC!U;#@H&-Cm3@}p&9%_XvR5(#;1hPEdD7ps}~AQApHlo7MgKy~c*P;0A`c58n+hS`{+9kWV(e_`cc@-}a{=UhZdzk1TAk6)R@2 z^;=`G=|T)K#qSzohAlP3taxsSS-sj2b8nI%w#80Etfr75Ax~37g3l~NLfBP9TvWW_ zMc2B97e~?zFD~pg#Fw0Ih}S^c9Bac2^8f z7Tq>HNZx3;`=*uQw(YPXN_E8$wYiKT>e~v#?Y*T8ck;J0+)>>(+{svAxZ5$)aCc{C z!`+uH4EO4V8}1F48}1$LZMf&O8Sb}kX1Kq8n&JNap@s*=yBHocU0`_7DAVv@;5WmA z-TMs>&V(5rWHmKB9M-|`P~XAus9a~mqvah9kE~q`k4qFXJnq}l@c3bv;Yq);h9{ZT z3=c-tGCYXyWq8)~tl`;3gW-7vli|7l1Hpx-+FkhHNk5Gu`w840W;yX z{UyBi-hr3054=tsfR|wzyly;zmvIxkZgqgy^<;S6EeZsjZgWBWts55FO>dcQvowa39XZv#0+2=$ZE&G=yq0ZyhsFSb^b>1(4 zcm6r>uC0T2m+SBz*A3qD@56iXTX-*c1n=pa;oXO_T4U6Ca|w0Ao}tdXyr|R20I&Cq z^Ez`3xd5-34dKJXb2<+4LJcvKGPP`Yd?JBjB-SB|H`mg2$5m@L11p z$6G(aBlrY7_9ehW=_IzWeFGjrui&xV43Bv);4ytZJcb{IM_0C+U4%#dGw^6!79Opq zz@y_&c=VvX0Y2~;(Gwn1PttZlc&t88+Z;dq3m&>d@VJ}|k30QPvjW#AWx(T~Dex?i z2+t}H;Mrg{Jlhq8XaBMAoaln*oLcZ)R34s7lHj@F5Ipy9hv%`m@I0Rh&r9C$ylR2x zgP!n=^Mq&WcX)m{1iv$J@VnRtepe>L?`k3VVJiF#4$4XJyY>@)H_O8B)(iOErhM=K zekL3I;=|yV5DGv0U--RT#(op{z4$^rW4k)#0?r%Zm(mn|=^KfOZ0CXBYYXR9@Qck0 zKU*n+{$drh6Ar(mTkuOZQ+D85JHhccw4K%rey;xT%VG>4web7V3I2Izz`y7Y_?K$| z|LTq5U++5nTU3L8TaI^V2>(t8;NNQ-{6|cN|G0YaA9J4VBJdwt1O7w82!FOO!he*2 z|D>hxpL+=Y%XRQyJ&t$>|CJ%|U)&D(WM3`1ebM zfB#YNAKZ`qTI@UF-<>ga2!emhNccCPe$^iEFGyZK7H0nB#xe|kPgcP19(gh5g5UL8 zL{H)+`6mws&R-bJx}1WaZZ!N(c7xxE+VBgP6MFa^+d*(lum+(q#3T3}F9E-ZNcf#& zea?Av9r-@Kl=|c}+(~_Zf_{#@qpfuW{h!F7-gpF5@Iyewxd^DV1p$@!BcO^I0o8^h zphhkP)MUH10|A~d5Kw0b0=%appzcNl)EiAW5m3K0v5t687&!Kuct|WG1{3S4pT>3l z5#UvrwrH;=^{Uf;)w1-rh<>}!{|O>TmS;R)i08y#_D$4bj8$V8PdL|pM?kGxv@;0- z9<3QybJ{&Z`^V|u8v%`0A)p!KYds$U9s3}l8~yfahJZn35HS2S0>->Wz@$A0m{AM? zbD|Njh_Nl-f`B#D3pj*;9n}!9mpHTl0ZPWBxq^T&&K)O5x{nArOB)xc7kQX?i2yxq zT;$j#+P_G<`ca(UhJbU`5fJeK0Vm!eKomfLnm*)e+OC0s?RNUQkASrtUpW#1OUcXp zDG2achkz*$5HOzkj9QI=q0D7K3k394q%%+ROmr?ATX>40>wlGhJQog@fZZ2 zJdQx!CIp_|i@>vC2t0oUfsr2&crgWmhT#aj+=4J65PgWP#44gc$Ml@vLA`7Q>KAZL zO$1)3#PxY;hkBQqBe!;uJT-U~r(&LOB%IRtg< zhoIgK5!8P%f(92s(9rh?8qotmqY4wh5i~jvWh0L5L(s7D2pW_~-7*O3W#Zaev{8w+ zGze-Ih@kq95#+H6K~>TaRJJ~XimpLW-Y*FJr9WCzI8BY<^c4tpUO;eW2?S^RAo#;31b-TX;D34{`0G{#e?NxcpN)t| zL=5qe_=4bHvxvP!J_P?bkKk`R5&WeP*RuT~1;Ov{(%wx3X9XiT!$BVz2zE@Szu5?W zcapyUAtn-j>^I?feu8m1dJ#l@U@P z1tCSYAfyCk+1&`KR0bi{C~HR{r0!OPG>$+>%S#Ape;y%S?;@o4QG^WAB4p$mgiP>9 z$n>=cnL7_5OBx|$#WRGgrkxF$2nlM4kR1aMvb!uo_N_+9!668dI}mat4I#>W2vPiC zJDg)<5OM$r*?SftAs$>)ookOFWK&UutZj^tl{aX64t=ym$g~Xz8JB>N;dKzwpD}l5 z9vz+{q{T0UG(3b5uT+H8V7@7x5u89@#U>%7AY;$l9l?KFFn2+oS)+f*(Z{0XHj&)s zC12z>{T|^-E-SLOz6eeRg5T~T))R5;S0hFcTWG%)h|NB_gou78+&Ps zWBrQL_Gid@zN4S&^v5-_L6CPHOaF1iC$_;cj(1xOdH0vpoe6obUy%2)(neRv2Th{g zAliRO|MWLL0P@KNA)nR+@|lcv_C3hw4}g5pAIN>V-me+tYr8?d=^x0qo`*bG4f*~< zkRQ1N`B6LM$M!=WF$D4pM#!%)|ENlkKU@g;bLM6Vg*+}d;x(MUn2zm9bkXO76`BeIvaGHG2g1kT~d3BLj)?h$> z@>-buQ1)b;-IK{(WpY@IT#~b{T-UWcqL04_-P-`6d&VMkH}ykKBQ%&}!JpZF#qnL#DS^=a zhpEr@;f7q_oOTjuZz64vq|f{ajmU@4vmpq*&<>&cg9tS+o-2BU-du~&JM{fvIzpdt zESj;|Dj_sJH$oG3BlOKTgubhe(9BK<{d5zd-}@u<_aB7iJ%g~Kl@M084#KKzLs-po z2=kbTusS0UR=*;`8dXJDlTd^;e~Yk|$p~w8mpDy;=s|E@%M_vz!kT6ytO3{8IfSs< zqY+l^GQuhjM_6fZgcWH=Uwi4Z4?@2%r}udgn$A4lkiYnc2sJZ@2bU0fgE?O$x{`|y zeDE?t4-F(LvmMB>^W@E+{4v)(afFs2w|nN2+qdMHH4g4e zbVlf|CJ5c}lDN)xVdyqCgf74iU62{NZ5HTuRe^4IZRqwcgl>N^=nia#?odzY4tqm) zgkwhz6CK!h64XE3m-FwTJ6ICBeIKCPa{#)K`q1r+hHm>N=(cigP%G#*%jxG4v6EOu zyyjRz%5Rh{pbM-8-If|$GY`7p`Sib(zQbu-Lm!i&Q%-^I=yB-68$fr0YjwMzJN*Q@ zbDf}z+z;L5;n4BBn$8#p-CY%Q58Ffcv?z4Z_o1^et{6RZ@jB>+HHI#wG;~hJl)V_b z&&=t^LPX@QgNT9^5mB@UB8q)RM2S*}DA^eirOqLu+*CwV?t_S`w~70RsPdiiDI%)r z5K+AuWf=R^=Nmer>?K4LpF&%p&3n-OJPF-Dx9E#;zk3B;(h0_~mobfp?r{`!cM33W z#&n7Ko!JJRNQO>P3A+3Zaq4eAa?+X{k;@G|pj$T^y7e~p`;e~`&gCYzsl)>I{}8mX znf2N1MV(JXBj`41DOnG`Z*>99SZ~I(p$_d-BX&XGbu;wcWqOLqD(t^h0h! zKP(OU;b#cS5lf*T)gAgVL!lqb`7sehN$5wrh|`>tQ)e~wL!HnMo<$qCpdX-szJF`z z`;~^ij|uwTH=ys;82TQEh+ISzeZJ*bRm!iFA=J}xT_o*Lchn5}S_J*LkI+vrLqFv% z^fRA9Ker$B3%sFUZ!(Ts^lhw8%!U5;b?EOCPdN81 zoR~oTf8)!JhyV5>mOu5e|HFN*r7*R=r3%5UiXx~7DBJhNBZUt2RR7L z3;jCUTA9E+lAxcrANmxOZseXxhzza^eJd@&TwCvjz8#?NJdHJDj6G@* z&4{tYV(7a+Ac_%%InS77HiF!DFV1<^u17Dn^_*k;vNNo`FO15MFlzKLh9NHQGBz+Y9MO2IB|0+kaXM z#&iP4IOb$|$GF=um%Pkt1B}xKtz-3*6uizkdbsk`?# zjE9rSUuW`1eQg5yWBy^Zsg2?I3!)Nj6lHD6mvSqNs$cAHCWaBLp_*~13UIC+%mMvj z4!Q*M_Qx>qdJl83fpQ1?N|?8F=3FyM3(PyV!5nf2=DkLk_qT!hz-pKeorU>u1(@ZJ zU_Mxeb91=ClzoJIlhH(FbPNMYh|+oPG#qexEU?G9P|-GrwvAb4(h2 zU4i*AIlg-V=9|P7wl6Y|^VgUIb34iS!^zdrJj@PtF^#HiKc_$r@~BT!D?q@onlc{*Bj(XqY#o5$_2~I>*qa4A$(QuzuJB>qiZ&pSZR#WBaIr^`9%Se!U0lHzlm!Q(^tdv0q(b zE%Ast)Xkp+Yd*&O`v7c(l3^<~1h#U`V5^u3Th)WGRZoYl#&tr?_Ep%b9i>#jR;wFq zwcAtfAUN;23bwj0VB>qv)?g`Y^{7{`3ftkZHO&oMi%8g7JtPEdt;fLDwgl0E?LgRC za;#}8b+xe7DNkQL=--R6e1ol60miiwwtPQf%f*=g7>W7JBL>#bQ(#T|0;_2Rtoq%s zM09}Vqy?5Me_@Hf1<)nT=%g{V5N!T2u&uZZ+cI*qj9e}!*DJOX_ld>q2XnkS zY`$PSmim?Gb0*ggBmTm+E)KR0?`gLzY(Yg}+d6=tzpbp_W;3h-S7BXqlc);oI{My_ zm;PVDx{G`T&xUpTQCK(rh1LH8tSiUEx}hO#`@>*6`UkdSzhKkVg6+&k*v{HvJ9i$o z^SvlF?E6u2Z$wsv?LsoggD9E9#hS2Pro6lWwoCiiXT1!}{aOlaQA1(7Ct!QXz42HN z+mkHVp4WiQRFW9awl{2#nE#y-u-zc1mq*bq=TEhOO?-q+%Q>|iHsxx{BkWhEWPJN4 zGgrnv^D1l|l3}gzi%5gDB|kGG$HBUAB&;hwz`ATbYzxWjstT}eDh=Bf#f^qc2+fx9q!^o6C@ELf)0;5s`jOCn)e7|ONWtJ9CeGUFU9GhecAqfNjv zzYx)j2!mzbO=26{*<3S`cujnPWzlll@rPxZ6_yoTw^B~|hy7n1r|$d*Fpt;(Gv71j z_HvjzeulY637GrwEFIbm<}vkPp41iQ@fMhe6ot8GMVNc{;P^sV*7t#BM-nW1Q(!q1 z2#X>YEb2zEXs5yw76S{ve_KLR*>4JqW*#gm#-$iWoMrnnb$EgGK&?ZSu6WH1n?n%djD^RQmu^67K=~7nu0FHnUe#=G+qIVUu9)wG-w(Mp&lZ zfW>z(ENk?zY}f?LI_Bb+7nY^uZ3$ysJeOF_b~5#tyKi%1EHQ&RIdhX?de9f9JIi6Z z=?l}P6qrN?Ogld#dP*ds`?W!I+xv)az8}$@jv;ztahMMHz;t0HOhye%_ou>i_c%A=?It}v;AlVWowRw!epi`i=bXK z`+CmxqwOfhrGN>|5gk+((Y{j4Ev1Jk#`1n<)C^~h;yVln&l@%bioZorhuSTdMT)-iYg zLpJ9a$9Jyd?+VO$^TS-I8q7slkK$EeF6{+#xi>IZ=nZrETrij6y-?lArNMStU9`g)hK)Fs-`cVXHlpZ{E?-(sxyH@;m$U^+1druo$nUAYvZ z3)M%oC-Wci0@1t4g|arHw-!Wn;BD3+iCmGpxS=p5(9Y}fFugd>x6l*v$sA2Z0EoS}syNu{+B@tcf8luYw@g>m1p-u=@M-DNYP+vP`e%f*Op zww-(PG@|=5w$ZHlL~r6d_sDWYZ+Xpq#{HZ~pM0B{-jmZ$i(vZb%^DAf$<>5>b8lqq zfXNX6Q)(hi$&F!p+Y6?YwlJlW(+t)qyEjZ9KC@PhVEVC#HS}XGceAFEJP%pp!hBN| zKg@aq=2CpGmgHNqIM2=^b73yXcWgfHi(KUT_XC18`c;W`_!;}jJ^7QqzBOb_+>g!{ ztY6M&g>^DA|ARkZ8vg{(zuv&}7rpV^5sc{i8Ho1#!o58X(aW+obMxoE{U_uKOv4CB}p@*l=J4TR-uV^}hIHoxK78+`|s zi$h?Zz&Cvdo~!jjiL)@bU@pD+)}37o<~@F}yqpb7T1i+^3c;+a$#-l|Sdz7{WN)Fp zr?Bu%X8FAomcM*=iAwe11tBpHJ?8$f8t^Jb`qA) z8;A<9eBv3M&9gFzXXZn`ZRehXc?$U)^%mwq88G)h2lH^g#i#qi>^mFgjZOI7;Va+O zc35+jhP5o;R5f11>X{E#j~}qs<~z6Mb6Be~rpn!6twO!(ONl^^alYz#w)r0Os0l0o zW@&ApB<{l6@Gz`R`Ce*GJMFi?+POX7M5kbt4S}^gZT04Rqwfad2&{b`!`kyZ=XghV z;e1EH+U^#?H*Nb6Sle#mU7iSQbH0}wI$^CvZYt^->u~bGyXt2hSklSw&AqVfN`u8G zf@cx;&R^ce{5^{K=w+C%-hx@{1M}{1m=$BmJ@;M|-$4cWhAvZ=;QPDO0&>ha3%BL? zK-%t2koQ8XVJ%P@);!%#3pS(?ArjgA0zhl zd&F8#GPWg%ja-aa%?iY{VhsF!YfP&s#B^37X5b6NOeuz#RpStIWGekmr=OyT%Qpsb z1s#MJ;__+`m+Jvy|1L-DA0L9_d0dDql!dtB9S~RI8{$gOMjX%ZxJnU-tJ()~RRR%L z@jd&c5m#**;ym6Uu5KLS8n#7Tqwz!`#5EpVa3v0ug`_8;=YzuUyd*bsYd2V(beZrTahELHf1XwDi2vgY`=~dt4{XJDTiU4xyN?g-ODDtbw;y(YJM1fe!M>c5?+*Lo%Y=dRO<`Yj z4E8y9V4vBSYv#f}w+8!DVV{vqROA@f`7DNg(s9^F_Jn-^+dT%t-hpUjMSRpC#6Mhz z_}Jly&la$cEeQLp3V7K#FYN2HVGkM)`&PyiG#2(v`(WQtn|Ka;z){$@T!DQ@3hW_l z@3X>wXgTbM55TUh0K29=>__XtuC4>SqBQJs+B)LFK4TAyhW$ik*iU|d{S+HP{XNeAZ<~T_j#ChQwpfk$B_+ z61S+3IIkQMVm~0^qzegKw<2NlY9wsgiG(#1i55sWm4bv{XYjgLJtP`hAn`&Z67~6z zXpBPQ^Q}na_sc{pE0%nTw>{wEH)Wn2w}ejo2?m+$Ic2 z{PP5fzm5^q`SX}#-o!`ZCD)!rQeGWVm+gH>;@_T=@=c`P4vnYb#!N zn1a_cWAS=UbLKo3uSYGx>ydfzdXSEFxPsRs0`Yoj9VEWlKwlM+_)keBeobJ$YZ+5D zjuk~>az`X4)u(Jpz4l0q%_8s1k(PT3Qa{y3s$&UK?Qf9!pb%0T{y;K+7m<9|jO6>v zkQ_4=$+w>(IieSmAK8&oeLhka7e(r~;Yj`26=`|;A+6+Yq*a`cwDO0LRz!)^9}SS2 zQ2?pQuaWxpBvP|3BCWs{&fiB`sZmJ%`x2=iXCXE7AyVJxLF(sANd2WoTK>04E2u?U zF3?7G+OLPSyf=|nLQne_kXCIe(rVX6T7wTrYkMAPT@sPj;wsXrY)4wLc%)UxMv8@b z^`ibjKcw+HQ(E5>NE>nhY2%DYo5poB80*B&NE_Y*X#>KkvznmpfB{@%r%WLRA#IR1 z(gt59)^YqO$0|{eKKhJ7TK6|d>v|k%-E{QT9BDoCB8~so^>ByTZ|Y zJ{*OLA>Epc^dv3&o8aWz&gEGXuBKYJx^#eRKpnV-^n+_~AGmsk!qt2pT-CE_?+9E4 zn!!~n9a;R zI9HW|b3!_tBWP##OU~VgbKe&@4>F(fZQ=9}fphXII0p@ab6k5k1LNpJ2dBvfXYzD7 z9UI^@cZb7s6w+Pq7*8zHGq1v_Y6a)h2sl&9!1-n(_3tsBXgKeY-}oJHW<7%Q^AI@S zMZy_70nW#VImfu~o`LiJNjM)YqNLAzci?TsaAvO~SBbQ{7tRy?i6A&H(g*)tvGZ(Q=Kq5GZ7iHS*D`Ld!x1nvI&f}Uex7y5cs7>dS;xJ0qZu6jSKye~3l3jTID%Z1)ws`g^Df{y z)QMv=JHoM#FH~pci4eb#2~l6Iv;hKz;$U3T-Un6b(694cb%?l@8CL5 zh_+nEb5`p`JH8xS!~CP+Qk;fsKl`c(#=;mcGv?C`sC%8ZPmmkdV$TN38*m+FOh>pr zi0Av|LU3sKxxa7#4tok5S*PJhtW1h9H zfV%W+gzMroxXk1`ekEKn!wLF&Sr9G@Ik@P<`>Ygew;CBy1(1<=85xdo$aqwOX9LfW zaDJwhJfAP{%(~6<@=|GjzRL2m%5!56Kb!JBJo96DX7b#~^h1UijEsni$S}M|hL!T& zQ)Hx!MaGL-$gpJ~!(I*YA@J7am+sMovjEql1kP(xNj4NtnME*hs?}N-ddSrSvBr=d$sy%JS zA>(H;+FgsxlIM_FWI8f_U*p>2$jI)2jI0f8cSgqF&B!demH3FvT#VyONQRlWfBVd=-58lo zY9h16L1eZagUmKYWO^iW?H<;8H2I6=|A&%}$Q&Dv%%KWo)+tIm`H^v=1~S5{BjbEq zWay_N<8%Nrj>U649T}&{$-We1tm%M^T~^j<8SV8$*7ftqGEYHP8fBIpSsC+@^+|^8 zJa%LkHd4Msma`eM?mt9k)%(b@)>2iXnoAUl@NsRqnS|^bw>g%I?B-jM-Q*y$Yj;F;gZKZ0X6-(78a&KKm}o~}oFMikQdo=msQ zM*7|6Naye4(<7sie&PkvLyjPQ#Z07c_=EJb*+~CZ4~`BR-ZR{*Rrp4(`4X-MUHHzO z3D>BbaP`~9H+@mg^9{pu-}#{mziSlZcLoz&b4T*6`~|MI72zt&ycf!Zn=d=#qN8Upn&L-3zXPZQ$zWpuc&1YxBOD z&oglGVtz;9oiVpCTywg@H9L{*mK;0AwfTAfH0GV;!S#%H@*B9OP;Z)+eZJ>s@?4&A zpKpBX&pXFAa|+-7pLo|XhNZP>!;|;aWrB7VOr)&^9H*Umx9BUD_DAzhJO$TU-e((_ zdmXb`|K-o>jTGX-t&PU@qe)N)p2bs@7J|k z_gi=OMsM13>n_liy1Tnu>h2T?AqkKah8mI@-JG_R3Y5B1g#;2paq6z`x$pP)$N9|8 z%rlS8?9R;2>}>AgrN}+SXa6Pk)pOP%fAJsocO@t*hPs(KCnCK!{WI|sW$vUdIkbax z6W%hcOM57L=BEGIKDj1x{o;^2js5B%_O)Z%BU8H!nal4Zb7&dP>;FP#iDt+w*&3N; z8pHk80QXIQxUIk7KH?2`Fs~1egnMs)xP4r32Oh(#B~|b$g!5K6=hJncaZYYVR&UOI zM{u4yX9lw7-{N%#&Oeioxhfo)Ggonr9*eAhFL55e6EK2e3` zmXb$pWNnorYa8dwyWC9g$B^$Qw(~u})@=meQ*eHmwatL6-2v1^MLv9gu(JKliUemR~RGewHEU@7pM2CFz&ZH*3P_hv)R!0H(d7-;N_|K^fLTU;1VjeaLza zorcWy?T|V4C^847uI!a8D`^ z_r#xY#ob`qYPc7%{cRe@cBqGYyMTM?F}SC_AU43gIGlOy;NDde?wuC6x1|tv#@oQX zzZ~3tRp1VK&9+>N;R(3?tFtY4VY@y+^*Dis38Q+~PYD>LF!M$cP^LZa_ngjQ$!*I`_9uw)8 zxt-v;rhqG{B3yTDa6Kyy_r65xNS^!2?*Q#Q5W(vZrU%j10D^WM7)~7r!_D^N-rS1( z7sI_`8)ed;i#x+T?<3r^7*6qodu%Vbhdb$C>ePKLT=T!d)%O8hZSKO=;t5=;32>x+ zgyZ^szVCetN7xL$-|Y{FbtxR-@8LL55{}goa4ar>L%WXeftkN=1RRHqNS|~Kjw#`A z{JR&9?#tlV7!AiAmix7n?_>WUJ?kyfKOE${Q7<^ljzQXh0;DxskF*jNIDYdz^TnBb z2U?%+MES0_pc2xWtVLSs4@j#;+MaBKRr~Quf31>gP?+v%YIrJ`^W7#&R z9f6a+bdF!kwAOGAbig@yFPx)O;hfBOyi+zY>_|LFUcM+V61f0t(a4w$= z=aMIIE@atx58#}?lGjU_pTY1MoO7taU zkgKqJa0Li z147{JF`6=8!Kvc?*g6=_#?+}UZLC&>@2V4#)~pkK-w=-M8gOK?-m|lis+x?{zxyKP z<3*(8okt4y-KBi&&CmlWIjKm=T8Wf`IY{~12r1kroAP5TQa(o_#dR4e*ES+0au`ys zbB&?gIHdODTFIzWa5T(1~>wK@jduKI8L{P!*&RcOGEzme*DI0IBqQD z^#M4pwj(GrF%FKKmEpMM!EiXk(rl9_;kf$-j(Zj0ctC#172$Z&gCX^NY=q+>%Rk`z z_6K~|{)n=kaNXl+M_#AEar-*?_lD#04a$v%+Pg>a>@-ZTbYq+D*KFY2)~NNZl}izJ1U2v_!6_4MA#;sz}i^LrVSMc=@L%*V3-y z<=Y~7`Nf2kBE7iYR1_%{*CM43*QA<+GxSADnHNaqUbmNzPT}SG>v(zV5ZC`cAjLz8 zl%a2s(&r&ke{zlN*HYL^orJx%9x2mrB4v>cDI2qpvb{L3TOeh|bfj!x`ras{Olpgi zzkcB5?Xq}z{~J;&97M`gCsLxxH=_FZBfc?%C*tyTie)Ssc+{TWy13KSP>mAJ~u3FP2`g8_vS+8%RBpVc%E@ z_9aDOpAd(X^_7q^`6Tt|gp{t8kJC3NL|{A z>xKQeR>}0$)sV_@acXZPQtPsAzW?Ggf$ga&*JlU#z|m$q99*BV|H?;di91MbAV+HF zv0VQpeJpuAyu|P&@BglRc63Kd;C$9UpE=uhAZ1=3r1WCF_b5*t$0H^GD%U6Z%y{|` z_NsQ++aHGAV-_5rxCT~q64J^JL|Ub6q*eKaw5o5BR_!U$YE)x55^0rhAgw|yX|jm1 zNGtyYX%%xwUxwE#SEd2dN-_Lv57LTt;94f11z&2z@tODJ7p`r7{|AmAbJ>s8BY5v| zU#{cJPB`9ggro2S9J!C+c(oc1&OaTgayXt}d9w=U>`b>_o1H8fci*1 zL;6Kr58k?w>%5`72Nxq1Re4|Y{!e51EIwOutHAM^b>z_!j)KQ<3gdjy3kv)^-^E%kjeZ()11eT9a#%<+;WxQTKzid1)Pv8T#_sHlEMU|NMP{Yn2&D8`}sD z2X(E#683`)VIReN+WREzBl(;dIt=!Sw_x`l0sB|h#ZJ~$efE(i)~mZa(tdNDxYT>3 z|NRxtWCNUO>`QZoz?uJ;;F@g#`})_V7_yJf-UDYQLw;B1Ow+-cJ_pW>T0~XGU2vvu zf%8QeoKM)--phyc#w9otio+Q{2+oVm;EbC{yoU2~IXJKMAlRQ?nF{A6zDa&{THOE{kVSn z25Cp>+l=3E#%kd_%W=!O8RX6N;UagaLrXX%j$vY-GyN8vHtHq~B^~?wvmAe%<2c~l zF=7L0Xy+;Nh~{`8DhAFI)H{NGzj+FE;`k$w<8&X6dA1DaSfUfhNd4h#dkkrvJ950W z9BFm-!!C+alZ|m2I^LpYtu_FBE92nr2oq? z;@DqEU(U7YZKaVO%>4j&=+j05&ZB?B)p-qE@eh!}XMTptf{d;TWb~iHbku;QbmP2rD>559klC~+GMlvFb$Mhq{KRyLXn@RmoYU4i%`zR3S()pp<@KE7 z(yl)-aDOv#t$P^Uugk)n-3;zm9Pj6FT%F5!&LOz-IsPuZN_>U8U_ab>+eycW+I3_wgjoiAlSuB-~Rs!?l`z8s`sJA01qt(e&E^xWdQ6m45~qHD4p6 z!8BylxrdCpUy#vMgN(Lg;eO9C`{%lFe<1I7R~W`I?F-y*#!$zO43~3WUXLO5d`%nj zISWVW7w%r^8_3s3rQKK)MlugbuEdmUVE0oRKra6L-roU8#{ zfluLDRu`@*m*JXpob6&1TVTZk+GHWO*a|d zLdLGev~wmhgg-KuQLoiCs3Y6-D)L{!cPGob(EdvV^;jH$OwLg<7gQkl&Sp_Pg8Kh= zJu{Vh&Y>;A$Q-vCnWJg%Fzz4drDT1+;5|?p?jMzrQDHS)za}9gPKk{B1CWtHAHQe( zM;2TiyCCCxcVv9ZfvZ_(WTbh*`JHw4WICMIX?)fN!nt-Y=OLW;c>m$^ZWS`d4o1d0 z(%2Tjc_9)`uCF`0RpDNjLr6cseJG>%A$`y)q>qe6`qb4(pU6EneYxl8U+x9!u@>pW z{gJ+q&j!OVq`%q%XBR!6BY(kpvkY9tjBvGU16TJ?e0RY1F`y6Je|LhLYh>`evozZkq$QsvO)MYQwD=0JnNRQHSpuI44&;hPw^#;by1dZhR8%hMZ&9<-J^+ z@mhAcYc_zpdM?9$aC3jCyKYrp52vh3aC2{^yHXz9RX#J$d!eeF@#)lu^jIo|Z_YEz(?6Zz{Pi5JHhdQR!)JqH0bJ}GTnfGe z=sFv&q12_xG`N*Az7yfRzRO*>y9L4h?;6^3leSSe`5Cy|uYi!{S5nfP42IWHlk4 z5*3M#%rAzAP4e*2n|s8UY52VX_su7C<39Wcc%02|+vMZ%sLN6ABOXS5w{VYd5clj- zw`ttl+KT&Z+i)*!JMMLr&Eme@>I}~i0?7&sl2!bsM!Avc-b~}(*!D&wx4w(y7Kf4C zjNe8zWw~bL-HdyITXMf{yEUZYKH#>$$(yuIxsSPF4$+^WoCg2x9e&6?xBPzN5BE?z zesM4PMed_7hlkVYn~rnvu%!_XW!zitbrTPF(|3=p+$+w#)m`~*pf|rQ9CD0%nfcw} zAnv6eewKD{A2Yw7O71e9dx;+)c?iGb8mS=WA=#Vzc)gw=c_8;0_qP$$WdQg7_W!|s z$lOcZKNHD=GMLZrhWeLb`D7&b=J!#(K5@^ooVN3uuX^8DM+;e33-G9RJ0xD+jKsTj zkfeKnq=4N>rd35E$7_k{uaKDWHxi#^A<@zqiM=}^q3{D8p)MT}+|c}8_1-8IsIVR7;ZWx6Nn%%o8=yn!U{B}xh^ zp#9~MwCDyBx6eT0AT<)J8jz6R4+*Kuk>D(jgx`TkR9cX@s5lZsULr9r8HtY~kobf; z1YSepQwfRoHuN8Pyq|)^pLLORwKtOPUqsULB1n2kJzxBTq$kIb^zbQ??#*G|03;>X zLekTAOrMRU$1jM#iFkr_`lvCI9&|_2-54a@DS@N}*6W#}NKEx+y&jKkA4jZ4T*3py-F=HVdke%B0`Xw_??l93G2&7Y55#?SA?`zW#C_b3 z_&62fe(puwYj4E8n1{G~i&!=XahHoC&ZI;f_m#$JVi8y3Brc{1T)d>g#i(Jp7#4|( zrgU5kc))NlE?Otx;sfAf;U&Zsy@t4|ml4-?JL0^5Aa2!i#Q9Q(qlJh&R~m8WPg2jn zsrLfJ`;|s~&;i7UR6zXUc*L{bb#o!8De)vV08bqlX~=D8tZki1#JkhBb)O z-$vYQ1L8V$MqKryh$};%R9}TSk0ywlR|Rna%McesKb-GP|1$q}7~&qcr;qC*E@L?2 zvL@1pI}rCd2yuU`xK#BbE;S#Gc=s{H=X^wbZgs@J5{SK)n7;K-j#^mPP*;a5W6!RvFj=zcEoSQR_D1053D%rAB?kHpE&Ey^BBB*adv!N zoL!rVvzA&sd!Yl)rtHGmH%7#k8i3fQQxV&(6JkfTLF|HIh~@sg*p11E-AtW#U53Qx zgEX-_q**@PM?hND9@6TXkk&PVv@sOYwyu!4J|pen^}ec*c2$71-$I7FVi4=ohZv7o|DA~CyeRf4)6SMctVDmuHbd-5 z`Y)R0V%&(m^bN7kDA&9iQuJ}0EB=db0>g;ckZjc9%vng*Z;(zpAaP%`baFkUm^qN5 z9zr@1Pqc*;xr1q&m`~eJ9wML0kdC#06sAU0zs`t$Z9-Jr8i;DV8c_{?FbqLdlTe1< zh-#FBsHOpkYFPqNEsiifjOk?&)y#~jcC8Vm9ET{`Xhbzn;aMoE{ud_BYR5;W;P;0}<)wIVD%h@T{5pi1hl2$htjw&dFAs;Mrl3C1r@L(G!u4 zJ0P;{Vnp(LnaI9-5jmRYoXo6`$W^0xKFoJS9_BeO(L8(R>>8e(!ZTfN-$Z1}FUpf6 z`qV8%OM?-O6hxn2kEkLY5PjhVqA!(2^tFA6z7>h+M|}|e+=b|u8yGf3bn?H5e*6Q` zPxmAGS!;&t5&gu-^Z=$85G?n&3Ck@(^aBl|?|nt|ZPHz)tmtQm_V+^c5}t{(gJ-i$ zh(cu7SVVGfKxCilh@7*W=jbr);v?#K3z7L35&8Zg&z)IIA3Z=+@r^udhG))H|ID*v zj`8f728jI4GhJT4LFB7IhV2nmVj!X#7DZGie?)&-hUnim#1y@Sm|`Un{ih6~zeOPW zlNHgQYa#mU-$Z@JZzKABSwz2Ef#{FjNmCQizqYe%8bjL6Z?vKYctfPkfOx5aNS+7r z@HxceTD*RS@J3~L7E%(!t~|@A1Hv1)5MI4I!b|duq3?EteSL#4?w<)S@d@EYc%IS6 zDhSKEhp>z$2z!;qbCwPw?A>pKz3$2Q352CzMi}=Fho#>_SaKo4F0A1BPCSojV=lt@ z{YThjp6xU)6=9<`B5V+8y7F8k+81HDU+Ci{JWI+B@nknd%1Q|K;r~B;ev9=*`ybgDwt#rX zAF{lx2JxIWJWgP`n`P3N))0c@3z3@7bFGH3>^+3HO@K+*VKTRbDaZzsZYtpm!>f)k z>6xZ8!jM@MhU`f&WTnE8-T;OtHDS2B8HPK5V0hdVhUe2@pj^ZCmoV7M!Eo|54CiLT zkXVW7OJPXu4#UfJyzULdgVQkFdke!26AVzp5K#jL|7S4lDu7{AIt*)W!my|-3=?x< z7{GkhRv79@FjNSGq09srD$jwTVPmGZgrScfhS6nUm_HhZjoo4DIuE9PqhK0Z8KyA` zm?rOqY34lYM7=JYg=uL+m{xy-X;WL6cGiJupD#>4;V}8qPTv|Z`MriI;0tYT2b2G0 zrkR+(5C;2sn060?Y0KX*tr`x~T>5F;BAEJDfJw0%rkb?5_X8NX&&tsJGYswBjPHVB z>fbPE7sGHM76$!e7((m75WXD-ju{L%4nxvY`tKj=Az<21+52b0w08h?qwVaY4YxYO zaJ3Z-mq~k5(9egdBW3DO!xR`o{Q_Y+aE4`!Fzso^a2xfFr~YP`4Ajql`qAV!eFn4CCb{FrGUHV^kjmeL9Aq4?KhKV|xUB842Tw*)T?4VO}?e z=MePu0BK4f=q=C3dvhK^?|J^;XP#U4X$XP}&mt(}8iG=85_!Bn$g}ZC|B+_^epM4A z$P>g@qCbM(jYLo;9+5#&jjo~F5EWaU35Hz1>0`8|Dj`JM83(UJ7ji4K}BZ=n) z-dRL@ZM2*Ay%|jZwuUi=VayZaG_j7p9!39*r(fyg&vjr7{{>@2OBjzjU_9~|M)M6A z{dm2zE{vPaw9f_Ork*hF-T|Y|i?+3)ZKY^iJKFXT#<)c=UM&G*!fhC@Z-enh1dLoC zHQs6tV^U3+y-LE|{Ugk(CJ0%ZkC5fF5VBN;kQL=&Ze9~18(tt}=Ou*flMu516he6B zOvpAfLN+Oxmq)mWZV1utK**ZuELXttI)p3+LY5RGPB7j8A~4(E}lxtqAcTuWtPj(vv)T&Ldc6Km>M{B}7vtW|$b&McjHGUqvAG9AeEx;G zXE4nD=fdoH2<9Q3iEX?#Atac#`FPQ;rwHLbr4Ucrr%p#mefp$SDuO@Wq^ubTeyK+A z^TpI51Hs>mBBbVE(j_8fAnhIHrp`yHb2dUo9%uL-A(M_EWKJLY;|cxt5B(krGrwau zcRmSoS1rS9jMsr#O}}i9N64-Yq+xnH511R2gSq5!gyuX(D9=*}DYu-ql4tSl2)#8C zp|14^{qYs%8Zy#9gSpE-+IE{Z6~N51+{~kDz&z#_%;RRmJk|p9$i4{P8i3%{pAkIo zC4wg=AlUmkjKP~=G%jNpi(oIdnd8f0Je>(6KEN0kgW%a#1kV@)qm)lPfH8InjAs?3 zc@N_$wm<7sw&g>_0fKFl`<{(afegcl6!O~xW5`z+_p<#ihOyU682j~t zk!SQ7rv@>v5{z5s^1dCwHvb()u1Oeo?c;r2ocAExct<{$8rWg1%(l-pC}Y)oFt+{y zW8ZKXCkq%?gz_Fa&-K#f zOZ^PglWRl9AR>^mxGrxDq+bo4XulRl@rnLX(}#8G)1L^Q^BBQ19SEMB!Su59<9-B> zTg!Vt5y5>fBDhrmOn=x%Kk$L!4g1`?O<>@8Y=%>FVYtG+#JCEk#2zrcdyn9{=MmhU zGGh3QxwaBUe%l|sZ6$(t??v!FAb3CTl|6nid4$6B5BtJ)U0`Z_7p5BbV5-m@hH;J9 zzZzi}Lpc+M!Z7U`408k6KZnDxl>H08!!y;}48yW>3^QS>dyIMGVVD;L!<@1(%sxk? z^Lh-^*nie#+ENSYt)#n7eji~N@(zYhePL*F6NaiRTV@CGhu0HfXhhrO)S*ik82WNd zFkm0UY{uKMj}3&O!FSSS!cd-LjPe)x%+|qB-x~%cb?xl|gLe%W#=eJP(g+wPa{MxG zFKwwu9A%v2kXh5HPY?Fdvtep^3Z^zG#A99$fvI63#})TkhW7PcPo8gJXw2&xhbZ$1 z4Am~dQ0ENwUI{~+S}-U`+uooRFTZ*dHj|dpQ(}(@NVf@K<{WA*2 zUu=JW8d67&FGj0jD&v4*{2_)MpOl;kQ}HS=3_nQBfWdPpj5%!cxA-g%Wt*7Idri*1 zsLN0UmEQw>t`+)}OVDT5Ls0Q`2vW08n$LS}66yUG^4ZV(D0dh6eut?}6zu~{O3LUM z2h+b9FnkVaVwZ1J5}#IOv-f$6_rI2 zMzlShbWYNz=djEs1hw8n8rt<}2*)vj2x?phLCvUB%kKzkTZC{Ss7(t5HGPYq+Uzqc zoIy~r4bXqD0R5-ygc16$41c#pPzm;vW!bM+hx4)qOaf9|yCIVsL9wNg9jsgB+yS>ErdXZ=HmZUGrY@tm&$G7Bi z`eP%@l&2iFw}TSzm$htHY_sK>@@!(-Yp@|Oa18=ysu9?&83NlZLSQ?F{njCHl?j2T zLlF2b9C}%G=;w52JElFm*e)Y^ZzNuY;pJ}nhrW5&1$wm-`cAu{=Xt>TF3$)bUNi3D zg1+NV=r!L+lM21kll<7$+MR)(b*OK*mhn!|%R`AryuJy2doy&a1oSQULf@=7^i4WL z-}nVVJsP!uzTrAzII)Hp2;JOq(9KGKZpKdNrgtD#LN{|dbh95b&j20Q=JX9EmYqb* z;Wg#*d9SZep7netIfw@@-d=ZMkD+>`AHWGpB+9Gh^0t8k&i-1a{;LrU`{zr`PU)vo1 zBVypscl-W5)bJnk9RAzaz&}15{(ssb;NOi1SbPJ4vttqHPkrN7A<(U*>>TLpZlWKi zA%OcD0+LQ6;O=Y$JY0i-WD5eGzC*xsZ=#TRPBdhi3jt4-GQTDQ9?KB$=ps=O0S}fV z;NBUQA4zl~Bm{8JUBKs( zM0A87_o4)t6A%!rhTl_v_&sa_zZ-J+U6_S{`Sjh$00ekcgWvFd@N27q@9#PAeJ~cj zhm+vD;vswob%t*%1$?Xi1K%=F;ak}S-!gjHhJ_%vHo@xZ$>etY1-2mO*x6t{{gf4Ii{9n9=zdaECX{qo}zYg7D6%kB0 zpbIO(@Dg;E!o(CPf43vfZ_9}b=0Nb+a856k%Q zJ~`BkVJ5+Q-)9K?@1_wYp*xU4?1pYX`R?xy-9A1K_D+CqPZ#KTwxwJ;D9y&uq;xj>B2I>~Xc;HFqmn9wb^NnSh5tPqo&%yEJafet! zxxa`uluJGL6sN5F&~0xC-IjXLZR`u3mig;uK(|&w|GD|h%Z6^{7wA^Jgl>5dbj#{M zw`31=i=5Cc_z2znPJEXA#b+92%_$(}L&tvH|1jVms(^pcVE7->!GFhc_^&z&{~3YM z@f~@E7FUPqjh7_aD&JtP7p88+83YK{wt+y*M7J%5i(G&m2c{%-O68jIDOU z*pA~=ImaogMKJyo3u6ZZ$GXj6^zbBp!uYR@A;v%BFe*7_Rh%P_ zOEAg+WBViIxf#Z`Q($bvF=?BU9H(a!!zphUj4dgn#V=wkjLjbq!>CIRp(ADx9LF@{ z*u6Q&P%YY!z9-9VW%(E6*@?WjQAQ+XQD-^F;0lht7387tV0;tPX;-^RlpDr!oFBBD zM_$yoSptkrKM|bYG;IoFlXPMbjEzUb*obyC@Q1P9OBn0W)>=KZj8$k$<;pNtnDH#7EO)=VA zgk#BH96NpE_~H}i`z`s-qanw3H8^G}$3F9q8Tt?L(C74o-o77(!!8&M?8AAMjNw2u z^w>r=A|H)+Ef8_w*xpiZoKY-&PKNw!K|NlA$#-j97 z+2Jt1x(IXjF_?2XSIl1wbDx7AM&N;ueVg7!iCu#GDRWN^kM3i9p2IQlH`2*+I?`OdL&Q5T?`%VyBiSCs5gE&hp zCp-ztecPW{NQ4tE(s5q)KWnd=(nr^ztC7!7j@^|S|2wgj__b?gp71ou=#WX~rqAR?y{F_k!lh}`1D z2qKJliwLf*MC2_d%tRvLAY8;F;w-V3m`O|^#uH8Hr9n0{)eCo)y8Iik|m_hU+YEZ`v;x=)Ua5Jq4F^hOWJS5%`e<-sw!;i#c z;xOSy2;x8e8xa55wSb5q?h)@1kzbqWK#U_c5?103A`0db@rWp-4uxsNBSgHeMCb_0 zdtE@-3DUd{CYBT3NJsgFyNSvyZy>6Z-yX`MFY=iG|9wI^xh-hRHDWgF<0oMxIwK;B zJhOHZ3Pil3-miQKFQO(QGV_T9g1*Wm@5}&#vNKN+$prcEjF1SPKM?VX`LCXlPeX#X zWStittHXuIs=dNvfBLyzrQ75gv2$h_1q8wqAJ5`a(<*9y9L~{e{Pj8-$nen4U_kAgT$EX%~oDL=)jL zwU8jq)KkPUB8s>~q%prTF`8iB^ag|=uP%g}yv>wDezPZ0_5|v&nELoruO#YOk^0V~ z-VbR*586T-S9GR*cWCQm;j!ke@L0E9c&x83JT}nwO)Z7TX7b)zLwIc4C_J`57alvC z3y)nhgvW0BWlxCk*c&Q5_U#uQ`{xLc0}3KXcpUT*9zM;8qr&45%N|-wJP;nfbqVU_ zt0S%u^p{^zf;4`l_iID6A{r7EnE!$JuN-yprQW_R$cH){qArK%FP~tdGUctL+}lJ| z>M@x5tfyYg+he7k67@Ypy=lh|>b!j}ZRtsyYSXqiwDGF&*yu-lY14X*@X(eO9&6q7 z#TB7BTU%&ORTY{Twa`Sa7aIO!TXSTG(1_kb6V^m%EcJyZR4FvUlZD1;5t^WPLZcrc zG`b|A2^c6ee$R!*ceT(QQVESuIiWf5L1^~B7MguOh&n>E*HdWr_zKPLx5R9r*~PqF zy9v_od_!C!LJ1$jpEyR`V19983CmO_+~o0z{Q8skOQAU!OL-|m<3Ctv0{<2ogM#`A z>bj76@1YF^LK7VzG^f`IjieWvi>X3$)kA2KP6*B2#zOOm_B~4wn$%}P9A z0d1!5Z$1&4M3%quUT7}S7w4ns8!w@JZxhNwgHYxa2xVrBP&ywArDKRtrc@Wo=QV`# zNk5@{6fTqx8VKcG=HFf*lsEqs%A~47nebUCujUKo8*&sQj4P8Z6wC)A%b`A>!N!)T$TFI6Rq3RQ)YLRDj`P&Irb zRIUAmO5IGTx;lhvKzE@UUP`DY{3TQ~stVQo%0ji|y-=-;7OJ&sq0;UWs`Y}<^Lp<8 z(pHgfS!Bf8%s!61o@s;u$3DweuLbd8T-R2vslry|rRh-Ldz7o|{*qaNP0yWd;d*n_s}gsLs|ZIUBYb+U!(?~X!M zE{lFC6w2Q-h4Sk#A$zkx$a3}zS;kHwb1?q0m5@D+5wgeEgzTZ0klk-5WOw@r+3hVt zcJrl>-RLM}*GZp{DP-6B3)z+HLUvgzWbw0w?9x^ti(|RC`a*VbqmW&AOjH!I^X-KU zorTQSTgar*LUzuFC@EyI%s*F*2xB^BqBdz$S=LT|e~~xKTw5b#iIR{deGsx+t%dB) zN+G*fNZm-E{9MSM`wCg=Pa$)WXHG-fb6&{4P8ITEWrV!EUdU_Og}mu8A(vkg^3It; zKENpCqq2p3`T`+elqKYATL}5q5<>}`gnTIZ zk60_@V*`bJvY(L8>L=t2eaL5!kZCzd5HR)5b}}pguMR;A@4#TD9;Od3+i5v zv{mN_dFccp`*T>xzBCrHN%Y^SsX{g+LdXWP?)$tGvYxEl?oWlR8|$-csE~DW619Y^ zGxPsty>w!IbYy*ZVEfSc5o3s0#-A{c?dacer2WdWl?k?~E6_a|DyaK^`IMDSn%v7mRzO|eOec>`|Ks)bI(f5g|532bE)(*< z!dX8M@^bz{Ua^>vSNbgERlJFGg6*g>+gOF=LSE*ake8r7zvEettiw0tmvw^oa8us9 zKX|Xw@0WH8+1ZIgcAS0)W!?DG@4M-v4XpQN1BGl(DWO==MJN^&7mAr0p_u$pC?+s% z>uyEMbBSC@$WgI_@|^$$V@`f#w`@B9E2y~V!Vw| zw7o18?Z-0zx=_gTghFv!D3t4jLRE!eS#=Q0=L$uKG@HXOvbbuu3Qry9&kqNTGPa@~$gFk=sou-uMf}$In9XEl?bt>;tP+63QCw zgtC@ODC<-a%DN|ovYuNg>mMV&5YGtqx%DP6|F=-qZb5nn%XMM-spK_}y!Q!Z=~T*O ze^_FqP!`h*Wzlv*@#lk3{H9I6Xxk5^P<%BD#iz+a@s55iI4u;}i)crHP^7jJiYE%8 zxJx^(vrPO%p^)hJm~lcO=y#J!D11q?tFcgQm?c!->_YXiw@|$qCsg^vgevD>p?bwW zB6EyTxqAs!2K!0pWuZ!AUzfsu?5R*!k_mO4pF-WDyihAT2(`y}q3*d{s0WyZdT5?dkG?9@liCRN zj9x-L_qR|lo-NcX<_YzhZ$iDHwoq?&3HA2rLcQy(Q12@$)CU&|wQq`0>v{@x&~~9V zX@xqJGAtv7`tS>(4!bSX;f)yP5dSc(tWcY033W(+p*98xwVrAI9^?}))O%J7^|l+7 z^Gv8$JB4~#L!n-9Q>f>Zp+1#_ddfwio=6+U*B9zBlZ1NYWT75bPN)YD5bA+-g__T8 zbq}Ho%XC^K)M|%N%bDN0K&YGT66yvlTPsVbt6mrC@|%UaWE-LSlPv{Ttu6)6i<5$N z8Y%c`f)rBox)jn|Cxz@9CxzVVA%&JNDus5mN}=m2NTKIDNujQ7lDT17$t*h}nU>d= zOokZA^t`SVT=S3=-21l_JkcZtPb)!wouuHjCQ?X?VCnFbmeS$F2c^R|K1w04Jfu*K zpA@?IofNuZjKnYGBy+QCl6k{O$$UIUGCy@omg1cy%Zfpg!m~2G}58ZDUxs7nUe363zBcZHOcp6 zl;oRnRPw9#R`Tn&N%C7WO!7N6PV#$pL-H$pC;1h)CI4bKB!79h3eeP*0#+@M0?w6?0_-|TSNyc38|NwMtZpf=UP;NnY6HoCzFG==|60;}_K@@k z7D)Pgi=?26Ev29#;Zl&#PbnyCh!pfyNCuBml40F4$q=3+8Saji3~6r3SbB+M99UE` z?ma0Pzc!Fevo}hnV_hZFqaL=r8g5%&i$k`&4&k=Eo`tr&0fTLM-g;Z!l*hKbl}&7U zJFeOC0{htV#B*ETxg)l`d%taYX%1Un)&g5z?mc3FEidnYE$@@umS428Ex)qNmfvii zEx)U$Eq~BpTmH~jw)`QFZTX(ZZTVxC*z%X2u;uTYYReCKW6KYpZ_AIqWy_BnWP7!9 zpDjOQlE&GSsmVKp} zE!+QzEqg;tTlUsNwrs;xTlTRyTlVScw(JYjY}vP9%f4ORmVK$6E&IefTefe1TlTib zw(Qj=TlSR6w(Q<3Y}xI#w(Ms9w(NE+)1{0pd%$a3_Gp(ad)7Bw_DZJjPqt+rtzqLE zKwI|xCbsPSCbpbX*|wbeCu}(i4_k(Jv@K_fWXsu7-z@zj=k=!Y#g%w@|xmTSvB)!3GsFwmB3uW8Hu zu*{bCmmQvG{=oCnLU>-EO}vHYwY~7XJPw{$p1?DyG(2y9fajgv@Z>iFp3lF)GiwJt zKOTqY_cVBZEdkF@!{PZs3D5VB3G(}JotOmApAX?xdKA1WC&R1i9C%fj53f?Q;8ipg zo_~4~YvEPwJiN-hhF7)w@M@@nS91rv8k43$0po6Xwd@5i`CfSOdtt9GZ{XGQJ-oag zz-vqzyk=a2*P>_eS`z}V9fk1nX$3Fe84O3jYkvuN`8|S{Wgl@BUOMvLPQI&1Gd~Aj z(+0q6+)a3m-2tzO#o;wO0$xje;KeCOn5#8JnwCx|7aT~!}D})cphsGPu){^ZY&GW>GY{*dwBMw zFMG5gR`EIso_z}8Ib z=fHEj3Z6SWh+Xj9-3^|;)G4GbJR?|tM^3}j@`~{yM0Z|C!BdQ6-Vt~nCEbax@I3yK z;VgK@ctKmqAKI!u&{kUuZFNs*Yk5FhZ#J}zen8uzA++tsLi^8tXghU>_TRsu?Oq#N zp68?;R1Dg|Tc91;4chVRp&d*5;g_Htu?^Z$4WOMc6WZA)pURZ^=%EUAg|*ap*>0-hb3rD zkf+I`ex2W4(~2<`St(DFe~}PC)y0 zJG8IJ(_RPKXSbkzJcIl+QzJh#$_3w4Q-o$pp|!qR>f;&0O_o(2L-fT1Jkyo&TSvFPL@F{qwQ^1 zLEHQ|^V{(HCh60l?Lhx2Z!oythxOC8En#~OX1ZLk^IhHs#)^Bmg0w?bR_zki1h zCit)Sw`2G(@mo0f@frWGZsI@2Kk{GW)8JF=EdM>elK)o!L}>V*^&I|h{Wrs-jQ`7j zlMDD%3WZN){zJV&8PYCi9{;^w`Z(kB;8W!fe5$@D4gaTJ?mElVgio0p3>T9o5%dm{CPc$J8kVjQmM%<>HMX-!t0n4;3Smv~V zWo}bgX6aeZO1ZCKS#E@7-8fj*e}ZMhQdl;Ag=I@gST;q&qTLJ28aXU$k1$>tmL0yZ z?0x{t9?~AX3QOQmSPplD<+wjAQQ5G>P~N#bSgw%vRzFyhd&ALBEFXfW<(+1$KqS=Qk`nH&b`|b^S!z7e@Oc5n81fLd&Pq zw&IKrVg4@Wm7*VL!`iyAtZD(v5(E81AI@4p|C7%G>amD4Gk3F&=)-{-tYgylrLX(9 zWH=L+Ue93ZS%#qu*0vt7s@B2Ey(ZSKPhjo7lBfl1w^p$BC=F}BAXtYUgmv_MSjSml z9eWhk;o-0jsswAV7qIp&gmv&%=C6g7XAW2=^yGCnScmz*>Ro{-OSC4<9+v4(exF(H zJb5#Xd)KU8Cs1A$^6_Tg8dwM3A?;hrd`=YMwSxFgneSm8^OF1@!8&O)tdm1wol2h5 z>cKjrFRZidu+EK#bwN#7`E8JOsV}Un=fS$EJghsL!g`SW^jcWWZdfDqupS)<>#>ip zT4}>2^1b&8))d;7UXI~Eus$bm``?U{PdfEZ&u5&pFWf z)@wPi+UTQ`gtZu~(LZ24UJBOm)v$(Eh1Dm%2@ZbxRhk>n6jx z_&Du2NqdGO>LdR_{)zl&b)Ze7n3ss?I@Eh9{kD|xCG_1wIjnPzkeuui1! zCo+BfBv?mvgq7b_T6=AvpWndR#$DP|?@n2{X|K*Rx|M$Hhgm^dcpHI6P z`^0~$n=nfJ=baaSzQN_caH#mtSwkLobEjPv|0&#e9p8^TvDxUOE*K#GLz{~KpbYUJ zc#hCt{D+&lAACssCzheldg4Es`^Hma#DC^<@t-{el!fL-XcwpN(&E4BmiVtsBAyKD z!_yP))k1gQ9`QfaUi=T+z;*FIzCrv?PvX99N%6lJ&YfH6-Lr5nd8PP=)Dr(7?iPnQ z#h?Ft`hS9s>y!9rG!g&gBHVYrDE^Lx#P@Q?`Jwoy&lG<*G}G^L-#LJ?1ny2FBYtxI zzx)vYq)go1JtO{4`a%Dw_}k&o}MTEC(+BXOYoXU{P)I+ z|1R$PZp$bBYc`4h67C8wm@C&Zhsw1aHRW3Vx8lE|J2F-Qk>bCaKGvbj^^C^`^u7^Y zZ=6ONzSdNxF5|N-5xpFO7J8daKa-(9o^iydmjCby;@_KgomBkWzZ82#o7k(b6notn zVsA7~>@D7iy?sfscU&#@POZh>d4kxx2Z_BO@j>sw7qJf*Aojis#oqg0vG*D+_8#x3 z(@gALZqu%?*xQkA6DIanvCt#k9#jKM1RB^rV@4+d*jYxZ`xDr zEz`vQ+jZ!jq>o!-?_ZU^a}iz;`;cv7ADJljF^|CSVjrDV?4zhR&LDPtq}Uha7W>jo zVqe|_91#18Sz=$gPVB4Zi2e6-V#goX{hBqzzN@;}_q`JP0q7sf5c{cKV*lGH_8X7I ze(R3d)g$)5pm`2n&Mp`GrC-E;bE(+xOcDD(T_}eS?I-r@@OASPdC`kc30?9-_;vA5XA`igzT53vt^ii~dbxF7tCE_WjLTy#E30^Ssn zfY|QnaFy6Qz*FaMj0ZaFR7&jliLm!ti%#HaU>svI7%(OSj)}d`QL*>T37zSThaYXr zpwFU=X>rD6pV*s`u2WC!l`4q69N$De(L^Fn`%1(g1^8Ym-*mmscTMl&^Wy>E@_WFy zMITg`h=<7%5q3->9#;mXB_hI2n(w8C4F&i3Zs5-n5y*E-@8SpJCf`oI$#(^B-Ij>k zdjQ|_`)3;8MeQOH`cfjU(B^V;zK_axRF&_QUi(!d1bSELTZ%}8KY0OsGc}ZN#Xc%7 z5mC^LyettBjh}OK>o&G$Y*?v9L6umDg21guBwZWcBQ7{4+Os=dHKfRe>8pJ z=RZD@R3Z|y!($x&K@8{)o`V)5qXQk)L@!44wLv28AkY0D=zAb(zWdktt%Nn{AYn}p zGd{Gt$hQ#x+z#5KKhig$9Z-fg&jI=hKu34b?ai6cB;TL!Giq;i$~OhCWg`z=-r!r0 z*PwBkvHOGXNuRvRx2pc&+f;A)P8EK6+kOEqRsRL|%R48(5xRpf%`7odyzXs16OM#Ad@sO)d-l?gaQ#%lnzW|B~2J z(fCv&zZAZKcz=#9aS8YW%FEEUTsQpg;76<+{8dDbYWRY!jX$3H_;zcs7N2a*@Jn=4 zVw>i|@03AeTi{=Vbu9giS2<8)t)8s>63yF*pB!; zYL1^)^fA0DzFG0#HYkMr;}Y972Q)&Vdl$d8@LH}Xd=17A>r4Ew#^KX8r^K3()BB&? zRnft@d|BE{hQWnk_`%!KqA6ki7Bk5-c8=`Ts&py6%<9sK1x@J;`Dt0X1%0WL{ODn#C&^jk-glJS?H zO5d(w^xuW{Wuapihb2NB*3;14Cyq?V#gUzK-UM;v`&AryHo->aa5j1Sdg?4zl!Nl*ArLu0;oQQQKM zt41;T%>MOrFS{f?y)8-cg6=**YsdfT+Tnw|KIqrYHvzt`MZZQ#djGQ|z2%z(Zy!n0 zyK;oJ!AA0zfM?`&W9)B9>PwN-_=1wkU71w=FP9qMMpC1#lFI+)Q|*}~_2Fho4gXtG zUqnl4;yX#r7%6GSd6M>XCy-Op%z={TJ3~^vgC#ZnsiZn8Nb09=)XglZAD>7n{wh=7 zyo7F9NqyyjP6&8FzjMJe@{UPrEctKf_x)Mgc#|KHv*-p}O znI$cEO~TpauadOPeIzZ*c}dH@OwzJtNLnTjaFAXuY59vv8gousA<}u*(FS_?x=C8G zXh|zGM$+)5m{wu0q*Xx=)#D_s+F8CaL&Z zOLbf!%u4=$dDA5Ia|v{FLsI>BO6s=C^7Y$5Nj7zo`XS`5g2= zLR=Yh#pP}-uCzVkN-8F}6|W)au7cyI^$v!Nd#F76w-zU&j1qmj5$nuF)! zO6~^Q(ykNrUQ@rgxW1Pc7rrIZeP>C!`HrNU8cN#!W0JNeQPNiYDrqZy;5C(g@=4mB zT9USNo220vC~Z*}NgGpI()zWOv{BXJVWgzzXe;RjQ{gwaq!$@M{FtN{noWEMVJAud zb(f^)&qKU9;S)*Eb3xMc1|r)eNiPuu9Mr8vyLsgOE$L-8NP5*^Nw3{Z((84T^oDyR zy~)3l-qI=QZ9jl&lHRN)<&Px2A~KY!4egSWo@E{8a4B#E<9b8)slCeFDLlAbh1+*yi>JAV}TN8Ft6r+bkxEiYiK zTqh+x13h?C#cfCtx6d(2|Bj3<@_}FJ7oC0eg4U!1X!Dz-|CmC%f|Bl0Nq1+Z zO?FA=Jkq(gs5r+I6K9X!;_RM9oCC_i->>2v|4^J$?eO(hoU>|*bI2rdw#|U=-^DpP zTGE|W#r>08+&Qm{``6##x0tvK`iZ+}5pkECE$&k7iDwpffn}7J5O=P<;{N4=xU-KH zckYMc&XbM2BeZWuyuP>#-4=IQ-d3xIiMw_tao0a4?k2Ou-2z&zvWUBNg1FlZB29Ut zJmRh~Lfqxoio4h@=>HAvC*uD3thlWzZp)wIwoO2n(9AYm+$L-;eTJmJ9V_X)8>i1$ zC~1-3B<;xu*m0EbzpZ`5 zx#=cxY~UR-_6ZmFye{HiQ%2nD+KGE#(JMa-21PJdl%!k8Q#`Z6Zdk)e#t6$K2E&Q|FDvA!_S6$8@96#IgV53 zEb%{&vTI=+|G{={KNI)O0C6khe-(O{e;4<~DdIl8S=@*6(|2v!V_i<9-; z-F~0A$D7dg3vn-=B<_WC(BU2QivGqU>p0{bn?f30Ok2UYAQNZR?wN;3=ODkmWE5E; z8F@-eMz$`J;fJj6?~41ORa}!Uh^yBqakZKxE`FbKP53CTqqW8LqOxR^+awvaDoRGR zkc?dWBqP&3aXVXxYre0zhBA+J*deYa@5EL2l(=dR z6j!4tarG!FuCc6RGuDV}PF8WPXDz$27M{jPMw3v<=(tESdbXB~KGBlVo4hW6kmn%u zlZ?)NC8O;f$!OXU-fF;e0eJlr9tV>qY;;XBc$d#;nyp9#J)`|c$!IfFGFnuVjOM%P z^9uSw&MM_4qhu@ezewC(qqyA>;(o$hwZSaT59#9kyCAk*OI%eZi0c4z_uX^iI#Z4K z5OG~UOgs204xd%(qRS_eQ307tpMoYjC<lN5zGeKuGqDEyU``p#oHA@bcEfYI*DCIy6v;475KsG&;^`hBo({;E zlm7Dd5>LyH;^}x)JU#yrPv1@;Q9S)Oh-Xkc@eG|Sp3$GhGcjB|Rl9159)S2C7xF0#M79vh9|{S_qlkg`HH7dANlsJhqxE;d`3 zXJ{ta*iYOic@FMJi#r{8y44cTg!kf^w^KZe7m8;Id@PA0ew(~x@ywbpo|y@tu6U;O z7tffN;u&5WK1DnuYKmtp^~cA;?=O`170;X&;#rDZYkm^X#@XW83hnL3#Iy5*cy`|> zTufdy%Ib(`brN-vYr$W%qo3&y#WVGhc&4ou&m?G#%8cy&(9<*Vw8%jBe&VU#NIWIx ziYE{JW&bSSlGez#KgP+o@vQ6rgh=}9Kyi0u4DM!#`%^A)hlGmz2J1{%BXNJWOGaks zb{-|30e_>@8RF@MZo4lRPp6IOb)|T^px>U*?z4=tZQ|*AhkW$V_JVlVG!-xBQ{HKX z#55J5DEZDg*7vce z#Jgjycn{gd%UPQD+%WN;-yz<=E{pf-Qt@7k5$}yBpgHOB;=K?m-ZPEFd$NyskJ0BL z=^&bE z$#=_p`Cio{--m~Zmvas8y~g4VpD*4g@c6X0cppC!@1qLhjjSf#keA{OY%5;-1MxoG z1ICK?K{4?LKtF)I2Zt%23+hp46y@*48#)nM=;t~7zv?I6*q_Dw@`QL_KBuFk+&wYw-T@974H>joxdX9!+FHJeUo@sp@TWo#XAdr=d2g+ z#6QJ5h(6mPcY##-5zROjT_c_>*j?KK;<>R&zT@BgdxPonz35%}o^`x@=RN#;=?n6m zZ!CFtJx8zC#CsWi{#6h?{4U1)|lk+(tS^4U4C#VPQM!QGf0YHOv`JyE9-6v_`05 z#zr+vr_I#WYMA6w!}y%=S3(WL%xdWGRYT9bYG|LJK7~N3^XgM0v--5{uZGD%&}*SS zv-_!GT|G5yzOROD-PN$Y2H_V_R}EW(!6wi`4O>=Fr>7b=Z2?`tK=7J8_}EOFO~cf% zk^BwN-f)RD{A{@nok!68iLfZ?ebj}&P5G(QR}E`6s9|Lv(D2U$ScgP?$3`X|8m(l-6HPkPphWzON0Cqo>aBFt)Cb`t7+!Xa` zjf`!2s86G6>Qk#X@^4n3+6BlrsL$Xx^b0@BI)P)1!8`Olp1j@wnHC_+oG5rXgpPkj zXQapft-djB)%TIF`i3l0-+S}b_wOp`5`G5fQ$x2^YUtBm4TI3Z0OG?2K&zG-PN%6M zoOX`g>ihK|7^S|ckJUH*H}!R=fE>i9s;?^^kmrK7i#*px@;8J0>gy~;n_}vlatnL| z3)DB6wBx+`eu0lq(EkAK_bt@-T?Uw~zHdms9St^6^;02(5=E2k{tbVp&_4}zU@mwGa@z3P_4GjmO z4>&>I1>Kn04eS?Oo@Abp#dpHL=$NrlAriS*rcs9zY;SBb6H%<)! zs)n03HC*>s!;PV8@GqqX*-6jhbqAP*b&+YASzCO@&1LK8#hrzya#_ca-``vijYwuYUfw)bIKZ z^$X9ee(t$yEdEJNzs0JlYi~7myQ-$1^VHP43SpRW|?P}^$OHG|=-=!UOMuQ(9fpjA^b$q9$4$$rRFKIvWhtTFdG)B?>Bz?V7Q~zyh z8j@K}gJ%-&s;0h`)ij{48vBK)Y3xEZO@z+GE1(;Ap{7ajJNXzmOZrbWjf00V&>x!( z1d*oRxLL$Ekq;kZkCSc<$`YslvDwLg27J{ts+5|Br>kjLBFGOP6Y1NlramW-|2aI9 z?lc1#ZmX%?0X4Nj?p8O{)S{P~8l#W;Q`J?bu=Sf<9>z17%ig_;Jwg`Zn$ z8W9DcJM3rj@}MuqwkJIIScWbcllC{&)cm1Zn*F1eMz7Vcs2VX^JjGV zvN~9=rdNS#3LU4WTj$jj)Im*2Gt^QwP%Ui^s-?>swe+Z=mR|8{={=0_rdoPlBCokx zx;;@#*M(~7ycUcgJ&LkrfONMA@@lH3;}*ilYUy5{{9|A-Wx3SctcjXio+kfy@+PRI z&jqy%fac)gY8f_2EyM4sW#|gE^gBxXA8KwGuI45+)ZA{lT82QIv9(MprIv|j)iNPB z$V7acHC=+K(Bpma+6P?jpEF+yf$sPgM(RkY(H_Lc{+r+oqP$yVWwX zFMK?Lhu`QAUIr(rWgs;A{wJdka<@=Rr+jMp?VVcM#i6UAYH8J5EiJmLrD-*_G(4x4 zdWY3g107ZPqLw0FwPb3grog3YvUgFFWC&*i_<~HTv13itT|Iycws}S+A@YfUF2^a@Z>J}M@0 z@ej2$N6$^MfktZqyfxWC-W(+Y<9JjdZCsQ{%R@yms*Np>xGfG zUiZ8^fG|4AMl1}Y&y?SteQ*yqULO!)%;60!oSt*hYfsUevC^|%g@75SuH`v$LKMU5!TGk&k&+uL1Nr@&pJbzLGXscwP#rWzcTgcBXAQ z`N!4VyM>y2@k|V_NBIghw~tqIhmkyIS=2nDms+MbRm%+eoQ!VzOvaXxx5q6)^gEpX zCT^q-&*9Ae)Kj%gt3_ROJqevmrtUOkpH91JOVl#8pIWAD1&r5}ATX3V+d)?9?*+xc zA@Y%D>c0Tno6>{48_Xwlk&!+}bcHwcJQSKkX)^@-9t@p{Uz&_0_7pS@GWHt9}jeatz zx!(pg4}s5d!_hO(*NFS*mgj2pdFCn;G9cfyIC#ePCN8I5EO-lcPzPHXirx2bs+PDM zYI%!|#9vlRG|$o%`q}HJ<_!nbysdC3x{D`?A*sbOe=Dnv?)spmHEveK`{z@BUdW+4xsjrp~-6(6XmNevcodcp4 z=UKqGr!}U26y<-Y-p+x{p*`rFFaH0 z;a_=XsNVzqb;^v5!f*F3#NoYl5$2O>ydzA&P7g8{((m9zwJbV7yZ356Sy!#+rYNqg z)hgvcakXA~q}Hpy60Zb)5O>p8HMRah{zYH4>M6C}nxob`)b~ffCz`3{P2;!%X~^BH}8t_|q-3p!ZdIv4BNJl?^2^X|5ScMEKE z%ud?;tCouiYJDF4V^(%5I*OtplJrp^V@*lfx!wYtZu z)%A~B(S!B#MYSg0CSHyF(rWeatbQA=*7S2~b@EIypIcLTPMCwOY3RiLn_7RgP^$InBj$<@VB^`roq59tlvm3#AIkfyWw(X-?3r5j zVn^Etk%w(<;+R$mPfj{lm0s1|Bo-}LpL2w`eF!%1q=J&jSUSar3 z`JWuy;s{%-Wee>#k+%u`ZGzUOw(yM}HWmlG2W&*_i{i zSev(?%gxYY_*cuS29QRcT9tVhLFcvZkv6HN$~3_I zRILH;CyeQT`?4CWpVcQ2W+J_qx~JK*!BJ`Dyi`sHeq?Arc9sJ^z%_QS3IW1Yg}kq>W<1xARns)=tUqUHT`zDJvw>K6BpD>;T@6 zPN8lmH6I+V<~>i;ypuifjvVZRCh>00o@*=fRny4wYMOpTP4k98r!ajW=ibq3-sj}q zzq*9nMyWa1Fg53`!5q;BY~=i>FKP6V7di9p2V2yf=d7A@GneG94}J$7 zK~eGxP|mqxUiMIVe7nLK==Y~xYPKFwv!ylj z7%;V0v(ZP*2A-84+tieSoUV3iN@Xrg-mj(<#y7PYXQd<5^ai_*dWU{4sBvExHC~8O z)4K_3N=LsRd#Ld%ec0$D(^%x4qDJppLU_#52@X|xWK98g05+awCuP{k zPtD*B-u-y?e3+xZ7gdwn4<2~-5*Mi{&ZEW~Bhc>`HJ%x$#uJO-=ZzY#=2qj`BWh${ zVcfMtjeBycku}$N__`X8&0sv+z%O=jHLn^kl_70X)2#;Bz%Fb9I;#d#Z!EMKv%5Le zcs*K8kCv$E*=aSsZp=CIYBj!Wtj2fi)aby@($=dnsQ@U9oO#q_=5~+|pr0DAKwHX|a!y*4 zdRM>=fSpxAPdJ7$R~if$(@M0hxLVB>u&;8<)LeEJA-XMx%~i+*qA9l!)}oBO(md}a z(l~2-qvoQect0Yq&>r5A-Z5`8$K@T(9FFZ|WiGIJuoIq9D`zEU8@9pR=5VR;pQY%V z^ZH{vL+nk>2Kvs8O&4T-EJ|3CzDhL!*j))>{X@1AN6E*YOO_N~WRnv}tK!}41@8TR;nDb)HoGG)`yx@0$AB?%F1p4GXD{n08C~IEHaMohh@Y3znQtH1r?*gTH zHz+lfwVAcI*l6mpP8RkMo&x@$A8Rq~3zlac=6x*RSLmHnORh$&kADN!=?H{^*Re<|fN)bh(g&TNmWC9_Q}Kap>vjwMnpeizjI{ji#y^LP)s!ddPM zHGe?9I{|7Af25{8=yvU5!dg5(jMw)t&Qf{D$T}68yl>=;SF4t>zop*1y>2ecq!Y(^=1veNsi%nu4r-A6R#@aAx%iUC_rT^l@&yn%geN z77H>*V0%HyYKmy8rYDE74W9L{|JgX>J?pcYLy`OEWX83IS}K*}ou{^1YO$|t!``(c zd(3X^;kxG}tf|(n>i^t=8Yzx3zPqwLSZ$_QBwH((Ikua4yp7pjumG zR%??JYHe6it@XC6wN6>J)(Ruct=77}YOPOx!!l}ZLFA8exL(!_RaO!z@4Y#gL$9;Wt!l};&>GH@PjIeWfcF&Uu$@n_Dc02= z=QuNHYYJKk32l`TpDQ?_lq<)WkTJ;CGGNYN^58oE@8sXHHPw z`HxQktlKA;Gfwf2cJ{TJuQpY4#2}t|Xbya+mOsaHw!KlU+3%{g)Je5gnWwfE_tn;n z=b`d>#_%EbP(!WNIoqh!m+_sUwvOx6)^)MkdN2+>Dygj-^|~TMm!;UlJ?vnY+S+c! z2693tQf-abs;wci)_;WV&~>dzYO6U?ZFLN!(OYZ+0M~hEykRza~)aV%BuObKl3|tZ9&epDji|nMBW0d zXSP9Vjya;{)6|{C_;z})=BBJgt@^0BGh^R%m73dRS99aP)!e{H*g?%5S!*U&W=}?b zRHB;GnUDTkfB)zJR%?kC`*-^_)=`SL+;jnQu_* zyh&=E+YR`Dt(j+V z0HNf2NC%MroW58CrZCq}K0|(pT8Fn{J<80S%vx<2z#1CMI*AOY53n!2sOCRTA=3=r z6*{t(@r<-)PH2mbv}Db1{)%vmS~_6=z0%b>y(0WKRO_NcYF!FnOXJnLxQSX997iUe z)48#r8gd2`((k+*)W@b4EGJw_z7tFV(~)_PT4y2m^uoxOsMaAFoZtGYr6RnPc%we^ z|5U>+w;B#_Q^WBmYBXV*``W9|@=NM7eT4c9o~ShM3mbLzt!=F_?!{qn!36-nx| zz^FdUPQdeA_~O5G2l=nt;WP9ZqlRNwk>RZxPV)b$ll;f-B>#6i$^Z0DUPGo2YBemLh ziKDZV{4eli5!ywlVebMptVM79AH=7zU44E<&a_U*-BG*`XN&iXm;X*~RiBQ8Q|7AA z?}yP9`q|*8K5NkJnp32Es?V-(>T`tuB3?YLhRdJSaE<@k-GKKSm-*jUV>M_WHC#o` zEBtTr3jc|``Vt`L)yc|NlBxe6rcjSDkrDtor32;s5AY(qp7^Lf>l&%gqLGZxB5`GI zDXu(A#LaiJ+!ZE^D_248#cUUMiH?+c#g)CPxbRQo%FEr4;+?n~lPu}p6_Wm`GIug= zi8FY%xXs7JUFL?kDn^T|dKln-QT-$0ss;S|ij(hLI=e3uXNfxE%;zJ{LQBM1^tCvP z{K5T(F2uRdlWU|nZ7Jg9E}pX}cYG?3;_gu`aW?5MPJHP&H-vJxV~#kVkoTTD5q!VY z<@* zhQ-8PryJ!zQE!~MtNlwnMBFvEPlo=8@4_bVfAd?L7q zS;alF8uw-3X*eQbN9FNt6`sGJpN0wtTfdm#=;q z^40uIQu39PuUU)BSD$m@NNg>Rx1r+Tn+cBiJ>t0S6vx&D;uy4qd$xDE+dD#%VzW!q zeeT3wog_(roRg$u36gYXnIzpPAW1K?a(DPIab(WV9oj_h({dlT-!XA4Jt&R?yTx(I zD~?;xx_ea|5B9_7a!HP_D9IllN%H3&lAK&vlAWQF>>eh`-!@6|k6=mhIV~xEizJ0Z z))cEqitUV~{9Hy-%;zM z)P}F>10r8t__slcd!8PW(4XsdiUVs`^nj9^96cDpf#3NvRkkDJ54*O5W#^ zl4T40cax8$W94I>yzgXbEp-NJ7j>Nr;Y;gy%&i;nfRCh&mw& z7ne)Ij0KWVxq*E6GFv|Qd*s9Z^YUS{%7@)a@x$!Z-69p%OeSV zQ#N7U97#B`LlUkwm4tw|lJMZ8Bs{Dn34B-j(~E-g>3xcPay*q!>5b)+XQ_Pp5hS0y z2PM%bOcMQKCDB?+5;MJ!MDD65<}4zKxh4=lCW$|nrp!SZ`B}D0V%E8Y319)ZNIFmw zf1zIXl7z*{`wCXmt{HR($fvY_%NqDjr{T!2o*MCXE*FZ_|wvdm06Xc_hLlQiDBq7}<3GRyW(e_b37W*KH z<-L+vcZq!LG*%LuUz5aU$kC*}B-V3EV(r3`SZ4=uV&t* z%S(A4$~Q}5t;6I`kwoqa#hvad@#|(u+=ibeZbdeUTRK_dmiGnU!FuxMNZjwoByQ~* ziCc44;#L)rxRv*aQ@*90#O-iM+*azXZ6L@V2SU4AJEsx zaS}gZwZsn^FY&|YNc^~35e8;-^g^URdI1_aY1hKS}(|&Y(1C z1{wl`#7`dvpfz=}#E%_JUo9oRS(L;VDk=C{L1Ith?_uw5i9L<)hkuVqtfv8d z21#7YMiSTki^TP<1g@fwrxM3~vA7w#ByJ&kTTo2mX8tR2V>eSCk33Tehf3V!J?J}1 z;*S0!@p}eI{K2OZ&v)434?d9ieOV=b$4ZIca!um5!teG<62Hw!dWXbs$N*<0e)BTw zAlsha62JEmVY0;Uo+oi<=Soa$n#9B$mzc**Bqpkz#5{T;F;6x~^yh^V{Uw`3C(V^; z?x{zA{v^@yF%tdAN1}sjNp#=>i4NZ>(NBCOII)RmZf8zkm_060iIR$>BPN{oHC!~_~Zd5L*=9;_q|?T3x2 z_g-S|L`n>EcZ}qhnDdwE>zu^k3ovF}e~IqVMWU`M^twzV!w@ z|G{`5!;5DU^SY+Q#I2W@SIF?Zn#4T)nQ>dg7;crA=*f&{L-Ou0p1ma|^1rdYAu-Rs z(DsMKyrh3;X5x$y-&u`}uPKo))=T7*E)p5l4FpSgT6KxUuUX`yq7t5Jm2lThP+h{4 zYD)N<))M}rCFm>R@pcLSkV@WB34iiY!UM-k_`M^9TP6JI1_^(+Q^G&DmGDnpiCd`m z8<;KOuQp3~czFrGxmv<6+>`JVZzcS&zl86Lk?^h1T?)PNBPG21UlQKnk%U)RF5#t9 zB)oEE32#W<&Z{JR=+6>9wTpz~UnzW7Rf+ulphRx^B9Xf%O61{I5_zh)ME)^L!o$i* zc*IOk!1OLE3WlsiBe*t5B zO(I{>M_f*ce7#m8jQAlk<1=XYS&7($FQAnNB}WDRnQ8;# zC-X?ysWN~zM}s76S5paFTSCH?43sd|jj-7>By2`q37bgXplK4;u7`xx7%E|9uS-}p zyM#3ild#tLB&@?532S=|-$>^=z2T9m1uY_G6EMXUK zO4zXl61HQQgsr+D5k3Esh^ENbYM+Ev_*=rVXOhq#`y|wTLc&b1C9KFn32U4oVf}JT z*o4y(HYYo>wUDs6Z6s{LNoaJGu&wLRO`1gfep4dW50;1x_}E(i3AB)iHMi&wnT|v6 zBs{NNPrZKdi;t@{t0iI?{+nh_k%%$)%^Dgc5d)8bxe_rlR3fI~gKOa^iCBdns14O5 zVhg^pwj=AVz5idvN62?nB97oc>Iiv<3P{wmJrWiANTTi};@b;-g~7|e|AlQO>heB` zh;ELrt$GsiK9fXz_!)m&gCrsrA7IbF`Yd(nzAB^GnprNQrv=QKH^%m#BBSBr5(D{=1ST!i8LKp&ho8_FE<5 z7&2|Zx7lpUMs$;i?(o>A5`Dve=M{{>VQi}{wsa0ZU7PTwwvIjn2+v4_=7&ZCKEFQU zE6ofKha}?FO^FIBCsAR)OH|ZMiFy<*QIWS8k3td^TvH+v8cW3IXQbmLDlkK$Zf2LL zlXi((vq7RJ{{bK9aU61t`$wXNUzMo2ZzO8#8i_jlL!xd`FQltPJy|4C&ynvH{ls2D zwh|KcCIy5_RNOKNuQyo2`L=U-$-5GsFSCSaEi9qsMMBHtHb{tjuY^3GAtCn`Nyz2b67uI#3Ayw^LhO$uewEaZ-h5rHu z_!%J~Nz0))5 zv|Mk(HKaF4Xmw<&-vPbRS4Y;fO5Ff!Z~0A(#SsZ@(FT1-NqC8;)G03ExtM=F8zc<> zdts;cNXRaK37HNbeRfMow=xnkh&gE#^XkT25~fiS7R!3+I|yBGlhCg9kmtX(^oxWy zVf}8_5$qveS3(x$laRS*C1lzl2^m{OLI!V=kltM+c;O=nUT}agRf3oGkl@vWC3s_3 z3EtXVg16_FkWpj6P6^&|R6@qIlaSG8D2t(vSAy5mew_<^BVIy+H$+hWD|P!y@UGuz z|6f@(37)q`f=6|h;O;+3aLd*b+;E@-H!3K>t;p+kUxJ6uli*29BzX2k37+3ag69^M z;K^Ggc*s);?&_xAAPH_3b1F#4{E+}%&RHTMQ~!lec=)}%1kYG2!NZ2bqd|h(!cWU) z^h>&3iUfCQFTuSYOYne1>Kv2c(WfPN{0n$oA;Hs-Y5F}0=~P-mI#iY5DgB5alaS5_ zC8TRJ3F)2*nO>sXoQ&N_cw0#*64HO8g!FAgKJ8chCBaKhNl3f?5>mOjgqZ8ezmZkt z-~BTsXwhK_>R4Wa$`q5JTpc7ReN zc%6`tDe=fyTS6B4NXTO3sGL)RE1|pcjCE0FnoPu=XdB(TLK32Yri{w4|RctQd@{Y^Mp z0^0;jVB>}o$k|R{+1(OYR3xzA4X{H3i)@y_k}84aze`}%XbG$tDS@?LgGm(6@`e-d>WRhb1NG;7<}X9v*5=mcVcL6pcD8_QQ?DKFc3J zqWJb4*j4N^%ZZ)$E&HpU@~}{pJnY*?9&Ujjd@KZoe2}0dcrf2%Y;&WJ>*%Mv*yrGP zbwLlYFKP;Ai+$l#vCqfX=)6{xT@(B4Ve~=08Qa7@?SKOPOGauY0-3RQK`v=;{{e8LRKK^*_@3w$}a({0( zxxfE{+{X#|{e$>=9gj@ogFqC(cj~ylw8h`)nAe00$%CKCr{Q5SeADL={KBrtkKfnP zge%3q3}35@@gcjQli24LLl!^y%pvwE_)wi(OYD=*0(hHT98ho4cKEMGKjZK%J5cOH z@1WbFV(){m+aAc8k_s{)R8r_FoC%s{?*?JEX|{8QJCjtS)kYHoj!LWB_>o zuT8fR_+*9m9&7R4Y7%?zeR3aPz4l%c@z=`O_4q&?7yfe-@cGKP*2ydO;>ct9OCG$w zBo8j1l?O|cC1Aun38*wf?tRNF_kyy^y`#D0-kLjd@An>ZZ+|Pf_rM?lKA$C^$!`)c z$0-5lO3Q<5mE^&jo?_2bSL}S3+ujkMv%{~WgKrYBb(jR~HA}$3*%EN1j|3c>Dgnno zNC11BfD>sT7EB=TZ_1KD2ap|9mVl!Pv{@|y|Bs~Wj%#XZ+KCOZ02(wB2qB~g0`}g! zNYiWY?b>_qiWThT+SeX?_i8euSSczNfeC~eygO4~Y|(!>->xmieQn^Gxl<5o)B0Gai$fpg;1HUJw-1Itk#L1~+( zL4FgZ+zF=o$@O-+)fgHwg7g$gO}p(hd1Q7vL7Io5P0$$a#UwEO43(ZVqr< z48RX-3tVqN`8}WnrR|hJcOK+-Lhl;reF{KMJi;~n+jGEIW>GIqLbooGQx_oOLUH(vkQa?K> zO%DGWAnv`wDQyz?EUtm^YNPZy!IZve0i`dKQ2MHHN=X~>-UjiOWs!@8|ImunFC3{3v~`9o%kE`noWrPk}}#!@9CY1DPC3Uzq@Iv1kL| zmUTt@eLxWGPlNp#&@$=;6KmY`XiA@^!W{n_bA1BlBo_0_nmSN}c)_r4DUJ z>A3f$R~bR+#a>csg-4X~tu3YG|DY7?@tbn!GNo)NqLf6G&q<+_S-mL*&mAetqbOx- zXN;{1l#*hil-$9T@^%!ZmTW@l((jZWG!gq>t0^7(Rji-!y zS1F?r?yF7nDWe7MU9Dc@o|;1$?Z#0?6zu2fcrQR=Dr~L|Pmv?(>cB1EH`xoiZBU1G?Z|eh27B8TBDo z-v+=|{hKIj1bl?dL&z7ReK_>s?~7)8U;~_5T!YO+uniuqptm)?P-wLoWw;l&0Iz1? z(sTsc{bvt+8a@S}qrpi4_u%@!fHsJO34d!qwss)y$!#g4dMsiD+Yy5)!;%Q#J{{2> zsD}I%0I>@P_i#P_W`Uf!BW0KtP=*nFj0(yytOTl1h8}Ho`1?njfWHYI;cp%<$_So} zzZ0NKo{lv?pRQKJ{n3Q`s|xoif7}B^16bme*$(ob{+&)ipX%F&g3@9D~OlzX1jwFme+;PnK$_7!#{_*&T2f}|^ORYpHf7@X^-Sz3nCXM(7~fNr>Gy~-E8uyl;ypZ*9iUAA zR={E4Hh|{~*>gaSdOUAcYyjEsz$U;KZTmsT7U;@^&ffr@Ysx(Z?gDrY`+qo*?gyp- zcqS@`XQKZ{9MYKp0Yzww=fd(kfj7_>28;xb0^eYxDS&4yUvTuR0t6s`9Z;Yi`SS3! z+$8w+1Xu*X7XCMvSr$CX{s8c7S+*nmpG}$G4B%O=^jyRRcD+^u&4FIPLzLs$u(TVv z4}5^E4nQ2e{cyb&@P^!c$}IH<=!y6uMkOxd88VeJi*2M#X=6Mq9;REr-_R}W<9zGK z1iDrDgl>JCOt-$^@BL4w>DEU)i@e7(?c33KHpa8T%QU+6!WYlTcs_UvxhLm<>A-nh zuSGf9zrf%6uTcNSO1IwO`eP!VrImE+>rv=`20Q7nhv(%#(Ud8HFU1y-)9(g3E0h3g z1H%CU+$3j32~Z8_3oHYU0{G3KB7Td&9%asouL1mSRPhKf1!w>iL)&LS9&jB<2KE6f zfiXaPAQXT-zpsD?!0#FrV5>qIKn4T=_+6@k5>NpCfH&ldfDZt4`<()kfT;j}-@yOy zoPHUQYXlqz&H>NCM+!h^MYON*9>8xR6)pjHQKkV#0?;AL3&$bwI<^G>op5{#6OPZ|^raYJ6^^fE0O*yx{m-{d_|+kiW1< zIL`GJ&PGFov+io)tg%=)Bh12STqT^LYlJhXwQyG2ESx?$!s&HXI3;S4UD!Zme<&@o zpRW?x_aBPv96ynL(nIe1k zb&);eq{yDuPh?Md1w;W;ME2BCB73?;WY63ovgcM5*$Zlm>?ALdy#l)a=>lC&k-hVR z$UXpj$8$tB`-tpIn?!cTB9Wb)Ewb;y#`BRP`{PcL{VNW9HwdS!gm4DW5Khxa;jGnF zIGfHA&i3%N^Go6ETO^#rRto3%g~B=gqHxZeDV$00d1ab#t}7{=o5~8O@D|Q(WrTBk z1ke}x6ye+oJzFA$bK`Kx1q^NOG zV>bzB)C%Ej-A6c^6^h)$vqbLBnj&|-M&vH-CUO&ch}NpW07lUD{_N2iQI}3k?S>IS5;T!8c&Pd>Jz|chRAIbByu}~W1sUPcVsD%JAJ>%P5LQv zH#tS_fwv;}Ohu8K;U#iilSS@}Mk4pCO5_#G6?s10L|&EMB2WENHhed9# zROC9qHRF%Sz2Yr$Q+J8nl-EezxJIn6gJ%}F$<_Nbm zShzmRh3jQi;d<~uxN=Fj90P>w+DhR{sViLPMhh4HC0wT(2-k`6!gahQa9p@f?h~%l zH-w9mAk$R1E`|x$rQHB@oganrGs5-n4BgtF$sbuYrV)&t1D8Y5hv76{jmt->u4!d_uh-beI`S=Z?+Qddp6;Iy+F9X4j1mC79zh?OOao$fyl40 zSL9dPE%K}M7Wn~BM82X4(heeDZV>shX(Hcew#YACP2@}d74E{7!u&ZreWLJ~2qRj~o&1LtB8s$ZLiBcsA;KK<1=yr)LRw)@bmlBiyeR zLuZD__Zkci9+9tGCGu+y5c$oEiTv15k>3s1{e4CLkTD{E#D0-K>Ku@RJYb_ODI@yG)1Ijw2 z3s;w(!qrbBTw}b2YhEki+E@quI8?YU*T=Z|DqMMeFox8^b>oL{ohcBmT~CB-F~-Tr zr^3}aPPkfNjMPK_)>$E3jY5R0?L6V?I$5}eOauon;Y#Q(T!|QOOD=$0yl^e-D_nE0 z2-hs|ocR*@n9wppS=$k%C!aZh@aL)=7?!_a7 zdqq{@Ub9=c*Bun@4VQ#_Q;cwLwhH%_*TTJ3CEVML0QzaW6e;?9E5@G~3Yo*gy)jw1 z*F)z&r=a5m##y*