From d3553c2e223427cdf76b729f004f68c8c02cd72c Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 24 Feb 2026 13:35:16 +0100 Subject: [PATCH 01/13] turn independent_column into a template argument, makes the kernel about 4% faster --- include_rt_kernels/raytracer_kernels_lw.h | 2 +- include_rt_kernels/raytracer_kernels_sw.h | 3 +- src_cuda_rt/raytracer_lw.cu | 53 ++++++++++----- src_cuda_rt/raytracer_sw.cu | 74 ++++++++++++++------- src_kernels_cuda_rt/raytracer_kernels_lw.cu | 33 ++++++++- src_kernels_cuda_rt/raytracer_kernels_sw.cu | 15 ++++- 6 files changed, 133 insertions(+), 47 deletions(-) diff --git a/include_rt_kernels/raytracer_kernels_lw.h b/include_rt_kernels/raytracer_kernels_lw.h index cc34dbc5..b25a51c5 100644 --- a/include_rt_kernels/raytracer_kernels_lw.h +++ b/include_rt_kernels/raytracer_kernels_lw.h @@ -18,10 +18,10 @@ constexpr int rt_lw_kernel_grid = 256; constexpr Float k_null_gas_min = Float(1.e-3); +template __global__ void ray_tracer_lw_kernel( const Int rng_offset, - const bool independent_column, const Int photons_to_shoot, const double* __restrict__ alias_prob, const int* __restrict__ alias_idx, diff --git a/include_rt_kernels/raytracer_kernels_sw.h b/include_rt_kernels/raytracer_kernels_sw.h index b1a88574..21243364 100644 --- a/include_rt_kernels/raytracer_kernels_sw.h +++ b/include_rt_kernels/raytracer_kernels_sw.h @@ -18,9 +18,8 @@ constexpr int rt_kernel_grid = 256; constexpr Float k_null_gas_min = Float(1.e-3); -__global__ +template __global__ void ray_tracer_kernel( - const bool independent_column, const Int photons_to_shoot, const Int qrng_grid_x, const Int qrng_grid_y, diff --git a/src_cuda_rt/raytracer_lw.cu b/src_cuda_rt/raytracer_lw.cu index 977031d1..c8346fca 100644 --- a/src_cuda_rt/raytracer_lw.cu +++ b/src_cuda_rt/raytracer_lw.cu @@ -317,23 +317,44 @@ void Raytracer_lw::trace_rays( const Int rng_offset = igpt*rt_lw_kernel_grid*rt_lw_kernel_block; - ray_tracer_lw_kernel<<>>( - rng_offset, - switch_independent_column, - photons_per_thread, - alias_prob.ptr(), - alias_idx.ptr(), - n_alias, - k_null_grid.ptr(), - tod_dn_count.ptr(), - tod_up_count.ptr(), - surface_dn_count.ptr(), - surface_up_count.ptr(), - atmos_count.ptr(), - k_ext.ptr(), scat_asy.ptr(), - emis_sfc.ptr(), - grid_size, grid_d, grid_cells, kn_grid); + if (switch_independent_column) + { + ray_tracer_lw_kernel<<>>( + rng_offset, + photons_per_thread, + alias_prob.ptr(), + alias_idx.ptr(), + n_alias, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_dn_count.ptr(), + surface_up_count.ptr(), + atmos_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + emis_sfc.ptr(), + grid_size, grid_d, grid_cells, kn_grid); + } + else + { + ray_tracer_lw_kernel<<>>( + rng_offset, + photons_per_thread, + alias_prob.ptr(), + alias_idx.ptr(), + n_alias, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_dn_count.ptr(), + surface_up_count.ptr(), + atmos_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + emis_sfc.ptr(), + grid_size, grid_d, grid_cells, kn_grid); + + } const Float power_per_photon = Float(total_power / (photons_per_thread * rt_lw_kernel_grid * rt_lw_kernel_block)); count_to_flux_2d<<>>( diff --git a/src_cuda_rt/raytracer_sw.cu b/src_cuda_rt/raytracer_sw.cu index e241f75f..ab40fadb 100644 --- a/src_cuda_rt/raytracer_sw.cu +++ b/src_cuda_rt/raytracer_sw.cu @@ -375,30 +375,58 @@ void Raytracer::trace_rays( const int mie_table_size = mie_cdf.size(); const int qrng_gpt_offset = (igpt-1) * rt_kernel_grid_size * rt_kernel_block_size * photons_per_thread; - ray_tracer_kernel<<>>( - switch_independent_column, - photons_per_thread, - qrng_grid_x, - qrng_grid_y, - qrng_gpt_offset, - k_null_grid.ptr(), - tod_dn_count.ptr(), - tod_up_count.ptr(), - surface_down_direct_count.ptr(), - surface_down_diffuse_count.ptr(), - surface_up_count.ptr(), - atmos_direct_count.ptr(), - atmos_diffuse_count.ptr(), - k_ext.ptr(), scat_asy.ptr(), - r_eff.ptr(), - tod_inc_direct, - tod_inc_diffuse, - surface_albedo.ptr(), - grid_size, grid_d, grid_cells, kn_grid, - sun_direction, - this->qrng_vectors_gpu, this->qrng_constants_gpu, - mie_cdf.ptr(), mie_ang.ptr(), mie_table_size); + if (switch_independent_column) + { + ray_tracer_kernel<<>>( + photons_per_thread, + qrng_grid_x, + qrng_grid_y, + qrng_gpt_offset, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_down_direct_count.ptr(), + surface_down_diffuse_count.ptr(), + surface_up_count.ptr(), + atmos_direct_count.ptr(), + atmos_diffuse_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + r_eff.ptr(), + tod_inc_direct, + tod_inc_diffuse, + surface_albedo.ptr(), + grid_size, grid_d, grid_cells, kn_grid, + sun_direction, + this->qrng_vectors_gpu, this->qrng_constants_gpu, + mie_cdf.ptr(), mie_ang.ptr(), mie_table_size); + } + else + { + ray_tracer_kernel<<>>( + photons_per_thread, + qrng_grid_x, + qrng_grid_y, + qrng_gpt_offset, + k_null_grid.ptr(), + tod_dn_count.ptr(), + tod_up_count.ptr(), + surface_down_direct_count.ptr(), + surface_down_diffuse_count.ptr(), + surface_up_count.ptr(), + atmos_direct_count.ptr(), + atmos_diffuse_count.ptr(), + k_ext.ptr(), scat_asy.ptr(), + r_eff.ptr(), + tod_inc_direct, + tod_inc_diffuse, + surface_albedo.ptr(), + grid_size, grid_d, grid_cells, kn_grid, + sun_direction, + this->qrng_vectors_gpu, this->qrng_constants_gpu, + mie_cdf.ptr(), mie_ang.ptr(), mie_table_size); + + } // convert counts to fluxes const Float toa_src = tod_inc_direct + tod_inc_diffuse; diff --git a/src_kernels_cuda_rt/raytracer_kernels_lw.cu b/src_kernels_cuda_rt/raytracer_kernels_lw.cu index 7dcd9686..5976d45f 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_lw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_lw.cu @@ -124,10 +124,9 @@ namespace } -__global__ -void ray_tracer_lw_kernel( +template +__global__ void ray_tracer_lw_kernel( const Int rng_offset, - const bool independent_column, const Int photons_to_shoot, const double* __restrict__ alias_prob, const int* __restrict__ alias_idx, @@ -442,3 +441,31 @@ void ray_tracer_lw_kernel( } } } + +template __global__ void ray_tracer_lw_kernel( + const Int, const Int, const double*, const int*, + int, const Float*, Float*, Float*, + Float*, + Float*, + Float*, + const Float*, + const Optics_scat*, + const Float*, + const Vector, + const Vector, + const Vector, + const Vector); + +template __global__ void ray_tracer_lw_kernel( + const Int, const Int, const double*, const int*, + int, const Float*, Float*, Float*, + Float*, + Float*, + Float*, + const Float*, + const Optics_scat*, + const Float*, + const Vector, + const Vector, + const Vector, + const Vector); diff --git a/src_kernels_cuda_rt/raytracer_kernels_sw.cu b/src_kernels_cuda_rt/raytracer_kernels_sw.cu index 18ed7d3a..2853c00f 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_sw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_sw.cu @@ -119,9 +119,8 @@ namespace } -__global__ +template __global__ void ray_tracer_kernel( - const bool independent_column, const Int photons_to_shoot, const Int qrng_grid_x, const Int qrng_grid_y, @@ -458,3 +457,15 @@ void ray_tracer_kernel( } } } + +template __global__ void ray_tracer_kernel( + const Int, const Int, const Int, const Int,const Float*, Float*, Float*, Float*, Float*, + Float*, Float*, Float*, const Float*, const Optics_scat*, const Float*, const Float, const Float, + const Float*, const Vector, const Vector, const Vector, const Vector, + const Vector, curandDirectionVectors32_t*, unsigned int*, const Float*, const Float*, const int); + +template __global__ void ray_tracer_kernel( + const Int, const Int, const Int, const Int,const Float*, Float*, Float*, Float*, Float*, + Float*, Float*, Float*, const Float*, const Optics_scat*, const Float*, const Float, const Float, + const Float*, const Vector, const Vector, const Vector, const Vector, + const Vector, curandDirectionVectors32_t*, unsigned int*, const Float*, const Float*, const int); From fdaa9782a9b7d47b18da31f506fe3985cabcda30 Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 24 Feb 2026 13:57:46 +0100 Subject: [PATCH 02/13] input camera settings through .ini file as well. Add options to provide custom sza and azi via [floats]. Disable writing cam settings to output file, this requires some extra effort to make it work --- src_test/test_rte_rrtmgp_bw.cu | 51 +++++++++++++++++++--------------- 1 file changed, 28 insertions(+), 23 deletions(-) diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index 55020210..c1801434 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -228,24 +228,23 @@ void solve_radiation(int argc, char** argv) const Vector kn_grid = {input_nc.get_variable("ngrid_x"), input_nc.get_variable("ngrid_y"), input_nc.get_variable("ngrid_z")}; - - // Reading camera data - Netcdf_group cam_in = input_nc.get_group("camera-settings"); Camera camera; - camera.fov = cam_in.get_variable("fov"); - camera.cam_type = int(cam_in.get_variable("cam_type")); - camera.position = {cam_in.get_variable("px"), - cam_in.get_variable("py"), - cam_in.get_variable("pz")}; - - camera.nx = int(cam_in.get_variable("nx")); - camera.ny = int(cam_in.get_variable("ny")); + camera.fov = get_ini_value(settings, "bw_camera", "field-of-view", 80.); + camera.cam_type = get_ini_value(settings, "bw_camera", "cam-type", 0); + camera.position = {get_ini_value(settings, "bw_camera", "px", 0.), + get_ini_value(settings, "bw_camera", "py", 0.), + get_ini_value(settings, "bw_camera", "pz", 100.)}; + camera.nx = get_ini_value(settings, "bw_camera", "nx", 0); + camera.ny = get_ini_value(settings, "bw_camera", "ny", 0); + camera.setup_rotation_matrix(get_ini_value(settings, "bw_camera", "yaw", 0.), + get_ini_value(settings, "bw_camera", "pitch", 0.), + get_ini_value(settings, "bw_camera", "roll", 0.)); + camera.setup_normal_camera(camera); + camera.npix = Int(camera.nx * camera.ny); - camera.setup_rotation_matrix(cam_in.get_variable("yaw"), - cam_in.get_variable("pitch"), - cam_in.get_variable("roll")); - camera.setup_normal_camera(camera); + const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); + const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); // Read the atmospheric fields. Array p_lay(input_nc.get_variable("p_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); @@ -582,6 +581,12 @@ void solve_radiation(int argc, char** argv) Array mu0(input_nc.get_variable("mu0", {n_col_y, n_col_x}), {n_col}); Array azi(input_nc.get_variable("azi", {n_col_y, n_col_x}), {n_col}); + if (input_sza < 0) + mu0.fill(cos(input_sza / Float(180.0) * M_PI)); + + if (input_sza < 0) + azi.fill(input_azi / Float(180.0) * M_PI); + Array sfc_alb(input_nc.get_variable("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col}); Array tsi_scaling({n_col}); @@ -877,14 +882,14 @@ void solve_radiation(int argc, char** argv) nc_azi.insert(azi({1})/M_PI * Float(180.), {0}); // camera position and direction - Netcdf_group output_cam = output_nc.add_group("camera-settings"); - - std::string cam_vars[] = {"yaw","pitch","roll","px","py","pz"}; - for (auto &&cam_var : cam_vars) - { - auto nc_cam_out = output_cam.add_variable(cam_var); - nc_cam_out.insert(cam_in.get_variable(cam_var), {0}); - } + // Netcdf_group output_cam = output_nc.add_group("camera-settings"); + + // std::string cam_vars[] = {"yaw","pitch","roll","px","py","pz"}; + // for (auto &&cam_var : cam_vars) + // { + // auto nc_cam_out = output_cam.add_variable(cam_var); + // nc_cam_out.insert(cam_in.get_variable(cam_var), {0}); + // } } Status::print_message("###### Finished RTE+RRTMGP solver ######"); From 5743462fd5475a84fd9ee17b5a089934ec20286c Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 24 Feb 2026 14:00:28 +0100 Subject: [PATCH 03/13] add option to provide sza/azi to forward version and fix typo in previous commit --- src_test/test_rte_rrtmgp_bw.cu | 2 +- src_test/test_rte_rrtmgp_rt.cu | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index c1801434..d712606d 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -584,7 +584,7 @@ void solve_radiation(int argc, char** argv) if (input_sza < 0) mu0.fill(cos(input_sza / Float(180.0) * M_PI)); - if (input_sza < 0) + if (input_azi < 0) azi.fill(input_azi / Float(180.0) * M_PI); Array sfc_alb(input_nc.get_variable("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col}); diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 37d1b98e..b616fc66 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -245,6 +245,9 @@ void solve_radiation(int argc, char** argv) Status::print_message("Longwave: using 2**"+std::to_string(lw_photon_power) + " ("+std::to_string(lw_photon_count) + ") rays per g-point"); + const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); + const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); + ////// READ THE ATMOSPHERIC DATA ////// Status::print_message("Reading atmospheric input data from NetCDF."); @@ -386,6 +389,12 @@ void solve_radiation(int argc, char** argv) mu0 = input_nc.get_variable("mu0", {n_col_y, n_col_x}); azi = input_nc.get_variable("azi", {n_col_y, n_col_x}); + if (input_sza < 0) + mu0.fill(cos(input_sza / Float(180.0) * M_PI)); + + if (input_azi < 0) + azi.fill(input_azi / Float(180.0) * M_PI); + if (switch_tica) { From 9886718cee24f686bdbaa72eba8e3bed83cd8782 Mon Sep 17 00:00:00 2001 From: Menno Date: Mon, 2 Mar 2026 08:59:29 +0100 Subject: [PATCH 04/13] remove old aerosol_optics file --- data/aerosol_optics.nc | Bin 111618 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 data/aerosol_optics.nc diff --git a/data/aerosol_optics.nc b/data/aerosol_optics.nc deleted file mode 100644 index 719889f3652c0b987ba1a08d3261a8d04f0e89e2..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 111618 zcmeEv1zc6j_WvdXq)`+^&?7CjaX=pI2jK2I<`1Q3W zDX&$MgM??{eY$Y-2iNG3fH-BHgBr@_as(CI8oYkdikk*8iYr@l)l_EW1vN0`Sz#Cz zO~$>h0%L#8^-;~0DbY~8_#3~~EEz`qZ@mBQm72-l##C$ax8KVs{aqTo4E_;r)QtWH zo?45)@$$F&e0LY}kCgcx{$IC$N#HLD{3U_EB=DC6{*pi~2}l#`Uk{24iD6{Uhw}H% zlHeKLuA|6OPvKKB-*B{RNYz}0*R=}x-6QV?$vh!>QutuDyo=;379f|$&F$xx5E8(1 zq|7h|JUH7nyN0+XBa{-+<>p2?R$UpEjsoi&Cj!>3UpjRTisxxD9~4C~DvDyhU27hB zmSULBjLe%;Al7qGym6{Q5U%Z}OP~5BL;S3|Y6u4b=U=n4xP*50bHqei0 zD68lk6?tK?Av64$9GR4p6%J2XLz5t8qs*7_;9V{L((CS&IL4o^4KrCuE}}IFD+`#0 zvMjuFDQkt;zWe$!mHZ7;r6j91bJQt)F4Is3@R9vun65o64lp}pS=^N5P5Pu`#Hl_3 z7eC40{`KqMA^~~(*Zo26KeYW1P4t5j{Ls8VH0=+X@PnTG(40Rg-Mw99Z~iS6`isav zB>}}ZV4{+eu`HXClEPvLUnRm-#ZYES44O5qF?S-yNP{uz#pkG_A$y~&tSrx{&wrrm zGTEyYRAwkHjPSjv&j@z7_Epo$ocI*<>gN)!HQ+G+jD&TV1__GVjjSbMz5hH3ySV;I2`j38v}2Rt zpG#P#IkF^N{}~C_Wrn>|F0-n}Q~Gbv;Xf;3g>7z-dHYhde@M7SJ(;Ge zD2<(}(AaeOGb3-kI)BQ>-yi$n$@i3Vwkjy=*Er|0(c_<|EXzDh?@9abDVOh1P-Yq^ zu78fQEN_;oQ*w~#cVnHa^PkjWnLBTt-Lv7(by#6VWMilz{WDfVg}Kp9VI`i_~Hx+$Nbf zr|363d#TF(JK;&@-!aTe`NyIB<4^YZ&A(Io<2Ff<;VChOybopk&4=ZRiP7R|^OS?~3kni53d(b2^mvc;Cv}v;1wU6I zKN$rf?d+x9OyTb-e9`iv{4yDLGGA1NC1AQ}x8@BA|EMW$@w_Y>@uX!c|YZOU-&tuC}7S@Z|+X0_l#z07OD1v(ag0Jp3K7x=NN1sCdN{ReQHT5zGhwzc2_o#wUR z0=``>xKJ;C6;qxBxe0Xa{{UVyz4*KKa{L2)r&@4f`F8?+0bVm5{0%PfBdJA4fU~vWLb==@ zz)3B*&|dTp;QX|%Jk{pYrxskG<68?Z)XT3Hyj(t>FpPgKxX@k!e*h1x1sCW9)q)HB z1pfg%q!wJ@Gqe_5pc7ULF5vU?;PNEMN8pEdHGY5#^@{uhcvLO8KqtBuT;M0>58$!2 z-~yjpT5tiseJ!|vUjtUG*g;B{3U_EB=DC6{@W$+ z{j8(PeMN11*L>7#^Mlymp7xM=sltrGw`q^U45`#1r*m0;QC693_=7K+^xr42Uo%AzWpg2s0qE?}hK5W#0+9nl_z&o147exn)zu99u^vC0qVA zqoi2a{Inm~vRobJ( z5c#~$(YKwef^VHoBxNKP#g^RKicn z<&>7?=SW4l#igRMyd06gZ)ChMMJJM$_2w zcPuV#=j7t-=;0!EcWPId?I7dd!NpzdC|jnNiL7muf6TY$t7!}i{wj#9Q%|{%NIumf7ujom5<$U6E0mEo(m>Awd=tBy!dC{z?9ya{dSTc_aHfRr)FOe~{nB zSt1HitQN_ACEgBZ6!B)4x7B4irS_ue7+-&p?GNioWumwc;_7VGT+@YDWWR!ptPY~A zjM7Y*!phJ!U3eY*5$spYj}5Z%!7rz@os_rlg)%XCYP#?m^&@6HKRT8w6k3MjrRl;e zMtm7>>c5Q$((-~_-tyPts<5!Qh#&bx!rKfPn>-0C2RCOcO(i}ER@OtaO8FtNxb)k? znrx(&Ew%_f)3+c;%A0YKd;t83iTbhs%MgwZ>p55!Oe*`w$A$64yU6>yh<_(g^_P_@ z8>&k4MZP7aJXa!KbRKRJ{)QiDyp{Nc8d-kcKgcf6Dr;LamdNSImPQ0Rp+WJA0aH|i zA#(RDb>W97Me$k1{6b(+nLOJ_K97)19k&1?Q(bvWzYL``FtbE$39sv^3ddM z@!!gsQB=ljUdbP3=R4ZpO+g`{-}(PCbl;k_9^aZZCEuF0Db9(8D;ln1?=8H}ZhX>!TRA6^ zYcQ=3m*_p7+thF=*C^`(H)UL*_wuotl5WSQNUk*8E?Lm^xn!K913N^$fPLI~4Ey$% z6|DKMSJ|e&*5?*23+IM;_U6pjjO9#*FXqM>pXWa17JDC4(U90!OpvVkxLHDuKayM& zH)oRu=CXGCMzEI$RIx?QSJ+;U4Y>6?Lbv_9z94V z{%Z}y2g7DalHKgsgWFoO#&JW~$sP8vU8+^N&MVEiwSA(w_N)4GWBrzLy`r9Q9a}XZ zRoX3xe^?MnncIV`ZD}Lk^JcB2(PLj$^>Sx++`VaRj`%vecDE77t`T#`_Z4!pHdS&S zANO$eOEkzY&+N$8!~P_=P8M;E8$&u-`-{hzo|25qNoOrOjAl0)u4nUmC~*;WoH_HL zB+jALP%glwiqqU*%~>yOLYl4#BB!Qg5ra?th_%gPQqR9g{LcQXr1geDEK*mn)1O>n zUs{-PR|d7@EV}pQOiLDVCmNpNY|Au>_X8jDU|dJ?D03(|GJF=9+V~iWzdTqxr)fQQ zsoqppy!RlxNcl7S=z$w|v0WM0A+VBLM|W{zpO>8XmL{aar42FYI+R?jnnAqnHxiE> zkBGI>B5}vFuI#kvZS1_+uh~QAOt=KwRBmN?C1+l?o-4lch}*2zfXwgdM}~JBM20O| zOh)zGP0ro9NpeRSP&vM_49nbmPS2eS&t;Hgz7-O@=DyGSL8UDw?2gv_0OtxZ-++2b8{eLv<{ND=)XO!@XM) z)u&HM>F0*jvCxM$4j4>r5H|zEnQZ#oooiThp{qp%T-F;dU9lMwzG{_nQ7P(^R5D6Mz4#2&XPz*>9 z$7z)a919G`v@=1tEbhn6tyiDw`OKp+9gonNN-yZk9)`#nYKad14rpKLjrJkH8}(2a z`$eGVrATP)j>LhX5s-`ug0b3UZtXdHs;0Ju?##YT&-GWuC!;1f*47>uhPvZ;R8Q>II2SWT27*Cx2QmcZHVbhtq_$f z!3pO8oSz$lr$eGJGcp!pX)9cL))EWV5)d8~i=b=an9cZLz)OGP@c0bf=WK!z5q8+{ z(hKjx!!UG03@&=c<5h7YMBYj0JEjc`>nCCP%Q&2#6NSSC!C1Zp_@p+B4E%bJYD(+D z?UpSJ8n|JNmIq4OIpTak3&i>bVL+%K%1?24x{W~dlLwmIb3@uER~+r_j54i}M7!@@ zTB+Xvo()@I;cHJU9KoXXX%DO~bw$bN5M=OIIe-DdC3&c=4>V?j&+;K-b znkb#VP4C(@gr0Q^=qGr?qOT9SjOH+Qgcurw!ol_kMbCHP=+C$J1o&b|Q{cof4x^fR zV?H;Q7{0$n?>jd{hq`v~@bJbi(*T_9>Id4=2VplN(PdEt&dlUxzB3Fh%Y$$vApl8^ zzOcE%B4dvc-BYR$BOP~Kb_~Md_OU38ZiRhraX2N4M&DUUI5#Z;=U1g-R#YU!Tri>? z{L#pf;B~GThNE0)moE)5TPlHe^DwM`n1D6LX((4o#(OS_$4f!do>q8~m5H(Gv3L?5 zf%H|u&^qITp`8e(-wLI(-kM@}F7WbVBtl*%fz!;uz$d8?pJ)wsOe*R%OGdBoY+SF3 z#mL?f=&un1i`hPKSNsY_!YD zz<|l=NEp}#fgkd*^IR;dX%s4p!!UVu0P3nmL7K6Q&gvt^(S$hEf0Y8e<5_rVBdE>+6})mb!|>rY zNSf}9`=K_Vflcw`p$R5cv_yraKSpP>P}|QU<%&Dvo;hRT9ea#CYK6<95tN;yg4AiI zP-$)hn=elA?PQG=UniMja%)5M7$1Yf-3cn&I3ab7E8=5XAV|#y6Q`O(SE(`l zjYiPPF)B!@XbdYaYp8s7#NuKrtdlgxKxYF?){MmT$6}a|vxVz;2dFo)#@)@$(8<>f z2WB?FEd3Gmt&s|rZ#F^4epVRX))Cc=6<%7I;O(81gL{%9KE3bMX2x_z*~ zZ{>&RdZ$BaS?UV@18f5 zbDQJW!@78@XM(`&`jC#)!@|>A__V4HCY>KnuSt|}`?VpwvzubsZCmtnG)I<-Av&y6 zMPes=1g&X`Q>|6-Y^XlsOtet9z6KTr)j`-f29YaaTP3ps8m^P-jy-i2FRD#RdFSMuh6E%7Mf&RMU9bHj1jP4)xoEGSR zrQ;T=$LUV@bbWXUMBbcR1w}RU}rLLkg-NkdKE=>Es0-bWB(v zy;_|}Cq66qkIz5!f-XBwN6Pj0)+NCU&4~7k5R%t%6lpo_FtJQaA=NlVuC-f5#&p%B zwu?Pz9lL0HYh7EanpH|02M?qT);3}h-j92I*^qo`VoiS85KCH39#1k0j+0BQ1p?XZ9S$olP?%tJgY{_-Co)aPCyX z*qW8H3sFNgWLL<6sT`*m2nn|0)OKEI;KYE!NLx;bYuxnq=<(U`Q=rfR-DSQ`cpK z=`v<4?H@mlRyr@G;b9}#-E*IEj!#;Veiq$G{F+Kqrh0_vj8~->4!tC;pV&~bhXY-s z7f1Ko4WO=b$5EHS>6BbqNP9f>mk>ULYl1VAk$eeK-Nd0|piP)(LeN|#X zZO28@HTBczb(?PV?Vf2g{`^wf>B%}e<=9Soxz7n2!WM|fwt6jT_iQNZ)O;sN@1gd|xutw~(w+8b}W8n?!~k*+TxW=Qn)2QM{|# zoxQ$zKRdp>4!1$WjvFzckh@bloin+-jr)4y6&JF}h%7euBd#m^lWv!06P2j-WJSVp z(z*HzG5C5_JiL2LR-}EO4K6m}S`GH+&P*B1HTZQMH=*(}mpW311h#P^%RSqY+kun# zKDLRdJDwnu+dm>hKj_o?r<5fJDmt@9x9f1ySI*qE=V_er&8gfxs}o#h!bk4L`(~tj zMhIDbvIjY|Y#CXff0DS=e@H&M)}d_=nNy!fMv`$KN3tWH8*{OSq1-js?%cGZmE2zA zM_m5_y2NRK7x}m+lXN>ZlJxt$i&PdqAS7Ic?j3GOr<`-APg5Nwv)|2O)5H#(bVeo@ zyJQ46f$ryK98@LUf<>f-aTM83dXkTsGfA6!=ZRQVnX)+sRMXIsIvWPlt}a26^P-Jx z(HcLF9omE2+GGY76nK@p5Nk;KKXxUJw&oC@iV?&;V4lm8l#T2_m(HvZZ(qu=!>^8E1FjxoKdn&VUWR&eW44uW_5(+Atifik^X|J` zyMzYhedkCLuhX4)myIM%dd?x8&mSQF*YlTWZ5L10lduOm9%0YK)a7PBbmTTj%eZjN z;Z7dh!#Pa)#4X)nN^&=akkq&#B)j`UVtaHeDcW(4*eR(}k#)6rb8<4P8vmG${9wk< zcLa0!TZeH=u$im6e}j8(VnFKOawEYFa>%aqX~d-Yb`tXCEb0CA1u1@NNUMD7NS-z9 z!ER}#!CjQNbIgNGPJQ_dZl2{CZn20V#>9%0OpPG1g6n#(7xzYg4X3^KDYxyQJ{dhoLQ*E@lKC-XNOJlPHA97>pt?LH|9G_I=q?3e#&>^Dl&4oCMHjp2-=EzRp?SX+YAd z-APWmSUz*Wv_nl}||3Io8TSyfFF#H~=|E~cs zYigpC(`wd1gpGdxzXM=ucJuRd#VNMz2KjL^vKW(}9w+m@hF2kAgeo@|>^w(XjlU>O@mC}m;QSQ|N(nhZ;i zr;u$+<|tTI`i6B#w|KweI2j?7NghDM(-f`(w~FIrMEs+IMs=0qjSwLrTq~Q$?vx)V zQfGQah&HpB$1-OOyzPvbU-(-i#%p6%Bjy&*k`vr0mX~5swGk12DfwAXVZN{W_7~;otXY~js{8MA6} z38T&vvtqgo53OKU^0)TP!wc6dn2mA@t_$M3FeBx+F6Pa<%D3P)Vb&z1bY=MEK&AzA zZ0p>vOmH2Bv0%!YCXHl@4BEtyT94jyA9a(0Y8E?i>T z%PIIH<}P!KFRcYL;PI$?jEwESAQzfUc&tuU;CQ@F+oML|=(45cv zFD7~W3);``4CTo+sZksLR?!AL--@JYD4q#LQZykiB1KZj?HetH4#~QrJSpla^tYP2 z))n3=8buxlCIcy|FY7Pzq-b<);ZdOuvUO#7QnZZFDQc8ZPk5`4x)AkOqhRZ8hJC{q z5dUvpWuc>rO~T4s6%r6OdHk;h7XJG;>M9nEg^jxD|L!nVnNsup^;_7gLaXoHjW_#E zx-XV1r2AZP`JV{e;w!4u<^`M|{UYPe*;wYXf` zQz%*(Gvr7=7TREvfFz70wJ6Do>(oMl#Q#=!*(rbLC3!Te|RQeEPPRBmLh!bcuNC-}BLoe;yu_lrIm&`QF%Ft=aC7 zvzEu;Q%6f#Hs;#_BD)N)X8n9GfKrH(+$m?QmU;`Sey^6f<_lb+NC-Xg1exBYR5cEgQc*gu6GYKWA6hkR)!c;zlew&gD63lTGJ6 zNd1?+$@mkqN!Q3tgS2WR?o^9{RjLmMsRd={yb2X<~m{qZrTWEp1dVrw>s0f#pyIe(w8oHJ&E>vwx8NRUQ1f1ccl9c527PyuckL_57R|1 zm+9?h^Jx3ZyY%1y1FT)q5R+3y$ZBhiab0XN^pG8<_HaO1H)s5(^B)|&N4)2~7yChH zJFAxPiG6v!2{*o333pmNo2!W1#OdoC;by#2Ar_aKktR#>iDS>P#B0?;qS|&hDP8!2 zG|%WS;V!;omHYJO7Ef5hnZG~BO^-4scGr?gOKCUqOOs(_!{oK(Oyo79Tt`GlZ3(3Y zj(N0che34Xkj3<*+7yy_F_N0rn?uK6UQZvn@1-NOu29jcv(&)*JuR>|z}8I-vGKYs z4(#>B+)RS1e!f_{JrLswAGPx@u9ELdTj9Ms&r_19e_S%AK#vViX4&J}{aAOiWh~=# zoV}~6#tm;_!`YW6b7jvaaXCY`b9pB&aswWJ4V*Ci&;&) zaN~*ueFFICp_0a2_sWNS1fDx#Z&i@K2WIjSlV`{#*Jd=KDvtJ^zMZC*=_C4TbG%Cx zV}E~&vjZf!@~#o*F!lj=Lfe(JOe!ah0%wqW8%~j16PnUNAu&|v))sndjXqYLvO{<= z3)P3d=r$4S{f1c4#s%@JyiOet z#?`a{%pc_q7yI>G(;36aif*TflC=RHx{RYnS9(yx9Xsh>v*)x#hmYzQ?S)<>2-Ene zHCNjZ+%e&>yy^t^+He+We(*U_A7M?eUT;OUzmB0s)|Y7YGZh?~WRBg29I_1fs5L&) ztJ9=l#2#SL_uva|K)20gO^^=N;3NH(WcHxtE=#D*v)6QWjxH)(>@fEiAN0Qyj?r#W za9JCK5M2%~8|#yjfz{+@vIBK5Eu^=9okBPCJxbU1)5NI9##ry=hVC>7Yc56M;M6E^ z%7GZcM<_XUX-*cORi^qRgf1I5n66u}inbBoqt@^1;Zf(N7?rY!q51czkHwM%}g}oW_H~gk^IN6q>bS^MiHvr3w z%`uOU%9&RPRE`Y93x7UR@Rc`4G<3lVsVn?VJdx9w#frNW%`W?(V688Drw77N#{$dR z5H!vV#_;-KNSPCW1JfvevG;)9Yj^Z`=#9LK1nt9p!7cWOpI1e!fhId=I0_X$UXpmiD?+>6Nk>OA-HW6 z0J~RwJ#7NP?h1wGstCRfqcG){D2&aG!QRgC=>IAPb#45y!#f=fR>i`}CIm<32ExfX z2$a^NMmH9H7ezukF9P?TN1;kR2Ad4ypt2?o zk=NpJ>1iUyc4~#ux~-9VBpAsJoZ$Z43G*9z@=Utm=^$fwm24NWP z6^_TQk+^Rhg=vY=c<&H}LvP&pXsaNcxAVq!12L*~eXw081iJ5o(OM@AD|Uy&pL{BWjxf9 z6L9rRItKdq;=toV=spj|xkeFK9TtsMucPthAg>b*<1oiO0ZXQ|MBt`GeAP}uPhOr= zI^)82Pd<544o4Xjk5{Zw74Y>UPyxHqgE(A-5N8ElfjxK<6^Ho7~PG< z{K-;2Dl`^5dE5w50uqf{qEJ5(_QpvVx~3JP9<|2d&B-7ali|KA85P~zW6$MQsG3~{ zqe(4M6qkfyJ}r^>DG{YF5^;WHE3{wJ8l9?>acyY|rtM3?vrZ{kII{z;XtY6rZ8?^! zwn9q(R+wd#2otSVC`nI(l74G+nUak1IVp(oOoe}3Dmr~g!BV|KG=H1{wTS0u|4}OFt5jqROhvch#aOM?4#)F4p|Wur4!=*u zvhZY7j!MR^CCOMlCIwY1Q(-)|4eowvgZW3=;CNyiSVsT@eXU^r#29Jo>Z9Pg9(0>) zz^%Ft#+7SeU=&-XC|^H3ix)OGo~YG963 z9dG(*Lb+Z&m>HO$*HRIBR$F0rk{u>Da)awc4=j;3!l6(D%$ukOr!XA^Y|y~yhILWc zSsRw!4Do2cDF#+pAb+7PURB%Uva$yPnmEI~sy+to)1}FvL=SM(7|>`0~jis zAnm9*9}(Xi=2IMD6YGj_bukX!wSum%COU^`Amf=fzAn^8|Cd@2N`P{WD$b+D^K8>(lu;n`RRInf62b`nl(jlJ6q&G zc7*vBH>^|j!qN;5x9^%j(pd%TnyKPhCZ8`#2W<*;(ITKe%&!^YOKvliSJ>6;Bh@s({A!5kLNi1@Ge<>&6&@#8|HtGz zl5I$@rVXakE=;HQo^PWE8y~0ht#8rmXFpP{o^{|{Y={sib2ucoz#Yy7`deLbd#@{I zw{A@BM~|Q_$IYR1$sX!eeTJ&Px<{w|!eAdCExOaz2un6wV$KE!e9d!%?i6=)N^!^V zUo7Z})#GTJm__tZhr={->1C?C;R&@dQ-$V4ExZ$(pygm2Xa+jt&{_|84)Fwe=K+;t zj&xn#R66Rw3Tp7{DSD*IZMs|kE%oO6{EPm2xPQS6n|9bC&dd#rvp0N3dt<;fFKC2= z7Njnq?5)l8)2nOr?)#^-Nvslb7V1Ft7eh=LY0gLXJ3^(OCyJCpEwwo$FxJM^HdPmpnuLDh;df45p5zcM4hR0f0=uin#`ccFjC0P5A#d@DM^!}(V zbW69RRMg}(-S|=&`$XEQ(8I1sV|-4s;B97eM7nt3-b0GSNdXvU&+lRx6o6TLfB3Y)m}ur5CjBz? z>6*oac2~`!1qDNB-Ht1$-_?V3Sq6iSai;K@=mg~*99o<3J72Uz{_FEO$3>I3UJJ?T_Aa88|R@;SQvwJAgx)hR-4yLqX zV-R(D-=Fr*n?`*WZ=?+d+@g!8exmzk7$da4Gxp4<@XrXw(34^CNeo5Fw1MPIA168` zESYZGI)ETO$u?xllf-=fn*YWz-=dN7P`j!FFYp54_T1lvYnXTJ!{Y8s3YrJIS{#XP!2 zX&CJ?Vm-at=p^mv@RaVZtBtN_jL=}d17gg8{+B{v$@_yPdm^y;Xb=`(J3_{;??#8b z9Z&6=?xuZST%qL72dY+P0RN+=(5i63@hV@WRfJ*f$4IDch`jEJ+&6Sa&+ zboPu8n)h1bUZC)E8ey0sqpLXGQC3!-viT8hfzRgjIQ&EAAU@)2I(X1i#6&xd0<~Jf}X>DpsoepN0O4_B_gVv z98MDt4Wh(uCS7W|g)XnIrZ=8_qM8v#`1;fyx3nba*3B1&1O3tFs1FWy9zqtoyHVTb zX*BWcIGS0yjOsfbqJ=&$XhOd_SitZ4==<6QDUB#@bO_*gI0RzeV_&SiKAB8?7)W)y zcBF?J%%bCxH__qC&Qa-_uk`F!E!?!WKwP*d+D7?e=K3JKZxD=9Pd`j}wuEe6+KP7X z+?Ph`FQZzo_R*ewt7(<78s8?NO7lq*Dn_#SU^TV8xTM6~1kX&7rVa9EgXL%>jEC%zf#G%jewDZlo)m z7$?MT;P(slZ-&7qEn${IaKwTi2@67DA&!FTzBv4h&yUU2vEXt8kc4JfwZs`sLOsy; zS37(d(+q1{I%D+|ub=XHezz0idCzX5zXYvYxWHw;4stWBvDU!{SFQwsDUA9lpYA7} zu)JP##PSw6>#8Njeqv$v(ih9VhT*zx67DTeg@&sOX?I$Mo_-ZdM=TseU!U4W6DEG9 z3-&fb(J2d@aC3#P*b`+^2Q-;)iy(UrAz^K>+vFt)KD2?}_o{;%CRR|gbH!_Z2FA3l zKTL}GJ)gHdz$O^L{-`}HXM`hpVJn2^6=Q1~2e@{OgZ-gYOebmh@-P*XmM5S)zgKB0 zZ#7^0wm|b{e%Q7-94bhL)ssB9-|zURzj_Qt7jh{!4kR_JFB#ln6){_Vj$B@?M3-1t zPzNnB9Z}kW+9!;s4acpbopO%T15fT zlG$8p^k%O6mm%bSpJeX+)g12W@}Xp7(|u&D-FrGhI}h_sI--7NIo$M$5Z>gjc-xhq8HO6Q!LQ*r5)72bbyCh52QWn z0iY*-z1R)^>HNQY+0pym>p>E=k)f>ms-bLP^GU4tfmm*R=m2g9vyyA0+Kfv+G@W}n zzZI7><}1-(l||DJccR%{2GLnbEvbdYG_k{_W$f#=-MHiTmvCvTHgRJ+OUZ3XCGoP^ zNyZHIAo;x(5_HWb;!oRX3tx5Qx7R{HPkk)y@`Q?~_U9Z;m1xnP_VlE2ANpu#Ax*uv zoeti&m}1ir>U6h1edD#9n%WMaPPbX~njC}8`O!%49fd1XBrv$@O+1FG(vOigG-Trw zl61s}Usae&jb?SDL-&d4`nWpuRNn%+Y`q6Hd~%&4f+( z@zOhCQ~ z6AnC)s4H8u{SJn))z7=J6O~r8Eqi`)pFfm%UGp*bHcjd8E#0)w`%=24IBb3=@v8T8 z#KG#%#OI@IBx!nGC0Z88B}Lcu*sgP(*y{Sp?9-={*>|GO-dg&Wk`ZQ_oUO7eCp~LS zj#mee2AXHN>Ptsi<%%4xaKL7EWWivv`Pxp}c>711D}6#^)?B8mQXS+cpC3JPgd9B* zL8nZfO!<2m3{xgIf}5KgcsPIG8o9Z-{IisQ{>O3o01EyDb&)T_tnlam7^gqe+bn<& z(B*$3`JZzBXHWpsdt3msb2)#%iT}CD|0oACv$O-5orZx-Z}Y%^DX%D;&(D`7qa(lX z8py`dR$M)?nl`tFx0M?P`$;_0J{|HNeD99aP8UXbvB&m z0>Z62)tT$iOoz;kILa-*;7^C-&7tk)o}nlGbwGZtj}KNYJTrY)d7G}AEpAXbT>Q>q zr�YOUdg;*V*c%`kd}+j$3cEg&V2Sn>*DX2#X$r(-a zUel*UY~CwOY~oQS#>1hKKKm!IAIfg9`*rGb0Url&Q!hAhtxlXK3z{{h?&|I6(Nn+B zbG|2O&+H4H&C4Uj$9=UV6ZXFn$CcKRkX|n&+T#qlzFKA6h+~Vn1Frf+rD`=tPDaqi zX2a>-7CY(G3r1LKVS>rdpFQIy=7`6X=}K;G(vYZzXi1)tbrSXIui3)fRBoQbQmzzQ zWWU8cuHVG~y2rI2owjH_^$%!>J}nzTHMWtLC}W8D^nM%3HFYbA2dpLO!}m%&cc^m@ z6FPGvX71<4U1~~BuG!2zIgvsiG#E`6HrP+c?`wiyH=04oi+b5aEE4x#7AASpJwoDL zKSDy(FG(hwHsVgt8p)}5ddLkaV2RV`3*3q|o#{2rS+w5h3)I8P7VTYHfUM~5Rj~iE zxb=X(l4Rw<5~oK4Bm?X|OJ3}7=g!Yu%1w@_PbPZgk+=X=;&*T~O^jPl2eo`orNi8C zFv}CYt*3bPn64^m)peR=toj0pn%P{5xw$T@ofyHjy|ah=s_#Ii-W*6ApBs})MziV8 z;sf--AtkIOED|>pe13Y|OS;QhVj6l}BEEZFQaSLlWd0K;)@@M_PXFmE&MhL1#LQYk z8ou`-Gp}x?Gnw0T6!al+2*Uag!ElyHyq9K7l`xCl*`R^J?1e% zkv}*cH+r_h%!my1yjTixdOl`t?S|>cd!UO;PkeOmi5#V#SboikRZp5j3SIiq zbBPP-w#hAUH838H8pmSBq$C*8d|cU)gc+m@_C|C=T5NZqV-FnJ*&S6go3R(GN0ET; z!8Ce|HGNkdOdtQE0(}=gPNJHhYiSjU$*OI!Yepb0cIbqseLCZbZ&xgw+8tZUyJ6y7 zL-ykMUL@*NV|q;Vp0sYGPH!~WMayqF!7_ovnsI)p@00+pu^WOsq*#<(j;hxcn7g7I zD$Ke8m&ps-z=iG8HDyZ=tf`MJwRuKd!e^x5TYJOVDQT%R6KM+pIfCUn=eKC zR-I6KsVn?=I)n2cu+uNrqupZDX;`9!POG%0%K4M&vG>+UXvbko?_l)n9gENwR#5gU zf$zK$Z22Tb_1P|XuH6Y$PCDGqDoeV$ZwGotJBa?0k)s$@F;Wy4cR`DK<%oFVz;zldrVZaz(1yMFqnKC;6|EajCoHo<-e3-W z5sHN!VxhE#pW$d<%-5+H1ItTs;z?&z?v>(wSvIF;$A6)`?MH8{&Y^=%Lg|8laWsZs z>+E}o!^D^{G%St91G8p0v$6=2wiM!4^NyHh(HW=qb;K>rpdsjq>@|luO~cVAJ{BXIi16IB2yLngFt}L>wv|<&qB(#4yq@ds+n4U1FqTFi z>`OCu=F(BuCQ$Pe*05wL$~s0Mq9PWXNmHzq6yhn<0f*WYVORh!y9tHRur(zeZcL?f zH_oL_SI5)J=l$u^-E-*Ps1{gHebJtc!u!K9aKCAavlH4wTEL&uvnU_$yrl3uork)v zm89*mCp68@9u1vhV0AwWyL^g}CN6;8!8Ej)Rm7io)CsZHxv**86_YP@LqpTZ{PXpd*tQ*dNYE^daj!OT1PI8>(u<@{QH&Zr6$t9FL( z_09-r-4$=bD-hGU7x}XLBn{N5gF&%xsLdS}Oqm&s_A}B^?!&Kn>*t`OdV9FP2t%g^ z<>)o84DCObqq}7njO6d52Nn>&Cwu7K{0DTU|1WgN;bXLOAwTQn-5RY^TftR79nYI) zBla3_XhkV<^E%>AiWC>>c1C@ZQht8(1^Hrgjjn0_fWEtNfOEerTie5NO){2LThQ%euG8{C zSLuM*qx9zWU9@se9ZZh$gSk-%>cm81=&MkKbT)=h$21t8YKhBjT4UgkO_b zYMfI|_l8}kBXv&E&tAgl5BU zH6LSR8$kzNyF(2R-JmBU&(ajP!*q4euhcA?L-1xl=tc$MdJ7+{&eKD|rvwx`hvCN6 zD41T%z|oPhkZ#PO(OZ9|C5x(Q>#$37t?>ye3HnUa(!{u*091|k#gSh8{RmBz7DXYz z&L0QeLUH9>8rtuU;$s$u(ZLO#QD@Iz=}7+UC+!QDXlmYPdcVILR&Vfvccv#C?$}|> z6J;!#;EMrEys>sTA4@hS8RJ(4qf@I)azEXQ9$uYC-3MAx5;c)7nX;MM4;@M$nZ;12 zZhG|1&F*x1N-nKDbCWJ@`kJnKp@bf7b#Zm0A&Py%$h`aYsCwG~dc2tyeK)3rIz}&{ zRW19`);;`bz$WMkSrVl{a&vETlT_>3r{(T%g{uloGSC=Wl9wR6&_XO13ywfPz~c>R*DPOO7P zj(luJi3r<1)*)``>&gEk?>*q6c$P(BqKITg34$UZA}C@YExS`YT|JASV#I(LNCJWc z0f`b611gFbKtMrJFeemI5e%47F@cI%5Hl)@Vt{YtaRnAM-_ln}OLk(3RQItjD|)G>cx1o4~L( zO``0&hfLoGpO_o>mB4guZ&0l@gCfg1;e~yxMYD&k6y5fiDq2yyL}XZF!SuO9nQg}1 z87K4iqMhbr7;WFjqPK;Qm`uYKrgN+^gtr+(6F=9-*c}#X&7UK3JTy%-xx`8|`hctG z@qQI%x8pG8k*hj0>C-DwS-dmz-MLznx#byiUFAEoYL^-;mss%otOf%&GE-P@&xxG7 z4-;u9YlnK60?lo{FEFGNY_Ic8k;B@r|BC1ZJC282mnp>@e%XwI_( zr{@Vm8^6AyyNk3%^70>rv7wzr!4uDjiaP5vaoL@i+4xp8*pFib6V8jO3!4}lK3Zg` zl@7c%v2VU0M=ah=>QVZf!mjx#68h*r+u!~9S>&z#Kp$dnCNfLHEG%)m~qcB?HvJBSA8 zJEpvQFGw$(>Ub$QK$xk2O!zf@m1D1>3gN_npdZpdet$aT}iavDvfyRe6{2_LGF6S;?N`q4EqG>R>bPTnC&$qz#yKj7TuWqOdWNx2y~R8aQUdqRJ)!jYe#W`l#@=DbaDg5lVMM|jIyQaxa+DmOAYA0~ zO4z*Rprf~mk!Z84ws3u5U(q7PTTEWqQ${~Q6D0X&z*ST+8J(T%?-=?EW~i9ZNBqGnn!kpOtqL~hNnA`GCnYFjO!JyxVp41@Wu8S1e%3?om+q%nlZ=3$>4vIV1=-+2B~Z>8s=7jrKyBiRr5IoCZ91Fgo@I&W9;U)`hQAFVkQJMB-VMM(*GyO|6<9V!vnV$QR z(K1knU5e}NK3yqraCa97+H4A( zN6ZXi%f!pVq*Adczi_o^=k6auzo0y($$BdjQ`E#v>SqKKdh&bO`m5Wg+o}jWnl}pi z3J(iR&5sK51I3R1N>M_OH&2A&Lnn$xHE$NFX{(6Lh5MLCdIuTzTp7sFGY171dwA<) zX>ZAIeD{t!CMdXnOOQ0YMsRn8tK)3NrNYnFtwK}@7e$OH6`ft6FRF1bXIf%UF|MVZ z;n`wqxWGp~Y#d!~e|fUHW9#z|j*-2naDn$~VYv96u(4o}Xv!2HQAqa!(NecM(USyy z#@cWela#fZNw6wpUVJ*iy!_h0{7fE}xupU2NJ_L`p zFERuxM=x;v+9z4)&`=|6Q51rp|ER%K-Eu3+%>! zKkAStJ6

F;TF>VXYwYvY+Fm`SXQu`rZ?EaY9i|PKGFt%80x!tYfqi^BJqP&zSw6 zwV-#BDV)i;Xcypd$HD(#gy8LobU{u|p5Vs#>5fVh<_lkpcq+WoI8x-@o-2CJ>m?(X zZeh%F_b{7Jeq=&i`oObxOW6PZvt3={PX}waRe~gqBEfvm0|GOPX^ttODZ(CYZ-ni& z(?z|07K&~cYKWZP?PL7v4l%_Fvhd!>9A1{$!4xNL`+o7-g4UNsf*#LL32f^w2t0~> z99`zG6zc7g5e*p>CK`49km%m7o}w3FCzy9Yh>R$#SE1Y$bao~jQNgg2iN z44&~=@G7TXkTovY@vdC9P}fjZbeNCYSLAm|6i54uHjTQ>WZb&KOjpr>3Eyo%xt`xw zGMtW+Q}g~0EG-=_u77K3>3IIHu8q6?4(ss0;qv-#v9u`gE5rZyEG=IDk)=hz7q4TH z{?FTc{u6O0{>`QqgWu)q@EHI9vzI=RmH$Qj2zz~A3Tgct=CiCnW4&Dh9%_LwqaX^L4#dN*aS6blNQAe)>Av!*FAr*MiQv^M zRp`8TA2U=-fr&Xy;j~TwXn%@=wzGVG;R%p@CK2BKrfcQXL>}hybf&DW&}L97V?AV- zs8baIS~Cj@eWD?zJ|0GHOn?cO6QSuh-FH6C=AoLW6YCxm)^&T#Tzk@$`MK5`+#Mod zu*`h8w=MxTcqPKY=ZVnro9+jn=JAlr)7hQdB^0PXVD^laV;<=Gfz6``$T~3}lIjv* z%KAj!aVQa_^^%&-fAZ-`9<+I|cz#AWJoX(kKuMh$K5;f&?>Y~<=p?|d_5>KWBN5Ji zOay7YI_N6-&n7%1kE$1ji)0}z-GGswF$Z1^je^0m5}>tDBGevBgn6G6;oEPzHa@+` zgB=eA@+zX+-_*eWjWyFHkzeGS6%EEE2@pO!5pGr`!j!LxAmxWNKPg=`4{kh6_+}(( zxYz>*C5f4xW9C9{ofy!3mH_6{6Cu1N5q`G*HC+u)6U4(37$Mrmd$rBx9WO&;BS80C z3?+&=b5l95ey$C zf>gGo`8D!zg{ND~^K?Gb7iONH0yRrJft?fIj+Su{eIo%LUPyqSYxp+c+qi@76;Ef( z_t|isH#ZdfLTDrJ#5PC)3`=4l`dBO&olSrz#}nXdP9m)1`@6K9()^_3+BrVI@x095 zYBGV7UPc&T1iGgGMxUzN)--ercyf=J*`b{StXWsGoRq#C3unvP^?$N+0 zP6H3CU_Mm-JkW}Z0pa0Tc#s?i3jCOMlgE|Hb2Fc|@cA1a;N`Q67sNV+0B>r6y4xWz z_vt*485|2SF!sMz6aALk6sAP2qhwLyoZx{s| z-o!%f%Q%pK91ndO6F|u0OUGVY9s=L-ZF)WtcAOp#PfBM(M4>7q1%-nFKgJvK2B>jy z@z6PKK3p(P1QCz_o(E|^M}C|-eT$dv%{H*#K?0O_Ci-d+0$=!Xu$y-r>^;mIAtuj< zPA0$Rgbq3#o^IFeL|Bz>0WaGHFx>Dv10^2@sIQkG(v;w48@|xWvn5AD*W#Ch5ZVumN!5 zZX2UFBp4Em!l3CW-#ZeBJ^n)y0YWq9e|8^&bH!lI%8$lVYG)A{X@ zTbIQF<@=%9=|m9l_|kTlj`1qITzBQ`I4LciSvyW0)Y6-oz$4vZWGCw3e;mflAu!462shMg6 zoGyU;!Lg7x+#8Hr{lFwC0QwK|g?@G;VXQ$Q!@bl1qcGk%d_`~Y7-0q7XAOd=iw=JWE?~lO@VO} zeSupZ0=v7<1#)LDRER=h6`cvQi9AztNDsILQ<=hz2n%$_!Pw9VP(|ILzhEYu{uuzi z9&;hyKLU<6MZm2B-mtpf2gIXVMQtg)K$JR}scswrvJj^Ojkn*$%R zA|SwV9_X)|2kCJUaQuQFgiLP|^*>?EI~=$%$F0W0z~C9M+hlyf5>S z_fE-t90kd}`MGR5KPN9qgll(sJbn(6^6+oz`kWXm3e#E0Z@=xz{0NPPt*Oy4!#feG z79~PmEk8H^rjzO+Qjq5Nfd~5uexg~~+05*>T1>ia3?%l80Z|a2$KphI&Fd2P|AtOl zKPlaS$6g}qr$-q1wx1#Jw6s=J&@O^LUQv(&T#P zsDGnKMJ^Vyc$tz9OoY+%c-~y$b%)<{KY4s9Na>{WwP(e8(T;pw7{`7Pc^bxn_o`S} z7t8B6alB5%(=Fubq-B%p$WoAgmg)zoZ|g;O2OneX6lGw(uqRyC901m?Lg=gJ3Tj>x zq4AR!d}^Kz+|^*-!8{aF4urxQ67t`ku3(wBD7Wb@v!Ex(^f&~(>&7T}&pW*f+`Yle z%pdA{g@8wVIC#ZILN>p7VZif9nCTP&Aq(9_dndeR5;vJLd#gskiE9$l4Yl+M%ihwymnnnt{_nzS0Mr6+%FEb2x0w ziG*0+XmD4H1)1Zq@a<$QR1S!Rv+JW_gVKCa{3T7;v8XdMud^=<+7Jj2`h~)T84>(O zkZ9P>>%b@U;vp|8o_9`-2fgQUP_ro(t_)ow%9!09q82MLtlMmOAqFE9StiF z#K3O;9cGhdJZQ{|he2cF;pm|_kb_w0RlG*D+13y&A9QAp-VFrbun>@B&x6(SF%Ztn zxiB~m${pk3(_CK1<8_(=iE&_G8v}DVOQy2r6QfdY2#btu;f*2%`TQ|Zx^x;$kra08k2> z0~goLg`gdgFoBnm8(pG6*>)aW{tynIoYF*?$s0gb#&&K(d9y@y1@diN;!Y#a#%mr_MC-d(_Nu_N;`b~?ySm<893g8(DL;O6mg zi2cHE9D5uMPgci(0zck;+!YNAvgX0C`71QHCAB%)?{$1WD zEe;r;cxWh$2miQuVEFeni-=fQS+zlQZ+=gRx!jAHR1pNB2BDxN8wL5JV<6Z$7G@U2 zK|jBEsLhFoEz{${Y*!p?=@tu96t{`)35>zk%YcbLGY1skhr*dTQBb)q1}4PBLfLU% z{sZG7Y*jpL67%oB;c>8YU(6p)gTGd&+JF99dg@qURsFZR68@W`?b+~kkuDr|EH~Qz zQ8&+m^re3t4*y}t+Vu+mp^^!ua#H$|2UTA_W%8b{V~!K|FZ#MrLGdvt0R>% z<-hBq{=d1GneeRp%U)(G<*Q6)b(Jvz#$O9Ma$N6zu`ojX`s>5m-|Rx#ibmS%qJ2PJe|kghv?sJJi3dg8nZ|8NEMoT^fx2FT<)Q z{qEti3dX;2HKl>sJ5Z#|865U|uKCp)0{;Vha!339t%0xO`M=sbrnLi|IPg`wQcMhrK7*Z z&hQ@`GS$%+J4T_7`+rr;R%ZU<8~dAm9b)#saOoSvSKbZI~{&;>)-A0xUpGU*Z(38kEwK;^dN{oO&oueTt?=R?IgKBzb8y-`}MeE zSW$KT#rH1!NnM>~f5UI2#9XlW%%=qK+pQHB=`)5xhncpyG>aXA8|L>U8E_#$S2PCZ_ zAmo>HCcKS@WZq}YV{QVp6(+#>hY7I6Isv-vjRPc+;atEIieZgn;-^-Y+ExJRkal*wY3It$AOr%4puS z5)T#K=R@stepyH(9x_M7LGbz*_<24O+HVKKLf&8H2_KE*wA>{4d*zh*G#?i5X2+|% ziTWLH`k$l7n>7>X!?3SWpq?2Bnm>I0ZaGI3C%}=g>r z+TCG;rQ6>vr+feT5W77ds(7=N=g>H)w}^!^G#p+X@CIc**4M_)j__GnEG!v^%=2-J zm}+Z#h!30rvz`V++cbW=2cO6HhFG{@7Y846BcYM^^D$iN43_0;(3E*tl#_Lz@f-1t z2^{1B4KqWbZV+!W5X3_1a-JqL8omsg3$0Im!GH5;Xlk#Kkekq$@dN~pt2Al8t!NW*DNbBncq04Q-kUyE-?Um9$ZF)LmoFdq9e~Y6C zdI*EG!i9LVNEoZ}PWa)FP?TyoPh`F5im02D2@~9ZCc}oNF|OaYGNOQEOx&0P!D;zR zjxk0K!nW^ggh@)5gq)F@=+!89(WGz7MJ{h2i%dc+neeebOukkcGyc#PrvIH1=FZhL z!OEh&j;1O`!kj@#!c$@UgrmE*2v;5l(KDAQk>37`qN_&?n4=RWF!TA-aQBVKVJa;4 zGykl~BP6mXyDRy<*x>GQ@uprT;_HFu;K$7(Sh%$o;yY;bsh&> zE1ke$<#gElWeMC{dLCxfYqJin^6aQcMK-VUy4ZTdZt<*L1I1f-f3CW$wysfjz! zse?(5Tf|C_(m;h<3(;fN!Oc}`p-%ldOrERCHrA=Lmgd$h#9FdDfA(OTuPhe7K0ZWT zI3rL@(+7ypUeFUyp72V1=Wa2m%{s^XWtM@@xy`U-o1A#0sT{lFtv-9E|8zDuKa@>= z>CP^?byQrte~LJ&GEuy@WTZIcjzqjYQ=gsr=^hM~e1o_R4Pdx&A9Tv=DjwfuKReLk z8(Yw41R1m}fSj+LOf>R~S<`9A=8hi6&Qh~s6S|1l@k!cbi2O70k)3Mn>a6lEp*HNZEQD?0;)2j(t1^ z1FE)@d)Lf~;Dr%+6{knk)vZX4unT&|KV{WN>JY#g3JD=EgBRl+^=`%hr+Yg)7L6o2jHhA)jQuh{J?Yeq>c?9B`Qwqu+dY8hF-9zXTJ=ha;MVpz*aTf{ayiY`*p(Ui;VD?_7{||P{N-Y z1aE~;!InkA7=N!0btm*@dxZC3CDVkg(CVak&%si0ceG~lpa=1|u#@#a*~X5Y z8bm(KUrKh(IZE!!R}*IDGqS-=4$DUBvZ`bIuz40f?A4xT?B>0CtVwSxcK)ku_TAAo zcJHb!Y@h1)?AhF@#ML;C+?~FQgf2csLj7xqsqznkz8dWOjRRO8jU{YlWDu)dIf)%L zaS(f7yo5cMQNteVoWrWVf5yIhGm=z&_9e$-Hj*Z#ViMNvBJnV5CVjRkuyft*S@k{} zS--U#*z3aO?1pdEzC{V=;1U~F~>evVj( z7Llv)gn2$5E4$1ttOz2p>=Ux@;X~rk?=-QNn?Zc1Eh9ghPLa>))#SxU4SX@u3ETJi zV&Q@`EIyKhH|lrb7KbwYzP_1t%FiOhbTrUKMh%T4KN0r)a&q*;Zo(A3Cykn&P|m>? zJ>G}ox+N>oENC}YFFB5ynwPQ3=^oZJYZDXYOJvPsdmK8=25%-Bpi#*|GGyojvf!mI zZnv?(SWS?C!D6*h6XQKvhlbkhBzkR$nnv5Ro)GipEMKC zvpQI{Wi-mW#NcTCZ766wkDhw(QSe2XdJgPLKR(u_{z6^a^=tr<9q54i(>7w^&<&_L zY$5Ku+5<1#72@UoNf>Ck38NQX#k1}2(Xw5OHor8Xm*Xty$2ZpW`*#4>d{T&5SDmM5mS){UWuZ@ST~fs<&B%t)HP zUll`_v}4c^2fAptH8r%*p^Yz&;{8#rSQXosR<@6&`I8c;*W)Bwe?6Nv`RCG8`?9HS zzf5X4V-4M~4cO2{huH>!4!NIkf;CM&&#tlRO)?E8lCE+=#CtZ7`~&VJIy{urb77=- zZy-5U=S?Q8nMzLo96?kU4P&Q{ILOW}Rwu`HSFpn#A7}m2yO2kw96^{#j+Wbz)NA93 z5%`mH`XOYK{T#A)fFHrI8DxsaX!5?aANy$WK9+r?M8?!rvdegRY%hDu`n?!T42O>; z78eH*&zO;Xur{fby(H%e1&pU`%=FcW?&-;>D-zJdEpHJEHpj2XJwujhN zuOjVZ6Ub=YB_#KJCHb`IA(@m{L%Ih)C${tJ$d_#|Nam6|r2g$y^6uIhqOpEIsX5oq zF7L9PTx+~Q8rt`gexI^Pr_`0?xZ`b-HNAy|d~YQhpA_(9?=QsPN)cDS|3Z3oeM4;B z9ufBQMdIJ0Pf90WCVNz0lg{6+l6}sHiFIWm35fYhF3ncQFE_PN>xnK-s#C$F&)PU{ zlrrLeIs9tULc-#oks?-1d@i+;V?HW)Rrf1ds&&E~b`LrjRXDH*1yI%>Y+2I0@QJ8Kr z0oQJG!)9Lw2MZASJ+@et-5+KBjj?TSAG|amha@!v?g$-@5s&ye#J3N2U9W}nu1!Gq z^+8A;g=3_D7~Wkp5$m-lqjm|0+Ve!XbE6%~4Hsgl{0(wYb{;-n7L8L~z41-d2)tqF zfUXIPu=K`i3`tsx#A_AKI+=(km&W2#^vC;NUg*4TChm3(L=~|J4b)HJ*R@CSy4r5E zSdxo^$_!lJtqONfs=@Uk*Rl7v(O=F%m+r^O z^V@RdmCPcdE=VU8D)-2-OX1|?;{;-xT|o?PX`u7#jbxf+KXJT}OwK+?B@-9qkzq1L zByY?P60!LLnP9F#uB}}|vhMF7Uf0{m!o(~xIba1ztGq}eZF}Oa$gQL$xs=pA&LZO7 ztB8|z5y|+lm!!PhPll`9Ci74CCPv~S@}{kdZ2iiQSxTqL(e9gwflV!0x2Heq-q}MM zzMmm4ijI&qcQZ+;;z3fz93w6b$I0iX&q;IRU~;Ytbxsz_@9|;#> zjmK=Gu?G=g^^xZ!W4R1|888MFDlQZIxv$8?FM2rq+F`OX`W|^)Sx@k51IgFX!GiN) zME#Z$u2}7gW=gA2ecW`^eJhL426V=^C#K?)idy1*xRtaQ3_+~{r-^XHYch6!3ps84 zm7L2mLo>GwGPJWfI;Jhb{GuaRZ?+C+b?SqA_UK_wN+g!rekRL)s-s%2AHI`$Ok{5; z<4P4ROkCI%zZi_fF0uQ_wZf6usInV(r8nToa~0Tk0L7V|tT1L&DrUr};=OaGnEY)u zPMRu*n@<^IyTbsK3>}0)3*#_nOEqcV8imFy?%@5N>h!beCmfj=jNW;OS|#gH(8CCO z4i}-r>JlvLtBYILixFE#;!}rlcttG-lQs14N5wIWlo?1T6-=X-k~mtEbqF&?EWj?4 zk6}pVG(4uc2-%fdv@mHVy4+ih&zEgPuhHA^!__tnej(4cUd?5ve%Z}VF!p10LVB~` z2N$z;!_;wE!Yx`}RzkGItEpJ8mQ8YA$VTkl$L>ERLu`5tCxg`I{FkS18PuDkKPx8R zwkhI(3r1+vQw4WD+(wGRrsBt{d^&W@Lt^dVPZx^c5U(%=QaF+&S=T2KvCax&_B92o zr%>`@LLGS;*$@rD&5$@>;9QYociL!woTYiw4ZL;Xo$f8x zKRmE+9#QCWlyn+!hBWtTBJpodAXj*wY&>`Gh_1zcb-p+&qZ0j|Rno~t z6r0o6QPtJ?m^*PVX}K$pE3~?xv0gu{c+`xME!udZ`$jxjr$o0-H>GD|d(cJpIk^A& zb=)#*KRpmQmFG(;-Tme^&UJiDK3umzP1k`q<(3Ox07YsTF$KqXJ;rgnTxj5paN26^ zMN?Xy;#sHX=-iY;>8BK&sz09!PU%x_Ssxr49EHg_kvQn?a$Kp~hw4eLW8nJnw0Ylo zdeiwGT@&|+s+v!wU7qODwSzL~vLSb|;N1f1qPLzJW$i-YNO}5$%F{=$wW%vF*V?55 z82Po6_|wFZ;>Y1(V%KfU#6L>6iKkqw5HmTq#3?IU#K-TbvXz7Tu{*wuXYcQiV;!%r zVAo8^XRm8qWdEnr`!dz+p!MU)LYFFXzSA3`BL9HQY-LHYd>y+gP9q6g9=F9285@Wdni0PQL$cGz zhC*kWZvw&C?l28pcW(tXd14ZMCz4QK2%#JyMh(Rh!y z#KPAOD@`pi<*p5`@?4KW2M0pf^zV>%gb}Z~>?t<+nk>%fo-5Ygc~~rnyd>^o_*T3` zL7Da1-si;>eQaQ6J^MeMzDvzHR;|>9MDaRSX2VmW)O3s79_dI5 z7Ce3mRA*$IE2n=WA@&hgzg1O_f+})hD4>rjfw5p5#0&AU>0fuzI=> zB?ks$V%j>A>u5|`mGw#CXhQ;LEy;`pCaC49!RvmmBuy2Fk(vwf@Lx~bCQK!FM6Zc@ zOaQuCEW^a61!#ZbE}6F0i#!gVK;)FCkZ+xWiLJpLtkn?`*UVJX*EW$bx0aB%{j15# zq}9Yn+#79E^Kjjx6Zk#%5c;HR=-0Q2hPX_TNBz1d5L&Z>#T+ht#f@N^?l4SGrD2FPJ< z&ra02?F=~-tAo$KsG{qx9+(;xfs-uM@N0QKru8$SS9Vg`E7OrWL}uYwWnJ9eNdupI ztD(*XLo91Iquv3Z$W^u86Jg)rw7n4cM6z7qUXN zt|W6v5P`?qSQB^+yS{LxX4?vAu-R2AORmyy3Sr{o{RWAbJn~^(l`6%`qcvEWLm})# z9YykE=|mFwz7yVlaSY39gfudH18sD?KwH$RDP40+Jny8p_~V;mHn?`b*zTYlyCh=@ zyKejk_U=`dn8Y@d5-39bPZl))+e$j<{7L%iZZ%!v-6qyuxI!#*?h8A4{3fxB!&~ti zSIRC}ahEk!wj#bs_eg5$T8s(NqZ=MAq_K5JsOp$%s+0AVb$eC9F5j$-Vy{-V)1;9^ zLw7sT8flKHvazUEt4RYE=F#<&l)2ThEEhht%^uY3X+~X62GjTAcT}L}$(37g=V~=Ba2A8hxYV{?+^vZ{uwtAcmVgbV1`Zf? zB@4UFdxt*8UNjQ(s9WEj-1bvx+|`+vxuh=*-1=@exP^Nwx#lSo@zm!e{Ptxbt*u;# zJJUYnvhHF!TI~Q$&TOVp+Qb>UUhc{|(GM)|cMaSL9OO z+0*z@8|cNHx3olUBzNQF60X7Y6?d(1sAO!w7|9!bkSuc?BI%-KCb_c5njD#?j0R)+ z;P8>XaZH;b-W${duV2zcE01nC<%I@*bJfB;{#{evyEFRtltq{NALQ1$e*9e539shp zq42B$=AZ0^ioNgG4e4vHZ>ooZFVomJtu7P#ymGOs~9FBiRO?-Az1MeN|f)lsNV~EoaQaxN1?UlCR^CPQq zKxrO+R>(r%%^7Hvy$)3uuEkpZ{r^S#8hl;363Z_x#ueQb;G?I>_`}~4&z~zp`*~Z@ zbWbUo&EAKzKkmTf!+Wrzd>1ZF+>N?3cj4PP+ix8zz(R4>C~ zo$Gk?$yv;_IDmj(|9q?-Tnn--0$I< zkMGbg^cl+3*5mz;4cNo|DLP)O$7u2tZP(nv?PXP%IqM?Ea~E-jMlo89?M7WP6{vfS z8om8pnU>#BrMISbp+Q%>PzNU!y4pdR`kTqoPm{mkaLpFv7pqcc6}q_l2y!6G&KX)*uXV^yf%axFT3?r>^4RzTYZ zv$Oy>I_epvr_UnIo61oiBT9!?38~f8A#}I$KsqIAAl19kpNdi=XzN1o2p05t8_^G+#W<}(Nx}kAR$Ae29VFd6as3>eB%U~<1%j4hU7UmJ7v`wU{HkVj<<&X(WvWaKLXL2~;JejX| zk(dXcBt=#i$#kR3B*V0VbPKB>3lAM9xo3`&hT)~;PGkuYeLqBETZ_rck1A;7Sx>gF zY#??v56OUn7euS~b28lW9x>|qfFzBtB~NbMCYu9q5Shwq!fv`o7WrQ#KCXQ*_?04- zX(-}Y6L~cLri3YLe`Ic(<9SVUjG49oOI(BS!s(g# zsmu?b4B@}3QtslEYSLb7n&NG~$r#987Qr;0)x zh%UurF+^oEmXDFAbFP`tQ_E6mWJ@kXJO5o186&{D?RRHM~806rju=zIQ%+}`}Qe|KT!4p7n#w*?WJO-U;_W9}b_PsRP`NK#0;o^Aih*tqO+5Z(c`jd`ihm}C`d7CjF(;tVs zkLu9;x#`$&BSS&s>+yqv}Oh8Qc%M+ywH)@pCFPg7$n*JI2jMV zE5yR0WZLQGb<{ntMNLAdQo~(`>CB=R)c2EwJ4JSKD^|%!F25ftdAY+!k~cL~(tY4& zys`H>=1km52Md0nq0N%#nq>CTbrL$URPy#x zjl?A@jlCUTf$w?y%C3N`r1v=^eA*O?vxdLK=12YLm%)3f%LQxBHf%jNJingXET=0u zd{XdVhHm+e3q*C~ZH)M`7K5A&(R*n$8my~8b1N$vd}JPtt^Go$u880)gQ~b@4SmV- z@2--S>k=eug1!>_Q7=$fSc>~&9I!1l19R-2V4T|sdUxz9syI)bGx@fRTcPuut828B zoZR3m`Cz?5QVs@K=qyW}a$ewrB2VmTbr2&0JJV-lBIr}tMGqO7bN1IZbILOnBx!q4 zGPf{UB6yrD8P(!~U9WYgHq&&d_+dJppH+`z)U0TUIGbK7uc6B_Msr8qDme>_UXs~i zvm`zr*GqIK9g@d!U+0Gqp)s=iV z4VK*Pxm|+FwUV>)-z0|?snDf4QM6@C6>aY{jD{7Y(bIJ|>6y*LIFI%S?s(cGE~`i) z*^-+nDR;ghS$SN+=~|19lYHe;+CBLWUATzj=K7zcMR&V%Ot=R(Q>latJ^6w2nK@3< zrguOxWuJ`GkUhD~-2QDpxZI~UlKXWtByR^tN}i@FV{n`k?k;V?>64e?o>jN8pOZ12 z!s~Z#i>qkfMi*|f^%?G?nxW)bZx6{SixkO=2P-7oI$I;UFGBf~y{VG^A)KM8LRape zNcZe1r?d7baPB>mxtQ`#-0IgLxjbp1MC05BN%_7!iKKicP7)o$Z=5q-QS<`m478y) z^_J0HJwDM*;r+Q*`2ucfLl23ZmY?JX$(HmVyI-Pl?L*?9YaIs^q|V> zv+2;p1GM>RPj0;4ByPIHW$xAmJ4t9zswCv*LCL;*=OtG9=OoM;Yj*bDpKN0MG-6R# zLVh&LenM+8aA?NcT06VoK3+k?>V zZ5qDg`^~eUtGMV*FKQGqk2Wp8NB4T@aCtXZa}700l9aMBlAguOB+p)Llw3&8lMJbv ziC2}%af!DpO;(eqx5Xo9``v8no7|auaZ1EZ=k>qznWmB<3!)?+>+&TnYl|g2H4aKH zmM+6p(>~*YD+zR8vKifOmP{QR&(ISqExB^bAa25x2i%ykoMg-UbjioEQi-n3CCM41 zvyxw1nbrU1?3NBcf{yfG-RNEM<<|z+s^9xu{}UU%|IeS@!aHF7{j*!7E;t< z7Yv-Set#%wNeU!<=wEV6CxL8UT|_G6ZV~tOia5f=2?rkaW%DK{k@sH|ac96j(&>}} z4mxa(F4jQ`ckZ|Sj&EnMqM z3NJa5J&q#?*}IqwuU<|b-O$Easwc_I)L;@KC?)NxV@bB^LAKiD0IPhjmi0_lAVDRz zMDNBj)@^biNqUZCRq58w{pj%r0y?1IaH?`~5}i{XPKTt%5M{+-=zCxZ>Q&}p>#he=t2|Ee{g?6)WP z>yk*FwFY0_C@d~3MGr-7+UZgp-G8Hqj{K^~$yN5_GCS*WJ96c@N3UhMJC;rKuH#xV zy7w&H9NB;t9)swx34ZkUvkY1lUPe7OUZ(Trme4s_vuNG>Na_?FO{d)Vrg3qz=YNx52`f@6U2Atjvd9IVO9Je+7-QO)|yxwZsaJH0Q`g(yb6qM1q z%R}g%reJzoFP7$d%%Y7J0e`og^j3Geqd|c_>EDD7%Pyl&pS$>|(`9s0y^OwsN<1HR zir;Z$Nls?CkbaTjq?2Ge52a*k@KnsZ+LPuUj-;WZzEZbxQ?Apvo?MEC0(YakEN3O( zL?3JfQmT?oj+GuG7TK+&tH}(UcDtR|Fqctz--A@IFaPFyDwIa*h0ti57`igTk1~P2 zH2LmQqEy@+@4d9d8s`+uURQ>0MT&IYb4~iPi#+YNy$YwEyo1K0F5|wKOZaeRB`%(K z3jf(v3*notv`8zFUX{zIZ?B%B{nkCC8BgRm|Ji2T)iT1(yRw`M8(Yg|m-LowJMAY~ zQ=BVl>~~G_#b7km-g$_!i=NTrle=&Q-?Rh9$orA?1Wt>n=@C{ZAuEd z*=U*$v{Ab^)Mj19M4K^VXW2{{;&0RRd@29^N}F!Mg*L5*_icKAP`BNEwx8YjKd#ZE z#=DIsqgh!w`QLxC(u=?Qb;-%e$;-cn44ou!^M|0Awn+LG#j5Pr48J??M4 zGdrIDtJ+kU*)MHn(Ql>YKcP+Wu1``OM_WdfcZ!KozhGiw-zgv|8`2=^n{(HpZ z=YDt7-Tn6Nf1kVOd5WH4YTl`?uC5N#U6c0b`*T$cW`+FJjnUQPd7dg2alYo`5aR1#)|^}$G3fBv#) zk~`?#PVIi%WVrU9+GH4(nfY(oWcWRvvs_W4%H)@4 z^8e+{r?TzS;(fntJ~jL;*nnYTMz*22Z&m!4VbA#Vgg<`!srY$V7ICM!cj2G+5Hde1 zOaEC9apA-F2g?0#@uxIpZr|UJ)nx%7>T@C|XY4S$?QfRQ=by5KUp&PBuq8Ns&(5tK zGn0Sai2Uux&+SKMmCZz(7B@iWvbH{$3N)Dd6Z&HVv=*yyVq20BpBtakJNL&ue>stD zluk%W>z$sHl#%{}^>20g{!(8;cAqTO_z_=2<_8M=>mML%8D-y`m0kaRXeL=q6WR5@ zWeKv#AF>=jPr0V_n`P|!y=4UbEgvBZsUYJ06(5n)HzPauU+5&9zvpN5vnM@%HI3L_ z_1Tm~>Ho*S1H2{SO5c`_&`c{CKoz{tcNas<4q;lb<6uJ}E6JfpPp7FZ<&z zQ7?W&lRy57obrbu2MGD{zamfjf$a9<>HnC_V%gY4pF8m@@})lvdE)NQrN1KUmDLOT zvp)PodKO2^LmM73{}nmp4?`vi&1cs4`M&+`{Y=XrnwF*iT%S)_Zf;n8P4fRyp=7;( z6*-yT_Sj_>~ai?n~vwyTo@ejlSRI@9sy% zjQsvCv7Joqo!%$6Zx1n)*7wFm?5Y3v>cIEW1RLp#!C;!PnGZJaAta@yhR3Jo^%kEh zC#P;`@N>RHyNR5(wEo}h>pJ!F%W@%0<>Yos#ZWmX0*>uxoSztezvHX>&O`#&+KRQ(W(Z%?)(@+2I;9#uqMf>{8V@RLx85grJC-)D; z1T&1!HOyg^p;vD2gzuw``g)P0oLO&%|Gg+jzt^N#j&yWr*rx)u3XfdLGE8? z9C2$=8{<}997OB3`q3Q%Pal19rY^n~bj$TNG@(?DM)y>pOEf<4wG|uryzxtT`@Eif zKwme0>24Q((}d^5uI@Bxn_59WSzqSL0{7! zqN_(*Qfa)beYc;zW?5%{D?&(_|-LOypK;3FF&yr-}Z_kpQyf^G>A|qI}Rsu z?JOp{Z?_&rO}>qwcC|w&oiLC-cgUgMHNEKcal!QbICuJ~bywPQkTUJ6r$DFfS!i_r;_(s>&;)g#1?uzdPF);^B}q_W^+vqcDkG1oku@SpGPA)&Y}LHGwAfLljyGZ zBdOi-f%M3!96F$L5}kX;fiA4EpgN0Q@Edj?=Pzg<1O^NDdYWJuaYP`XYwh;mNPB5M)yv?N#$zq(B<-XX?~X*^q&7k zdeiCzmA|u(HtbeIFX>m)cM6lJ!lM{!UhPi#(GvQtks>wMc*CD;v6!C|+=_R1=}!`8 z#c`$^6$H(K#)1y7B=k2`7G_>i7N)&b5_Xm=3jKZM1tX8IblU9mRDNC^RWX=ByLKN* zyYv`N<5aw8^8GHfPQ3yBvLu(E_v#oiD`T7PvK<7R?INu9bQM(eIN?mWw-7PXUsyUL zP|$xGAShM%3)y2Ng1oDv5SZ6W_%KdMaQW0k7<}mj<->NCm6QxxXb* z>D7>&zSxMo{i;Y(2Q?+`mdZqKa&w}Ws!9TzsgVFPHS)5h2Ju;^PE=-SkQgiUTio;d^S-dsmki)v$6(JjME@AQ)s%Gjszl}tZB8zYRUvIps}ScWs$}F+H4^hpl??J%BO2_yHT(bh8WqxgoI3f| zQG-})QzMTmHOZD1`sDnHKJM4870KCsO-R#TN@U;QX5{mM=A@yX3i-x@4fSVo?6FrP z4p&sk`V1x~Q6WR@R7jq{0v`-iCxiQ`lM5=U#C)+fF^bkBuNnlnPhUHb*jKSl|JfC! z$<-mmd-QO!q+%plXE&OxOdmtK^%=`HG>#|zFoBHkGlA5xjZ2eSPbBL06G%YHcv5_I z0_m+Zf$TMzM)uwc<7Qv0B!|{cC4+afb5s*^=`3<+-W)P@?>y2^c>$4hSV-i~E+%<1 zmXJn`mJ+Yrr9^+kQu0Z%l$bd$Ay->1CDkvNkar)~kQ4HG+z_lFh3+d!Lc`T0Vb2=U zGi?L$i`qn%wB1UEIcz6;o^B^qbT?@(xR=-c6o{ z?IY`w50dt?rgD?C50UDML!>SOLeg%MCW-gR63JaM`sH5EcIz9`_5B-icE?+CW!ih< zWB!>ep7)h-ehv76%jNm*gBAFFf+AlN(u6-aM2VkwxEVj~QZs&?e>1-EzGnP$uV(y8 zrRMzICQ7{9GDWiGoiaa(ohybl=Lalr&VL9~<(F%!^E0FxyxO9cyysF)-p5*(e{e>Z zZ_nxTNpS{z?QjF$VUGdtSYyDOKQ-j5*BSCH-nHSo+%qG!{vG&{Y7OP-wGs+Lo`5lwX`Ia-wd5y$!ezIdZpVGOUKRdgeAK}=a&pXkduX;a-S8h_lo2iW8 zFDH)V%|?vm?`|8zKe<1KZ@+gO|8U7Te!Aj#{^F{uWY3UR+&)IzTk3GJQ94}lD;;jp z2wkoTuggtz*5gbb>v0d0T5}qiZMc@fZMd|D`kcPK0T-&tHdpa&Io;5RkL)tM#uFy?LB&&abz1#|6Lo-wzduT z*3pprkZ;I6S!BrBUo_y9PaAP7=UZ?y9vg8Rx_04wJX>)kcUy6%7PaEaq&i%SSvuST zbzM$9RF6C6sK=4px}3$zHe9z6ZMaBv0~X|>EoWWQmg@)uZgzVEu93V6x3G&DCuAFQ z?L3URX46}7`b<7`YaMR?cpWa6k?sXuZvJy!uE8i>Zu=))E~2hAH>pc&E~|$jceH(5 zZi;4Gu3(owcVU46x6j*@o8oQ4ZLc!nj5JlaBf8VMoRI0<+VRu5kdxE7S);4D-TBpA z=7nl5aQ_T$?VuT)Q>$5A(fyg6j{6)==i+S6=jd#%^x|A@iSK;w#>IKu65VB7?t{f# zP5MU8W7>M|iqZydq`?Mm$#NDvsQpGRs9+;kIDR8%vw9O}ICT>@Fn9}>)@}!UR(>CM(Cav-5O$nfed##2q0I>{H2Va1 zIQt~GpKbP<_wpondDkg!`S#PC)9f?crn0kKj}hm&UX#vqJ7%5Z#&*2GHJN;cTio?J z*Sz))SL%0{Gjn>$y}bC6JFNSPiyioi)6aO#XnkjUB=#E@+S8>q%MMwhfC z$3L|s%|Ev!$s;sLt0S7E={7A=IaG^Co@o>RWNq>#TbrzZsY5LK>X0WlbV=8zx@2>V z2{B8xBw8uvq-n7^>AKRK3~z5iG$vXQqp=p`J+&m9rX?AZVoAEiS&`5N9Y~>FM{+;8 zBiR|kLPNybkb`G!$fBk!XwVEhq9Ao4gYU9M?#ly6+f@PN?dt%tA|;S)eIH1oz6Fv^ zc|qjnwIGt?A57-5&1D~+1e32HLWo?iFrqvyjJRzKBl|yxljGkaNb2xda`t#U84y%T z4t1YO4h)-0o^PE>X1$q8dfQi#C*@VdbZ`}kWE;&YPgRlY$EK0TXQq*ysOeRg> zSWVVz%^+H}Gsu*~Gl);sEOI+`Ht`s{hgff%L2IN_=-tqf^qbBg+IB@UZFM|?Hrqg{ zvZe$5I>d^$4pOAtq_6yykDK|O{pa&F{9IoCMINsn-;1|?VZrx&x|eiI;z{U~+1$?V zt7wP3MYOs7bedr@f#Sgcx=A66wo-_tlPbOGu-5K$(0*;IazTZLc09pXm~Y{e=56F( z=nv=5dP@1V@!$tUUL;RPhm#>AmT`>_Zl@!B*3)+;i|Dbm8T5YRQMBQe{#1~p)76<# z^tE#!y}sX=x_8i}FEnrRFDwu8mRbAx#tM`87K?}RFFyG3BPKl{<-Q4IdEI)hi~n(& zqH=&%Z{0)(7_OjU=F{k&_oL}{$AR?Qu1q>}Y;Rgv*NHx&7IeGQN50R(t9<*}m-w=k z^Z3|SmAqD~UVP%gZ=_2^A?eg%H+RV@jFVR`Z>{t0ZnOx@Yno%`4{}9^mTl*$_4xYzL0;=q#3_$P6bJ^KE@pld`l}c z9?|*JF4A!QI?CVNNZm&-rrS16qlynk)0W;t=|hDGYBj~1_RiF#r=B&WOAHjK@O}@! z;nq5S?XU{o`lBXK`-~%X3omocDjEyxmdOdz9iOp2?G4(f;voI@b~BZKw}kFzebX!F zCQ+@ssq|({5Bji+39Y)RLG!PwQRjs7{DqnW{Q824e7D!_c-{Up$;ly)xV~;0g4d^J z!WK7qVYI?q+Hu__+U43IdSclY8b4zx?R8=rtO98H+6*;BQ=Cup$dXc>sQpk={mKotfSk~ z_fh*@(`kv~MB3FhmR{e^)39b7l{QzW+aEQer)u`|R;kf^O~z@`#8j2!^zaqBw&ex$ zm@b0u-VVb4lGZ}g4(h@a!^VR0o_F+mizoDBdum5!TYOSd&| zN8inR#c#-~;NJ(!^U@+?7T79HP+=SKmUIjk7Q6WfeJ^$qjz6{#zO~aARN^&+K~BoT zSi2k4-{lPL`f3(EF?$Rx*)WP~ybq&x`zbAM*n(akF^4}|q0bxKbtD!&7Yd8}%@$l; zDuu=ta5?Jt~8aSU*h}mI9^K1vKk4;jK0!63J+<{g|)!l@S z;hsW(#92t{(q33p)=KEzzM1f*_ZMng`kc0FwTZf@Eu($<=F!qcp>(qb8`IBJr#F-u z(tGdL^P3u3^1&mfl2$9^h|8|QLP(fYIA4<@^skE%q_26w*ULeWyf+cFrfUl|Mry*5 zp;zg*KF6t2pKVl^vv2Ph|qLxo3J>#hH zsW$Y$uBp6BbtWfCGq_=Qpa+L7k ziKU>mt*tQs>=Sx(>2aE!bCfnnoj~`ktDs@eC2YgtHh!XQM}A>SA7WkiOc)+;S9mhx ztk5s)kT4~8qp*C>5~0#)h7hSWUP#|ILKtM-T?iT(Bow^Y6I_Qj6-YltVZ({N^iI+` zsvp~*zD|~>DmQX?<4xnp(#NhIvp2IZJX>_|*fX!4N9r_fk89IZJe-Rfd01M15GJg8 zCKPVoBkbC4_QcG~KchEk1XP&s%zpkDoG`@6P7;sDEevW9~w7{)rBGse8jc-OWtc z->_>+lvY|R(w5g@HCcWO_ zFkj=dk$>E+kUyk$hp6uiAg0Z&IGvlmLaz~2c%0xU2z{-Ef}A$Of-7pmr@Ds1j8!ja ziSZq(l(Ld)te;C0?fTF+y*+5=N~03>NwaN`=v*`wC+}L<)zddkVMvbrmMx zvKIEJC<$c~8VL5Aw$M@g=1^iZie!*O! zWAEw0uFjK%H+_Z+H>UI#Hk`~92Hfo~yjdM4?DH}gTw)D`Hv8|>gqMe@;_h0y-fc83 z^y^RcH`r4phh@A{xgI}nt2{X`zf)*xw@s-1vQF?>wnAv>J4etzQz>X48!c$>86*s< zEE9Su`w9K0x(hqJRfWCJKhcSY-_aM^YiLU79O`hm7rj(;ng15&#ZMe)Oy<5mMlZFf zrF*_?p^vw$qK7(G)BGl5>Dk5uX*8SjYq}whUSk`)Bj4!J!m(#~mu2huLtWPJ)@Et^ z;Sr&{+)8DBfln!k8+M-?F?*%kfAIb}{oM%D6_VL@i_sy>CJb`>G`hy-6q5o~2_9=F_q9oyClWxMZ zBroChOgF)0a4TVNg{rWn)k!+CVj~UgxsD#tFQcztCem>Qh7`x9^6Iuqd{Z`OP%~XfB7*#h;7)FN)ou-!yVLBPYHJ9GPrCD8sG`@o%NstpR$K9Yk=3SwZ zv(4aal)n&tpG6lw|vF-Y;|Avc4#@nN;HM-@0Jo z?|Y@$U{T4Jvw!+6K=zJ3?(12_1`WTzJZ9PW^Bu!#013O+?FWT(7Hd7t{$DWdjb+RRf1O!GerqaYHvGZd zIyEUFKD$RYTl$jmU$Dnp$k^NdEp`$8@}_#O1`UmWu>Z+t_e2fMUnE0p_k%^1eM3s7 zBN60SrHsSm&v1wsrVC{p;>W7QpL!bmx?k$AD7Ti#c-Z!a@4r9g5pfQzs^|H^z>;6x z(FvRStFlEm|AK`tZU3&9?c&OlzrN&!Eq(gF=0VIzPEvZG)ZT_Uz2b9oduOu~!}!#m zy%REiS(zc`{co?#U{4j7UHx(EUA%rsDO>;I{kwB<;ZB5G?uLKQ3XH~$>Wf?WVH1Cc z%ldb&z*x&K+52-n3scOS?E`xmybf3&TN9^r4>YTQ3!tFpB( zvQ*TK%r^bl@vrCN9*Gxg!cW)Rrrp1}c^y++UylEctz-YX^?*iy$%%=2t(TK|rEI$w zob21N?EmaamTQD)n?F9498%BO@W%~*yt9>)D>f;<`Tfy<`A^fe3bb{R0!{?RKkFWWo^X~A8%U1C% zdOGq|yJiue-Rrn!tupzb_#8gKynx>^teDs6%3k*xl=J-qOZcZVv-tNPQ~1?AV)(_w zd+<7&NAlb6&)}1G&Egxgo%2&%Yx#)aBm7IXx%{G&!F+UaV{-h#UgB_YKiTSgh=k5P zOd=Hyksya!5)o2Mb{QTfE(ed2t0fmmNZCaa9(tOvFHsQBGj(K!ZykBn>p1D_4w+On@Pr5dHSVVp628#Q151q=;?vY=-a#Ml+0~K&yUfklF+tP=aDI`jxnVf z^L6NuPpY)~o(escd!M&!dy!wac_BZxS*M@*m2Z;Ghdj^aCq@+VE-g#>XYI=Qm}jN@ zoVB_9q_||hY(oMcSJj=rIHNc3c4;C%s%!yYV!VhqnzfIg?{bL0D>=j~b(r-tzbP&E zlij42+L}>Rl&8SGVcIUw>(Ii*c~PM3df0J?FBM&=LM3Yb%K;_IYO3>IznD< zJ4Wu@K1r(5AV6Vl6SUV%(q#+gx|G$AFuwIm2d1O z-qo%%-}JaWdEXo4L+AbE;+q5Hljk8)Irk75GqRQh-#$Rb3^+(O_#P!YSp3IG;R5mY zx?j^hH1j- z-0p&IM=G=$;~mTPw3H}H|V5<)ii(TY$`cClTMvmPVZOc(19~N z=>tVo>f>;l*I~<7%Aao(*56twR1KLc%pO=J#0?!Tm}Zs>lb5Cm@ocDmQY}On_SsSx zHP2Wm%x@rgzIjCVr9Gfe4{oIvYgg0um&egK_|vO(+Vteg=WMCVav^!)d_kdMl~6Zw zoIocJ5ORXE1)r2S!LC`LuxpWrushsXkejS0pzqh(Pwu2hSB4%Qj|THV%$?HRHRZI50%0< zy$p|^6k(rAPt;CIb2-8I9x-~$=x^bDWq4CxNOA0$6wil9;TOkbsVKuPy<+xl(KuYU zRdYGdaGo-H_Y0m(E?FwdCqAzQqYp2oh+>#Z6QtOmF2&^q?3<_2#Xv@PJPP0Be4X*U z5cwF1^2soMO=Wo2?E1>rQh2#>4wEIgDaJ&o}>x%KwHwH(hW;)+tJfiMC{sk-% z_sfs#<5=3`qZDgCNTGU4iqU(em|ZMI?~DE5IX?riRSHt+4e34=W z!#sIkii#VoY_>{aGP?|SuV=$?eVwyRcVeDo^2u;ysfZ^&PqeMM3@@6|>?cglds65e zmSXJrGMEg>#ftp<&aW7cs5=>+4CB{ShIfYXQD!q5j8?I7p3dyBisgBnh7@vc`FQEj z)J3#YQ8(f}qI|#N$d1Kp;`1&un#-_C8I6AD=s6>?jAU4%E<_tW%_xoW)?Jrkh76qOIl z@WH$k0UdPFf$2sfo-oX1JEbsMCq>8z_Wk9jWk|kJ zg0YD@pp0(@qjQWzJW)28?qqnP9&R%d{k;LJ7w?$=+q_>2ek&{EN$mJ}8Qu*mL3mv& zSTVkpjKsPr>P;q}3`@k8$*0VY+rE>c8LJo9j!IFyPm1vQ%-_5ygQ-pl`%ZK#6fvHi zj70s2c)yoT#1-+x>!PlWSUpQ+m=jJjy))l3dKI$^R=$ddit(j)D-2~k2N>ODB$G|d zvrIlQ747*NBV$JK%r{3e%pz7#cAaFhZDwuW?J_J3FUB}yT`XigqVB{t?)T*?Vu|H^ zk&$RG4vZRnlA?gs+uqEVUb`g4umjBB94>>Ce-RGY8elbh)*@Cf!vO2t^AEIsgA$(N{6y%Drwk!woFVnuaEk3LsZDM%4**(V? zU0@{YOIDx7{}(XQV|w$L%F1Ls%hShx7(TThKE|-`f3Ge^bF*A*=-3BanCv?k-b;oj z`g1YwqOL^Q#BvbpnrQ1{xmhsXOvx>hSv^)yUzXuUTp5mJl;Y~~ zY;e~(3}txN7+xvkzsyMV>7u>8W+eI=5pO%Axr|OTyyh2}4KU2dR#IG1lp=jl86KpS zqHA(CQViX2oZ%@kUePanucIs#WfSp4ESY>W*l_~W&8N$(kHs(*-B|y?fcefvW!Rcj z3ZI)<$Z&T>M<$zChbJ%+^D4?F>P&_urlS7lGrGp`o-sR~a+Q@~cPSdXGycx@_)0>i6dzWE93XLM)vBq&Xa+}2WrWtW1^ zO2>^cozR5w4P+$x7|}*#c@=$)OulXGShPQLhPQ|HO&73!#)g?ppCg(5crkq%vhwgv zLH1y4G+dX6qs4vjKCBO(*!01sGl^(3BoT?`iAaByfT8XQI2;>~cKhOxvN;lmi$gK` zO$au6d7=3n0i~0iG5ww?8eGwbYQxEr9e4U7WquNN4NF4%)Fc$OO~Qz4eKB@SUku6a zgT-};$hJ;E@q%8s?;3;8o{^|H5`nCJemL*J0_vW0!>98d@a}dyIGB%?*o{iYWu)Ln z`xNYGmI6s_GVYH^MsF?|I~|h{F}p9?+9kqmZEskg=nnhV(MVhvg)be0(D}4KTHoX` zb%GsMmYXB;+z^RUaw-(n(y-`LDvEBW;@Xl__~)l$lmk2GQ_yR3GVToQi$3aoK@Z1c zbbNQ5d=m?O%P>s&6pZGpykO+t7318kaU-`}Vrr6(?OW4fwj>>cC#6HKS2~6?yC(^Ao2k_s2*->%0cK)*=>;+Dt2`&<+r zy$FWfQi7py!feyt5)HE~+?$()X_Z;nJ}e7*J+pAYH47WmvaqmcCN_P^z=HU6WU8j) zL0J+~pCsb`v_$w{?*VQ5?zr$e9Feu2ST+<$k?$tCnUxLxRyG!2%7!QN5s#N;GryRP z4~f~xJ(7jp*;#n1kcnAaGf=rO1qs1PuxpZp!pA)^r%ya|Peh^59zV<*!(&>$pQLnD zE*?J5#REni9_M1|xm={~$wlzuT&Q%+h3DNI?B2rq?)|gT#vudpUTN6ZEDd>I`e5yr zMC|Syk5r>@oSNc?gUXH)Lx+5XkIBcwiTNmIBsV%Aegha~=3_~-e7ruMhp-#D$XS<* zjqkFM9+S!JEE9*`rl85gWGwei#PU}$NbeMenmeYF(F+RT_ND-yUkX?sumFYfh0HG% z;PJBpe41GRyQ~7ltjovdarsb<%tf;%IWW7LjaK;?h^R`(#-qtdZPp8#8=^7xU@M8U zN-=zVS=$~}4Ch|O7?fU&`B|(loLLNsdoj#Sieb2^2vzfoz;T7xpIm@@0R?ELn+J!n zxgeI=_-dL4-CK#cnfTdNtn+`3L|-E#(N~CVf!MxOGHS=jocY9zF3dk%$;GGPX8-5# zM87GL*tZh;Y(`TU`HJlh!;>)oZxO)85RnBKJ^S0A;E8=nkz`o1^jExT?D|+nHVjXP z;caJa$qUvlmAx&%>=yE_vV4j6%J6=l%Gwf9zKM*6F>1x|J}_GuahLTidb0l4l0qC- zR{V}9!}`qlMEY}hqKqP*DBB=LO&H#NhWC!OpS?0!J2tEc%GoV`;)yu2{-_LBmdbP} zwtoW{k=fF zy>;1mC==~=W+GvFCK_jFV%)b3w5`m*4xJ2y4^2mfO*%}CQZZvl3iinNMWG}SXBs6U zdQ}|yHtqqHEfMfv#=ZhSu02A+>;8lHPwMZj6>B>>=c9Aid~E2D568y&u)LUuL38p@ z8J>sMR(WWwlMDTWIoQGapy@9&;e9$2p|{g8>O(3*@At)Z#? z(y$~r0lgGMk&|WU@=UG>uQZFW+^Pu1_C;93D7_s^Rf=%&VIkU_D#SOT5Y6=pVW3uk ziXr(pmy-{i%0+;0F1o(V!q(@>D2?rok8Sl`ibobBeNHjfZ!N~0+G41jVq>vm#W=FN z7#r9)Drs&pUeaP*G%H4BZV`xf5mq!ULgJ$W_>V5Y#z%SBq?`r0?A|aiZ0FLpb14RS zm7;TeDT-LXI;@oSt)->-Qc{ZPtRL%H<;jy02uh{6-?0?&d?~v3WU{dF$xbHQ zafbbl$-Fzb6h;hl#>*1y-d6&-ni71mDZyvC5-3>|WAMBJC@#xHO;mu3LqRrnP0Pl{ z-Pss)A{*zAWMk!;Y+M^ETJDAg?V+|2BR^$}N2PQ%fhG}vh;gKwCGWi|;o zs2z^O6M^#Nwyyud`{(p$uqz)suIIy{VF7luD8M+q0wl0{5d9_}pHAju=+=B>s^_Eb zavp5PBJk?W1deqr#?%yQ{;oAvP;79!uG5bmyp2=^<* zY_CGJ?plb``h~DiDn#jER+pm-U=^5;?7evyxGIm$HRRxMb`FleXTE)Q6211Q@Fbpih(gHSrGO`FZvst~JTZ9JFig0FV5&ZiVp{0Hi0>2eP*SZh~_7`Bn<^n88 zV*S}Kd5F=@Ww_~hESG=|>fSCPtbXp?S&ZB}#pw637=1KKP{!t>s|)gZfqf*ITV1`&p{`rB!rqXJPE5~;iXuYTZ)zyr5I3M z3e)+eP+VAwE-dXbyA;L!Se;?@?uKqDLLW2yJ0)o0UxI;cO3-0#5jL*O!McH|c)B#m z86B@=ILt>pXLYZeNh$7pEWs$&Hh$?_0!Y54n5dPFRQ*(!p6~O~HM{_Z%%5e&7GgzI5xNEy z!+`t@s9iI!Q@p-WSkcV*>3lV>&2tg}LA?e=_v6K7ZGn;o= z&E|7vvbk6__Rc5yV?Qhz$xA{$$HCR2Cw8;)o~hp(&&m_=CATkTOi9MHo2fY6B^`^I zt=(Yl5_wX9PR|SB)S?J^pNn8`Q;dPeB^bY;6w2pI@lf-nq;u7!ztFDOCk{1QYno3w1)?`MARZw?3dItn)pdSdIdJ_z2LhPMH1EW9%pwz2tWXIFsg z4#g;DayGol+CRk-)HE%@Yj(bLaxr?c`ncIQ2i2Fl!F))DOLn;{-gy#CyyuI-n?e!v zBpQ|T;!xh8H@=Kc#055wvR5u0oxC$~B{&ZY@8v^YvjE2;i?F|jjq8RLWBH3hyuFYH zt%xD6<&nmyQ#Z$kJbSF}>w@hkc~pn_ph_VGFI6KD9oP+(y?Q}@NN*g??G2^4M9eHp zL?5dpT-}z8FEcYSW=cMErFlqdH|r0dFIh8ps^tBGNsibj4`N4gKMgH7WfibMNL$BDW<=J{;*Tc@`BZX%OZp`TWXO;_R*_ zd046-Ax-2Z!H*BSUL3o_HAoomYHw-k8vA;^OYQ9bjz>q2m25vTQnL1Fg`}dx0LjLq z*^-zuiINX%qa+=!`b$*$drDNzEF{75#**t>|7R!TXrssxhp%DMdFs752~uim7VMK|M8$7$?v1}|E+L%KX(@V zxc}#FfFI}SQ#%^}S{f6c)?fI(E}9*R5s#bye_I+&e&j-y%CsR{4*YXV#}7W);%eFT zpPMv(T$gD}cKzp8iyzlz+LK-Xxmn`J_2Tv}nRz2Sqx)kAgna9S_Ydt6y4o7UE3B}m zuQ?9QH^D-EV+5*c;;^;`HlJ*O7_+Amhu}w&(chLz0wb#>HqCoU6g9uNdS^~@-B^|A z99-jw29LVpz=SS%q|+JWN7^Enc0|$@3%IOnk5eCw5xGSZN*gqAr(py5raY32lHQlZ z4p=P7FqkHpCW({e&j09Y_I!eCD-!SAti}b4H#_6@D@SA*b;Y0ndnDem#$ih;}jz>D=BhMfJ?STWj zOtHstdlM|^r;nQLZSW#Z1(NPfaW3GxIC0e^LINxvY1B)F#IM%`!L+|*Y z(+Y2-xAsCcH4rrUgJk~AiIM?H zR+4c;-aBue>4VqlzNmEeLq?t-=nP-19Or|6Z@ln8;)(f|9@uN@hMUbLSa;tJ2Q@mN zF2o9zOIqU^TmMltL=oq&Es-dx3KHel&0P{E_~6Yb{Kq%`5!f`_Yv>yB5)LL&` z*C1HY5s+Ue!QxZ)$O^SV=BBnNjB9}PL+41&oKA3=5bA+35-+T1=8Kc6{-AMz(DMmG z!QmiO-w1-*z94k%9)LH&{_u$R#0CwD%Edf(nLEQesVjsdmiQvChPsCPB}3ngaJgbe zP@2!f^S*!zc`r0MF{h_)e5KRn&F=0|Lj%kFzX>TwlZx4WNr62fkKYXka&}|lv z5o?{X;JGo9Iyc7iiaM7W;S%IIOJL(ALBwDQcHMMEd6FwMpSoaG8y7@wa>kX#POzEg zh%X1az+2K8hu=A%*4G}xp4*|#IU5YTWQ{SyT=3GmM)L838`O8YLF1+y#>l&4eor@~ zJq1W0FfU&M(p7@?#jNhEaz^`8j!<~t1;^HOfltfM*t6OJ58v3~<=IZSea{_{%{NOt zrx09^C)hrZz@~s;t^$GKI1Xj&+;P2wJJw!tgUvDtKDtPdvCbJw&O5?sjU(=y?246% zT~QF|fYJ~r41Opex$j!ZtWE;fcnNsePQVQv0oJ=HPTt_riQ}>4IKj*a0^JmM^m^)s zH8Wjd#<{@G-UWTCS)FO@j1hcS7#X<1a)b|BOIAoy9(kkja&P2}@P_9oZ;T7@#yyEQ zB#PdMo#}-#W>3ABdm>}7CtI&3zq`du}CK(-&6C;A~Ce(c?-A4Wa%MJF|1q#g3XjmKWNY2*dp7GCgB@<1Oo0lSoV zY@W7TzCxXMgn7^T)a8ei+%)51l{xB5slo z?hN(D#YAs}sCvQ8+7sskDc&*NM8^3beR*5dObmn)ABeU$0pyy4>D4QnkAyn0D7{JuBzCYfWsdobvWAZ*PJ z!V~8pte+o<%c+6T2oJwmH6PEjuTYF!_faiD0+h-g%y| zd+mX~4?OTz%M-rcS(`JRqSlnhPOb}Df3?8?%>YOn2jE?^0N7viN1q&joMpC{P|R%5 z%n!wU06(eZv4^p|IKMKcH~je<~d zI1s|+Kt!4b;^Qu6FCzo6`vvpaRsPuD%^%%Y`@ta37bgC`XeIc0qqWV7g-V^C>I{$-4vLH#ZPp>zGXn0VrPQkHOR* zPfYx=vAZv*hA$$%_~1h}fl`hO#+MW591F+TZQ-!pAC94$!ZFx092**jqd{#LmXw5H zt#=pT2voBC zJl`3PzFuK?#DyWYVHj+5*?J--r=eYUjNKarqXSX6w=@bn<}>UbQFy7t=I@$BpA=-D(J&6|Z|RFhCt>V~3O>v*_rjfUHbXjrrL z7yIT%LpLlMJ6)rZ@`mAFi^8SLQOI|Tg2fjmbI)$rvxlwASs97DYLQs7E&^BghQnq+ z7)Eu7!1I6v?2e7WGyfR&jx7c&17qOPI0h~282;vH1Xo1EBqtg%ccQT7coYtrMM1K; z8y?MNx>4zdlUpN^8W;h(H5}_MM*4K=|WhQs8aF&|3(`LF)eYW zL<4d}1G^VA$D2vY_$1d9V=WuQ*hn204j7>AOH=sUnPP%fdt8e)VV__%#=GH0XzJS* zkJ}i+CD#D%m-TS$u`YZ%YJ*5Lk*uJJFJ08&&{7quhRrbVUSlkHslw#d$FaEwy}yU$xe?zSZ&3@lMuVU9=L%@Ezr6dA8g5aMNmzO&k*+Y&?E ztZt1{<8?7JR2L7R3GHt!P$X4FhoGh?Y1$kehqgw$J~n85z#9ISJL0B(N92sM!nY(# zDBD|Lxq>;CMVi5SX*<-dF-Fw@1N2F4gJV{0aQdS*W-Zr(+e>xyyr%@M8Ok`fNDqts z>|wRp4$V@9nX|!yfBS+G6ll8_XSPjr_Qdc)XzlW++== z6Pu%cwzeI*R~q44gb^k*V!B$W2b)P+IBTbZ;;G7b@mU8}nO%@|zcZ$@{!**H4)ATl z)+~OoL*F`E44-EMhtJmN(9H^$Vl5H*wmoKFZ-;BM+QG=%0LgdTVAsZ0IB-f8*KR6f zgt`t|40FN)O()z9bi|;QUGY4j3tsf>3^(cktM;8xJkTEbS~eJJY>g9}EugT}4DVCS z@RE(|8f|S0^^(?X{k#S|Pc%n^GA)d|5-gicJ zY-coYX%BgGJFL3h0RcBHv2dOx_9wK*_C_Y?b=Uy&FKFV1yDAPjw?J77pmdJ}Dz*}Q zTJMU(3_D?(GiI!E!nY-kP}Fk7lP3$ne+36%{Q=r>SDVKBo6{R+bu#5T~qwEhTRydncTgBrIAl3Aj+}`hN{?y{#Q) zKek1eL|bfrWP_Jl*2r7m0S)b}U{z&-&#TSxct<-(l8n*2qBRVpItVb*!DVlC_!_HX z*M-Iy_W6-yr}sfg?u@OjxelEW{Dr;i8EKC}M)r6SZ40G_Hi#PD5uZL;v3F-yXarLP zKQlqKsR5iUT4P9^9$b`LV%!6DB&}=;o!S=?UiFxy$Bvz@9|JpMmsw|Iop3;NHmagrH!=KGDV!at=wlu;v69cS0tdBHvZS+>rM29EJh#c@p`Aiv!C3UwdBxP}R1rji9KA=&?m56f6)4Y4)0X ze0u|uida|(NGT0A28xuVfQXbJNC~2tgaL|zT}O{`Y&~{&y#=_>^_&~;zwe&+zx)3G zOk%BX#awfa8EcMMWBjbg>H>6GWR4Ecua95{Z)h;H#p-Nt-auwyGmtI0tjG=>>%+R7 z@5*MmNiiFTpSe`8GlJ4iraaGW%9J|t{K*p&Hay&zWy~>R(t3u>|BgP(Q`cux-f6SK zGb7o#6Z~A;e+YZnIGAyj%4|>70Jf`3FXqr)mdzJF;QFn(EjT}xV-tNiHYJQ>Ny|AF ztR`ST^fF`felcO51;)&tm(5JSJc`NJ>M@tgBiKj-4VHI%IQwva5F7DMm5JOGnO?XY zGcbI~-K}~k2)RNmraRcA$zbOffvq88DKF|WzK)UX2GVkwPbSxEZL0|3-;lCU(B#}TZwqb9qNrOIM1E3ty&u54MD9J_m{BRja@ z758T3BhG!(an8^vjr&qPfs@XfWcFu%{~GwEqW1UB`C0jL_$~YDF;;G9V_gk{~A`5nj=XJ`oz*28B_I8zkS*e>b4Q>=G>7~mSWesM|vsKyEMiu7Q zi+}fw?!r2jOR*o#k8on+6z*(~QnQ^sZ5S07ZF?WM#!et z@cnuXFTYGTVEtq?S+?d7wkmxvuS@UCRND1opB}eoFZC~QQT4gpo}7BKbLzJ2_&poe z!PbVUuC-?NCr2~yNtSF@xjD0{5;2#2FilR^vtxSi?AU2(d-md%9g}@*%N{D(vU9y{*w`CZEO(F<8#hhNK70n79c0RG zoZ;ob?uJZ_m#@YJ3}^9;O6=Llc1+>ISp z`HURKrpzDAy4_P_-p>75n;wd+&kqXh(EbjrX+Q^-eBu=sKIJi2GwC$fQ;@}7_i*7v zJ>1Ox%IR_a6o9<&S2J_~ zcB-)tbJps^UOkp$q12v@pZ0L8wO6-nuH|==@^S?KYU6{hl+>fjA za*{r*(}_qf6j&PXMWUdzf}1?&0I!Sieu^tg>3ncu(@o%MQ%?gBknkZ_1LM>oU_fTI|^pE%rKX5SuVx zh24|w!Ge0dOnK$lf+R{+rJ;Vutg#GM~14B<)14rY%bP6iUik(zpAyYe00{o4oart z#ybs@lG0Un)|8TxX!H16T9QP2C&4#fuheRHN}a9H9mHzxD6>xjC3dYwiHUEjum{m< z?5B(&EaBELHm_8hWe4f7!|}Ro&0QTf_p}aMq@u~H5{9#PR|hhsrJdQ#4kx*jTOOJz zG%=>oUC3ewGj>|YG0U?iY;3Lpn$O>yrXG7zXUO8N>oVOo>dd7oPcZd`HLIOr&H5j(V*8(uW)Zy*kxy4ySQfrYwS3Jsj8{6^l)W%IFpwWl<=Iw20d2bw1*2>GlWIF=VyG` zfo#GkRp!HUe>v;aS?QEvZ0iC|rk<|FqGpX^sZU3-?aDmIXfTR}bQ{I)){o#ha7}hA zU=TCsB`cEpUfU~TgYSa9ecuk-tyI-ckK6M+&Y)w_!yhwhP-Bs(?`fA+NVfQ9$u~h}b zn6=3emd~%4t89m{LESak&~sYs!d7io!8nUfF7_f|6`mAQnNH%4m7K?bM z&VEws%=*tL;~ngzncZ`6WS2+s+Nqbkl=m60wU7laDH7%+Zx8y?17|SyNeY;v}XM zK8}^Ck74UAZCM*18|HM$n(2gEu-jwBtTVr|Q+3s4d1X3GG-Mdd>olC#e(CeGsXg10 zFr6s`HF1BmG|KhX70$r*A~(n8GwU(Vvl}&;d!Z%Zq-7PjzFQgv>SlSS7TVRES>IjU#F`4uM0zu) z?YfccuqTVl8nli(acT|s{7V8?etrcfh+D?ByWqw}i=4P<`N`a+cf6tMJ7Z38K#fbN z{VKQ^yI0Ujx7hUSuzK#w&wIHd&E4FGk(HcO*Ai}zLOvIvoXs_ctmC4CQ@IJAkzCxs za89{kHs?Iand{nV2KRcRn3H=h;54TX;TDXM;!>Sz1QjwHP5)@=f@H4wbNeR+{-nUa z9|gX*b+TX5y0ZAy`p&j>HS!0>J$C!g+ETfW7#&Yx8QfIaYJ9 zdkDYxk5g~{baT$LynWsMg1!CyzFjVP+TXkTufn~53sfyxI!m@XfvpW8X^vmr+@{O7!_N}l%pRzh>n<|iNdF}bw>8Z!s_1UI=M+{6p! zCnI3?B@vG;qp`v%7Lzuu#l5>}h;3Mdn<}0-YB>v0TWugeZYb(54?^wYYSPqQdJG;SpE;7B;h!kUC>Miaw=b<&Y#C8RZEbfjhdm zkH-F8gYas!GIV?HCcA+Jbnf#!DpWHgBhHZ6Q-ATeF%j5U5(EDM8Q5Ty0nfNZEW5lG z=09h_x+V=CF$u_@;)Os3Ph5N@28~rlZ$9o@;VC-^vTM!C?PnTBM&c$6m0#_p}&c|hs+RO|t}@iCNaIF4?Ht`+<4D#U@JBBV)IV9c{EP|)6n zxeIsT9yJdo|kva70yfbZg zs1lD?*$S@>+cEfHCG5ZKhVP8+sN;4ay-yYHpWTW+#TA%*a~pQeDZ!cXTk!EvF4~*t zpkmMh+@3HO(N)uNk!t8(_W9&f8azFh>2`z*uJJwZ4id`NEo;dIn0O)M-qj$U*U4l~c< zbLLr`l{pHV%%hn1_Aok@)#2{EgLtA|g$*-n@GiF;GryL?Y~N-aY?F_iZfg;~a3P}G ztD$0hCdFJ>F1ArWiO!R+L$mriQmk)-HeQ8(mvaaUJBxEOPNE^a346CS;PldccD4B8@xi0}u!JsSr~7>*Sii=up)aBH`aUk6 zx`U|XD@cEH1#OB?;L`3>&^mVzZ~c#;`-H>zpi_a_y?3DSVJ3nuhw$+bUy$H|iMU_6 z4cwQ9L6M43=6Vc+AAW}HvJF1{qLbt^&$==p24;1Q#kYL zI3f-oLE*hhL}Wx_P;Xc4>Gnbt`qGvTPq3%q;&HTdiZiV|?L<#MIFn(lGuf`1MW;8< zrk$d2vi>QImj9ee)=G)=)FqzoEKVT1l6aDBn?QFKU8Mag49CZPq>ZLCDLci5t||nP zpHU!HPF_y?x&%`H2SKE}Hi+7UE}}~m(ma2BWd&M2pV}koHp!-@-*d?_w> z1?^uFM#ZIZHUL2F4iERUq-Zy-G#z8waXkV<$FE$mlJUrjfW_Q?XO zjQ&igv-iWv=M%cAE~4S(%jl77B#qHXrqjt$bOiAfV6~3kcgZE?j7&;X%cW}rH;{6- zGOC@lnVd^D(V|2+^e`lchIA^RDy>4g!!~AE8I!uV5gYCj8Z*dY;|FD{7M)5rqZB(vAJG6d09A zs)OUHWMnKI>ybz)m5H>tB$kEVofq4sZfkZ(vO-Ji|xAKXQ`?+%cl zeh29-kcU}ix#(qQUur0yPQ&y4No|flT|2Xwn#2pJ*Bx(i+7w7t!G6>!-GT0On@khk zEJVxWdi^SlvX$d#jv$p5HfEByI2^AASW!)#J&klVr-=pTv@^+ycKtGj zW`1#^E~e8d*>(o$Z?K{p53H$g8Alzu+t9CBwnwEy5&ERDD;bnw3;G_!dkeBJG}u)d!RDXr{C6IaTS+Pkha zSf@9I-_Ro8Ri-4DF^S&DO{L8aGpNvF8ttAlnNl2Q(Zu4Vl$x}kqFvNPR*mO`N7r2u z&JOu1+@I5tyox%Kd$`uOq}I{) zgnWvueJb)(2&Auy0W>mvDS0jnqpke))y`4$WNSD%+N`3kccN%mN)in{l|VbH6KJ1l z0=;QU8)=&Ni{er!%W*v!4qQib-lx#JxFnkOJb^xT%%;$yd{VJ1q(cL?(wC+KlxEGxKOT}r zzdTq^Q{A$O?6PQc*Y(uael5A&$fN_=S!8UOMSXXq(6eP}l-MDgt}RO;4ZgqK7?Dl> zKV_4%P6@?_RnWWk$7%jWbFpf6E_r_~pl4UI$i|)z@YS%9_I%2sF)uSn?L#_cIb>3& zrx_H;&xx)#vZ;1^E^XGzq_xSZv}R{2owF;Vxp7tWMBzNinuLqR_p|8yf^2%Fl|fe0 znKW(E29n>po@VaMps!{bG)60(mi1mou{^yZbJo%eqg0x5JdHkd-av1Ib7`VNDkYq# zqk890^mSC7c)U-4L=Ntc0f~J4h|)7B6csYY=6Y>3@XW=Qy8it%?$QCDOLITHgMB0V^t&_{x08aDi>CK`A`2d{SN zV0kBf=#DW&)G8BnS2Tlxh5*%v85Exiaqz1U;}wKh>Hzf7B0Q8PxY(K^Kg<~WyXm9O zX%s#WNJL=A!6MBXhLbxPo?ZpYUyIQp!xCxxM#I+03gI(s;G$!Tc16~(o@9j~v#f9` z#1h+dEivICLvpnMNri?O)=v+bXCpyTLqt0dje+>c7>q0)2d}X4NcWzIbyAbyBs&R; zj*~H9|71QE^8_qNam1sJ4#*m1k4+ctu=k`GHhMt202A0b=wkfCOYvo{-2Vw*sNNkviz{DA7JKPC0Wdi!&9FMRkRv7!U1z!gN zez>B84es->?4^dt^WGc;9-f1T>VCMf-X9t}0w7EY!0SAJ_>}uYFUTKT+j(PetOpb- zo#8%zDk5G?!RQP2$eCw@Hyed`)z$xREc=gjtOg_LvblA2s5;O@haK})3T;x z`f?|nUSID3Mp&RX7o|3K{Dn@hC3}O&-y(h>6Cn>CxbHqfz-b3fC8{ z!jZP&@QGLg<16!FH=O_D`LCJfuGnqoi1seVNR+ci+tJ3N6Hg*>end3<3s)oe4bQ9S z$HQ-7JW9648R0OU8V2o8%dsUP00vLzV8NKF=+c+NksL8j z6mp`t+GtE!xf+KL#6#qggr&VwpuIa4GW*w}lJE0#ebR6&J{5*zQgGKO0S|R!k#}Y_ zn&elaU-ELC&Gy5kZnnr)1G0>)MfGNHsYbC49=GWTFOx2)zSUh*|2yVL$M%xl4jBBrmzCnGVK2#nLgFE1;&ucPR*FYY*MU?5b^#8m3?&$(S{P#B4=1)OJ+Hk}jRG ze&01JDc?g80h>e}@7dr}!PFR0-`1C<8li~5gr zM*D{z&`%76+zfuaE)2!z=p|4}T7qppp->H2j7pgxj5_QG->%cp)6WrCW{tzt1!i#F zWq^==12N~qS9*Nr6v>R)AW9DO#5i4l9QtuF#{L)vX{~UyyBm%#a;sotq)gp0&oOvtq33I}$sP%@{QgDydCzS^q_7Dc4I1Hm;V^1V>QSe*58}^zFf)HUFk&kX^0K=T zvFk9%Effi(Jz;Zg7G?=gyUEia10F} zk6?QJA$)Q8}NQ?OxWG`iSE z;pn^R5N~k6u2(}r)pyB1r<@#qaTjw9O{mv7j!n$MVzpx2Y>Km}_ z_Cd5O+=VymE73(g3ybfq!HKFQ36|vGC=3tT~*H*H4$@mop(yI%bJ4 zW80&_ppGmQ9mSuD4ax9 zWhkAV4`abb%()eV!tvo4?c{`kkGsP9_%S*$lEld;>u`pb7)-Z5iTVC#(RIps_{#l= ztrh37eBn8)UUC)#mo}m(v;l3xccA9hW@tPrMWAj98h6BDviTf5DjxtB#Vd4eioW>N z=OAooTm*}AD^O$*fp?)Xm|DCVGp4LYyUb|*x+wVSEI>+T0G|Fd8jDWxYu`@>IFZ`} z7jMer{l)uK8N8LISgxnF+op>C%5YXILVZ)w-z9W>NkU)+20 zCJfwMg4yCNh-}z~a;+WcxqK&1$m~MM=N;&pumf#-ZGmm(X2?9qMC{EJh+ifnuumwu z3<*HaQ*-3^ZjTT14p7JhZLxCu62#}1<7LHm45{CRo(_}Lpg~9LkL|i)76|Yl{k!<8(u{fd>eitgB-exym7FEMH zVLv{{)?s2!Jr*U@BVv0UE<4piQ@0u&ua{wQ;wI$P6rd(Q5l^d*h7;Bn48L7K@UQN_*V_Ydt!ncHfeXVIoj zvtgQ4n{TQ_i3(7%m6310xmoepO5zSuzqQB<`-??hSgS?eiO|+7=QU@`Cy7O#%-^@j zYYrdF=P$=AMLS6O`1^S}`?>h?#&9jE*%GihjA?7Gf7^C0*n5$?3Lm*so)7WbY=zgH zPX{SCe_vj|8q#7*DDmggoTb=3*gYsjW!Bu_5EbuWm08~Nyl1)3Qt|c+@mC4)3UcS; zResAa*hEEJOIvS*wwAuWCcmq%&+i(mxcJTDf$<1kZ53Z{x1eVIZQzzqGn(^mP7+@l ztw>5LyjVJLtNC5w(dCbYN`2bUjLw~Dhgh3>zb5kcT}Z|j@w8{-M#}d2Oh?vsz}TMV zNNt*ngC`dwZ~00%dBmd3F%wRe%TTE#4dsZDqOk#esCrTdDzdpPOqMw&G<>j4xXz|b z7?-Vgvk}v2 z9Oj0er1x2;Xr|pMnlPq``WhCJ3tz4Vl>pkY!;%`5I2x|do@U&7Cv0PQKzLJTk1%f3 zZegu?r4S>kg$>dVgt>;3Da30p%^Wfrew(_{ge%?X%cK6}E7GDB^Cwe}x;a$PJ(YxK zN+~+3j=X%j;)}unh;=6+erq^te~84@m67QFFcAS`*5mV*ztp!R#){NB=u-6g&b;d3 znb171Ntow;K**Mq32ojK37Z~l6prWi2|t^i5;n$IQ|aVj${xR%ZdW`f)dXqWt7g0i zA>!Bie(Kd7ajYvkTC_(_ySKFW%t5L%*hoDZqNwhoE1mZoOLBp|Y0hm0TDt0Ip{CwX zLM6G2!YABM!XWvZ!c+%2np*1hYkeOs4y5cg9(4GLGj(6%PK)ItY3PG^YTRB%y&4YD z$xB!1Rhbf2J7_>=ZU|OBN`{TpT673a#`F6*xN|ZcC(UEt8`< zu7%LOOD2>!petRDlO|u`RiV=2bHZGw*$1bAEB?Cp1tK>yjz0Tr5@&;V#?TCyC{qSSsT`KV0LDIAGXeG0u zx{oF#A2WdV%gfS}*Be>=~zPMk4(DWVLOxlVTi`$7jGL=PpMMFiN?qfyc z<%&d!KX(y_=BinC?U0EH1sRwhkpa7l8QA8Z1=rCz7<48VufF8NMy3!ihHgix#coud zJdG(H*U+GP1E(BiMe#nJMc!&!q7O%6M7rZ1ieg{1v;5^lDf+~eV$AGP=pHY{)$%fo zzO)7V&Thlnn4NfYXE#2X9zjz3<9Jy601?q2a3=N>9<1&yayX(aDy}gX2@YqAQn~ix zw`coW{?Xze$z1d2_D>4@Nr8Vq3ShjIy_D<^l7hwardq?LDw|g&eCjI5Z(ZH=GU~ch zruEpL|Noz&0LI%|NXd4S)Mbm6!@onCww{)fj_KSD9BsucGVf#J1x_W7AiJQ5a<|NN^4$bD4f6cXU{8O9-970>i1uL{hqIs znE8Gqt<*>5oi^V>!b`Uq+I+N)bSupyRF|~l@A3cXH+LOcCNP;cx1O{`Gc$vi)bquB6AvK5ZqVBo9fC`3F*=mc`LOk_rv~R*UAG|JE0kLGB(w zE^Z;-^W1sa&Mfx@D!zPU`?N}_UT$e^*K@7KlKSIK-uWp*O1o$L{O@Oe$!7&dlBX_c zHnRVB^$^_oTPY>m`S#xqeB0w48ny2v{;CH`l`7NGQ`CDwY?YCIs<5qI?Kl=*y zKbsFvu9J-2)*Nq@{1#gnc1~zL^ZP3oPqSX}n|6EBdYxq6{NsCz7yp3LjHPwIdFr?B z^6NU`f2BX^R?12J>Q7pNCH|x(SmIAwf+hZ>C0OE5T7o4Wk|kK;5m|yI9*`wi;_+C5 zB_G8|@Q}ja_KvfcLPMKdgEt`HN9f;Loo(W8mE5Bv5YKh>pXWVC#XUsD#YaU~TU%Gf z!hdeyTo+!{+*MOW}Q595siH3c z{~BpoEx Date: Tue, 3 Mar 2026 19:59:55 +0100 Subject: [PATCH 05/13] add bw solver to _rt executable. All raytracing work in one executable is a lot easier --- include_test/radiation_solver_bw.h | 57 +---- src_test/CMakeLists.txt | 6 +- src_test/radiation_solver_bw.cu | 19 +- src_test/test_rte_rrtmgp_bw.cu | 4 +- src_test/test_rte_rrtmgp_rt.cu | 377 ++++++++++++++++++++++++++++- 5 files changed, 384 insertions(+), 79 deletions(-) diff --git a/include_test/radiation_solver_bw.h b/include_test/radiation_solver_bw.h index 6a31a76c..a57677dc 100644 --- a/include_test/radiation_solver_bw.h +++ b/include_test/radiation_solver_bw.h @@ -34,35 +34,14 @@ -class Radiation_solver_longwave +class Radiation_solver_bw_longwave { public: - Radiation_solver_longwave( - const Gas_concs& gas_concs, - const std::string& file_name_gas, - const std::string& file_name_cloud); - - Radiation_solver_longwave( + Radiation_solver_bw_longwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud); - void solve( - const bool switch_fluxes, - const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, - const Gas_concs& gas_concs, - const Array& p_lay, const Array& p_lev, - const Array& t_lay, const Array& t_lev, - const Array& col_dry, - const Array& t_sfc, const Array& emis_sfc, - const Array& lwp, const Array& iwp, - const Array& rel, const Array& dei, - Array& tau, Array& lay_source, - Array& lev_source_inc, Array& lev_source_dec, Array& sfc_source, - Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net, - Array& lw_bnd_flux_up, Array& lw_bnd_flux_dn, Array& lw_bnd_flux_net) const; #ifdef __CUDACC__ void solve_gpu( @@ -108,42 +87,15 @@ class Radiation_solver_longwave }; -class Radiation_solver_shortwave +class Radiation_solver_bw_shortwave { public: - Radiation_solver_shortwave( - const Gas_concs& gas_concs, - const std::string& file_name_gas, - const std::string& file_name_cloud, - const std::string& file_name_aerosol); - Radiation_solver_shortwave( + Radiation_solver_bw_shortwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud, const std::string& file_name_aerosol); - void solve( - const bool switch_fluxes, - const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, - const bool switch_delta_cloud, - const bool switch_delta_aerosol, - const Gas_concs& gas_concs, - const Array& p_lay, const Array& p_lev, - const Array& t_lay, const Array& t_lev, - const Array& col_dry, - const Array& sfc_alb, - const Array& tsi_scaling, const Array& mu0, - const Array& lwp, const Array& iwp, - const Array& rel, const Array& dei, - Array& tau, Array& ssa, Array& g, - Array& toa_src, - Array& sw_flux_up, Array& sw_flux_dn, - Array& sw_flux_dn_dir, Array& sw_flux_net, - Array& sw_bnd_flux_up, Array& sw_bnd_flux_dn, - Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const; - void load_mie_tables( const std::string& file_name_mie_bb, const std::string& file_name_mie_im, @@ -152,7 +104,6 @@ class Radiation_solver_shortwave #ifdef __CUDACC__ void solve_gpu( - const bool tune_step, const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, diff --git a/src_test/CMakeLists.txt b/src_test/CMakeLists.txt index 16342d04..6ad4a779 100644 --- a/src_test/CMakeLists.txt +++ b/src_test/CMakeLists.txt @@ -28,11 +28,11 @@ if(USECUDA) add_executable(test_rte_rrtmgp radiation_solver.cu test_rte_rrtmgp.cu) target_link_libraries(test_rte_rrtmgp rte_rrtmgp rte_rrtmgp_cuda ${LIBS} m) - add_executable(test_rte_rrtmgp_rt radiation_solver_rt.cu test_rte_rrtmgp_rt.cu ../src_tilt/tilt_utils.cpp) + add_executable(test_rte_rrtmgp_rt radiation_solver_rt.cu test_rte_rrtmgp_rt.cu radiation_solver_bw.cu ../src_tilt/tilt_utils.cpp) target_link_libraries(test_rte_rrtmgp_rt rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) - add_executable(test_rte_rrtmgp_bw radiation_solver_bw.cu test_rte_rrtmgp_bw.cu) - target_link_libraries(test_rte_rrtmgp_bw rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) + #add_executable(test_rte_rrtmgp_bw radiation_solver_bw.cu test_rte_rrtmgp_bw.cu) + #target_link_libraries(test_rte_rrtmgp_bw rte_rrtmgp rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) add_executable(test_rt_lite test_rt_lite.cu) target_link_libraries(test_rt_lite rte_rrtmgp_cuda rte_rrtmgp_cuda_rt curand ${LIBS} m) diff --git a/src_test/radiation_solver_bw.cu b/src_test/radiation_solver_bw.cu index 94a4be71..2f6570b0 100644 --- a/src_test/radiation_solver_bw.cu +++ b/src_test/radiation_solver_bw.cu @@ -477,7 +477,7 @@ namespace } -Radiation_solver_longwave::Radiation_solver_longwave( +Radiation_solver_bw_longwave::Radiation_solver_bw_longwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud) @@ -491,7 +491,7 @@ Radiation_solver_longwave::Radiation_solver_longwave( } -void Radiation_solver_longwave::solve_gpu( +void Radiation_solver_bw_longwave::solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_output_optical, @@ -508,7 +508,7 @@ void Radiation_solver_longwave::solve_gpu( Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, Array_gpu& lw_bnd_flux_up, Array_gpu& lw_bnd_flux_dn, Array_gpu& lw_bnd_flux_net) { - throw std::runtime_error("Longwave raytracing is not implemented"); + throw std::runtime_error("Longwave backward raytracing is not implemented"); /* const int n_col = p_lay.dim(1); @@ -720,7 +720,7 @@ Float xyz_irradiance( } -Radiation_solver_shortwave::Radiation_solver_shortwave( +Radiation_solver_bw_shortwave::Radiation_solver_bw_shortwave( const Gas_concs_gpu& gas_concs, const std::string& file_name_gas, const std::string& file_name_cloud, @@ -737,7 +737,7 @@ Radiation_solver_shortwave::Radiation_solver_shortwave( load_and_init_aerosol_optics(file_name_aerosol)); } -void Radiation_solver_shortwave::load_mie_tables( +void Radiation_solver_bw_shortwave::load_mie_tables( const std::string& file_name_mie_bb, const std::string& file_name_mie_im, const bool switch_broadband, @@ -790,8 +790,7 @@ void Radiation_solver_shortwave::load_mie_tables( } -void Radiation_solver_shortwave::solve_gpu( - const bool tune_step, +void Radiation_solver_bw_shortwave::solve_gpu( const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, @@ -886,7 +885,7 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu albedo; if (!switch_lu_albedo) albedo = sfc_alb.subset({{ {band, band}, {1, n_col}}}); - if (!tune_step && (! (band == 10 || band == 11 || band ==12))) continue; + if (! (band == 10 || band == 11 || band ==12)) continue; const Float solar_source_band = kdist_gpu->band_source(band_limits_gpt({1,band}), band_limits_gpt({2,band})); @@ -929,8 +928,6 @@ void Radiation_solver_shortwave::solve_gpu( gas_optics_subset(col_s, n_col_residual); } - if (tune_step) return; - constexpr bool do_scattering = true; toa_src.fill(toa_src_temp({1}) * tsi_scaling({1})); @@ -1112,7 +1109,7 @@ void Radiation_solver_shortwave::solve_gpu( } } -void Radiation_solver_shortwave::solve_gpu_bb( +void Radiation_solver_bw_shortwave::solve_gpu_bb( const bool switch_cloud_optics, const bool switch_cloud_mie, const bool switch_aerosol_optics, diff --git a/src_test/test_rte_rrtmgp_bw.cu b/src_test/test_rte_rrtmgp_bw.cu index d712606d..ccd8e346 100644 --- a/src_test/test_rte_rrtmgp_bw.cu +++ b/src_test/test_rte_rrtmgp_bw.cu @@ -385,7 +385,7 @@ void solve_radiation(int argc, char** argv) Gas_concs_gpu gas_concs_gpu(gas_concs); - Radiation_solver_longwave rad_lw(gas_concs_gpu, "coefficients_lw.nc", "cloud_coefficients_lw.nc"); + Radiation_solver_bw_longwave rad_lw(gas_concs_gpu, "coefficients_lw.nc", "cloud_coefficients_lw.nc"); // Read the boundary conditions. const int n_bnd_lw = rad_lw.get_n_bnd_gpu(); @@ -572,7 +572,7 @@ void solve_radiation(int argc, char** argv) Gas_concs_gpu gas_concs_gpu(gas_concs); - Radiation_solver_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc","aerosol_optics_sw.nc"); + Radiation_solver_bw_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc","aerosol_optics_sw.nc"); // Read the boundary conditions. const int n_bnd_sw = rad_sw.get_n_bnd_gpu(); diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index b616fc66..9d6a3863 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -28,6 +28,8 @@ #include "raytracer_kernels_sw.h" #include "radiation_solver_rt.h" #include "aerosol_optics_rt.h" +#include "raytracer_kernels_bw.h" +#include "radiation_solver_bw.h" #include "gas_concs.h" #include "types.h" #include "mem_pool_gpu.h" @@ -185,6 +187,13 @@ void solve_radiation(int argc, char** argv) const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); const bool switch_tica = get_ini_value(settings, "switches", "tica", false); + const bool switch_bw_raytracing = get_ini_value(settings, "switches", "bw-raytracing", true); + const bool switch_lu_albedo = get_ini_value(settings, "switches", "lu-albedo", false); + const bool switch_image = get_ini_value(settings, "switches", "image", true); + const bool switch_broadband = get_ini_value(settings, "switches", "broadband", false); + const bool switch_cloud_cam = get_ini_value(settings, "switches", "cloud-cam", false); + + const int photons_per_pixel = get_ini_value(settings, "ints", "bw-raytracing", 1); if (!switch_shortwave) switch_sw_raytracing = false; @@ -209,6 +218,7 @@ void solve_radiation(int argc, char** argv) } } + Int lw_photon_power; Int lw_photon_count; if (switch_lw_raytracing) @@ -217,6 +227,11 @@ void solve_radiation(int argc, char** argv) lw_photon_count = 1 << lw_photon_power; } + if (switch_longwave && switch_bw_raytracing) + { + Status::print_warning("No longwave radiation implemented in the backward ray tracer"); + } + if (switch_cloud_optics) { switch_liq_cloud_optics = true; @@ -228,6 +243,9 @@ void solve_radiation(int argc, char** argv) if (switch_tica) switch_independent_column = true; + if (switch_cloud_cam && !(switch_cloud_optics)) + throw std::runtime_error("Enable cloud optics (liquid and/or ice) when cloud-cam is switch on!"); + if (switch_cloud_mie && switch_ice_cloud_optics) { std::string error = "Thou shall not use mie tables as long as ice optics are enabled"; @@ -272,20 +290,40 @@ void solve_radiation(int argc, char** argv) Array grid_y(input_nc.get_variable("y", {n_col_y}), {n_col_y}); Array grid_yh(input_nc.get_variable("yh", {n_col_y+1}), {n_col_y+1}); Array grid_z(input_nc.get_variable("z", {n_z_in}), {n_z_in}); + Array z_lev(input_nc.get_variable("z_lev", {n_lev}), {n_lev}); const Vector grid_d = {grid_xh({2}) - grid_xh({1}), grid_yh({2}) - grid_yh({1}), grid_z({2}) - grid_z({1})}; const Vector kn_grid = {input_nc.get_variable("ngrid_x"), input_nc.get_variable("ngrid_y"), input_nc.get_variable("ngrid_z")}; + Camera camera; + if (switch_bw_raytracing) + { + camera.fov = get_ini_value(settings, "bw_camera", "field-of-view", 80.); + camera.cam_type = get_ini_value(settings, "bw_camera", "cam-type", 0); + camera.position = {get_ini_value(settings, "bw_camera", "px", 0.), + get_ini_value(settings, "bw_camera", "py", 0.), + get_ini_value(settings, "bw_camera", "pz", 100.)}; + camera.nx = get_ini_value(settings, "bw_camera", "nx", 0); + camera.ny = get_ini_value(settings, "bw_camera", "ny", 0); + camera.setup_rotation_matrix(get_ini_value(settings, "bw_camera", "yaw", 0.), + get_ini_value(settings, "bw_camera", "pitch", 0.), + get_ini_value(settings, "bw_camera", "roll", 0.)); + camera.setup_normal_camera(camera); + + camera.npix = Int(camera.nx * camera.ny); + + const Float input_sza = get_ini_value(settings, "floats", "sza", -1.0); + const Float input_azi = get_ini_value(settings, "floats", "azi", -1.0); + } + // Read the atmospheric fields. Array p_lay(input_nc.get_variable("p_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); Array t_lay(input_nc.get_variable("t_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); Array p_lev(input_nc.get_variable("p_lev", {n_lev, n_col_y, n_col_x}), {n_col, n_lev}); Array t_lev(input_nc.get_variable("t_lev", {n_lev, n_col_y, n_col_x}), {n_col, n_lev}); - - if (input_nc.variable_exists("col_dry") && switch_tica) { std::string error = "col_dry is not supported in tica mode"; @@ -300,6 +338,19 @@ void solve_radiation(int argc, char** argv) } + // read land use map if present, used for choosing between spectral and lambertian reflection and for spectral albedo + // 0: water, 1: "grass", 2: "soil", 3: "concrete". Interpolating between 1 and 2 is currently possible + Array land_use_map({n_col}); + if (input_nc.variable_exists("land_use_map") && switch_lu_albedo) + { + land_use_map = std::move(input_nc.get_variable("land_use_map", {n_col_y, n_col_x})); + } + else + { + // default to grass with some soil + land_use_map.fill(Float(1.3)); + } + // Fetch the col_dry in case present. Array col_dry; if (input_nc.variable_exists("col_dry")) @@ -358,6 +409,11 @@ void solve_radiation(int argc, char** argv) dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); } } + else + { + rel.set_dims({n_col, n_lay}); + rel.fill(Float(0.)); + } Array rh; Aerosol_concs aerosol_concs; @@ -507,6 +563,7 @@ void solve_radiation(int argc, char** argv) Netcdf_file output_nc(case_name + "_output.nc", Netcdf_mode::Create); output_nc.add_dimension("col", n_col); + output_nc.add_dimension("x", n_col_x); output_nc.add_dimension("y", n_col_y); output_nc.add_dimension("lay", n_lay); @@ -534,13 +591,27 @@ void solve_radiation(int argc, char** argv) Netcdf_file coef_nc_lw("coefficients_lw.nc", Netcdf_mode::Read); nbnds = std::max(coef_nc_lw.get_dimension_size("bnd"), nbnds); ngpts = std::max(coef_nc_lw.get_dimension_size("gpt"), ngpts); + + output_nc.add_dimension("gpt_lw", ngpts); + output_nc.add_dimension("band_lw", nbnds); + } - if (switch_shortwave) + if (switch_shortwave || switch_bw_raytracing) { Netcdf_file coef_nc_sw("coefficients_sw.nc", Netcdf_mode::Read); nbnds = std::max(coef_nc_sw.get_dimension_size("bnd"), nbnds); ngpts = std::max(coef_nc_sw.get_dimension_size("gpt"), ngpts); + + output_nc.add_dimension("gpt_sw", ngpts); + output_nc.add_dimension("band_sw", nbnds); } + + if (switch_bw_raytracing) + { + output_nc.add_dimension("px", camera.nx); + output_nc.add_dimension("py", camera.ny); + } + configure_memory_pool(n_lay, n_col, 1024, ngpts, nbnds); @@ -750,9 +821,6 @@ void solve_radiation(int argc, char** argv) Array rt_flux_sfc_dn_cpu(rt_flux_sfc_dn); Array rt_flux_abs_cpu(rt_flux_abs); - output_nc.add_dimension("gpt_lw", n_gpt_lw); - output_nc.add_dimension("band_lw", n_bnd_lw); - auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber_gpu().v(), {0, 0}); @@ -844,7 +912,6 @@ void solve_radiation(int argc, char** argv) // Initialize the solver. Status::print_message("Initializing the shortwave solver."); - Gas_concs_gpu gas_concs_gpu(gas_concs); Radiation_solver_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc", "aerosol_optics_sw.nc"); @@ -1081,9 +1148,6 @@ void solve_radiation(int argc, char** argv) Array rt_flux_abs_dir_cpu(rt_flux_abs_dir); Array rt_flux_abs_dif_cpu(rt_flux_abs_dif); - output_nc.add_dimension("gpt_sw", n_gpt_sw); - output_nc.add_dimension("band_sw", n_bnd_sw); - auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); @@ -1226,6 +1290,299 @@ void solve_radiation(int argc, char** argv) } } + ////// RUN THE BACKWARD SOLVER ////// + if (switch_bw_raytracing || switch_cloud_cam) + { + // Initialize the solver. + Status::print_message("Initializing the shortwave backward solver."); + + Gas_concs_gpu gas_concs_gpu(gas_concs); + Radiation_solver_bw_shortwave rad_sw(gas_concs_gpu, "coefficients_sw.nc", "cloud_coefficients_sw.nc","aerosol_optics_sw.nc"); + + // Read the boundary conditions. + const int n_bnd_sw = rad_sw.get_n_bnd_gpu(); + const int n_gpt_sw = rad_sw.get_n_gpt_gpu(); + + Array mu0(input_nc.get_variable("mu0", {n_col_y, n_col_x}), {n_col}); + Array azi(input_nc.get_variable("azi", {n_col_y, n_col_x}), {n_col}); + + if (input_sza < 0) + mu0.fill(cos(input_sza / Float(180.0) * M_PI)); + + if (input_azi < 0) + azi.fill(input_azi / Float(180.0) * M_PI); + + Array sfc_alb(input_nc.get_variable("sfc_alb_dir", {n_col_y, n_col_x, n_bnd_sw}), {n_bnd_sw, n_col}); + + Array tsi_scaling({n_col}); + if (input_nc.variable_exists("tsi")) + { + Array tsi(input_nc.get_variable("tsi", {n_col_y, n_col_x}), {n_col}); + const Float tsi_ref = rad_sw.get_tsi_gpu(); + for (int icol=1; icol<=n_col; ++icol) + tsi_scaling({icol}) = tsi({icol}) / tsi_ref; + } + else if (input_nc.variable_exists("tsi_scaling")) + { + Float tsi_scaling_in = input_nc.get_variable("tsi_scaling"); + for (int icol=1; icol<=n_col; ++icol) + tsi_scaling({icol}) = tsi_scaling_in; + } + else + { + for (int icol=1; icol<=n_col; ++icol) + tsi_scaling({icol}) = Float(1.); + } + + Array_gpu XYZ; + Array_gpu radiance; + + if (switch_broadband) + { + radiance.set_dims({camera.nx, camera.ny}); + } + if (switch_image) + { + XYZ.set_dims({camera.nx, camera.ny, 3}); + } + + if (switch_cloud_mie) + rad_sw.load_mie_tables("mie_lut_broadband.nc", "mie_lut_visualisation.nc", switch_broadband, switch_image); + + Array_gpu liwp_cam; + Array_gpu tauc_cam; + Array_gpu dist_cam; + Array_gpu zen_cam; + + if (switch_cloud_cam) + { + liwp_cam.set_dims({camera.nx, camera.ny}); + tauc_cam.set_dims({camera.nx, camera.ny}); + dist_cam.set_dims({camera.nx, camera.ny}); + zen_cam.set_dims({camera.nx, camera.ny}); + } + + const Vector grid_cells = {n_col_x, n_col_y, n_z_in}; + + auto run_solver_bb = [&]() + { + Array_gpu p_lay_gpu(p_lay); + Array_gpu p_lev_gpu(p_lev); + Array_gpu t_lay_gpu(t_lay); + Array_gpu t_lev_gpu(t_lev); + Array_gpu z_lev_gpu(z_lev); + Array_gpu col_dry_gpu(col_dry); + Array_gpu sfc_alb_gpu(sfc_alb); + Array_gpu tsi_scaling_gpu(tsi_scaling); + Array_gpu mu0_gpu(mu0); + Array_gpu azi_gpu(azi); + Array_gpu lwp_gpu(lwp); + Array_gpu iwp_gpu(iwp); + Array_gpu rel_gpu(rel); + Array_gpu dei_gpu(dei); + Array_gpu rh_gpu(rh); + Aerosol_concs_gpu aerosol_concs_gpu(aerosol_concs); + + Array_gpu land_use_map_gpu(land_use_map); + + cudaDeviceSynchronize(); + cudaEvent_t start; + cudaEvent_t stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start, 0); + + rad_sw.solve_gpu_bb( + switch_cloud_optics, + switch_cloud_mie, + switch_aerosol_optics, + switch_lu_albedo, + switch_delta_cloud, + switch_delta_aerosol, + switch_cloud_cam, + switch_bw_raytracing, + grid_cells, + grid_d, + kn_grid, + photons_per_pixel, + gas_concs_gpu, + p_lay_gpu, p_lev_gpu, + t_lay_gpu, t_lev_gpu, + z_lev_gpu, + col_dry_gpu, + sfc_alb_gpu, + tsi_scaling_gpu, + mu0_gpu, azi_gpu, + lwp_gpu, iwp_gpu, + rel_gpu, dei_gpu, + land_use_map_gpu, + rh_gpu, + aerosol_concs, + camera, + radiance, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); + + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float duration = 0.f; + cudaEventElapsedTime(&duration, start, stop); + + cudaEventDestroy(start); + cudaEventDestroy(stop); + + Status::print_message("Duration shortwave backward solver (broadband version): " + std::to_string(duration) + " (ms)"); + }; + + auto run_solver = [&]() + { + Array_gpu p_lay_gpu(p_lay); + Array_gpu p_lev_gpu(p_lev); + Array_gpu t_lay_gpu(t_lay); + Array_gpu t_lev_gpu(t_lev); + Array_gpu z_lev_gpu(z_lev); + Array_gpu col_dry_gpu(col_dry); + Array_gpu sfc_alb_gpu(sfc_alb); + Array_gpu tsi_scaling_gpu(tsi_scaling); + Array_gpu mu0_gpu(mu0); + Array_gpu azi_gpu(azi); + Array_gpu lwp_gpu(lwp); + Array_gpu iwp_gpu(iwp); + Array_gpu rel_gpu(rel); + Array_gpu dei_gpu(dei); + + Array_gpu rh_gpu(rh); + Aerosol_concs_gpu aerosol_concs_gpu(aerosol_concs); + + Array_gpu land_use_map_gpu(land_use_map); + + cudaDeviceSynchronize(); + cudaEvent_t start; + cudaEvent_t stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start, 0); + + rad_sw.solve_gpu( + switch_cloud_optics, + switch_cloud_mie, + switch_aerosol_optics, + switch_lu_albedo, + switch_delta_cloud, + switch_delta_aerosol, + switch_cloud_cam, + switch_bw_raytracing, + grid_cells, + grid_d, + kn_grid, + photons_per_pixel, + gas_concs_gpu, + p_lay_gpu, p_lev_gpu, + t_lay_gpu, t_lev_gpu, + z_lev_gpu, + col_dry_gpu, + sfc_alb_gpu, + tsi_scaling_gpu, + mu0_gpu, azi_gpu, + lwp_gpu, iwp_gpu, + rel_gpu, dei_gpu, + land_use_map_gpu, + rh_gpu, + aerosol_concs, + camera, + XYZ, + liwp_cam, + tauc_cam, + dist_cam, + zen_cam); + + cudaEventRecord(stop, 0); + cudaEventSynchronize(stop); + float duration = 0.f; + cudaEventElapsedTime(&duration, start, stop); + + cudaEventDestroy(start); + cudaEventDestroy(stop); + + + Status::print_message("Duration shortwave backward solver (image version): " + std::to_string(duration) + " (ms)"); + }; + + if (switch_broadband) + { + // Profiling step; + run_solver_bb(); + } + if (switch_image) + { + // actual solve + run_solver(); + } + + // Store the output. + Status::print_message("Storing the backward output."); + + if (switch_bw_raytracing) + { + if (switch_broadband) + { + Array radiance_cpu(radiance); + + auto nc_var = output_nc.add_variable("radiance", {"py","px"}); + nc_var.insert(radiance_cpu.v(), {0, 0}); + nc_var.add_attribute("long_name", "shortwave radiance"); + nc_var.add_attribute("units", "W m-2 sr-1"); + } + if (switch_image) + { + Array xyz_cpu(XYZ); + output_nc.add_dimension("n",3); + + auto nc_xyz = output_nc.add_variable("XYZ", {"n","py","px"}); + nc_xyz.insert(xyz_cpu.v(), {0, 0, 0}); + + nc_xyz.add_attribute("long_name", "X Y Z tristimulus values"); + } + } + + if (switch_cloud_cam) + { + + Array liwp_cam_cpu(liwp_cam); + Array tauc_cam_cpu(tauc_cam); + Array dist_cam_cpu(dist_cam); + Array zen_cam_cpu(zen_cam); + + auto nc_var_liwp = output_nc.add_variable("liq_ice_wp_cam", {"py","px"}); + nc_var_liwp.insert(liwp_cam_cpu.v(), {0, 0}); + nc_var_liwp.add_attribute("long_name", "accumulated liquid+ice water path"); + + auto nc_var_tauc = output_nc.add_variable("tau_cld_cam", {"py","px"}); + nc_var_tauc.insert(tauc_cam_cpu.v(), {0, 0}); + nc_var_tauc.add_attribute("long_name", "accumulated cloud optical depth (441-615nm band)"); + + auto nc_var_dist = output_nc.add_variable("dist_cld_cam", {"py","px"}); + nc_var_dist.insert(dist_cam_cpu.v(), {0, 0}); + nc_var_dist.add_attribute("long_name", "distance to first cloudy cell"); + + auto nc_var_csza = output_nc.add_variable("zen_cam", {"py","px"}); + nc_var_csza.insert(zen_cam_cpu.v(), {0, 0}); + nc_var_csza.add_attribute("long_name", "zenith angle of camera pixel"); + } + + auto nc_mu0 = output_nc.add_variable("sza"); + nc_mu0.insert(acos(mu0({1}))/M_PI * Float(180.), {0}); + + auto nc_azi = output_nc.add_variable("azi"); + nc_azi.insert(azi({1})/M_PI * Float(180.), {0}); + + + } + Status::print_message("###### Finished RTE+RRTMGP solver ######"); } From 178b056ab812df78b1bde4227354726da59deaf2 Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 3 Mar 2026 20:27:14 +0100 Subject: [PATCH 06/13] remove profiling, output-optical and output-bnd-fluxes switches. Also remove a few comment lines containing old code --- include_test/radiation_solver.h | 30 +---- src_test/radiation_solver.cpp | 200 +++++--------------------------- src_test/radiation_solver.cu | 187 ++--------------------------- src_test/test_rte_rrtmgp.cpp | 139 +--------------------- src_test/test_rte_rrtmgp.cu | 182 +---------------------------- src_test/test_rte_rrtmgp_rt.cu | 17 --- 6 files changed, 50 insertions(+), 705 deletions(-) diff --git a/include_test/radiation_solver.h b/include_test/radiation_solver.h index 7d311270..ded98f77 100644 --- a/include_test/radiation_solver.h +++ b/include_test/radiation_solver.h @@ -46,8 +46,6 @@ class Radiation_solver_longwave void solve( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs& gas_concs, const Array& p_lay, const Array& p_lev, const Array& t_lay, const Array& t_lev, @@ -55,10 +53,7 @@ class Radiation_solver_longwave const Array& t_sfc, const Array& emis_sfc, const Array& lwp, const Array& iwp, const Array& rel, const Array& dei, - Array& tau, Array& lay_source, - Array& lev_source, Array& sfc_source, - Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net, - Array& lw_bnd_flux_up, Array& lw_bnd_flux_dn, Array& lw_bnd_flux_net) const; + Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net) const; int get_n_gpt() const { return this->kdist->get_ngpt(); }; int get_n_bnd() const { return this->kdist->get_nband(); }; @@ -73,8 +68,6 @@ class Radiation_solver_longwave void solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs_gpu& gas_concs, const Array_gpu& p_lay, const Array_gpu& p_lev, const Array_gpu& t_lay, const Array_gpu& t_lev, @@ -82,10 +75,7 @@ class Radiation_solver_longwave const Array_gpu& t_sfc, const Array_gpu& emis_sfc, const Array_gpu& lwp, const Array_gpu& iwp, const Array_gpu& rel, const Array_gpu& dei, - Array_gpu& tau, Array_gpu& lay_source, - Array_gpu& lev_source, Array_gpu& sfc_source, - Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, - Array_gpu& lw_bnd_flux_up, Array_gpu& lw_bnd_flux_dn, Array_gpu& lw_bnd_flux_net); + Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net); int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); }; int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); }; @@ -140,8 +130,6 @@ class Radiation_solver_shortwave const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs& gas_concs, @@ -154,12 +142,8 @@ class Radiation_solver_shortwave const Array& rel, const Array& dei, const Array& rh, const Aerosol_concs& aerosol_concs, - Array& tau, Array& ssa, Array& g, - Array& toa_src, Array& sw_flux_up, Array& sw_flux_dn, - Array& sw_flux_dn_dir, Array& sw_flux_net, - Array& sw_bnd_flux_up, Array& sw_bnd_flux_dn, - Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const; + Array& sw_flux_dn_dir, Array& sw_flux_net) const; int get_n_gpt() const { return this->kdist->get_ngpt(); }; int get_n_bnd() const { return this->kdist->get_nband(); }; @@ -177,8 +161,6 @@ class Radiation_solver_shortwave const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs_gpu& gas_concs, @@ -191,12 +173,8 @@ class Radiation_solver_shortwave const Array_gpu& rel, const Array_gpu& dei, const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, - Array_gpu& tau, Array_gpu& ssa, Array_gpu& g, - Array_gpu& toa_src, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, - Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_bnd_flux_up, Array_gpu& sw_bnd_flux_dn, - Array_gpu& sw_bnd_flux_dn_dir, Array_gpu& sw_bnd_flux_net); + Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net); int get_n_gpt_gpu() const { return this->kdist_gpu->get_ngpt(); }; int get_n_bnd_gpu() const { return this->kdist_gpu->get_nband(); }; diff --git a/src_test/radiation_solver.cpp b/src_test/radiation_solver.cpp index 7ae67f91..2d2ebc69 100644 --- a/src_test/radiation_solver.cpp +++ b/src_test/radiation_solver.cpp @@ -384,8 +384,6 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs& gas_concs, const Array& p_lay, const Array& p_lev, const Array& t_lay, const Array& t_lev, @@ -393,10 +391,7 @@ void Radiation_solver_longwave::solve( const Array& t_sfc, const Array& emis_sfc, const Array& lwp, const Array& iwp, const Array& rel, const Array& dei, - Array& tau, Array& lay_source, - Array& lev_source, Array& sfc_source, - Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net, - Array& lw_bnd_flux_up, Array& lw_bnd_flux_dn, Array& lw_bnd_flux_net) const + Array& lw_flux_up, Array& lw_flux_dn, Array& lw_flux_net) const { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -445,8 +440,7 @@ void Radiation_solver_longwave::solve( std::unique_ptr& cloud_optical_props_subset_in, Source_func_lw& sources_subset_in, const Array& emis_sfc_subset_in, - Fluxes_broadband& fluxes, - Fluxes_broadband& bnd_fluxes) + Fluxes_broadband& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -487,44 +481,14 @@ void Radiation_solver_longwave::solve( dynamic_cast(*cloud_optical_props_subset_in)); } - // Store the optical properties, if desired. - if (switch_output_optical) - { - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int ilay=1; ilay<=n_lay; ++ilay) - for (int icol=1; icol<=n_col_in; ++icol) - { - tau ({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_tau()({icol, ilay, igpt}); - lay_source({icol+col_s_in-1, ilay, igpt}) = sources_subset_in.get_lay_source()({icol, ilay, igpt}); - lev_source({icol+col_s_in-1, ilay, igpt}) = sources_subset_in.get_lev_source()({icol, ilay, igpt}); - } - - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int icol=1; icol<=n_col_in; ++icol) - lev_source({icol+col_s_in-1, n_lev, igpt}) = sources_subset_in.get_lev_source()({icol, n_lev, igpt}); - - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int icol=1; icol<=n_col_in; ++icol) - sfc_source({icol+col_s_in-1, igpt}) = sources_subset_in.get_sfc_source()({icol, igpt}); - } - if (!switch_fluxes) return; Array gpt_flux_up; Array gpt_flux_dn; - // Save the output per gpt if postprocessing is desired. - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - } + gpt_flux_up.set_dims({n_col_in, n_lev, 1}); + gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); constexpr int n_ang = 1; @@ -537,42 +501,14 @@ void Radiation_solver_longwave::solve( gpt_flux_up, gpt_flux_dn, n_ang); - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_flux_up ({icol+col_s_in-1, ilev}) = fluxes.get_flux_up ()({icol, ilev}); - lw_flux_dn ({icol+col_s_in-1, ilev}) = fluxes.get_flux_dn ()({icol, ilev}); - lw_flux_net({icol+col_s_in-1, ilev}) = fluxes.get_flux_net()({icol, ilev}); - } - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - for (int ibnd=1; ibnd<=n_bnd; ++ibnd) - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_bnd_flux_up ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_up ()({icol, ilev, ibnd}); - lw_bnd_flux_dn ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_dn ()({icol, ilev, ibnd}); - lw_bnd_flux_net({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_net()({icol, ilev, ibnd}); - } - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up({icol, ilev, 1}); - lw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}); - lw_flux_net({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } + // Copy the data to the output. + for (int ilev=1; ilev<=n_lev; ++ilev) + for (int icol=1; icol<=n_col_in; ++icol) + { + lw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up({icol, ilev, 1}); + lw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}); + lw_flux_net({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); + } }; for (int b=1; b<=n_blocks; ++b) @@ -584,8 +520,6 @@ void Radiation_solver_longwave::solve( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -593,8 +527,7 @@ void Radiation_solver_longwave::solve( cloud_optical_props_subset, *sources_subset, emis_sfc_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -605,8 +538,6 @@ void Radiation_solver_longwave::solve( Array emis_sfc_residual = emis_sfc.subset({{ {1, n_bnd}, {col_s, col_e} }}); std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -614,8 +545,7 @@ void Radiation_solver_longwave::solve( cloud_optical_props_residual, *sources_residual, emis_sfc_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } @@ -646,8 +576,6 @@ void Radiation_solver_shortwave::solve( const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs& gas_concs, @@ -660,12 +588,8 @@ void Radiation_solver_shortwave::solve( const Array& rel, const Array& dei, const Array& rh, const Aerosol_concs& aerosol_concs, - Array& tau, Array& ssa, Array& g, - Array& toa_src, Array& sw_flux_up, Array& sw_flux_dn, - Array& sw_flux_dn_dir, Array& sw_flux_net, - Array& sw_bnd_flux_up, Array& sw_bnd_flux_dn, - Array& sw_bnd_flux_dn_dir, Array& sw_bnd_flux_net) const + Array& sw_flux_dn_dir, Array& sw_flux_net) const { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -714,8 +638,7 @@ void Radiation_solver_shortwave::solve( std::unique_ptr& optical_props_subset_in, std::unique_ptr& cloud_optical_props_subset_in, std::unique_ptr& aerosol_optical_props_subset_in, - Fluxes_broadband& fluxes, - Fluxes_broadband& bnd_fluxes) + Fluxes_broadband& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -786,23 +709,6 @@ void Radiation_solver_shortwave::solve( } - // Store the optical properties, if desired. - if (switch_output_optical) - { - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int ilay=1; ilay<=n_lay; ++ilay) - for (int icol=1; icol<=n_col_in; ++icol) - { - tau({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_tau()({icol, ilay, igpt}); - ssa({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_ssa()({icol, ilay, igpt}); - g ({icol+col_s_in-1, ilay, igpt}) = optical_props_subset_in->get_g ()({icol, ilay, igpt}); - } - - for (int igpt=1; igpt<=n_gpt; ++igpt) - for (int icol=1; icol<=n_col_in; ++icol) - toa_src({icol+col_s_in-1, igpt}) = toa_src_subset({icol, igpt}); - } - if (!switch_fluxes) return; @@ -810,19 +716,9 @@ void Radiation_solver_shortwave::solve( Array gpt_flux_dn; Array gpt_flux_dn_dir; - // Save the output per gpt if postprocessing is desired. - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, 1}); - } + gpt_flux_up.set_dims({n_col_in, n_lev, 1}); + gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); + gpt_flux_dn_dir.set_dims({n_col_in, n_lev, 1}); Rte_sw::rte_sw( optical_props_subset_in, @@ -836,44 +732,16 @@ void Radiation_solver_shortwave::solve( gpt_flux_dn, gpt_flux_dn_dir); - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, gpt_flux_dn_dir, optical_props_subset_in, top_at_1); - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_flux_up ({icol+col_s_in-1, ilev}) = fluxes.get_flux_up ()({icol, ilev}); - sw_flux_dn ({icol+col_s_in-1, ilev}) = fluxes.get_flux_dn ()({icol, ilev}); - sw_flux_dn_dir({icol+col_s_in-1, ilev}) = fluxes.get_flux_dn_dir()({icol, ilev}); - sw_flux_net ({icol+col_s_in-1, ilev}) = fluxes.get_flux_net ()({icol, ilev}); - } - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, gpt_flux_dn_dir, optical_props_subset_in, top_at_1); - - for (int ibnd=1; ibnd<=n_bnd; ++ibnd) - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_bnd_flux_up ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_up ()({icol, ilev, ibnd}); - sw_bnd_flux_dn ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_dn ()({icol, ilev, ibnd}); - sw_bnd_flux_dn_dir({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_dn_dir()({icol, ilev, ibnd}); - sw_bnd_flux_net ({icol+col_s_in-1, ilev, ibnd}) = bnd_fluxes.get_bnd_flux_net ()({icol, ilev, ibnd}); - } - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up ({icol, ilev, 1}); - sw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}); - sw_flux_dn_dir({icol+col_s_in-1, ilev}) = gpt_flux_dn_dir({icol, ilev, 1}); - sw_flux_net ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } + // Copy the data to the output. + for (int ilev=1; ilev<=n_lev; ++ilev) + for (int icol=1; icol<=n_col_in; ++icol) + { + sw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up ({icol, ilev, 1}); + sw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}); + sw_flux_dn_dir({icol+col_s_in-1, ilev}) = gpt_flux_dn_dir({icol, ilev, 1}); + sw_flux_net ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); + } + }; for (int b=1; b<=n_blocks; ++b) @@ -883,16 +751,13 @@ void Radiation_solver_shortwave::solve( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_subset, cloud_optical_props_subset, aerosol_optical_props_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -902,15 +767,12 @@ void Radiation_solver_shortwave::solve( std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_residual, cloud_optical_props_residual, aerosol_optical_props_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } diff --git a/src_test/radiation_solver.cu b/src_test/radiation_solver.cu index 3e240be8..a0e6f475 100644 --- a/src_test/radiation_solver.cu +++ b/src_test/radiation_solver.cu @@ -419,8 +419,6 @@ Radiation_solver_longwave::Radiation_solver_longwave( void Radiation_solver_longwave::solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const Gas_concs_gpu& gas_concs, const Array_gpu& p_lay, const Array_gpu& p_lev, const Array_gpu& t_lay, const Array_gpu& t_lev, @@ -428,10 +426,7 @@ void Radiation_solver_longwave::solve_gpu( const Array_gpu& t_sfc, const Array_gpu& emis_sfc, const Array_gpu& lwp, const Array_gpu& iwp, const Array_gpu& rel, const Array_gpu& dei, - Array_gpu& tau, Array_gpu& lay_source, - Array_gpu& lev_source, Array_gpu& sfc_source, - Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net, - Array_gpu& lw_bnd_flux_up, Array_gpu& lw_bnd_flux_dn, Array_gpu& lw_bnd_flux_net) + Array_gpu& lw_flux_up, Array_gpu& lw_flux_dn, Array_gpu& lw_flux_net) { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -472,8 +467,7 @@ void Radiation_solver_longwave::solve_gpu( std::unique_ptr& cloud_optical_props_subset_in, Source_func_lw_gpu& sources_subset_in, const Array_gpu& emis_sfc_subset_in, - Fluxes_broadband_gpu& fluxes, - Fluxes_broadband_gpu& bnd_fluxes) + Fluxes_broadband_gpu& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs_gpu gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -514,43 +508,12 @@ void Radiation_solver_longwave::solve_gpu( dynamic_cast(*cloud_optical_props_subset_in)); } - // Store the optical properties, if desired. - if (switch_output_optical) - { - Subset_kernels_cuda::get_from_subset( - n_col, n_lay, n_gpt, n_col_in, col_s_in, tau.ptr(), lay_source.ptr(), lev_source.ptr(), - optical_props_subset_in->get_tau().ptr(), sources_subset_in.get_lay_source().ptr(), - sources_subset_in.get_lev_source().ptr()); - - Subset_kernels_cuda::get_from_subset( - n_col, n_gpt, n_col_in, col_s_in, sfc_source.ptr(), sources_subset_in.get_sfc_source().ptr()); - } - if (!switch_fluxes) return; Array_gpu gpt_flux_up({n_col_in, n_lev, n_gpt}); Array_gpu gpt_flux_dn({n_col_in, n_lev, n_gpt}); - - // CvH The structure below is valid if broadband flux solvers are implemented - /* - Array_gpu gpt_flux_up; - Array_gpu gpt_flux_dn; - - // Save the output per gpt if postprocessing is desired. - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - } - */ - constexpr int n_ang = 1; rte_lw.rte_lw( @@ -570,47 +533,6 @@ void Radiation_solver_longwave::solve_gpu( fluxes.get_flux_up().ptr(), fluxes.get_flux_dn().ptr(), fluxes.get_flux_net().ptr()); - if (switch_output_bnd_fluxes) - { - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, lw_bnd_flux_up.ptr(), lw_bnd_flux_dn.ptr(), lw_bnd_flux_net.ptr(), - bnd_fluxes.get_bnd_flux_up().ptr(), bnd_fluxes.get_bnd_flux_dn().ptr(), bnd_fluxes.get_bnd_flux_net().ptr()); - - } - - // CvH The structure below is valid if broadband flux solvers are implemented - /* - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - // Copy the data to the output. - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_col_in, col_s_in, lw_flux_up.ptr(), lw_flux_dn.ptr(), lw_flux_net.ptr(), - fluxes.get_flux_up().ptr(), fluxes.get_flux_dn().ptr(), fluxes.get_flux_net().ptr()); - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, lw_bnd_flux_up.ptr(), lw_bnd_flux_dn.ptr(), lw_bnd_flux_net.ptr(), - bnd_fluxes.get_bnd_flux_up().ptr(), bnd_fluxes.get_bnd_flux_dn().ptr(), bnd_fluxes.get_bnd_flux_net().ptr()); - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - lw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up({icol, ilev, 1}); - lw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}); - lw_flux_net({icol+col_s_in-1, ilev}) = gpt_flux_dn({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } - */ }; for (int b=1; b<=n_blocks; ++b) @@ -622,8 +544,6 @@ void Radiation_solver_longwave::solve_gpu( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -631,8 +551,7 @@ void Radiation_solver_longwave::solve_gpu( cloud_optical_props_subset, *sources_subset, emis_sfc_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -643,8 +562,6 @@ void Radiation_solver_longwave::solve_gpu( Array_gpu emis_sfc_residual = emis_sfc.subset({{ {1, n_bnd}, {col_s, col_e} }}); std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, @@ -652,8 +569,7 @@ void Radiation_solver_longwave::solve_gpu( cloud_optical_props_residual, *sources_residual, emis_sfc_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } @@ -684,8 +600,6 @@ void Radiation_solver_shortwave::solve_gpu( const bool switch_fluxes, const bool switch_cloud_optics, const bool switch_aerosol_optics, - const bool switch_output_optical, - const bool switch_output_bnd_fluxes, const bool switch_delta_cloud, const bool switch_delta_aerosol, const Gas_concs_gpu& gas_concs, @@ -698,12 +612,8 @@ void Radiation_solver_shortwave::solve_gpu( const Array_gpu& rel, const Array_gpu& dei, const Array_gpu& rh, const Aerosol_concs_gpu& aerosol_concs, - Array_gpu& tau, Array_gpu& ssa, Array_gpu& g, - Array_gpu& toa_src, Array_gpu& sw_flux_up, Array_gpu& sw_flux_dn, - Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net, - Array_gpu& sw_bnd_flux_up, Array_gpu& sw_bnd_flux_dn, - Array_gpu& sw_bnd_flux_dn_dir, Array_gpu& sw_bnd_flux_net) + Array_gpu& sw_flux_dn_dir, Array_gpu& sw_flux_net) { const int n_col = p_lay.dim(1); const int n_lay = p_lay.dim(2); @@ -745,8 +655,7 @@ void Radiation_solver_shortwave::solve_gpu( std::unique_ptr& optical_props_subset_in, std::unique_ptr& cloud_optical_props_subset_in, std::unique_ptr& aerosol_optical_props_subset_in, - Fluxes_broadband_gpu& fluxes, - Fluxes_broadband_gpu& bnd_fluxes) + Fluxes_broadband_gpu& fluxes) { const int n_col_in = col_e_in - col_s_in + 1; Gas_concs_gpu gas_concs_subset(gas_concs, col_s_in, n_col_in); @@ -768,8 +677,10 @@ void Radiation_solver_shortwave::solve_gpu( optical_props_subset_in, toa_src_subset, col_dry_subset); + auto tsi_scaling_subset = tsi_scaling.subset({{ {col_s_in, col_e_in} }}); scaling_to_subset(n_col_in, n_gpt, toa_src_subset, tsi_scaling_subset); + if (switch_cloud_optics) { Array cld_mask_liq({n_col_in, n_lay}); @@ -809,16 +720,6 @@ void Radiation_solver_shortwave::solve_gpu( dynamic_cast(*aerosol_optical_props_subset_in)); } - // Store the optical properties, if desired. - if (switch_output_optical) - { - Subset_kernels_cuda::get_from_subset( - n_col, n_lay, n_gpt, n_col_in, col_s_in, tau.ptr(), ssa.ptr(), g.ptr(), optical_props_subset_in->get_tau().ptr(), - optical_props_subset_in->get_ssa().ptr(), optical_props_subset_in->get_g().ptr()); - - Subset_kernels_cuda::get_from_subset( - n_col, n_gpt, n_col_in, col_s_in, toa_src.ptr(), toa_src_subset.ptr()); - } if (!switch_fluxes) return; @@ -827,26 +728,6 @@ void Radiation_solver_shortwave::solve_gpu( Array_gpu gpt_flux_dn({n_col_in, n_lev, n_gpt}); Array_gpu gpt_flux_dn_dir({n_col_in, n_lev, n_gpt}); - // CvH The structure below is valid if broadband flux solvers are implemented - /* - Array_gpu gpt_flux_up; - Array_gpu gpt_flux_dn; - Array_gpu gpt_flux_dn_dir; - - if (switch_output_bnd_fluxes) - { - gpt_flux_up.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn.set_dims({n_col_in, n_lev, n_gpt}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, n_gpt}); - } - else - { - gpt_flux_up.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn.set_dims({n_col_in, n_lev, 1}); - gpt_flux_dn_dir.set_dims({n_col_in, n_lev, 1}); - } - */ - rte_sw.rte_sw( optical_props_subset_in, top_at_1, @@ -866,50 +747,8 @@ void Radiation_solver_shortwave::solve_gpu( n_col, n_lev, n_col_in, col_s_in, sw_flux_up.ptr(), sw_flux_dn.ptr(), sw_flux_dn_dir.ptr(), sw_flux_net.ptr(), fluxes.get_flux_up().ptr(), fluxes.get_flux_dn().ptr(), fluxes.get_flux_dn_dir().ptr(), fluxes.get_flux_net().ptr()); - if (switch_output_bnd_fluxes) - { - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, sw_bnd_flux_up.ptr(), sw_bnd_flux_dn.ptr(), sw_bnd_flux_dn_dir.ptr(), sw_bnd_flux_net.ptr(), - bnd_fluxes.get_bnd_flux_up().ptr(), bnd_fluxes.get_bnd_flux_dn().ptr(), bnd_fluxes.get_bnd_flux_dn_dir().ptr(), bnd_fluxes.get_bnd_flux_net().ptr()); - } - - // CvH The structure below is valid if broadband flux solvers are implemented - /* - if (switch_output_bnd_fluxes) - { - // Aggegated fluxes. - fluxes.reduce(gpt_flux_up, gpt_flux_dn, gpt_flux_dn_dir, optical_props_subset_in, top_at_1); - - // Copy the data to the output. - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_col_in, col_s_in, sw_flux_up, sw_flux_dn, sw_flux_dn_dir, sw_flux_net, - fluxes.get_flux_up(), fluxes.get_flux_dn(), fluxes.get_flux_dn_dir(), fluxes.get_flux_net()); - - // Aggegated fluxes per band - bnd_fluxes.reduce(gpt_flux_up, gpt_flux_dn, optical_props_subset_in, top_at_1); - - Subset_kernels_cuda::get_from_subset( - n_col, n_lev, n_bnd, n_col_in, col_s_in, sw_bnd_flux_up, sw_bnd_flux_dn, sw_bnd_flux_dn_dir, sw_bnd_flux_net, - bnd_fluxes.get_bnd_flux_up(), bnd_fluxes.get_bnd_flux_dn(), bnd_fluxes.get_bnd_flux_dn_dir(), bnd_fluxes.get_bnd_flux_net()); - } - else - { - // Copy the data to the output. - for (int ilev=1; ilev<=n_lev; ++ilev) - for (int icol=1; icol<=n_col_in; ++icol) - { - sw_flux_up ({icol+col_s_in-1, ilev}) = gpt_flux_up ({icol, ilev, 1}); - sw_flux_dn ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}); - sw_flux_dn_dir({icol+col_s_in-1, ilev}) = gpt_flux_dn_dir({icol, ilev, 1}); - sw_flux_net ({icol+col_s_in-1, ilev}) = gpt_flux_dn ({icol, ilev, 1}) - gpt_flux_up({icol, ilev, 1}); - } - } - */ }; - for (int b=1; b<=n_blocks; ++b) { const int col_s = (b-1) * n_col_block + 1; @@ -917,16 +756,13 @@ void Radiation_solver_shortwave::solve_gpu( std::unique_ptr fluxes_subset = std::make_unique(n_col_block, n_lev); - std::unique_ptr bnd_fluxes_subset = - std::make_unique(n_col_block, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_subset, cloud_optical_props_subset, aerosol_optical_props_subset, - *fluxes_subset, - *bnd_fluxes_subset); + *fluxes_subset); } if (n_col_block_residual > 0) @@ -936,15 +772,12 @@ void Radiation_solver_shortwave::solve_gpu( std::unique_ptr fluxes_residual = std::make_unique(n_col_block_residual, n_lev); - std::unique_ptr bnd_fluxes_residual = - std::make_unique(n_col_block_residual, n_lev, n_bnd); call_kernels( col_s, col_e, optical_props_residual, cloud_optical_props_residual, aerosol_optical_props_residual, - *fluxes_residual, - *bnd_fluxes_residual); + *fluxes_residual); } } diff --git a/src_test/test_rte_rrtmgp.cpp b/src_test/test_rte_rrtmgp.cpp index 734a693d..1ee29fd8 100644 --- a/src_test/test_rte_rrtmgp.cpp +++ b/src_test/test_rte_rrtmgp.cpp @@ -142,8 +142,6 @@ void solve_radiation(int argc, char** argv) const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); const bool switch_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_output_optical = get_ini_value(settings, "switches", "output-optical", false); - const bool switch_output_bnd_fluxes = get_ini_value(settings, "switches", "output-bnd-fluxes", false); const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta-cloud", true); const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); @@ -270,20 +268,6 @@ void solve_radiation(int argc, char** argv) Array emis_sfc(input_nc.get_variable("emis_sfc", {n_col_y, n_col_x, n_bnd_lw}), {n_bnd_lw, n_col}); Array t_sfc(input_nc.get_variable("t_sfc", {n_col_y, n_col_x}), {n_col}); - // Create output arrays. - Array lw_tau; - Array lay_source; - Array lev_source; - Array sfc_source; - - if (switch_output_optical) - { - lw_tau .set_dims({n_col, n_lay, n_gpt_lw}); - lay_source.set_dims({n_col, n_lay, n_gpt_lw}); - lev_source.set_dims({n_col, n_lev, n_gpt_lw}); - sfc_source.set_dims({n_col, n_gpt_lw}); - } - Array lw_flux_up; Array lw_flux_dn; Array lw_flux_net; @@ -295,18 +279,6 @@ void solve_radiation(int argc, char** argv) lw_flux_net.set_dims({n_col, n_lev}); } - Array lw_bnd_flux_up; - Array lw_bnd_flux_dn; - Array lw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - lw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_net.set_dims({n_col, n_lev, n_bnd_lw}); - } - - // Solve the radiation. Status::print_message("Solving the longwave radiation."); @@ -315,8 +287,6 @@ void solve_radiation(int argc, char** argv) rad_lw.solve( switch_fluxes, switch_cloud_optics, - switch_output_optical, - switch_output_bnd_fluxes, gas_concs, p_lay, p_lev, t_lay, t_lev, @@ -324,9 +294,7 @@ void solve_radiation(int argc, char** argv) t_sfc, emis_sfc, lwp, iwp, rel, dei, - lw_tau, lay_source, lev_source, sfc_source, - lw_flux_up, lw_flux_dn, lw_flux_net, - lw_bnd_flux_up, lw_bnd_flux_dn, lw_bnd_flux_net); + lw_flux_up, lw_flux_dn, lw_flux_net); auto time_end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration(time_end-time_start).count(); @@ -343,25 +311,6 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber().v(), {0, 0}); - if (switch_output_optical) - { - auto nc_lw_band_lims_gpt = output_nc.add_variable("lw_band_lims_gpt", {"band_lw", "pair"}); - nc_lw_band_lims_gpt.insert(rad_lw.get_band_lims_gpoint().v(), {0, 0}); - - auto nc_lw_tau = output_nc.add_variable("lw_tau", {"gpt_lw", "lay", "y", "x"}); - nc_lw_tau.insert(lw_tau.v(), {0, 0, 0, 0}); - - auto nc_lay_source = output_nc.add_variable("lay_source", {"gpt_lw", "lay", "y", "x"}); - auto nc_lev_source = output_nc.add_variable("lev_source", {"gpt_lw", "lev", "y", "x"}); - - auto nc_sfc_source = output_nc.add_variable("sfc_source", {"gpt_lw", "y", "x"}); - - nc_lay_source.insert(lay_source.v(), {0, 0, 0, 0}); - nc_lev_source.insert(lev_source.v(), {0, 0, 0, 0}); - - nc_sfc_source.insert(sfc_source.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); @@ -372,16 +321,6 @@ void solve_radiation(int argc, char** argv) nc_lw_flux_dn .insert(lw_flux_dn .v(), {0, 0, 0}); nc_lw_flux_net.insert(lw_flux_net.v(), {0, 0, 0}); - if (switch_output_bnd_fluxes) - { - auto nc_lw_bnd_flux_up = output_nc.add_variable("lw_bnd_flux_up" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_dn = output_nc.add_variable("lw_bnd_flux_dn" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_net = output_nc.add_variable("lw_bnd_flux_net", {"band_lw", "lev", "y", "x"}); - - nc_lw_bnd_flux_up .insert(lw_bnd_flux_up .v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_dn .insert(lw_bnd_flux_dn .v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_net.insert(lw_bnd_flux_net.v(), {0, 0, 0, 0}); - } } } @@ -422,20 +361,6 @@ void solve_radiation(int argc, char** argv) tsi_scaling({icol}) = Float(1.); } - // Create output arrays. - Array sw_tau; - Array ssa; - Array g; - Array toa_source; - - if (switch_output_optical) - { - sw_tau .set_dims({n_col, n_lay, n_gpt_sw}); - ssa .set_dims({n_col, n_lay, n_gpt_sw}); - g .set_dims({n_col, n_lay, n_gpt_sw}); - toa_source.set_dims({n_col, n_gpt_sw}); - } - Array sw_flux_up; Array sw_flux_dn; Array sw_flux_dn_dir; @@ -449,20 +374,6 @@ void solve_radiation(int argc, char** argv) sw_flux_net .set_dims({n_col, n_lev}); } - Array sw_bnd_flux_up; - Array sw_bnd_flux_dn; - Array sw_bnd_flux_dn_dir; - Array sw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - sw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn_dir.set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_net .set_dims({n_col, n_lev, n_bnd_sw}); - } - - // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -486,12 +397,8 @@ void solve_radiation(int argc, char** argv) rel, dei, rh, aerosol_concs, - sw_tau, ssa, g, - toa_source, sw_flux_up, sw_flux_dn, - sw_flux_dn_dir, sw_flux_net, - sw_bnd_flux_up, sw_bnd_flux_dn, - sw_bnd_flux_dn_dir, sw_bnd_flux_net); + sw_flux_dn_dir, sw_flux_net); auto time_end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration(time_end-time_start).count(); @@ -508,23 +415,6 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber().v(), {0, 0}); - if (switch_output_optical) - { - auto nc_sw_band_lims_gpt = output_nc.add_variable("sw_band_lims_gpt", {"band_sw", "pair"}); - nc_sw_band_lims_gpt.insert(rad_sw.get_band_lims_gpoint().v(), {0, 0}); - - auto nc_sw_tau = output_nc.add_variable("sw_tau", {"gpt_sw", "lay", "y", "x"}); - auto nc_ssa = output_nc.add_variable("ssa" , {"gpt_sw", "lay", "y", "x"}); - auto nc_g = output_nc.add_variable("g" , {"gpt_sw", "lay", "y", "x"}); - - nc_sw_tau.insert(sw_tau.v(), {0, 0, 0, 0}); - nc_ssa .insert(ssa .v(), {0, 0, 0, 0}); - nc_g .insert(g .v(), {0, 0, 0, 0}); - - auto nc_toa_source = output_nc.add_variable("toa_source", {"gpt_sw", "y", "x"}); - nc_toa_source.insert(toa_source.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); @@ -548,31 +438,6 @@ void solve_radiation(int argc, char** argv) nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); nc_sw_flux_net.add_attribute("units","W m-2"); - - if (switch_output_bnd_fluxes) - { - auto nc_sw_bnd_flux_up = output_nc.add_variable("sw_bnd_flux_up" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn = output_nc.add_variable("sw_bnd_flux_dn" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn_dir = output_nc.add_variable("sw_bnd_flux_dn_dir", {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_net = output_nc.add_variable("sw_bnd_flux_net" , {"band_sw", "lev", "y", "x"}); - - nc_sw_bnd_flux_up .insert(sw_bnd_flux_up .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn .insert(sw_bnd_flux_dn .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn_dir.insert(sw_bnd_flux_dn_dir.v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_net .insert(sw_bnd_flux_net .v(), {0, 0, 0, 0}); - - nc_sw_bnd_flux_up.add_attribute("long_name","Upwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_up.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_net.add_attribute("long_name","Net shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_net.add_attribute("units","W m-2"); - } } } diff --git a/src_test/test_rte_rrtmgp.cu b/src_test/test_rte_rrtmgp.cu index b385bf65..a63ffdd7 100644 --- a/src_test/test_rte_rrtmgp.cu +++ b/src_test/test_rte_rrtmgp.cu @@ -168,9 +168,6 @@ void solve_radiation(int argc, char** argv) const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); const bool switch_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_output_optical = get_ini_value(settings, "switches", "output-optical", false); - const bool switch_output_bnd_fluxes = get_ini_value(settings, "switches", "output-bnd-fluxes", false); - const bool switch_timings = get_ini_value(settings, "switches", "timings", false); const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta-cloud", true); const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); @@ -185,6 +182,7 @@ void solve_radiation(int argc, char** argv) const int n_lay = input_nc.get_dimension_size("lay"); const int n_lev = input_nc.get_dimension_size("lev"); const int n_col = n_col_x * n_col_y; + // Read the atmospheric fields. Array p_lay(input_nc.get_variable("p_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); Array t_lay(input_nc.get_variable("t_lay", {n_lay, n_col_y, n_col_x}), {n_col, n_lay}); @@ -317,20 +315,6 @@ void solve_radiation(int argc, char** argv) Array emis_sfc(input_nc.get_variable("emis_sfc", {n_col_y, n_col_x, n_bnd_lw}), {n_bnd_lw, n_col}); Array t_sfc(input_nc.get_variable("t_sfc", {n_col_y, n_col_x}), {n_col}); - // Create output arrays. - Array_gpu lw_tau; - Array_gpu lay_source; - Array_gpu lev_source; - Array_gpu sfc_source; - - if (switch_output_optical) - { - lw_tau .set_dims({n_col, n_lay, n_gpt_lw}); - lay_source.set_dims({n_col, n_lay, n_gpt_lw}); - lev_source.set_dims({n_col, n_lay, n_gpt_lw}); - sfc_source.set_dims({n_col, n_gpt_lw}); - } - Array_gpu lw_flux_up; Array_gpu lw_flux_dn; Array_gpu lw_flux_net; @@ -342,18 +326,6 @@ void solve_radiation(int argc, char** argv) lw_flux_net.set_dims({n_col, n_lev}); } - Array_gpu lw_bnd_flux_up; - Array_gpu lw_bnd_flux_dn; - Array_gpu lw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - lw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_lw}); - lw_bnd_flux_net.set_dims({n_col, n_lev, n_bnd_lw}); - } - - // Solve the radiation. Status::print_message("Solving the longwave radiation."); @@ -383,8 +355,6 @@ void solve_radiation(int argc, char** argv) rad_lw.solve_gpu( switch_fluxes, switch_cloud_optics, - switch_output_optical, - switch_output_bnd_fluxes, gas_concs_gpu, p_lay_gpu, p_lev_gpu, t_lay_gpu, t_lev_gpu, @@ -392,9 +362,7 @@ void solve_radiation(int argc, char** argv) t_sfc_gpu, emis_sfc_gpu, lwp_gpu, iwp_gpu, rel_gpu, dei_gpu, - lw_tau, lay_source, lev_source, sfc_source, - lw_flux_up, lw_flux_dn, lw_flux_net, - lw_bnd_flux_up, lw_bnd_flux_dn, lw_bnd_flux_net); + lw_flux_up, lw_flux_dn, lw_flux_net); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); @@ -410,30 +378,11 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - - if (switch_timings) - { - constexpr int n_measures=10; - for (int n=0; n lw_tau_cpu(lw_tau); - Array lay_source_cpu(lay_source); - Array sfc_source_cpu(sfc_source); - Array lev_source_cpu(lev_source); Array lw_flux_up_cpu(lw_flux_up); Array lw_flux_dn_cpu(lw_flux_dn); Array lw_flux_net_cpu(lw_flux_net); - Array lw_bnd_flux_up_cpu(lw_bnd_flux_up); - Array lw_bnd_flux_dn_cpu(lw_bnd_flux_dn); - Array lw_bnd_flux_net_cpu(lw_bnd_flux_net); output_nc.add_dimension("gpt_lw", n_gpt_lw); output_nc.add_dimension("band_lw", n_bnd_lw); @@ -441,25 +390,6 @@ void solve_radiation(int argc, char** argv) auto nc_lw_band_lims_wvn = output_nc.add_variable("lw_band_lims_wvn", {"band_lw", "pair"}); nc_lw_band_lims_wvn.insert(rad_lw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_output_optical) - { - auto nc_lw_band_lims_gpt = output_nc.add_variable("lw_band_lims_gpt", {"band_lw", "pair"}); - nc_lw_band_lims_gpt.insert(rad_lw.get_band_lims_gpoint_gpu().v(), {0, 0}); - - auto nc_lw_tau = output_nc.add_variable("lw_tau", {"gpt_lw", "lay", "y", "x"}); - nc_lw_tau.insert(lw_tau_cpu.v(), {0, 0, 0, 0}); - - auto nc_lay_source = output_nc.add_variable("lay_source", {"gpt_lw", "lay", "y", "x"}); - auto nc_lev_source = output_nc.add_variable("lev_source", {"gpt_lw", "lev", "y", "x"}); - - auto nc_sfc_source = output_nc.add_variable("sfc_source", {"gpt_lw", "y", "x"}); - - nc_lay_source.insert(lay_source_cpu.v(), {0, 0, 0, 0}); - nc_lev_source.insert(lev_source_cpu.v(), {0, 0, 0, 0}); - - nc_sfc_source.insert(sfc_source_cpu.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_lw_flux_up = output_nc.add_variable("lw_flux_up" , {"lev", "y", "x"}); @@ -469,17 +399,6 @@ void solve_radiation(int argc, char** argv) nc_lw_flux_up .insert(lw_flux_up_cpu .v(), {0, 0, 0}); nc_lw_flux_dn .insert(lw_flux_dn_cpu .v(), {0, 0, 0}); nc_lw_flux_net.insert(lw_flux_net_cpu.v(), {0, 0, 0}); - - if (switch_output_bnd_fluxes) - { - auto nc_lw_bnd_flux_up = output_nc.add_variable("lw_bnd_flux_up" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_dn = output_nc.add_variable("lw_bnd_flux_dn" , {"band_lw", "lev", "y", "x"}); - auto nc_lw_bnd_flux_net = output_nc.add_variable("lw_bnd_flux_net", {"band_lw", "lev", "y", "x"}); - - nc_lw_bnd_flux_up .insert(lw_bnd_flux_up_cpu.v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_dn .insert(lw_bnd_flux_dn_cpu.v(), {0, 0, 0, 0}); - nc_lw_bnd_flux_net.insert(lw_bnd_flux_net_cpu.v(), {0, 0, 0, 0}); - } } } @@ -522,20 +441,6 @@ void solve_radiation(int argc, char** argv) tsi_scaling({icol}) = Float(1.); } - // Create output arrays. - Array_gpu sw_tau; - Array_gpu ssa; - Array_gpu g; - Array_gpu toa_source; - - if (switch_output_optical) - { - sw_tau .set_dims({n_col, n_lay, n_gpt_sw}); - ssa .set_dims({n_col, n_lay, n_gpt_sw}); - g .set_dims({n_col, n_lay, n_gpt_sw}); - toa_source.set_dims({n_col, n_gpt_sw}); - } - Array_gpu sw_flux_up; Array_gpu sw_flux_dn; Array_gpu sw_flux_dn_dir; @@ -549,20 +454,6 @@ void solve_radiation(int argc, char** argv) sw_flux_net .set_dims({n_col, n_lev}); } - Array_gpu sw_bnd_flux_up; - Array_gpu sw_bnd_flux_dn; - Array_gpu sw_bnd_flux_dn_dir; - Array_gpu sw_bnd_flux_net; - - if (switch_output_bnd_fluxes) - { - sw_bnd_flux_up .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn .set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_dn_dir.set_dims({n_col, n_lev, n_bnd_sw}); - sw_bnd_flux_net .set_dims({n_col, n_lev, n_bnd_sw}); - } - - // Solve the radiation. Status::print_message("Solving the shortwave radiation."); @@ -597,8 +488,6 @@ void solve_radiation(int argc, char** argv) switch_fluxes, switch_cloud_optics, switch_aerosol_optics, - switch_output_optical, - switch_output_bnd_fluxes, switch_delta_cloud, switch_delta_aerosol, gas_concs_gpu, @@ -611,12 +500,8 @@ void solve_radiation(int argc, char** argv) rel_gpu, dei_gpu, rh_gpu, aerosol_concs_gpu, - sw_tau, ssa, g, - toa_source, sw_flux_up, sw_flux_dn, - sw_flux_dn_dir, sw_flux_net, - sw_bnd_flux_up, sw_bnd_flux_dn, - sw_bnd_flux_dn_dir, sw_bnd_flux_net); + sw_flux_dn_dir, sw_flux_net); cudaEventRecord(stop, 0); cudaEventSynchronize(stop); @@ -632,32 +517,12 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - - if (switch_timings) - { - constexpr int n_measures=10; - for (int n=0; n sw_tau_cpu(sw_tau); - Array ssa_cpu(ssa); - Array g_cpu(g); - Array toa_source_cpu(toa_source); Array sw_flux_up_cpu(sw_flux_up); Array sw_flux_dn_cpu(sw_flux_dn); Array sw_flux_dn_dir_cpu(sw_flux_dn_dir); Array sw_flux_net_cpu(sw_flux_net); - Array sw_bnd_flux_up_cpu(sw_bnd_flux_up); - Array sw_bnd_flux_dn_cpu(sw_bnd_flux_dn); - Array sw_bnd_flux_dn_dir_cpu(sw_bnd_flux_dn_dir); - Array sw_bnd_flux_net_cpu(sw_bnd_flux_net); output_nc.add_dimension("gpt_sw", n_gpt_sw); output_nc.add_dimension("band_sw", n_bnd_sw); @@ -665,23 +530,6 @@ void solve_radiation(int argc, char** argv) auto nc_sw_band_lims_wvn = output_nc.add_variable("sw_band_lims_wvn", {"band_sw", "pair"}); nc_sw_band_lims_wvn.insert(rad_sw.get_band_lims_wavenumber_gpu().v(), {0, 0}); - if (switch_output_optical) - { - auto nc_sw_band_lims_gpt = output_nc.add_variable("sw_band_lims_gpt", {"band_sw", "pair"}); - nc_sw_band_lims_gpt.insert(rad_sw.get_band_lims_gpoint_gpu().v(), {0, 0}); - - auto nc_sw_tau = output_nc.add_variable("sw_tau", {"gpt_sw", "lay", "y", "x"}); - auto nc_ssa = output_nc.add_variable("ssa" , {"gpt_sw", "lay", "y", "x"}); - auto nc_g = output_nc.add_variable("g" , {"gpt_sw", "lay", "y", "x"}); - - nc_sw_tau.insert(sw_tau_cpu.v(), {0, 0, 0, 0}); - nc_ssa .insert(ssa_cpu .v(), {0, 0, 0, 0}); - nc_g .insert(g_cpu .v(), {0, 0, 0, 0}); - - auto nc_toa_source = output_nc.add_variable("toa_source", {"gpt_sw", "y", "x"}); - nc_toa_source.insert(toa_source_cpu.v(), {0, 0, 0}); - } - if (switch_fluxes) { auto nc_sw_flux_up = output_nc.add_variable("sw_flux_up" , {"lev", "y", "x"}); @@ -706,30 +554,6 @@ void solve_radiation(int argc, char** argv) nc_sw_flux_net.add_attribute("long_name","Net shortwave fluxes (TwoStream solver)"); nc_sw_flux_net.add_attribute("units","W m-2"); - if (switch_output_bnd_fluxes) - { - auto nc_sw_bnd_flux_up = output_nc.add_variable("sw_bnd_flux_up" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn = output_nc.add_variable("sw_bnd_flux_dn" , {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_dn_dir = output_nc.add_variable("sw_bnd_flux_dn_dir", {"band_sw", "lev", "y", "x"}); - auto nc_sw_bnd_flux_net = output_nc.add_variable("sw_bnd_flux_net" , {"band_sw", "lev", "y", "x"}); - - nc_sw_bnd_flux_up .insert(sw_bnd_flux_up_cpu .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn .insert(sw_bnd_flux_dn_cpu .v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_dn_dir.insert(sw_bnd_flux_dn_dir_cpu.v(), {0, 0, 0, 0}); - nc_sw_bnd_flux_net .insert(sw_bnd_flux_net_cpu .v(), {0, 0, 0, 0}); - - nc_sw_bnd_flux_up.add_attribute("long_name","Upwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_up.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn.add_attribute("long_name","Downwelling shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_dn_dir.add_attribute("long_name","Downwelling direct shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_dn_dir.add_attribute("units","W m-2"); - - nc_sw_bnd_flux_net.add_attribute("long_name","Net shortwave fluxes per spectral band (TwoStream solver)"); - nc_sw_bnd_flux_net.add_attribute("units","W m-2"); - } } } diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 9d6a3863..4f74cd73 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -182,7 +182,6 @@ void solve_radiation(int argc, char** argv) const bool switch_cloud_mie = get_ini_value(settings, "switches", "cloud-mie", false); const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); const bool switch_single_gpt = get_ini_value(settings, "switches", "single-gpt", false); - const bool switch_profiling = get_ini_value(settings, "switches", "profiling", false); const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta-cloud", false); const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); const bool switch_tica = get_ini_value(settings, "switches", "tica", false); @@ -779,14 +778,6 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - if (switch_profiling) - { - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - } - //// Store the output. Status::print_message("Storing the longwave output."); Array lw_tau_tot_cpu(lw_tau_tot); @@ -1113,14 +1104,6 @@ void solve_radiation(int argc, char** argv) // Tuning step; run_solver(); - // Profiling step; - if (switch_profiling) - { - cudaProfilerStart(); - run_solver(); - cudaProfilerStop(); - } - // Store the output. Status::print_message("Storing the shortwave output."); Array sw_tot_tau_cpu(sw_tot_tau); From 550269959c61c6abb0bb0ce6a91c05a33b054a13 Mon Sep 17 00:00:00 2001 From: Menno Date: Tue, 3 Mar 2026 20:31:08 +0100 Subject: [PATCH 07/13] forgotting in previous commit --- src_test/test_rte_rrtmgp.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src_test/test_rte_rrtmgp.cpp b/src_test/test_rte_rrtmgp.cpp index 1ee29fd8..79c91bd7 100644 --- a/src_test/test_rte_rrtmgp.cpp +++ b/src_test/test_rte_rrtmgp.cpp @@ -383,8 +383,6 @@ void solve_radiation(int argc, char** argv) switch_fluxes, switch_cloud_optics, switch_aerosol_optics, - switch_output_optical, - switch_output_bnd_fluxes, switch_delta_cloud, switch_delta_aerosol, gas_concs, From c66f173a7ca5edb1014ed270b9cb0532c5d5ed2a Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 13:15:19 +0100 Subject: [PATCH 08/13] enable lw delta scaling in RT executable --- include_test/radiation_solver_rt.h | 2 ++ src_test/radiation_solver.cpp | 2 -- src_test/radiation_solver.cu | 2 -- src_test/radiation_solver_rt.cu | 10 ++++++++-- src_test/test_rte_rrtmgp_rt.cu | 4 ++-- 5 files changed, 12 insertions(+), 8 deletions(-) diff --git a/include_test/radiation_solver_rt.h b/include_test/radiation_solver_rt.h index 75d0601e..9df1af22 100644 --- a/include_test/radiation_solver_rt.h +++ b/include_test/radiation_solver_rt.h @@ -50,6 +50,8 @@ class Radiation_solver_longwave const bool switch_raytracing, const bool switch_cloud_optics, const bool switch_aerosol_optics, + const bool switch_delta_cloud, + const bool switch_delta_aerosol, const bool switch_single_gpt, const bool switch_lw_scattering, const bool switch_independent_column, diff --git a/src_test/radiation_solver.cpp b/src_test/radiation_solver.cpp index 2d2ebc69..1fa0bcfa 100644 --- a/src_test/radiation_solver.cpp +++ b/src_test/radiation_solver.cpp @@ -473,8 +473,6 @@ void Radiation_solver_longwave::solve( dei.subset({{ {col_s_in, col_e_in}, {1, n_lay} }}), *cloud_optical_props_subset_in); - // cloud->delta_scale(); - // Add the cloud optical props to the gas optical properties. add_to( dynamic_cast(*optical_props_subset_in), diff --git a/src_test/radiation_solver.cu b/src_test/radiation_solver.cu index a0e6f475..75c3fdd1 100644 --- a/src_test/radiation_solver.cu +++ b/src_test/radiation_solver.cu @@ -500,8 +500,6 @@ void Radiation_solver_longwave::solve_gpu( dei.subset({{ {col_s_in, col_e_in}, {1, n_lay} }}), *cloud_optical_props_subset_in); - //cloud->delta_scale(); - // Add the cloud optical props to the gas optical properties. add_to( dynamic_cast(*optical_props_subset_in), diff --git a/src_test/radiation_solver_rt.cu b/src_test/radiation_solver_rt.cu index 9f77d26e..7278181a 100644 --- a/src_test/radiation_solver_rt.cu +++ b/src_test/radiation_solver_rt.cu @@ -483,6 +483,8 @@ void Radiation_solver_longwave::solve_gpu( const bool switch_raytracing, const bool switch_cloud_optics, const bool switch_aerosol_optics, + const bool switch_delta_cloud, + const bool switch_delta_aerosol, const bool switch_single_gpt, const bool switch_lw_scattering, const bool switch_independent_column, @@ -643,7 +645,9 @@ void Radiation_solver_longwave::solve_gpu( dei, switch_lw_scattering, *cloud_optical_props); - // cloud->delta_scale(); + + if (switch_lw_scattering && switch_delta_cloud) + cloud_optical_props->delta_scale(); } @@ -670,7 +674,9 @@ void Radiation_solver_longwave::solve_gpu( rh, p_lev, switch_lw_scattering, *aerosol_optical_props); - // aerosol->delta_scale(); + + if (switch_lw_scattering && switch_delta_aerosol) + aerosol_optical_props->delta_scale(); } // Add the cloud optical props to the gas optical properties. diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index 4f74cd73..e884173c 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -738,6 +738,8 @@ void solve_radiation(int argc, char** argv) switch_lw_raytracing, switch_cloud_optics, switch_aerosol_optics, + switch_delta_cloud, + switch_delta_aerosol, switch_single_gpt, switch_lw_scattering, switch_independent_column, @@ -1497,12 +1499,10 @@ void solve_radiation(int argc, char** argv) if (switch_broadband) { - // Profiling step; run_solver_bb(); } if (switch_image) { - // actual solve run_solver(); } From 8f03fd6e567aecc354649e616c2c743a32d2806f Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 13:55:14 +0100 Subject: [PATCH 09/13] pass reference to r_eff and index ijk to mie sampling functions to avoid out of bounds read if we don't initalize rel --- include_rt/raytracer_functions.h | 12 ++++++------ src_kernels_cuda_rt/raytracer_kernels_bw.cu | 14 +++++++------- src_kernels_cuda_rt/raytracer_kernels_sw.cu | 2 +- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/include_rt/raytracer_functions.h b/include_rt/raytracer_functions.h index dc5d40e2..07eb7517 100644 --- a/include_rt/raytracer_functions.h +++ b/include_rt/raytracer_functions.h @@ -106,11 +106,11 @@ namespace Raytracer_functions } __device__ - inline Float mie_sample_angle(const Float* mie_cdf, const Float* mie_lut, const Float random_number, const Float r_eff, const int n_mie) + inline Float mie_sample_angle(const Float* mie_cdf, const Float* mie_lut, const Float random_number, const Float* r_eff, const int ijk, const int n_mie) { // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP) - const int r_idx = min(max(int(r_eff-2.5), 0), 18); - const Float r_rest = fmod(r_eff-Float(2.5),Float(1.)); + const int r_idx = min(max(int(r_eff[ijk]-2.5), 0), 18); + const Float r_rest = fmod(r_eff[ijk]-Float(2.5),Float(1.)); const int i = min(max(0, find_index(mie_cdf, n_mie, random_number)), n_mie - 2); @@ -125,11 +125,11 @@ namespace Raytracer_functions } __device__ - inline Float mie_interpolate_phase_table(const Float* mie_phase, const Float* mie_lut, const Float scat_ang, const Float r_eff, const int n_mie) + inline Float mie_interpolate_phase_table(const Float* mie_phase, const Float* mie_lut, const Float scat_ang, const Float* r_eff, const int ijk, const int n_mie) { // interpolation over effective radius. Currently, r_eff should range between 2.5 and 21.5 (similar to RRTMGP) - const int r_idx = min(max(int(r_eff-2.5), 0), 18); - const Float r_rest = fmod(r_eff-Float(2.5),Float(1.)); + const int r_idx = min(max(int(r_eff[ijk]-2.5), 0), 18); + const Float r_rest = fmod(r_eff[ijk]-Float(2.5),Float(1.)); // interpolation between 1800 equally spaced scattering angles between 0 and PI (both inclusive). constexpr Float d_pi = Float(1.74629942e-03); diff --git a/src_kernels_cuda_rt/raytracer_kernels_bw.cu b/src_kernels_cuda_rt/raytracer_kernels_bw.cu index 053739d6..dc60442a 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_bw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_bw.cu @@ -275,7 +275,8 @@ namespace const Float sun_solid_angle, const Float g, const Float* __restrict__ mie_phase_ang, const Float* __restrict__ mie_phase, - const Float r_eff, + const Float* __restrict__ r_eff, + const int ijk, const int mie_table_size, const Vector& normal, const Phase_kind kind) @@ -289,7 +290,7 @@ namespace else if (kind == Phase_kind::Mie) { // return interpolate_mie_phase_table(mie_phase_ang, mie_phase, max(0.05, acos(cos_angle)), r_eff, mie_table_size) * sun_solid_angle; - return mie_interpolate_phase_table(mie_phase_ang, mie_phase, acos(cos_angle), r_eff, mie_table_size) * sun_solid_angle; + return mie_interpolate_phase_table(mie_phase_ang, mie_phase, acos(cos_angle), r_eff, ijk, mie_table_size) * sun_solid_angle; } else if (kind == Phase_kind::Rayleigh) { @@ -497,7 +498,7 @@ void ray_tracer_kernel_bw( // direct contribution const Phase_kind kind = (scatter_type==0) ? Phase_kind::Rayleigh : Phase_kind::HG; - const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, Float(0.), 0, surface_normal, kind); + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff, 0, 0, surface_normal, kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, bg_tau_cum, z_lev_bg, bg_idx, @@ -585,7 +586,7 @@ void ray_tracer_kernel_bw( : Phase_kind::Lambertian; // SUN SCATTERING GOES HERE - const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, Float(0.), mie_phase_ang_shared, mie_phase, Float(0.), 0, + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, Float(0.), mie_phase_ang_shared, mie_phase, r_eff, 0, 0, surface_normal, surface_kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, @@ -714,7 +715,7 @@ void ray_tracer_kernel_bw( } const Float cos_scat = scatter_type == 0 ? rayleigh(rng()) : // gases -> rayleigh, 1 ? ( (mie_cdf_table_size > 0) //clouds: Mie or HG - ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff[ijk], mie_cdf_table_size) ) + ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff, ijk, mie_cdf_table_size) ) : henyey(g, rng())) : henyey(g, rng()); //aerosols const Float sin_scat = max(Float(0.), sqrt(Float(1.) - cos_scat*cos_scat + Float_epsilon)); @@ -725,8 +726,7 @@ void ray_tracer_kernel_bw( ? Phase_kind::Mie : Phase_kind::HG : Phase_kind::HG; - const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff[ijk], mie_phase_table_size, - surface_normal, kind); + const Float p_sun = probability_from_sun(photon, sun_direction, sun_solid_angle, g, mie_phase_ang_shared, mie_phase, r_eff, ijk, mie_phase_table_size, surface_normal, kind); const Float trans_sun = transmission_direct_sun(photon,n,rng,sun_direction, k_null_grid,k_ext, bg_tau_cum, z_lev_bg, bg_idx, diff --git a/src_kernels_cuda_rt/raytracer_kernels_sw.cu b/src_kernels_cuda_rt/raytracer_kernels_sw.cu index 2853c00f..9663499a 100644 --- a/src_kernels_cuda_rt/raytracer_kernels_sw.cu +++ b/src_kernels_cuda_rt/raytracer_kernels_sw.cu @@ -410,7 +410,7 @@ void ray_tracer_kernel( // 0 (gas): rayleigh, 1 (cloud): mie if mie_table_size>0 else HG, 2 (aerosols) HG const Float cos_scat = scatter_type == 0 ? rayleigh(rng()) : // gases -> rayleigh, 1 ? ( (mie_table_size > 0) //clouds: Mie or HG - ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff[ijk], mie_table_size) ) + ? cos( mie_sample_angle(mie_cdf_shared, mie_ang, rng(), r_eff, ijk, mie_table_size) ) : henyey(g, rng())) : henyey(g, rng()); //aerosols const Float sin_scat = max(Float(0.), sqrt(Float(1.) - cos_scat*cos_scat + Float_epsilon)); From d9416e8c8eacb4e04fd25a497a8962694815077c Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 13:56:13 +0100 Subject: [PATCH 10/13] remove switch_cloud_optics and split into switch_liquid_cloud_optics and switch_ice_cloud_optics in all executables --- include_tilt/tilt_utils.h | 2 +- src_test/test_rte_rrtmgp.cpp | 45 +++++++++++------ src_test/test_rte_rrtmgp.cu | 48 +++++++++++++------ src_test/test_rte_rrtmgp_rt.cu | 88 ++++++++++++++-------------------- src_tilt/tilt_utils.cpp | 20 ++++---- 5 files changed, 111 insertions(+), 92 deletions(-) diff --git a/include_tilt/tilt_utils.h b/include_tilt/tilt_utils.h index 2d49568f..0d9a81cc 100644 --- a/include_tilt/tilt_utils.h +++ b/include_tilt/tilt_utils.h @@ -121,4 +121,4 @@ void tica_tilt(const Float sza, const Float azi, Array& lwp_out, Array& iwp_out, Array& rel_out, Array& dei_out, Array& rh_out, Gas_concs& gas_concs_out, Aerosol_concs& aerosol_concs_out, std::vector gas_names, std::vector aerosol_names, - bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics); + bool switch_liquid_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics); diff --git a/src_test/test_rte_rrtmgp.cpp b/src_test/test_rte_rrtmgp.cpp index 79c91bd7..2f8a1aa4 100644 --- a/src_test/test_rte_rrtmgp.cpp +++ b/src_test/test_rte_rrtmgp.cpp @@ -137,13 +137,14 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - const bool switch_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta-cloud", true); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); + const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); + const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); + const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); + const bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); + const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", true); + const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); ////// READ THE ATMOSPHERIC DATA ////// @@ -200,19 +201,35 @@ void solve_radiation(int argc, char** argv) Array rel; Array dei; + const bool switch_cloud_optics = switch_liquid_cloud_optics || switch_ice_cloud_optics; if (switch_cloud_optics) { lwp.set_dims({n_col, n_lay}); - lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - + rel.set_dims({n_col, n_lay}); iwp.set_dims({n_col, n_lay}); - iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei.set_dims({n_col, n_lay}); - rel.set_dims({n_col, n_lay}); - rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + if (switch_liquid_cloud_optics) + { + lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); + rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + } + else + { + lwp.fill(Float(0.)); + rel.fill(Float(0.)); + } - dei.set_dims({n_col, n_lay}); - dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + if (switch_ice_cloud_optics) + { + iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + } + else + { + iwp.fill(Float(0.)); + dei.fill(Float(0.)); + } } Array rh; diff --git a/src_test/test_rte_rrtmgp.cu b/src_test/test_rte_rrtmgp.cu index a63ffdd7..0bb86ca1 100644 --- a/src_test/test_rte_rrtmgp.cu +++ b/src_test/test_rte_rrtmgp.cu @@ -163,13 +163,14 @@ void solve_radiation(int argc, char** argv) const auto settings = toml::parse(case_name + ".ini"); ////// FLOW CONTROL SWITCHES ////// - const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); - const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); - const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - const bool switch_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta-cloud", true); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); + const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); + const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); + const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); + const bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); + const bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); + const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); + const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", true); + const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); ////// READ THE ATMOSPHERIC DATA ////// @@ -226,19 +227,35 @@ void solve_radiation(int argc, char** argv) Array rel; Array dei; + const bool switch_cloud_optics = switch_liquid_cloud_optics || switch_ice_cloud_optics; if (switch_cloud_optics) { lwp.set_dims({n_col, n_lay}); - lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - + rel.set_dims({n_col, n_lay}); iwp.set_dims({n_col, n_lay}); - iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei.set_dims({n_col, n_lay}); - rel.set_dims({n_col, n_lay}); - rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + if (switch_liquid_cloud_optics) + { + lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); + rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); + } + else + { + lwp.fill(Float(0.)); + rel.fill(Float(0.)); + } - dei.set_dims({n_col, n_lay}); - dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + if (switch_ice_cloud_optics) + { + iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); + } + else + { + iwp.fill(Float(0.)); + dei.fill(Float(0.)); + } } Array rh; @@ -411,7 +428,8 @@ void solve_radiation(int argc, char** argv) Gas_concs_gpu gas_concs_gpu(gas_concs); - Radiation_solver_shortwave rad_sw(gas_concs_gpu, switch_cloud_optics, switch_aerosol_optics, "coefficients_sw.nc", "cloud_coefficients_sw.nc", "aerosol_optics_sw.nc"); + Radiation_solver_shortwave rad_sw(gas_concs_gpu, switch_cloud_optics, switch_aerosol_optics, + "coefficients_sw.nc", "cloud_coefficients_sw.nc", "aerosol_optics_sw.nc"); // Read the boundary conditions. const int n_bnd_sw = rad_sw.get_n_bnd_gpu(); diff --git a/src_test/test_rte_rrtmgp_rt.cu b/src_test/test_rte_rrtmgp_rt.cu index e884173c..4dd5911b 100644 --- a/src_test/test_rte_rrtmgp_rt.cu +++ b/src_test/test_rte_rrtmgp_rt.cu @@ -170,29 +170,28 @@ void solve_radiation(int argc, char** argv) const bool switch_shortwave = get_ini_value(settings, "switches", "shortwave", true); const bool switch_longwave = get_ini_value(settings, "switches", "longwave", true); const bool switch_fluxes = get_ini_value(settings, "switches", "fluxes", true); - bool switch_sw_twostream = get_ini_value(settings, "switches", "sw-two-stream", false); - bool switch_sw_raytracing = get_ini_value(settings, "switches", "sw-raytracing", true); - bool switch_lw_raytracing = get_ini_value(settings, "switches", "lw-raytracing", true); - bool switch_independent_column = get_ini_value(settings, "switches", "independent-column", false); - bool switch_cloud_optics = get_ini_value(settings, "switches", "cloud-optics", false); - bool switch_liq_cloud_optics = get_ini_value(settings, "switches", "liq-cloud-optics", false); - bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice-cloud-optics", false); - const bool switch_lw_scattering = get_ini_value(settings, "switches", "lw-scattering", false); - const bool switch_min_mfp_grid_ratio= get_ini_value(settings, "switches", "min-mfp-grid-ratio", true); - const bool switch_cloud_mie = get_ini_value(settings, "switches", "cloud-mie", false); - const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol-optics", false); - const bool switch_single_gpt = get_ini_value(settings, "switches", "single-gpt", false); - const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta-cloud", false); - const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta-aerosol", false); + bool switch_sw_twostream = get_ini_value(settings, "switches", "sw_two_stream", false); + bool switch_sw_raytracing = get_ini_value(settings, "switches", "sw_raytracing", true); + bool switch_lw_raytracing = get_ini_value(settings, "switches", "lw_raytracing", true); + bool switch_independent_column = get_ini_value(settings, "switches", "independent_column", false); + bool switch_liquid_cloud_optics = get_ini_value(settings, "switches", "liquid_cloud_optics", false); + bool switch_ice_cloud_optics = get_ini_value(settings, "switches", "ice_cloud_optics", false); + const bool switch_lw_scattering = get_ini_value(settings, "switches", "lw_scattering", false); + const bool switch_min_mfp_grid_ratio= get_ini_value(settings, "switches", "min_mfp_grid_ratio", true); + const bool switch_cloud_mie = get_ini_value(settings, "switches", "cloud_mie", false); + const bool switch_aerosol_optics = get_ini_value(settings, "switches", "aerosol_optics", false); + const bool switch_single_gpt = get_ini_value(settings, "switches", "single_gpt", false); + const bool switch_delta_cloud = get_ini_value(settings, "switches", "delta_cloud", false); + const bool switch_delta_aerosol = get_ini_value(settings, "switches", "delta_aerosol", false); const bool switch_tica = get_ini_value(settings, "switches", "tica", false); - const bool switch_bw_raytracing = get_ini_value(settings, "switches", "bw-raytracing", true); - const bool switch_lu_albedo = get_ini_value(settings, "switches", "lu-albedo", false); + const bool switch_bw_raytracing = get_ini_value(settings, "switches", "bw_raytracing", true); + const bool switch_lu_albedo = get_ini_value(settings, "switches", "lu_albedo", false); const bool switch_image = get_ini_value(settings, "switches", "image", true); const bool switch_broadband = get_ini_value(settings, "switches", "broadband", false); - const bool switch_cloud_cam = get_ini_value(settings, "switches", "cloud-cam", false); + const bool switch_cloud_cam = get_ini_value(settings, "switches", "cloud_cam", false); - const int photons_per_pixel = get_ini_value(settings, "ints", "bw-raytracing", 1); + const int photons_per_pixel = get_ini_value(settings, "ints", "bw_raytracing", 1); if (!switch_shortwave) switch_sw_raytracing = false; @@ -209,7 +208,7 @@ void solve_radiation(int argc, char** argv) Int sw_photons_per_pixel; if (switch_sw_raytracing) { - sw_photons_per_pixel = get_ini_value(settings, "ints", "sw-raytracing", Int(256)); + sw_photons_per_pixel = get_ini_value(settings, "ints", "sw_raytracing", Int(256)); if (Float(int(std::log2(Float(sw_photons_per_pixel)))) != std::log2(Float(sw_photons_per_pixel))) { std::string error = "number of photons per pixel should be a power of 2 "; @@ -222,7 +221,7 @@ void solve_radiation(int argc, char** argv) Int lw_photon_count; if (switch_lw_raytracing) { - lw_photon_power = get_ini_value(settings, "ints", "lw-raytracing", Int(22)); + lw_photon_power = get_ini_value(settings, "ints", "lw_raytracing", Int(22)); lw_photon_count = 1 << lw_photon_power; } @@ -231,18 +230,10 @@ void solve_radiation(int argc, char** argv) Status::print_warning("No longwave radiation implemented in the backward ray tracer"); } - if (switch_cloud_optics) - { - switch_liq_cloud_optics = true; - switch_ice_cloud_optics = true; - } - if (switch_liq_cloud_optics || switch_ice_cloud_optics) - switch_cloud_optics = true; - if (switch_tica) switch_independent_column = true; - if (switch_cloud_cam && !(switch_cloud_optics)) + if (switch_cloud_cam && !(switch_liquid_cloud_optics || switch_ice_cloud_optics)) throw std::runtime_error("Enable cloud optics (liquid and/or ice) when cloud-cam is switch on!"); if (switch_cloud_mie && switch_ice_cloud_optics) @@ -253,7 +244,7 @@ void solve_radiation(int argc, char** argv) int single_gpt = get_ini_value(settings, "ints", "single-gpt", 1); - const Float min_mfp_grid_ratio = switch_min_mfp_grid_ratio ? get_ini_value(settings, "floats", "min-mfp-grid-ratio", Float(1.)) : Float(0.); + const Float min_mfp_grid_ratio = switch_min_mfp_grid_ratio ? get_ini_value(settings, "floats", "min_mfp_grid_ratio", Float(1.)) : Float(0.); if (switch_sw_raytracing) Status::print_message("Shortwave: using "+ std::to_string(sw_photons_per_pixel) + " rays per pixel per g-point"); @@ -387,31 +378,24 @@ void solve_radiation(int argc, char** argv) Array rel; Array dei; - if (switch_cloud_optics) - { - - if (switch_liq_cloud_optics) - { - lwp.set_dims({n_col, n_lay}); - lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - - rel.set_dims({n_col, n_lay}); - rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); - } + const bool switch_cloud_optics = switch_liquid_cloud_optics || switch_ice_cloud_optics; - if (switch_ice_cloud_optics) - { - iwp.set_dims({n_col, n_lay}); - iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + if (switch_liquid_cloud_optics) + { + lwp.set_dims({n_col, n_lay}); + lwp = std::move(input_nc.get_variable("lwp", {n_lay, n_col_y, n_col_x})); - dei.set_dims({n_col, n_lay}); - dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); - } + rel.set_dims({n_col, n_lay}); + rel = std::move(input_nc.get_variable("rel", {n_lay, n_col_y, n_col_x})); } - else + + if (switch_ice_cloud_optics) { - rel.set_dims({n_col, n_lay}); - rel.fill(Float(0.)); + iwp.set_dims({n_col, n_lay}); + iwp = std::move(input_nc.get_variable("iwp", {n_lay, n_col_y, n_col_x})); + + dei.set_dims({n_col, n_lay}); + dei = std::move(input_nc.get_variable("dei", {n_lay, n_col_y, n_col_x})); } Array rh; @@ -527,7 +511,7 @@ void solve_radiation(int argc, char** argv) lwp_out, iwp_out, rel_out, dei_out, rh_out, gas_concs_out, aerosol_concs_out, gas_names, aerosol_names, - switch_cloud_optics, switch_liq_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics + switch_liquid_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics ); lwp_out.expand_dims({n_col, n_lay}); diff --git a/src_tilt/tilt_utils.cpp b/src_tilt/tilt_utils.cpp index fa78aacb..116349e4 100644 --- a/src_tilt/tilt_utils.cpp +++ b/src_tilt/tilt_utils.cpp @@ -227,12 +227,12 @@ void restore_bkg_profile_bundle(const int n_col_x, const int n_col_y, Array* lwp, Array* iwp, Array* rel, Array* dei, Array* rh, Gas_concs& gas_concs, Aerosol_concs& aerosol_concs, std::vector gas_names, std::vector aerosol_names, - bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics + bool switch_liquid_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics ) { const int n_col = n_col_x*n_col_y; - if (switch_liq_cloud_optics) + if (switch_liquid_cloud_optics) { restore_bkg_profile(n_col_x, n_col_y, n_lay, n_z_in, bkg_start_z, lwp_copy->v(), lwp->v()); lwp_copy->expand_dims({n_col, n_lay_tot}); @@ -333,7 +333,7 @@ void post_process_output(const std::vector& col_results, Array* iwp_out, Array* rel_out, Array* dei_out, - const bool switch_liq_cloud_optics, + const bool switch_liquid_cloud_optics, const bool switch_ice_cloud_optics) { const int total_cols = n_col_x * n_col_y; @@ -343,7 +343,7 @@ void post_process_output(const std::vector& col_results, const ColumnResult& col = col_results[col_idx]; for (int j = 0; j < n_z; ++j) { int out_idx = col_idx + j * total_cols; - if (switch_liq_cloud_optics) { + if (switch_liquid_cloud_optics) { lwp_out->v()[out_idx] = col.lwp.v()[j]; rel_out->v()[out_idx] = col.rel.v()[j]; } @@ -822,7 +822,7 @@ void tica_tilt( Array& lwp_out, Array& iwp_out, Array& rel_out, Array& dei_out, Array& rh_out, Gas_concs& gas_concs_out, Aerosol_concs& aerosol_concs_out, std::vector gas_names, std::vector aerosol_names, - bool switch_cloud_optics, bool switch_liq_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics + bool switch_liquid_cloud_optics, bool switch_ice_cloud_optics, bool switch_aerosol_optics ) { // if t lev all 0, interpolate from t lay @@ -869,7 +869,7 @@ void tica_tilt( gas_concs_out, gas_names, aerosol_concs_out, aerosol_names, switch_aerosol_optics); ////// SETUP FOR RANDOM START POINT TILTING ////// - if (switch_cloud_optics) + if (switch_liquid_cloud_optics || switch_ice_cloud_optics) { std::vector> x_y_start_arr(n_col); std::vector> by_col_paths(n_col); @@ -954,7 +954,7 @@ void tica_tilt( for (int ilay = 1; ilay <= n_zh_in; ++ilay) { Float dz = zh({ilay + 1}) - zh({ilay}); - if (switch_liq_cloud_optics) + if (switch_liquid_cloud_optics) { (lwp_norm_reshaped)({ilay, icol}) = (lwp)({icol, ilay})/dz; (rel_reshaped)({ilay, icol}) = (rel)({icol, ilay}); @@ -967,7 +967,7 @@ void tica_tilt( } } - if (switch_liq_cloud_optics){ + if (switch_liquid_cloud_optics){ Array lwp_compress; lwp_compress.set_dims({n_z_in}); Array rel_compress; @@ -1151,7 +1151,7 @@ void tica_tilt( post_process_output(col_results, n_col_x, n_col_y, n_z_in, n_zh_in, &lwp_out, &iwp_out, &rel_out, &dei_out, - switch_liq_cloud_optics, switch_ice_cloud_optics); + switch_liquid_cloud_optics, switch_ice_cloud_optics); } @@ -1165,7 +1165,7 @@ void tica_tilt( &lwp, &iwp, &rel, &dei, &rh, gas_concs, aerosol_concs, gas_names, aerosol_names, - switch_liq_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics + switch_liquid_cloud_optics, switch_ice_cloud_optics, switch_aerosol_optics ); } From c82a4af961321bffbb22bbc8d646797e26578e83 Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 14:09:20 +0100 Subject: [PATCH 11/13] update ini files of test cases --- allsky/allsky.ini | 41 ++++++++++++++++++++--------------------- rcemip/rcemip.ini | 41 ++++++++++++++++++++--------------------- rfmip/rfmip.ini | 41 ++++++++++++++++++++--------------------- 3 files changed, 60 insertions(+), 63 deletions(-) diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 85de65b4..68b4e10a 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -2,34 +2,33 @@ shortwave = true longwave = true fluxes = true -cloud-optics = true -aerosol-optics = false -output-optical = false -output-bnd-fluxes = false -delta-cloud = true -delta-aerosol = false +liquid_cloud-optics = true +ice_cloud-optics = true +aerosol_optics = false +delta_cloud = true +delta_aerosol = false -# gpu-only (test_rte_rrtmgp_gpu) +# gpu_only (test_rte_rrtmgp_gpu) timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw-two-stream = false -sw-raytracing = true -lw-raytracing = true -independent-column = false -liq-cloud-optics = false -ice-cloud-optics = false -lw-scattering = false -cloud-mie = false -single-gpt = false +sw_two-stream = false +sw_raytracing = true +lw_raytracing = true +independent_column = false +liq_cloud-optics = false +ice_cloud-optics = false +lw_scattering = false +cloud_mie = false +single_gpt = false profiling = false -min-mfp-grid-ratio = true +min_mfp-grid-ratio = true tica = false [ints] -sw-raytracing = 256 -lw-raytracing = 22 -single-gpt = 1 +sw_raytracing = 256 +lw_raytracing = 22 +single_gpt = 1 [floats] -min-mfp-grid-ratio = 1.0 +min_mfp-grid-ratio = 1.0 diff --git a/rcemip/rcemip.ini b/rcemip/rcemip.ini index d1957694..00363448 100644 --- a/rcemip/rcemip.ini +++ b/rcemip/rcemip.ini @@ -2,34 +2,33 @@ shortwave = true longwave = true fluxes = true -cloud-optics = false -aerosol-optics = false -output-optical = false -output-bnd-fluxes = false -delta-cloud = true -delta-aerosol = false +liquid_cloud_optics = false +ice_cloud_optics = false +aerosol_optics = false +delta_cloud = true +delta_aerosol = false -# gpu-only (test_rte_rrtmgp_gpu) +# gpu_only (test_rte_rrtmgp_gpu) timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw-two-stream = false -sw-raytracing = true -lw-raytracing = true -independent-column = false -liq-cloud-optics = false -ice-cloud-optics = false -lw-scattering = false -cloud-mie = false -single-gpt = false +sw_two-stream = false +sw_raytracing = true +lw_raytracing = true +independent_column = false +liq_cloud-optics = false +ice_cloud-optics = false +lw_scattering = false +cloud_mie = false +single_gpt = false profiling = false -min-mfp-grid-ratio = true +min_mfp-grid-ratio = true tica = false [ints] -sw-raytracing = 256 -lw-raytracing = 22 -single-gpt = 1 +sw_raytracing = 256 +lw_raytracing = 22 +single_gpt = 1 [floats] -min-mfp-grid-ratio = 1.0 +min_mfp-grid-ratio = 1.0 diff --git a/rfmip/rfmip.ini b/rfmip/rfmip.ini index d1957694..00363448 100644 --- a/rfmip/rfmip.ini +++ b/rfmip/rfmip.ini @@ -2,34 +2,33 @@ shortwave = true longwave = true fluxes = true -cloud-optics = false -aerosol-optics = false -output-optical = false -output-bnd-fluxes = false -delta-cloud = true -delta-aerosol = false +liquid_cloud_optics = false +ice_cloud_optics = false +aerosol_optics = false +delta_cloud = true +delta_aerosol = false -# gpu-only (test_rte_rrtmgp_gpu) +# gpu_only (test_rte_rrtmgp_gpu) timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw-two-stream = false -sw-raytracing = true -lw-raytracing = true -independent-column = false -liq-cloud-optics = false -ice-cloud-optics = false -lw-scattering = false -cloud-mie = false -single-gpt = false +sw_two-stream = false +sw_raytracing = true +lw_raytracing = true +independent_column = false +liq_cloud-optics = false +ice_cloud-optics = false +lw_scattering = false +cloud_mie = false +single_gpt = false profiling = false -min-mfp-grid-ratio = true +min_mfp-grid-ratio = true tica = false [ints] -sw-raytracing = 256 -lw-raytracing = 22 -single-gpt = 1 +sw_raytracing = 256 +lw_raytracing = 22 +single_gpt = 1 [floats] -min-mfp-grid-ratio = 1.0 +min_mfp-grid-ratio = 1.0 From c65ba900554b21354490702c584a16288f98df4f Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 14:13:11 +0100 Subject: [PATCH 12/13] didn't do all - to _ substitues correctly --- allsky/allsky.ini | 14 +++++++------- rcemip/rcemip.ini | 10 +++++----- rfmip/rfmip.ini | 10 +++++----- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 68b4e10a..02052c2a 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -2,8 +2,8 @@ shortwave = true longwave = true fluxes = true -liquid_cloud-optics = true -ice_cloud-optics = true +liquid_cloud_optics = true +ice_cloud_optics = true aerosol_optics = false delta_cloud = true delta_aerosol = false @@ -12,17 +12,17 @@ delta_aerosol = false timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw_two-stream = false +sw_two_stream = false sw_raytracing = true lw_raytracing = true independent_column = false -liq_cloud-optics = false -ice_cloud-optics = false +liq_cloud_optics = false +ice_cloud_optics = false lw_scattering = false cloud_mie = false single_gpt = false profiling = false -min_mfp-grid-ratio = true +min_mfp_grid_ratio = true tica = false [ints] @@ -31,4 +31,4 @@ lw_raytracing = 22 single_gpt = 1 [floats] -min_mfp-grid-ratio = 1.0 +min_mfp_grid_ratio = 1.0 diff --git a/rcemip/rcemip.ini b/rcemip/rcemip.ini index 00363448..b7db423f 100644 --- a/rcemip/rcemip.ini +++ b/rcemip/rcemip.ini @@ -12,17 +12,17 @@ delta_aerosol = false timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw_two-stream = false +sw_two_stream = false sw_raytracing = true lw_raytracing = true independent_column = false -liq_cloud-optics = false -ice_cloud-optics = false +liq_cloud_optics = false +ice_cloud_optics = false lw_scattering = false cloud_mie = false single_gpt = false profiling = false -min_mfp-grid-ratio = true +min_mfp_grid_ratio = true tica = false [ints] @@ -31,4 +31,4 @@ lw_raytracing = 22 single_gpt = 1 [floats] -min_mfp-grid-ratio = 1.0 +min_mfp_grid_ratio = 1.0 diff --git a/rfmip/rfmip.ini b/rfmip/rfmip.ini index 00363448..b7db423f 100644 --- a/rfmip/rfmip.ini +++ b/rfmip/rfmip.ini @@ -12,17 +12,17 @@ delta_aerosol = false timings = false # raytracer (test_rte_rrtmgp_rt_gpu) -sw_two-stream = false +sw_two_stream = false sw_raytracing = true lw_raytracing = true independent_column = false -liq_cloud-optics = false -ice_cloud-optics = false +liq_cloud_optics = false +ice_cloud_optics = false lw_scattering = false cloud_mie = false single_gpt = false profiling = false -min_mfp-grid-ratio = true +min_mfp_grid_ratio = true tica = false [ints] @@ -31,4 +31,4 @@ lw_raytracing = 22 single_gpt = 1 [floats] -min_mfp-grid-ratio = 1.0 +min_mfp_grid_ratio = 1.0 From d5e5a8b19b4a58e03369b8c6d4646d5c1ee503a9 Mon Sep 17 00:00:00 2001 From: Menno Date: Wed, 4 Mar 2026 14:17:19 +0100 Subject: [PATCH 13/13] third attempt, remove duplicate entries in .ini files --- allsky/allsky.ini | 2 -- rcemip/rcemip.ini | 2 -- rfmip/rfmip.ini | 2 -- 3 files changed, 6 deletions(-) diff --git a/allsky/allsky.ini b/allsky/allsky.ini index 02052c2a..8e1fea55 100644 --- a/allsky/allsky.ini +++ b/allsky/allsky.ini @@ -16,8 +16,6 @@ sw_two_stream = false sw_raytracing = true lw_raytracing = true independent_column = false -liq_cloud_optics = false -ice_cloud_optics = false lw_scattering = false cloud_mie = false single_gpt = false diff --git a/rcemip/rcemip.ini b/rcemip/rcemip.ini index b7db423f..c81c2e49 100644 --- a/rcemip/rcemip.ini +++ b/rcemip/rcemip.ini @@ -16,8 +16,6 @@ sw_two_stream = false sw_raytracing = true lw_raytracing = true independent_column = false -liq_cloud_optics = false -ice_cloud_optics = false lw_scattering = false cloud_mie = false single_gpt = false diff --git a/rfmip/rfmip.ini b/rfmip/rfmip.ini index b7db423f..c81c2e49 100644 --- a/rfmip/rfmip.ini +++ b/rfmip/rfmip.ini @@ -16,8 +16,6 @@ sw_two_stream = false sw_raytracing = true lw_raytracing = true independent_column = false -liq_cloud_optics = false -ice_cloud_optics = false lw_scattering = false cloud_mie = false single_gpt = false