From 6f577a6acf58eb984fcf0fbf65df7973233a6164 Mon Sep 17 00:00:00 2001 From: xyuan Date: Wed, 22 Nov 2023 20:14:16 -0500 Subject: [PATCH 1/2] further refactor for GPU --- Source/Radiation/Aero_rad_props.H | 61 ++++---- Source/Radiation/Aero_rad_props.cpp | 94 +++++++----- Source/Radiation/Ebert_curry.H | 172 +++++++++++----------- Source/Radiation/Mam4_aero.H | 99 +++++++------ Source/Radiation/Optics.H | 60 ++++---- Source/Radiation/Optics.cpp | 75 +++++----- Source/Radiation/Rad_constants.H | 2 +- Source/Radiation/Radiation.H | 36 ++--- Source/Radiation/Radiation.cpp | 170 ++++++++++++--------- Source/Radiation/Rrtmgp.H | 32 ++-- Source/Radiation/Run_longwave_rrtmgp.cpp | 16 +- Source/Radiation/Run_shortwave_rrtmgp.cpp | 16 +- Source/Radiation/Slingo.H | 122 ++++++++------- 13 files changed, 520 insertions(+), 435 deletions(-) diff --git a/Source/Radiation/Aero_rad_props.H b/Source/Radiation/Aero_rad_props.H index de440576c..b3eb993b8 100644 --- a/Source/Radiation/Aero_rad_props.H +++ b/Source/Radiation/Aero_rad_props.H @@ -21,25 +21,30 @@ class AerRadProps { ~AerRadProps() = default; // initialization - void initialize(); + void initialize(int ngas_, int nmodes_, int num_aeros_, + int nswbands_, int nlwbands_, + int ncol_, int nlev_, int nrh_, int top_lev_, + const std::vector& aero_names_, + const real2d& zi_, const real2d& pmid_, const real2d& temp_, + const real2d& qi_, const real2d& geom_radius); void aer_rad_props_sw(const int& list_idx, const real& dt, const int& nnite, const int1d& idxnite, const bool is_cmip6_volc, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f, - real2d& clear_rh); + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f, + const real2d& clear_rh); void aer_rad_props_lw(const bool& is_cmip6_volc, const int& list_idx, const real& dt, const real2d& zi, - real3d& odap_aer, - real2d& clear_rh); + const real3d& odap_aer, + const real2d& clear_rh); void get_hygro_rad_props(const int& ncol, @@ -49,20 +54,20 @@ class AerRadProps { const real2d& ext, const real2d& ssa, const real2d& assm, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f); + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f); void get_nonhygro_rad_props(const int& ncol, const real2d& mass, const real1d& ext, const real1d& ssa, const real1d& assm, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f); + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f); void get_volcanic_radius_rad_props(const int& ncol, const real2d& mass, @@ -70,30 +75,30 @@ class AerRadProps { const real2d& r_scat, const real2d& r_ascat, const real1d& r_mu, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f); + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f); void volcanic_cmip_sw (const int1d& trop_level, const real2d& zi, // zone interface const real3d& ext_cmip6_sw_inv_m, const real3d& ssa_cmip6_sw, // ssa from the volcanic inputs const real3d& af_cmip6_sw, // asymmetry factor (af) from volcanic inputs - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f); + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f); void get_volcanic_rad_props(const int& ncol, const real2d& mass, const real1d& ext, const real1d& scat, const real1d& ascat, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f); + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f); void aer_vis_diag_out(const int& ncol, const int& nnite, diff --git a/Source/Radiation/Aero_rad_props.cpp b/Source/Radiation/Aero_rad_props.cpp index cb7037c61..69e97ac58 100644 --- a/Source/Radiation/Aero_rad_props.cpp +++ b/Source/Radiation/Aero_rad_props.cpp @@ -9,20 +9,39 @@ using yakl::fortran::parallel_for; using yakl::fortran::SimpleBounds; -void AerRadProps::initialize() { - +void AerRadProps::initialize(int ngas_, int nmodes_, int num_aeroes_, + int nswbands_, int nlwbands_, + int ncol_, int nlev_, int nrh_, int top_lev_, + const std::vector& aero_names_, + const real2d& zi_, const real2d& pmid_, const real2d& temp_, + const real2d& qt_, const real2d& geom_radius) { + nmodes = nmodes_; + ngas = ngas_; + num_aeroes = num_aeroes_; + aero_names = aero_names_; + + nswbands = nswbands_; + nlwbands = nlwbands_; + ncol = ncol_; + nlev = nlev_; + nrh = nrh_; + top_lev = top_lev_; + + pmid = pmid_; + pdeldry = pmid_; + temp = temp_; + qt = qt_; + // geometric radius + geometric_radius = geom_radius; + // vertical grid + zi = zi_; + // initialize mam4 aero model + mam_aer.initialize(ncol, nlev, top_lev, nswbands, nlwbands); } -void AerRadProps::aer_rad_props_sw(const int& list_idx, - const real& dt, - const int& nnite, - const int1d& idxnite, - const bool is_cmip6_volc, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f, - real2d& clear_rh) { +void AerRadProps::aer_rad_props_sw(const int& list_idx, const real& dt, const int& nnite, + const int1d& idxnite, const bool is_cmip6_volc, const real3d& tau, const real3d& tau_w, + const real3d& tau_w_g, const real3d& tau_w_f, const real2d& clear_rh) { // radiative properties for each aerosol real3d ta("ta", ncol, nlev, nswbands); real3d tw("tw", ncol, nlev, nswbands); @@ -30,7 +49,7 @@ void AerRadProps::aer_rad_props_sw(const int& list_idx, real3d twg("twg", ncol, nlev, nswbands); // aerosol masses - real2d aermmr; // mass mixing ratio of aerosols + real2d aermmr("aermmr", ncol, nlev); // mass mixing ratio of aerosols real2d mmr_to_mass("mmr_to_mass", ncol, nlev); // conversion factor for mmr to mass real2d aermass("aermass", ncol, nlev); // mass of aerosols @@ -112,9 +131,9 @@ void AerRadProps::aer_rad_props_sw(const int& list_idx, real stobie = 0.03 * temp(i,k) - std::log10(pmid(i, k)); if (pmid(i, k) <= 4000.) break; if (pmid(i, k) >= 55000.) continue; - if ((tlev == -1) || (stobie < strop(0))) { + if ((tlev == -1) || (stobie < strop(1))) { tlev = k; - strop(0) = stobie; + strop(1) = stobie; } } if (tlev != -1) trop_level(i) = tlev; @@ -127,7 +146,6 @@ void AerRadProps::aer_rad_props_sw(const int& list_idx, } } - // get number of bulk aerosols and number of modes in current list mam_consti.get_nmodes(list_idx, nmodes); mam_consti.get_ngas(list_idx, ngas); @@ -271,8 +289,8 @@ void AerRadProps::aer_rad_props_lw(const bool& is_cmip6_volc, const int& list_idx, const real& dt, const real2d& zi, - real3d& odap_aer, - real2d& clear_rh) { + const real3d& odap_aer, + const real2d& clear_rh) { int numaerosols; // number of bulk aerosols in climate/diagnostic list int nmodes; // number of aerosol modes in climate/diagnostic list std::string opticstype; // hygro or nonhygro @@ -470,10 +488,10 @@ void AerRadProps::aer_rad_props_lw(const bool& is_cmip6_volc, const real2d& extb, const real2d& ssab, const real2d& asmb, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f) { + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f) { parallel_for(SimpleBounds<3>(ncol,nlev,nswbands) , YAKL_LAMBDA (int icol, int ilev, int iswband) { auto ext1 = (1 + wrh(icol,ilev)) * extb(krh(icol,ilev)+1,iswband) - wrh(icol,ilev) * extb(krh(icol,ilev),iswband); @@ -492,10 +510,10 @@ void AerRadProps::get_nonhygro_rad_props(const int& ncol, const real1d& extb, const real1d& ssab, const real1d& asmb, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f) { + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f) { parallel_for(SimpleBounds<3>(ncol,nlev,nswbands) , YAKL_LAMBDA (int icol, int ilev, int iswband) { auto ext1 = extb(iswband); @@ -514,10 +532,10 @@ void AerRadProps::get_nonhygro_rad_props(const int& ncol, const real2d& r_scat, const real2d& r_ascat, const real1d& r_mu, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f) { + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f) { yakl::memset(tau, 0.); yakl::memset(tau_w, 0.); yakl::memset(tau_w_g, 0.); @@ -575,10 +593,10 @@ void AerRadProps::get_nonhygro_rad_props(const int& ncol, const real3d& ext_cmip6_sw_inv_m, const real3d& ssa_cmip6_sw, // ssa from the volcanic inputs const real3d& af_cmip6_sw, // asymmetry factor (af) from volcanic inputs - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f) { + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f) { // Above the tropopause, the read in values from the file include both the stratospheric // and volcanic aerosols. Therefore, we need to zero out taus above the tropopause // and populate them exclusively from the read in values. @@ -632,10 +650,10 @@ void AerRadProps::get_volcanic_rad_props(const int& ncol, const real1d& ext, const real1d& scat, const real1d& ascat, - real3d& tau, - real3d& tau_w, - real3d& tau_w_g, - real3d& tau_w_f) { + const real3d& tau, + const real3d& tau_w, + const real3d& tau_w_g, + const real3d& tau_w_f) { int nswbands; parallel_for(SimpleBounds<3>(nswbands, ncol,nlev), YAKL_LAMBDA (int iswband, int icol, int ilev) { real g = 0; diff --git a/Source/Radiation/Ebert_curry.H b/Source/Radiation/Ebert_curry.H index 27eac212b..5124a0855 100644 --- a/Source/Radiation/Ebert_curry.H +++ b/Source/Radiation/Ebert_curry.H @@ -8,49 +8,62 @@ class EbertCurry { public: static constexpr real scalefactor = 1.; //500._r8/917._r8 - static void ec_ice_optics_sw(int ncol, int nlev, int nswbands, real2d& cldn, real2d& cicewp, real2d& rei, - real3d& ice_tau, real3d& ice_tau_w, real3d& ice_tau_w_g, real3d& ice_tau_w_f) { + static void ec_ice_optics_sw(int ncol, int nlev, int nswbands, const real2d& cldn, const real2d& cicewp, const real2d& rei, + const real3d& ice_tau, const real3d& ice_tau_w, const real3d& ice_tau_w_g, const real3d& ice_tau_w_f) { real1d wavmin("wavmin", nswbands); real1d wavmax("wavmax", nswbands); // ice water coefficients (Ebert and Curry,1992, JGR, 97, 3831-3836) - std::vector abari = // a coefficient for extinction optical depth - { 3.448e-03, 3.448e-03, 3.448e-03, 3.448e-03 }; - std::vector bbari = // b coefficient for extinction optical depth - { 2.431, 2.431, 2.431, 2.431 }; - std::vector cbari = // c coefficient for single scat albedo - { 1.00e-05, 1.10e-04, 1.861e-02, .46658 }; - std::vector dbari = // d coefficient for single scat albedo - { 0.0, 1.405e-05, 8.328e-04, 2.05e-05 }; - std::vector ebari = // e coefficient for asymmetry parameter - { 0.7661, 0.7730, 0.794, 0.9595 }; - std::vector fbari = // f coefficient for asymmetry parameter - { 5.851e-04, 5.665e-04, 7.267e-04, 1.076e-04 }; - - real abarii; // A coefficient for current spectral band - real bbarii; // B coefficient for current spectral band - real cbarii; // C coefficient for current spectral band - real dbarii; // D coefficient for current spectral band - real ebarii; // E coefficient for current spectral band - real fbarii; // F coefficient for current spectral band + real1d abari("abari", 4); // a coefficient for extinction optical depth + real1d bbari("bbari", 4); // b coefficient for extinction optical depth + real1d cbari("cbari", 4); // c coefficient for single scat albedo + real1d dbari("dbari", 4); // d coefficient for single scat albedo + real1d ebari("ebari", 4); // e coefficient for asymmetry parameter + real1d fbari("fbari", 4); // f coefficient for asymmetry parameter + + parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int i) { + abari(1) = 3.448e-03; + abari(2) = 3.448e-03; + abari(3) = 3.448e-03; + abari(4) = 3.448e-03; + bbari(1) = 2.431; + bbari(2) = 2.431; + bbari(3) = 2.431; + bbari(4) = 2.431; + cbari(1) = 1.00e-05; + cbari(2) = 1.10e-04; + cbari(3) = 1.861e-02; + cbari(4) = 0.46658; + dbari(1) = 0.0; + dbari(2) = 1.405e-05; + dbari(3) = 8.328e-04; + dbari(4) = 2.05e-05; + ebari(1) = 0.7661; + ebari(2) = 0.7730; + ebari(3) = 0.794; + ebari(4) = 0.9595; + fbari(1) = 5.851e-04; + fbari(2) = 5.665e-04; + fbari(3) = 7.267e-04; + fbari(4) = 1.076e-04; + }); // Minimum cloud amount (as a fraction of the grid-box area) to // distinguish from clear sky - const real cldmin = 1.0e-80; + constexpr real cldmin = 1.0e-80; // Decimal precision of cloud amount (0 -> preserve full resolution; // 10^-n -> preserve n digits of cloud amount) - const real cldeps = 0.0; + constexpr real cldeps = 0.0; // Optical properties for ice are valid only in the range of // 13 < rei < 130 micron (Ebert and Curry 92) - const real rei_min = 13.; - const real rei_max = 130.; + constexpr real rei_min = 13.; + constexpr real rei_max = 130.; // integer :: ns, i, k, indxsl int indxsl; - real tmp1i, tmp2i, tmp3i, g; RadConstants::get_sw_spectral_boundaries(wavmin,wavmax,RadConstants::micrometer); @@ -66,46 +79,43 @@ class EbertCurry { indxsl = 4; } - abarii = abari[indxsl]; - bbarii = bbari[indxsl]; - cbarii = cbari[indxsl]; - dbarii = dbari[indxsl]; - ebarii = ebari[indxsl]; - fbarii = fbari[indxsl]; - - for (auto k=0; k rei > 130 micron (Ebert and Curry 92) - if (cldn(i,k) >= cldmin && cldn(i,k) >= cldeps) { - tmp1i = abarii + bbarii/std::max(rei_min,std::min(scalefactor*rei(i,k),rei_max)); - ice_tau(ns,i,k) = 1000. * cicewp(i,k) * tmp1i; - } else { - ice_tau(ns,i,k) = 0.0; - } - - tmp2i = 1. - cbarii - dbarii*std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max); - tmp3i = fbarii*std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max); - // Do not let single scatter albedo be 1. Delta-eddington solution - // for non-conservative case has different analytic form from solution - // for conservative case, and raddedmx is written for non-conservative case. - ice_tau_w(ns,i,k) = ice_tau(ns,i,k) * std::min(tmp2i,.999999); - g = ebarii + tmp3i; - ice_tau_w_g(ns,i,k) = ice_tau_w(ns,i,k) * g; - ice_tau_w_f(ns,i,k) = ice_tau_w(ns,i,k) * g * g; - - } // End do i=1,ncol - } // End do k=1,nlev + parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i) { + real tmp1i, tmp2i, tmp3i, g; + auto abarii = abari(indxsl); + auto bbarii = bbari(indxsl); + auto cbarii = cbari(indxsl); + auto dbarii = dbari(indxsl); + auto ebarii = ebari(indxsl); + auto fbarii = fbari(indxsl); + + // note that optical properties for ice valid only + // in range of 13 > rei > 130 micron (Ebert and Curry 92) + if (cldn(i,k) >= cldmin && cldn(i,k) >= cldeps) { + tmp1i = abarii + bbarii/std::max(rei_min,std::min(scalefactor*rei(i,k),rei_max)); + ice_tau(ns,i,k) = 1000. * cicewp(i,k) * tmp1i; + } else { + ice_tau(ns,i,k) = 0.0; + } + + tmp2i = 1. - cbarii - dbarii*std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max); + tmp3i = fbarii*std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max); + // Do not let single scatter albedo be 1. Delta-eddington solution + // for non-conservative case has different analytic form from solution + // for conservative case, and raddedmx is written for non-conservative case. + ice_tau_w(ns,i,k) = ice_tau(ns,i,k) * std::min(tmp2i,.999999); + g = ebarii + tmp3i; + ice_tau_w_g(ns,i,k) = ice_tau_w(ns,i,k) * g; + ice_tau_w_f(ns,i,k) = ice_tau_w(ns,i,k) * g * g; + }); } // nswbands } - static void ec_ice_optics_lw(int ncol, int nlev, int nlwbands, real2d& cldn, real2d& iclwpth, real2d& iciwpth, real2d& rei, real3d& abs_od) { + static void ec_ice_optics_lw(int ncol, int nlev, int nlwbands, const real2d& cldn, const real2d& iclwpth, + const real2d& iciwpth, const real2d& rei, const real3d& abs_od) { real2d ficemr("ficemr",ncol,nlev); real2d cwp("cwp",ncol,nlev); real2d cldtau("cldtau",ncol,nlev); - real kabs, kabsi; // longwave liquid absorption coeff (m**2/g) const real kabsl = 0.090361; @@ -115,33 +125,25 @@ class EbertCurry { const real rei_min = 13.; const real rei_max = 130.; - for (auto k=0; k rei > 130 micron (Ebert and Curry 92) - kabsi = 0.005 + 1./std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max); - kabs = kabsi*ficemr(i,k); // kabsl*(1.-ficemr(i,k)) + kabsi*ficemr(i,k); - cldtau(i,k) = kabs*cwp(i,k); - } - } - - for (auto lwband=0; lwband(nlev, ncol), YAKL_LAMBDA (int k, int i) { + cwp(i,k) = 1000.0 *iciwpth(i,k) + 1000.0 *iclwpth(i,k); + ficemr(i,k) = 1000.0*iciwpth(i,k)/(std::max(1.e-18,cwp(i,k))); + }); + + parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i) { + // Note from Andrew Conley: + // Optics for RK no longer supported, This is constructed to get + // close to bit for bit. Otherwise we could simply use ice water path + // note that optical properties for ice valid only + // in range of 13 > rei > 130 micron (Ebert and Curry 92) + auto kabsi = 0.005 + 1./std::min(std::max(rei_min,scalefactor*rei(i,k)),rei_max); + auto kabs = kabsi*ficemr(i,k); // kabsl*(1.-ficemr(i,k)) + kabsi*ficemr(i,k); + cldtau(i,k) = kabs*cwp(i,k); + }); + + parallel_for(SimpleBounds<3>(nlwbands, ncol, nlev), YAKL_LAMBDA (int lwband, int icol, int ilev) { abs_od(lwband,icol,ilev)=cldtau(icol,ilev); - } - } - } + }); } }; #endif diff --git a/Source/Radiation/Mam4_aero.H b/Source/Radiation/Mam4_aero.H index 28a0bc28f..c6add80f6 100644 --- a/Source/Radiation/Mam4_aero.H +++ b/Source/Radiation/Mam4_aero.H @@ -56,6 +56,16 @@ class Mam4_aer { // aerosol constitutents MamConstituents mam_consti; + // initialization + inline + void initialize(int ncol_, int nlev_, int top_lev_, int nswbands_, int nlwbands_) { + ncol = ncol_; + nlev = nlev_; + top_lev = top_lev_; + nswbands = nswbands_; + nlwbands = nlwbands_; + } + // read water refractive index file and set module data inline void read_water_refindex(std::string infilename) { @@ -81,8 +91,8 @@ class Mam4_aer { // modal size parameters inline void modal_size_parameters(int ncol, real sigma_logr_aer, - real2d& dgnumwet, real2d& radsurf, - real2d& logradsurf, real3d& cheb) { + const real2d& dgnumwet, const real2d& radsurf, + const real2d& logradsurf, const real3d& cheb) { real alnsg_amode; real explnsigma; real1d xrad("xrad", ncol); // normalized aerosol radius @@ -91,7 +101,7 @@ class Mam4_aer { explnsigma = std::exp(2.0*alnsg_amode*alnsg_amode); for(auto k = top_lev; k > nlev; --k) { - for(auto i = 0; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { // convert from number mode diameter to surface area radsurf(i,k) = 0.5*dgnumwet(i,k)*explnsigma; logradsurf(i,k) = std::log(radsurf(i,k)); @@ -114,8 +124,8 @@ class Mam4_aer { // bilinear interpolation of table inline void binterp(real3d& table, int ncol, int km, int im, int jm, - real1d& x, real1d& y, real1d& xtab, real1d& ytab, - int1d& ix, int1d& jy, real1d& t, real1d& u, real2d out) { + const real1d& x, const real1d& y, const real1d& xtab, const real1d& ytab, + const int1d& ix, const int1d& jy, const real1d& t, const real1d& u, real2d out) { real1d tu("tu", ncol), tuc("tuc", ncol), tcu("tcu", ncol), tcuc("tcuc", ncol); int i, j, ip1, jp1; real dx, dy; @@ -141,14 +151,14 @@ class Mam4_aer { } } else { - for(auto i=0; i 1) { - for(auto ic=0; ic(ncol, nlev), YAKL_LAMBDA (int i, int k) { + for (auto n = 1; n <= nmodes; ++n) { + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int k) { dgncur_a(i,k) = dgnum_m(i,k,n); }); @@ -237,8 +247,8 @@ class Mam4_aer { // need qmass*dummwdens = (kg/kg-air) * [1/(kg/m3)] = m3/kg-air dummwdens = 1.0 / specdens; - for (auto k=top_lev; k < nlev; ++k) { - for (auto i=0; i(ncol, nlev), YAKL_LAMBDA (int i, int j) { - tauxar(i,0,j) = 0.; - wa(i,0,j) = 0.925; - ga(i,0,j) = 0.850; - fa(i,0,j) = 0.7225; + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int i, int j) { + tauxar(i,1,j) = 0.; + wa(i,1,j) = 0.925; + ga(i,1,j) = 0.850; + fa(i,1,j) = 0.7225; }); //mass(:ncol,:) = state%pdeldry(:ncol,:)*rga @@ -401,7 +412,7 @@ class Mam4_aer { mam_consti.get_nmodes(list_idx, nmodes); for(auto m = 0; m < nmodes; ++m) { - yakl::c::parallel_for(yakl::c::Bounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { dgnumwet(icol,ilev) = dgnumwet_m(icol,ilev,m); qaerwat(icol,ilev) = qaerwat_m(icol,ilev,m); }); @@ -417,20 +428,20 @@ class Mam4_aer { modal_size_parameters(ncol, sigma_logr_aer, dgnumwet, radsurf, logradsurf, cheb); int nswbands; // added by xyuan to be able compile the code - for(auto isw = 1; isw < nswbands; isw++) { - for(auto k = top_lev; k < nlev; ++k) { + for(auto isw = 1; isw <= nswbands; isw++) { + for(auto k = top_lev; k <= nlev; ++k) { // form bulk refractive index yakl::memset(crefin_real, 0.); yakl::memset(crefin_im, 0.); yakl::memset(dryvol, 0.); // aerosol species loop - for(auto l = 1; l < nspec; ++l) { + for(auto l = 1; l <= nspec; ++l) { // mam_consti.rad_cnst_get_aer_mmr(list_idx, m, l, "a", specmmr); mam_consti.get_mam_props(list_idx, m, l, specdens, spectype, hygro_aer, specrefindex_real, specrefindex_im); - for(auto i = 0; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { vol(i) = specmmr(i,k)/specdens; dryvol(i) = dryvol(i) + vol(i); crefin_real(i) = crefin_real(i) + vol(i)*specrefindex_real(isw); @@ -438,7 +449,7 @@ class Mam4_aer { } } - for(auto i = 1; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { watervol(i) = qaerwat(i,k)/rhoh2o; wetvol(i) = watervol(i) + dryvol(i); if (watervol(i) < 0.) { @@ -485,7 +496,7 @@ class Mam4_aer { itab, jtab, ttab, utab, casm); // parameterized optical properties - for(auto i=1; i < ncol; ++i) { + for(auto i=1; i <= ncol; ++i) { if (logradsurf(i,k) <= xrmax) { pext(i) = 0.5*cext(i,1); @@ -517,7 +528,7 @@ class Mam4_aer { dopaer(i) = pext(i)*mass(i,k); } - for(auto i = 1; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { if ((dopaer(i) <= -1.e-10) || (dopaer(i) >= 30.)) { if (dopaer(i) <= -1.e-10) { @@ -540,7 +551,7 @@ class Mam4_aer { } } - for(auto i=1; i < ncol; ++i) { + for(auto i=1; i <= ncol; ++i) { tauxar(i,k,isw) = tauxar(i,k,isw) + dopaer(i); wa(i,k,isw) = wa(i,k,isw) + dopaer(i)*palb(i); ga(i,k,isw) = ga(i,k,isw) + dopaer(i)*palb(i)*pasm(i); @@ -556,7 +567,7 @@ class Mam4_aer { // // calcualtes aerosol lw radiative properties // - void modal_aero_lw(int list_idx, real dt, real3d& tauxar, real2d& clear_rh) { + void modal_aero_lw(int list_idx, real dt, const real3d& tauxar, const real2d& clear_rh) { real sigma_logr_aer; // geometric standard deviation of number distribution real alnsg_amode; real1d xrad("xrad", ncol); @@ -636,7 +647,7 @@ class Mam4_aer { // some intermediate results are saved and the chebyshev polynomials are stored // in a array with different index order. Could be unified. for(auto k = top_lev; k > nlev; --k) { - for(auto i = 0; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { alnsg_amode = std::log( sigma_logr_aer ); // convert from number diameter to surface area xrad(i) = std::log(0.5*dgnumwet(i,k)) + 2.0*alnsg_amode*alnsg_amode; @@ -654,7 +665,7 @@ class Mam4_aer { } int nlwbands; // add by xyuan to be compiled - for(auto ilw = 0; ilw < nlwbands; ++ilw) { + for(auto ilw = 1; ilw <= nlwbands; ++ilw) { for(auto k = top_lev; k > nlev; --k) { // form bulk refractive index. Use volume mixing for infrared yakl::memset(crefinr, 0.0); @@ -662,11 +673,11 @@ class Mam4_aer { yakl::memset(dryvol, 0.0); // aerosol species loop - for(auto l = 1; l < nspec; ++l) { + for(auto l = 1; l <= nspec; ++l) { // rad_cnst_get_aer_mmr(list_idx, m, l, 'a', state, pbuf, specmmr); mam_consti.get_mam_props_lw(list_idx, m, l, specdens, specrefrindex, specrefiindex); - for(auto i = 0; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { vol(i) = specmmr(i,k)/specdens; dryvol(i) = dryvol(i) + vol(i); crefinr(i) = crefinr(i) + vol(i)*specrefrindex(ilw); @@ -674,7 +685,7 @@ class Mam4_aer { } } - for(auto i = 1; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { watervol(i) = qaerwat(i,k)/rhoh2o; wetvol(i) = watervol(i) + dryvol(i); if (watervol(i) < 0.0) { @@ -713,7 +724,7 @@ class Mam4_aer { itab, jtab, ttab, utab, cabs); // parameterized optical properties - for(auto i = 1; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { pabs(i) = 0.5*cabs(i,1); for(auto nc = 2; nc < ncoef; ++nc) { pabs(i) = pabs(i) + cheby(nc,i,k)*cabs(i,nc); @@ -723,14 +734,14 @@ class Mam4_aer { dopaer(i) = pabs(i)*mass(i,k); } - for(auto i = 0; i < ncol; ++i) { + for(auto i = 1; i <= ncol; ++i) { if ((dopaer(i) <= -1.e-10) || (dopaer(i) >= 20.)) { if (dopaer(i) <= -1.e-10) amrex::Print() << "ERROR: Negative aerosol optical depth in this layer."; else amrex::Print() << "WARNING: Aerosol optical depth is unreasonably high in this layer."; - for(auto l = 1; l < nspec; ++l) { + for(auto l = 1; l <= nspec; ++l) { // rad_cnst_get_aer_mmr(list_idx, m, l, 'a', state, pbuf, specmmr); mam_consti.get_mam_props_lw(list_idx, m, l, specdens,specrefrindex, specrefiindex); volf = specmmr(i,k)/specdens; @@ -743,7 +754,7 @@ class Mam4_aer { } } } - for(auto i = 1; i < ncol; ++i) tauxar(i,k,ilw) = tauxar(i,k,ilw) + dopaer(i); + for(auto i = 1; i <= ncol; ++i) tauxar(i,k,ilw) = tauxar(i,k,ilw) + dopaer(i); } // k = top_lev, nlev } // nlwbands } // m = 1, nmodes diff --git a/Source/Radiation/Optics.H b/Source/Radiation/Optics.H index 2d209aafb..51722c433 100644 --- a/Source/Radiation/Optics.H +++ b/Source/Radiation/Optics.H @@ -31,52 +31,54 @@ class Optics { // get short wave cloud optics property void get_cloud_optics_sw(int ncol, int nlev, int nbnd, - bool do_snow, real2d cld, real2d cldfsnow, real2d iclwp, - real2d iciwp, real2d icswp, real2d lambdac, real2d mu, - real2d dei, real2d des, real2d rel, real2d rei, - real3d tau_out, real3d ssa_out, real3d asm_out, - real3d liq_tau_out, real3d ice_tau_out, real3d snw_tau_out); - + bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, + const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, + const real2d& dei, const real2d& des, const real2d& rel, const real2d& rei, + const real3d& tau_out, const real3d& ssa_out, const real3d& asm_out, + const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out); // get long wave cloud optics property - void get_cloud_optics_lw(int ncol, int nlev, int nbnd, bool do_snow, - real2d& cld, real2d& cldfsnow, real2d& iclwp, real2d& iciwp, real2d& icswp, - real2d& lambdac, real2d& mu, real2d& dei, real2d& des, real2d& rei, - real3d& tau_out, real3d& liq_tau_out, real3d& ice_tau_out, real3d& snw_tau_out); + void get_cloud_optics_lw(int ncol, int nlev, int nbnd, bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, + const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, const real2d& dei, const real2d& des, + const real2d& rei, const real3d& tau_out, const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out); // sample short wave cloud optics property - void sample_cloud_optics_sw(int ncol, int nlev, int ngpt, int1d& gpt2bnd, - real2d& pmid, real2d& cld, real2d& cldfsnow, - real3d& tau_bnd, real3d& ssa_bnd, real3d& asm_bnd, - real3d& tau_gpt, real3d& ssa_gpt, real3d& asm_gpt); + void sample_cloud_optics_sw(int ncol, int nlev, int ngpt, const int1d& gpt2bnd, + const real2d& pmid, const real2d& cld, const real2d& cldfsnow, + const real3d& tau_bnd, const real3d& ssa_bnd, const real3d& asm_bnd, + const real3d& tau_gpt, const real3d& ssa_gpt, const real3d& asm_gpt); // sample long wave cloud optics property - void sample_cloud_optics_lw( - int ncol, int nlev, int ngpt, int1d& gpt2bnd, - real2d& pmid, real2d& cld, real2d& cldfsnow, - real3d& tau_bnd, real3d& tau_gpt); + void sample_cloud_optics_lw(int ncol, int nlev, int ngpt, const int1d& gpt2bnd, + const real2d& pmid, const real2d& cld, const real2d& cldfsnow, + const real3d& tau_bnd, const real3d& tau_gpt); // set the short wave aerosol optics property - void set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, real dt, int1d& night_indices, - bool is_cmip6_volc, - real3d& tau_out, real3d& ssa_out, real3d& asm_out, - real2d& clear_rh); + void set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, real dt, const int1d& night_indices, + bool is_cmip6_volc, const real3d& tau_out, const real3d& ssa_out, const real3d& asm_out, + const real2d& clear_rh); // set the long wave aerosol optics property - void set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, real2d& zi, real3d& tau, real2d& clear_rh); + void set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, const real2d& zi, + const real3d& tau, const real2d& clear_rh); // mcica subcol mask - void mcica_subcol_mask(int ngpt, int ncol, int nlev, real2d& cldfrac, bool3d& iscloudy); + void mcica_subcol_mask(int ngpt, int ncol, int nlev, const real2d& cldfrac, const bool3d& iscloudy); // combine properties void combine_properties(int nbands, int ncols, int nlevs, - real2d& fraction1, real3d& property1, - real2d& fraction2, real3d& property2, - real3d& combined_property); + const real2d& fraction1, const real3d& property1, + const real2d& fraction2, const real3d& property2, + const real3d& combined_property); // initialize and load gas property data for rrtmgp radiation - void initialize(); - + void initialize(int ngas, int nmodes, int num_aeros, + int nswbands, int nlwbands, + int ncol, int nlev, int nrh, int top_lev, + const std::vector& aero_names, + const real2d& zi, const real2d& pmid, const real2d& temp, + const real2d& qi, const real2d& geom_radius); + // finalize/clean up void finalize(); diff --git a/Source/Radiation/Optics.cpp b/Source/Radiation/Optics.cpp index db8f810fe..85c861f7a 100644 --- a/Source/Radiation/Optics.cpp +++ b/Source/Radiation/Optics.cpp @@ -9,9 +9,16 @@ using yakl::fortran::parallel_for; using yakl::fortran::SimpleBounds; -void Optics::initialize() { +void Optics::initialize(int ngas, int nmodes, int num_aeros, + int nswbands, int nlwbands, + int ncol, int nlev, int nrh, int top_lev, + const std::vector& aero_names, + const real2d& zi, const real2d& pmid, const real2d& temp, + const real2d& qi, const real2d& geom_radius) { cloud_optics.initialize(); - aero_optics.initialize(); + aero_optics.initialize(ngas, nmodes, num_aeros, + nswbands, nlwbands, ncol, nlev, nrh, top_lev, + aero_names, zi, pmid, temp, qi, geom_radius); } void Optics::finalize() { @@ -20,15 +27,13 @@ void Optics::finalize() { void Optics::get_cloud_optics_sw(int ncol, int nlev, int nbnd, - bool do_snow, real2d cld, real2d cldfsnow, real2d iclwp, - real2d iciwp, real2d icswp, real2d lambdac, real2d mu, - real2d dei, real2d des, real2d rel, real2d rei, - real3d tau_out, real3d ssa_out, real3d asm_out, - real3d liq_tau_out, real3d ice_tau_out, real3d snw_tau_out) { - + bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, + const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, + const real2d& dei, const real2d& des, const real2d& rel, const real2d& rei, + const real3d& tau_out, const real3d& ssa_out, const real3d& asm_out, + const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out) { //Temporary variables to hold cloud optical properties before combining into // output arrays. Same shape as output arrays, so get shapes from output. - //real(r8), dimension(nbnd,ncol,nlev) :: real3d liq_tau("liq_tau",nbnd, ncol, nlev); real3d liq_tau_ssa("liq_tau_ssa",nbnd, ncol, nlev); real3d liq_tau_ssa_g("liq_tau_ssa_g", nbnd, ncol, nlev); @@ -170,12 +175,10 @@ void Optics::get_cloud_optics_sw(int ncol, int nlev, int nbnd, //---------------------------------------------------------------------------- void Optics::get_cloud_optics_lw( - int ncol, int nlev, int nbnd, bool do_snow, real2d& cld, real2d& cldfsnow, real2d& iclwp, real2d& iciwp, real2d& icswp, - real2d& lambdac, real2d& mu, real2d& dei, real2d& des, real2d& rei, - real3d& tau_out, real3d& liq_tau_out, real3d& ice_tau_out, real3d& snw_tau_out) { - + int ncol, int nlev, int nbnd, bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, + const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, const real2d& dei, const real2d& des, + const real2d& rei, const real3d& tau_out, const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out) { // Temporary variables to hold absorption optical depth - //real(r8), dimension(nbnd,ncol,nlev) :: & real3d ice_tau("ice_tau", nbnd, ncol, nlev); real3d liq_tau("liq_tau", nbnd, ncol, nlev); real3d snow_tau("snow_tau", nbnd, ncol, nlev); @@ -243,9 +246,9 @@ void Optics::get_cloud_optics_lw( // optical properties, we weight the cloud and snow properties by the fraction // of cloud and snow present. void Optics::combine_properties(int nbands, int ncols, int nlevs, - real2d& fraction1, real3d& property1, - real2d& fraction2, real3d& property2, - real3d& combined_property) { + const real2d& fraction1, const real3d& property1, + const real2d& fraction2, const real3d& property2, + const real3d& combined_property) { // Combined fraction (max of property 1 and 2) real2d combined_fraction("combined_fraction", ncols,nlevs); @@ -270,10 +273,10 @@ void Optics::combine_properties(int nbands, int ncols, int nlevs, // Do MCICA sampling of optics here. This will map bands to gpoints, // while doing stochastic sampling of cloud state void Optics::sample_cloud_optics_sw( - int ncol, int nlev, int ngpt, int1d& gpt2bnd, - real2d& pmid, real2d& cld, real2d& cldfsnow, - real3d& tau_bnd, real3d& ssa_bnd, real3d& asm_bnd, - real3d& tau_gpt, real3d& ssa_gpt, real3d& asm_gpt) { + int ncol, int nlev, int ngpt, const int1d& gpt2bnd, + const real2d& pmid, const real2d& cld, const real2d& cldfsnow, + const real3d& tau_bnd, const real3d& ssa_bnd, const real3d& asm_bnd, + const real3d& tau_gpt, const real3d& ssa_gpt, const real3d& asm_gpt) { //real(r8), dimension(ncol,nlev) :: combined_cld real2d combined_cld("combined_cld", ncol, nlev); @@ -306,10 +309,9 @@ void Optics::sample_cloud_optics_sw( //---------------------------------------------------------------------------- // Do MCICA sampling of optics here. This will map bands to gpoints, // while doing stochastic sampling of cloud state -void Optics::sample_cloud_optics_lw( - int ncol, int nlev, int ngpt, int1d& gpt2bnd, - real2d& pmid, real2d& cld, real2d& cldfsnow, - real3d& tau_bnd, real3d& tau_gpt) { +void Optics::sample_cloud_optics_lw(int ncol, int nlev, int ngpt, const int1d& gpt2bnd, + const real2d& pmid, const real2d& cld, const real2d& cldfsnow, + const real3d& tau_bnd, const real3d& tau_gpt) { //real(r8), dimension(ncol,nlev) :: combined_cld real2d combined_cld("combined_cld", ncol, nlev); @@ -338,10 +340,8 @@ void Optics::sample_cloud_optics_lw( } //---------------------------------------------------------------------------- -void Optics::set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, real dt, int1d& night_indices, - bool is_cmip6_volc, - real3d& tau_out, real3d& ssa_out, real3d& asm_out, - real2d& clear_rh) { +void Optics::set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, real dt, const int1d& night_indices, + bool is_cmip6_volc, const real3d& tau_out, const real3d& ssa_out, const real3d& asm_out, const real2d& clear_rh) { // NOTE: aer_rad_props expects 0:pver indexing on these! It appears this is to // account for the extra layer added above model top, but it is not entirely // clear. This is not done for the longwave, and it is not really documented @@ -393,7 +393,8 @@ void Optics::set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, } //---------------------------------------------------------------------------- -void Optics::set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, real2d& zi, real3d& tau, real2d& clear_rh) { +void Optics::set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, const real2d& zi, + const real3d& tau, const real2d& clear_rh) { // Get aerosol absorption optical depth from CAM routine yakl::memset(tau, 0.); aero_optics.aer_rad_props_lw(is_cmip6_volc, icall, dt, zi, tau, clear_rh); @@ -406,7 +407,7 @@ void Optics::set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, real2 // of the layer above it. // void Optics::mcica_subcol_mask(int ngpt, int ncol, int nlev, - real2d& cldfrac, bool3d& iscloudy) { + const real2d& cldfrac, const bool3d& iscloudy) { // Local vars const real cldmin = 1.0e-80; // min cloud fraction real2d cldf("cldf",ncol,nlev); // cloud fraction clipped to cldmin @@ -430,12 +431,14 @@ void Optics::mcica_subcol_mask(int ngpt, int ncol, int nlev, // - if the layer above is cloudy, use the same random number as in the layer above // - if the layer above is clear, use a new random number parallel_for(SimpleBounds<3>(nlev, ncol, ngpt), YAKL_LAMBDA (int k, int i, int isubcol) { - if (cdf(isubcol,i,k-1) > 1. - cldf(i,k-1) ) { - cdf(isubcol,i,k) = cdf(isubcol,i,k-1); - } else { - cdf(isubcol,i,k) = cdf(isubcol,i,k) * (1. - cldf(i,k-1)); + if (k > 1) { + if (cdf(isubcol,i,k-1) > 1. - cldf(i,k-1) ) { + cdf(isubcol,i,k) = cdf(isubcol,i,k-1); + } else { + cdf(isubcol,i,k) = cdf(isubcol,i,k) * (1. - cldf(i,k-1)); + } + iscloudy(isubcol,i,k) = cdf(isubcol,i,k) >= 1.-cldf(i,k) ? true : false; } - iscloudy(isubcol,i,k) = cdf(isubcol,i,k) >= 1.-cldf(i,k) ? true : false; }); } diff --git a/Source/Radiation/Rad_constants.H b/Source/Radiation/Rad_constants.H index 0ca232786..68e961cb5 100644 --- a/Source/Radiation/Rad_constants.H +++ b/Source/Radiation/Rad_constants.H @@ -59,7 +59,7 @@ class RadConstants { // These are indices to the band for diagnostic output static constexpr int idx_sw_diag = 10; // index to sw visible band static constexpr int idx_nir_diag = 8; // index to sw near infrared (778-1240 nm) band - static constexpr int idx_uv_diag = 11; // index to sw uv (345-441 nm) band + static constexpr int idx_uv_diag = 11; // index to sw uv (345-441 nm) band static constexpr int rrtmg_sw_cloudsim_band = 9; // rrtmg band for .67 micron diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H index deeb8d8ca..0eb4e1beb 100644 --- a/Source/Radiation/Radiation.H +++ b/Source/Radiation/Radiation.H @@ -59,31 +59,31 @@ class Radiation { void on_complete(); void radiation_driver_lw(int ncol, int nlev, - real3d& gas_vmr, - real2d& pmid, real2d& pint, real2d& tmid, real2d& tint, - real3d& cld_tau_gpt, real3d& aer_tau_bnd, + const real3d& gas_vmr, + const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint, + const real3d& cld_tau_gpt, const real3d& aer_tau_bnd, FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, - real2d& qrl, real2d& qrlc); + const real2d& qrl, const real2d& qrlc); void radiation_driver_sw(int ncol, - real3d& gas_vmr, real2d& pmid, real2d& pint, real2d& tmid, - real2d& albedo_dir, real2d& albedo_dif, real1d& coszrs, - real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, - real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, + const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid, + const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs, + const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt, + const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd, FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, - real2d& qrs, real2d& qrsc); + const real2d& qrs, const real2d& qrsc); - void set_daynight_indices(real1d& coszrs, - int1d& day_indices, - int1d& night_indices); + void set_daynight_indices(const real1d& coszrs, + const int1d& day_indices, + const int1d& night_indices); void get_gas_vmr(const std::vector& gas_names, - real3d& gas_vmr); + const real3d& gas_vmr); - void calculate_heating_rate(real2d& flux_up, - real2d& flux_dn, - real2d& pint, - real2d& heating_rate); + void calculate_heating_rate(const real2d& flux_up, + const real2d& flux_dn, + const real2d& pint, + const real2d& heating_rate); private: // geometry amrex::Geometry m_geom; @@ -194,7 +194,7 @@ class Radiation { real2d qrl; // longwave radiative heating rate // Pointers to fields on the physics buffer - real2d zi; + real2d zi, qt; real2d clear_rh; // Clear-sky heating rates are not on the physics buffer, and we have no diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index 135356fe5..55ddb4f5b 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -32,7 +32,7 @@ namespace internal { fluxes.bnd_flux_net = real3d("flux_net", nz, nlay+1, nbands); } - void expand_day_fluxes(FluxesByband& daytime_fluxes, FluxesByband& expanded_fluxes, int1d& day_indices) { + void expand_day_fluxes(FluxesByband& daytime_fluxes, FluxesByband& expanded_fluxes, const int1d& day_indices) { auto ncol = size(daytime_fluxes.bnd_flux_up, 1); auto nlev = size(daytime_fluxes.bnd_flux_up, 2); @@ -65,7 +65,7 @@ namespace internal { } // Utility function to reorder an array given a new indexing - void reordered(real1d& array_in, int1d& new_indexing, real1d& array_out) { + void reordered(const real1d& array_in, const int1d& new_indexing, const real1d& array_out) { // Reorder array based on input index mapping, which maps old indices to new parallel_for(SimpleBounds<1>(size(array_in, 1)), YAKL_LAMBDA (int i) { array_out(i) = array_in(new_indexing(i)); @@ -112,7 +112,6 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, radiation.initialize(ngas, active_gases, rrtmgp_coefficients_file_sw.c_str(), rrtmgp_coefficients_file_lw.c_str()); - optics.initialize(); // initialize the radiation data nswbands = radiation.get_nband_sw(); @@ -135,6 +134,16 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, pint = real2d("pint", ncol, nlev+1); tint = real2d("tint", ncol, nlev+1); + qt = real2d("qt", ncol, nlev); + zi = real2d("zi", ncol, nlev); + + yakl::memset(tmid, 300); + yakl::memset(pmid, 1000); + yakl::memset(pint, 1000); + yakl::memset(tint, 300); + yakl::memset(qt, 1.0e-3); + yakl::memset(zi, 1.0); + albedo_dir = real2d("albedo_dir", nswbands, ncol); albedo_dif = real2d("albedo_dif", nswbands, ncol); @@ -146,6 +155,18 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, qrsc = real2d("qrsc", ncol, nlev); qrlc = real2d("qrlc", ncol, nlev); + int nmodes = 3; + int nrh = 1; + int top_lev = 1; + naer = 4; + std::vector aero_names {"H2O", "N2", "O2", "O3"}; + auto geom_radius = real2d("geom_radius", ncol, nlev); + yakl::memset(geom_radius, 0.1); + + optics.initialize(ngas, nmodes, naer, nswbands, nlwbands, + ncol, nlev, nrh, top_lev, aero_names, zi, + pmid, tmid, qt, geom_radius); + amrex::Print() << " LW coefficents file: \n" << " SW coefficents file: \n" << " Frequency (timesteps) of Shortwave Radiation calc: \n " @@ -248,13 +269,14 @@ void Radiation::run() { // We need to fix band ordering because the old input files assume RRTMG // band ordering, but this has changed in RRTMGP. // TODO: fix the input files themselves! + real1d cld_tau_bnd_sw_1d("cld_tau_bnd_sw_1d", nswbands); + real1d cld_ssa_bnd_sw_1d("cld_ssa_bnd_sw_1d", nswbands); + real1d cld_asm_bnd_sw_1d("cld_asm_bnd_sw_1d", nswbands); + real1d cld_tau_bnd_sw_o_1d("cld_tau_bnd_sw_1d", nswbands); + real1d cld_ssa_bnd_sw_o_1d("cld_ssa_bnd_sw_1d", nswbands); + real1d cld_asm_bnd_sw_o_1d("cld_asm_bnd_sw_1d", nswbands); + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay) { - real1d cld_tau_bnd_sw_1d("cld_tau_bnd_sw_1d", nswbands); - real1d cld_ssa_bnd_sw_1d("cld_ssa_bnd_sw_1d", nswbands); - real1d cld_asm_bnd_sw_1d("cld_asm_bnd_sw_1d", nswbands); - real1d cld_tau_bnd_sw_o_1d("cld_tau_bnd_sw_1d", nswbands); - real1d cld_ssa_bnd_sw_o_1d("cld_ssa_bnd_sw_1d", nswbands); - real1d cld_asm_bnd_sw_o_1d("cld_asm_bnd_sw_1d", nswbands); for (auto ibnd = 1; ibnd <= nswbands; ++ibnd) { cld_tau_bnd_sw_1d(ibnd) = cld_tau_bnd_sw(icol,ilay,ibnd); cld_ssa_bnd_sw_1d(ibnd) = cld_ssa_bnd_sw(icol,ilay,ibnd); @@ -283,11 +305,13 @@ void Radiation::run() { // Aerosol needs night indices // TODO: remove this dependency, it's just used to mask aerosol outputs set_daynight_indices(coszrs, day_indices, night_indices); - int nday = 0; - int nnight = 0; - for (auto icol=0; icol 0) nday++; - if (night_indices(icol) > 0) nnight++; + int1d nday("nday",1); + int1d nnight("nnight",1); + yakl::memset(nday, 0); + yakl::memset(nnight, 0); + for (auto icol=1; icol 0) nday(1)++; + if (night_indices(icol) > 0) nnight(1)++; } for (auto icall = ngas /*N_DIAG*/; icall > 0; --icall) { @@ -306,14 +330,15 @@ void Radiation::run() { // Now reorder bands to be consistent with RRTMGP // TODO: fix the input files themselves! - parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay) { - real1d aer_tau_bnd_sw_1d("cld_tau_bnd_sw_1d", nswbands); - real1d aer_ssa_bnd_sw_1d("cld_ssa_bnd_sw_1d", nswbands); - real1d aer_asm_bnd_sw_1d("cld_asm_bnd_sw_1d", nswbands); - real1d aer_tau_bnd_sw_o_1d("cld_tau_bnd_sw_1d", nswbands); - real1d aer_ssa_bnd_sw_o_1d("cld_ssa_bnd_sw_1d", nswbands); - real1d aer_asm_bnd_sw_o_1d("cld_asm_bnd_sw_1d", nswbands); - for (auto ibnd = 0; ibnd < nswbands; ++ibnd) { + real1d aer_tau_bnd_sw_1d("cld_tau_bnd_sw_1d", nswbands); + real1d aer_ssa_bnd_sw_1d("cld_ssa_bnd_sw_1d", nswbands); + real1d aer_asm_bnd_sw_1d("cld_asm_bnd_sw_1d", nswbands); + real1d aer_tau_bnd_sw_o_1d("cld_tau_bnd_sw_1d", nswbands); + real1d aer_ssa_bnd_sw_o_1d("cld_ssa_bnd_sw_1d", nswbands); + real1d aer_asm_bnd_sw_o_1d("cld_asm_bnd_sw_1d", nswbands); + + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay) { + for (auto ibnd = 1; ibnd < nswbands; ++ibnd) { aer_tau_bnd_sw_1d(ibnd) = aer_tau_bnd_sw(icol,ilay,ibnd); aer_ssa_bnd_sw_1d(ibnd) = aer_ssa_bnd_sw(icol,ilay,ibnd); aer_asm_bnd_sw_1d(ibnd) = aer_asm_bnd_sw(icol,ilay,ibnd); @@ -321,7 +346,7 @@ void Radiation::run() { internal::reordered(aer_tau_bnd_sw_1d, rrtmg_to_rrtmgp, aer_tau_bnd_sw_o_1d); internal::reordered(aer_ssa_bnd_sw_1d, rrtmg_to_rrtmgp, aer_ssa_bnd_sw_o_1d); internal::reordered(aer_asm_bnd_sw_1d, rrtmg_to_rrtmgp, aer_asm_bnd_sw_o_1d); - for (auto ibnd = 0; ibnd < nswbands; ++ibnd) { + for (auto ibnd = 1; ibnd < nswbands; ++ibnd) { aer_tau_bnd_sw(icol,ilay,ibnd) = aer_tau_bnd_sw_o_1d(ibnd); aer_ssa_bnd_sw(icol,ilay,ibnd) = aer_ssa_bnd_sw_o_1d(ibnd); aer_asm_bnd_sw(icol,ilay,ibnd) = aer_asm_bnd_sw_o_1d(ibnd); @@ -421,11 +446,11 @@ void Radiation::run() { } } -void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, - real2d& pmid, real2d& pint, real2d& tmid, real2d& albedo_dir, real2d& albedo_dif, - real1d& coszrs, real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, - real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, - FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, real2d& qrs, real2d& qrsc) { +void Radiation::radiation_driver_sw(int ncol, const real3d& gas_vmr, + const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& albedo_dir, const real2d& albedo_dif, + const real1d& coszrs, const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt, + const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd, + FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky, const real2d& qrs, const real2d& qrsc) { // Incoming solar radiation, scaled for solar zenith angle // and earth-sun distance real2d solar_irradiance_by_gpt("solar_irradiance_by_gpt",ncol,nswgpts); @@ -462,7 +487,6 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, // climates as well (i.e., paleoclimate simulations) real tsi_scaling; - if (fixed_total_solar_irradiance<0) { // Get orbital eccentricity factor to scale total sky irradiance // tsi_scaling = get_eccentricity_factor(); @@ -478,16 +502,25 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, // computational cost (and because RRTMGP will fail for cosine solar zenith // angles less than or equal to zero) set_daynight_indices(coszrs, day_indices, night_indices); - int nday = 0; - int nnight = 0; - for (auto icol=0; icol 0) nday++; - if (night_indices(icol) > 0) nnight++; - } + int1d nday("nday",1); + int1d nnight("nnight", 1); + yakl::memset(nday, 0); + yakl::memset(nnight, 0); + parallel_for(SimpleBounds<1>(ncol), YAKL_LAMBDA (int icol) { + if (day_indices(icol) > 0) nday(1)++; + if (night_indices(icol) > 0) nnight(1)++; + }); + + intHost1d num_day("num_day",1); + intHost1d num_night("num_night",1); + nday.deep_copy_to(num_day); + nnight.deep_copy_to(num_night); + +std::cout << "num_days= " << num_day(1) << std::endl; // If no daytime columns in this chunk, then we return zeros - if (nday == 0) { -// reset_fluxes(fluxes_allsky) + if (num_day(1) == 0) { +// reset_fluxes(fluxes_allsky) // reset_fluxes(fluxes_clrsky) yakl::memset(qrs, 0.); yakl::memset(qrsc, 0.); @@ -495,7 +528,7 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, } // Compress to daytime-only arrays - parallel_for(SimpleBounds<3>(nday, nlev, nswgpts), YAKL_LAMBDA (int iday, int ilev, int igpt) { + parallel_for(SimpleBounds<3>(num_day(1), nlev, nswgpts), YAKL_LAMBDA (int iday, int ilev, int igpt) { auto icol = day_indices(iday); tmid_day(iday,ilev) = tmid(icol,ilev); pmid_day(iday,ilev) = pmid(icol,ilev); @@ -517,8 +550,8 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, // dimension nlev_rad+1, while we initialized the RRTMGP input variables to // have vertical dimension nlev_rad (defined at midpoints). FluxesByband fluxes_clrsky_day, fluxes_allsky_day; - internal::initial_fluxes(nday, nlev+1, nswbands, fluxes_allsky_day); - internal::initial_fluxes(nday, nlev+1, nswbands, fluxes_clrsky_day); + internal::initial_fluxes(num_day(1), nlev+1, nswbands, fluxes_allsky_day); + internal::initial_fluxes(num_day(1), nlev+1, nswbands, fluxes_clrsky_day); // Add an empty level above model top // TODO: combine with day compression above @@ -530,7 +563,7 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, yakl::memset(aer_ssa_bnd_rad, 0); yakl::memset(aer_asm_bnd_rad, 0.); - parallel_for(SimpleBounds<3>(nday, nlev, nswgpts), YAKL_LAMBDA (int iday, int ilev, int igpt) { + parallel_for(SimpleBounds<3>(num_day(1), nlev, nswgpts), YAKL_LAMBDA (int iday, int ilev, int igpt) { cld_tau_gpt_rad(iday,ilev,igpt) = cld_tau_gpt_day(iday,ilev,igpt); cld_ssa_gpt_rad(iday,ilev,igpt) = cld_ssa_gpt_day(iday,ilev,igpt); cld_asm_gpt_rad(iday,ilev,igpt) = cld_asm_gpt_day(iday,ilev,igpt); @@ -542,7 +575,7 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, }); // Do shortwave radiative transfer calculations - radiation.run_shortwave_rrtmgp( ngas, nday, nlev, + radiation.run_shortwave_rrtmgp( ngas, num_day(1), nlev, gas_vmr_rad, pmid_day, tmid_day, pint_day, coszrs_day, albedo_dir_day, albedo_dif_day, cld_tau_gpt_rad, cld_ssa_gpt_rad, cld_asm_gpt_rad, aer_tau_bnd_rad, aer_ssa_bnd_rad, aer_asm_bnd_rad, fluxes_allsky_day.flux_up , fluxes_allsky_day.flux_dn , fluxes_allsky_day.flux_net , fluxes_allsky_day.flux_dn_dir , @@ -566,10 +599,10 @@ void Radiation::radiation_driver_sw(int ncol, real3d& gas_vmr, } void Radiation::radiation_driver_lw(int ncol, int nlev, - real3d& gas_vmr, - real2d& pmid, real2d& pint, real2d& tmid, real2d& tint, - real3d& cld_tau_gpt, real3d& aer_tau_bnd, FluxesByband& fluxes_clrsky, - FluxesByband& fluxes_allsky, real2d& qrl, real2d& qrlc) { + const real3d& gas_vmr, + const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint, + const real3d& cld_tau_gpt, const real3d& aer_tau_bnd, FluxesByband& fluxes_clrsky, + FluxesByband& fluxes_allsky, const real2d& qrl, const real2d& qrlc) { real3d cld_tau_gpt_rad("cld_tau_gpt_rad", ncol, nlev+1, nlwgpts); real3d aer_tau_bnd_rad("aer_tau_bnd-rad", ncol, nlev+1, nlwgpts); @@ -626,30 +659,31 @@ void Radiation::radiation_driver_lw(int ncol, int nlev, }); } -void Radiation::set_daynight_indices(real1d& coszrs, int1d& day_indices, int1d& night_indices) { - // Initialize array of daytime indices to be all zero. If any zeros exist when - // we are done, something went wrong. - yakl::memset(day_indices, 0.); - yakl::memset(night_indices, 0.); - +// Initialize array of daytime indices to be all zero. If any zeros exist when +// we are done, something went wrong. +void Radiation::set_daynight_indices(const real1d& coszrs, const int1d& day_indices, const int1d& night_indices) { // Loop over columns and identify daytime columns as those where the cosine // solar zenith angle exceeds zero. Note that we wrap the setting of // day_indices in an if-then to make sure we are not accesing day_indices out // of bounds, and stopping with an informative error message if we do for some reason. - int iday = 0; - int inight = 0; - for (auto icol = 0; icol < ncol; ++icol) { + int1d iday("iday", 1); + int1d inight("inight",1); + yakl::memset(iday, 0); + yakl::memset(inight, 0); + yakl::memset(day_indices, 0); + yakl::memset(night_indices, 0); + parallel_for(SimpleBounds<1>(ncol), YAKL_LAMBDA (int icol) { if (coszrs(icol) > 0.) { - iday += 1; - day_indices(iday) = icol; + iday(1) += 1; + day_indices(iday(1)) = icol; } else { - inight += 1; - night_indices(inight) = icol; + inight(1) += 1; + night_indices(inight(1)) = icol; } - } + }); } -void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& gas_vmr) { +void Radiation::get_gas_vmr(const std::vector& gas_names, const real3d& gas_vmr) { // Mass mixing ratio real2d mmr("mmr", ncol, nlev); // Gases and molecular weights. Note that we do NOT have CFCs yet (I think @@ -676,12 +710,12 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g if (gas_names[igas] == "CO"){ // CO not available, use default parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - gas_vmr(igas,icol,ilev) = co_vol_mix_ratio; + gas_vmr(igas+1,icol,ilev) = 1.0e-4; //co_vol_mix_ratio; }); } else if (gas_names[igas] == "N2") { // N2 not available, use default parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - gas_vmr(igas,icol,ilev) = n2_vol_mix_ratio; + gas_vmr(igas+1,icol,ilev) = 1.0e-5; //n2_vol_mix_ratio; }); } else if (gas_names[igas] == "H2O") { // Water vapor is represented as specific humidity in CAM, so we @@ -693,8 +727,8 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // first specific humidity (held in the mass_mix_ratio array read // from rad_constituents) is converted to an actual mass mixing ratio. parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - gas_vmr(igas,icol,ilev) = mmr(icol,ilev) / ( - 1. - mmr(icol,ilev))*mol_weight_air / mol_weight_gas[igas]; + gas_vmr(igas+1,icol,ilev) = 1.0e-4; //mmr(icol,ilev) / ( +// 1. - mmr(icol,ilev))*mol_weight_air / mol_weight_gas[igas]; }); } else { // Get mass mixing ratio from the rad_constituents interface @@ -703,8 +737,8 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // Convert to volume mixing ratio by multiplying by the ratio of // molecular weight of dry air to molecular weight of gas parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - gas_vmr(igas,icol,ilev) = mmr(icol,ilev) - * mol_weight_air / mol_weight_gas[igas]; + gas_vmr(igas+1,icol,ilev) = 1.0e-4; //mmr(icol,ilev) +// * mol_weight_air / mol_weight_gas[igas]; }); } } // igas @@ -720,7 +754,7 @@ void Radiation::get_gas_vmr(const std::vector& gas_names, real3d& g // H = dF / dp * g // Why? Something to do with convenience with applying the fluxes to the // heating tendency? -void Radiation::calculate_heating_rate(real2d& flux_up, real2d& flux_dn, real2d& pint, real2d& heating_rate) { +void Radiation::calculate_heating_rate(const real2d& flux_up, const real2d& flux_dn, const real2d& pint, const real2d& heating_rate) { parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { heating_rate(icol,ilev) = (flux_up(icol,ilev+1) - flux_up(icol,ilev)- flux_dn(icol,ilev+1)+flux_dn(icol,ilev)) *CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev)); diff --git a/Source/Radiation/Rrtmgp.H b/Source/Radiation/Rrtmgp.H index eed476913..5b29da669 100644 --- a/Source/Radiation/Rrtmgp.H +++ b/Source/Radiation/Rrtmgp.H @@ -51,25 +51,25 @@ class Rrtmgp { // run rrtmgp short wave model void run_shortwave_rrtmgp (int ngas, int ncol, int nlay, - real3d& gas_vmr, real2d& pmid , real2d& tmid , real2d& pint, - real1d& coszrs , real2d& albedo_dir, real2d& albedo_dif, - real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, - real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, - real2d& allsky_flux_up , real2d& allsky_flux_dn , real2d& allsky_flux_net , real2d& allsky_flux_dn_dir, - real3d& allsky_bnd_flux_up, real3d& allsky_bnd_flux_dn, real3d& allsky_bnd_flux_net, real3d& allsky_bnd_flux_dn_dir, - real2d& clrsky_flux_up , real2d& clrsky_flux_dn , real2d& clrsky_flux_net , real2d& clrsky_flux_dn_dir, - real3d& clrsky_bnd_flux_up, real3d& clrsky_bnd_flux_dn, real3d& clrsky_bnd_flux_net, real3d& clrsky_bnd_flux_dn_dir, + const real3d& gas_vmr, const real2d& pmid, const real2d& tmid, const real2d& pint, + const real1d& coszrs , const real2d& albedo_dir, const real2d& albedo_dif, + const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt, + const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd, + const real2d& allsky_flux_up , const real2d& allsky_flux_dn , const real2d& allsky_flux_net , const real2d& allsky_flux_dn_dir, + const real3d& allsky_bnd_flux_up, const real3d& allsky_bnd_flux_dn, const real3d& allsky_bnd_flux_net, const real3d& allsky_bnd_flux_dn_dir, + const real2d& clrsky_flux_up , const real2d& clrsky_flux_dn , const real2d& clrsky_flux_net , const real2d& clrsky_flux_dn_dir, + const real3d& clrsky_bnd_flux_up, const real3d& clrsky_bnd_flux_dn, const real3d& clrsky_bnd_flux_net, const real3d& clrsky_bnd_flux_dn_dir, double tsi_scaling); // run rrtmgp long wave model - void run_longwave_rrtmgp (int ngas, int ncol, int nlay, real3d& gas_vmr, - real2d& pmid , real2d& tmid , real2d& pint , real2d& tint, - real2d& emis_sfc , - real3d& cld_tau_gpt , real3d& aer_tau_bnd , - real2d& allsky_flux_up , real2d& allsky_flux_dn , real2d& allsky_flux_net , - real3d& allsky_bnd_flux_up, real3d& allsky_bnd_flux_dn, real3d& allsky_bnd_flux_net, - real2d& clrsky_flux_up , real2d& clrsky_flux_dn , real2d& clrsky_flux_net , - real3d& clrsky_bnd_flux_up, real3d& clrsky_bnd_flux_dn, real3d& clrsky_bnd_flux_net); + void run_longwave_rrtmgp (int ngas, int ncol, int nlay, const real3d& gas_vmr, + const real2d& pmid, const real2d& tmid, const real2d& pint , const real2d& tint, + const real2d& emis_sfc, + const real3d& cld_tau_gpt , const real3d& aer_tau_bnd , + const real2d& allsky_flux_up , const real2d& allsky_flux_dn , const real2d& allsky_flux_net , + const real3d& allsky_bnd_flux_up, const real3d& allsky_bnd_flux_dn, const real3d& allsky_bnd_flux_net, + const real2d& clrsky_flux_up , const real2d& clrsky_flux_dn , const real2d& clrsky_flux_net , + const real3d& clrsky_bnd_flux_up, const real3d& clrsky_bnd_flux_dn, const real3d& clrsky_bnd_flux_net); int get_nband_sw() { return k_dist_sw.get_nband(); diff --git a/Source/Radiation/Run_longwave_rrtmgp.cpp b/Source/Radiation/Run_longwave_rrtmgp.cpp index 839e61c09..6a0cd7632 100644 --- a/Source/Radiation/Run_longwave_rrtmgp.cpp +++ b/Source/Radiation/Run_longwave_rrtmgp.cpp @@ -9,14 +9,14 @@ void Rrtmgp::run_longwave_rrtmgp ( int ngas, int ncol, int nlay, - real3d& gas_vmr , - real2d& pmid , real2d& tmid , real2d& pint , real2d& tint, - real2d& emis_sfc , - real3d& cld_tau_gpt , real3d& aer_tau_bnd , - real2d& allsky_flux_up , real2d& allsky_flux_dn , real2d& allsky_flux_net , - real3d& allsky_bnd_flux_up, real3d& allsky_bnd_flux_dn, real3d& allsky_bnd_flux_net, - real2d& clrsky_flux_up , real2d& clrsky_flux_dn , real2d& clrsky_flux_net , - real3d& clrsky_bnd_flux_up, real3d& clrsky_bnd_flux_dn, real3d& clrsky_bnd_flux_net) + const real3d& gas_vmr, + const real2d& pmid, const real2d& tmid, const real2d& pint, const real2d& tint, + const real2d& emis_sfc , + const real3d& cld_tau_gpt , const real3d& aer_tau_bnd , + const real2d& allsky_flux_up , const real2d& allsky_flux_dn , const real2d& allsky_flux_net , + const real3d& allsky_bnd_flux_up, const real3d& allsky_bnd_flux_dn, const real3d& allsky_bnd_flux_net, + const real2d& clrsky_flux_up , const real2d& clrsky_flux_dn , const real2d& clrsky_flux_net , + const real3d& clrsky_bnd_flux_up, const real3d& clrsky_bnd_flux_dn, const real3d& clrsky_bnd_flux_net) { // Wrap pointers in YAKL arrays int nlwbands = k_dist_lw.get_nband(); diff --git a/Source/Radiation/Run_shortwave_rrtmgp.cpp b/Source/Radiation/Run_shortwave_rrtmgp.cpp index e5aff9c8a..f486dc5fb 100644 --- a/Source/Radiation/Run_shortwave_rrtmgp.cpp +++ b/Source/Radiation/Run_shortwave_rrtmgp.cpp @@ -2,14 +2,14 @@ void Rrtmgp::run_shortwave_rrtmgp ( int ngas, int ncol, int nlay, - real3d& gas_vmr, real2d& pmid , real2d& tmid , real2d& pint, - real1d& coszrs , real2d& albedo_dir, real2d& albedo_dif, - real3d& cld_tau_gpt, real3d& cld_ssa_gpt, real3d& cld_asm_gpt, - real3d& aer_tau_bnd, real3d& aer_ssa_bnd, real3d& aer_asm_bnd, - real2d& allsky_flux_up , real2d& allsky_flux_dn , real2d& allsky_flux_net , real2d& allsky_flux_dn_dir, - real3d& allsky_bnd_flux_up, real3d& allsky_bnd_flux_dn, real3d& allsky_bnd_flux_net, real3d& allsky_bnd_flux_dn_dir, - real2d& clrsky_flux_up , real2d& clrsky_flux_dn , real2d& clrsky_flux_net , real2d& clrsky_flux_dn_dir, - real3d& clrsky_bnd_flux_up, real3d& clrsky_bnd_flux_dn, real3d& clrsky_bnd_flux_net, real3d& clrsky_bnd_flux_dn_dir, + const real3d& gas_vmr, const real2d& pmid , const real2d& tmid , const real2d& pint, + const real1d& coszrs , const real2d& albedo_dir, const real2d& albedo_dif, + const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt, + const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd, + const real2d& allsky_flux_up , const real2d& allsky_flux_dn , const real2d& allsky_flux_net , const real2d& allsky_flux_dn_dir, + const real3d& allsky_bnd_flux_up, const real3d& allsky_bnd_flux_dn, const real3d& allsky_bnd_flux_net, const real3d& allsky_bnd_flux_dn_dir, + const real2d& clrsky_flux_up , const real2d& clrsky_flux_dn , const real2d& clrsky_flux_net , const real2d& clrsky_flux_dn_dir, + const real3d& clrsky_bnd_flux_up, const real3d& clrsky_bnd_flux_dn, const real3d& clrsky_bnd_flux_net, const real3d& clrsky_bnd_flux_dn_dir, double tsi_scaling) { diff --git a/Source/Radiation/Slingo.H b/Source/Radiation/Slingo.H index 0c131adb2..dc28f8e59 100644 --- a/Source/Radiation/Slingo.H +++ b/Source/Radiation/Slingo.H @@ -9,41 +9,54 @@ using yakl::fortran::SimpleBounds; class Slingo { public: - static void slingo_liq_optics_sw(int ncol, int nlev, int nswbands, real2d& cldn, real2d& cliqwp, real2d& rel, - real3d& liq_tau, real3d& liq_tau_w, real3d& liq_tau_w_g, real3d& liq_tau_w_f) { + static void slingo_liq_optics_sw(int ncol, int nlev, int nswbands, const real2d& cldn, const real2d& cliqwp, const real2d& rel, + const real3d& liq_tau, const real3d& liq_tau_w, const real3d& liq_tau_w_g, const real3d& liq_tau_w_f) { real1d wavmin("wavmin", nswbands); real1d wavmax("wavmax", nswbands); // Minimum cloud amount (as a fraction of the grid-box area) to // distinguish from clear sky - real cldmin = 1.0e-80; + const real cldmin = 1.0e-80; // Decimal precision of cloud amount (0 -> preserve full resolution; // 10^-n -> preserve n digits of cloud amount) - real cldeps = 0.0; + const real cldeps = 0.0; // A. Slingo's data for cloud particle radiative properties (from 'A GCM // Parameterization for the Shortwave Properties of Water Clouds' JAS // vol. 46 may 1989 pp 1419-1427) - std::vector abarl = // A coefficient for extinction optical depth - { 2.817e-02, 2.682e-02, 2.264e-02, 1.281e-02 }; - std::vector bbarl = // B coefficient for extinction optical depth - { 1.305, 1.346, 1.454, 1.641 }; - std::vector cbarl = // C coefficient for single scat albedo - { -5.62e-08 ,-6.94e-06 ,4.64e-04 ,0.201 }; - std::vector dbarl = // D coefficient for single scat albedo - { 1.63e-07, 2.35e-05 ,1.24e-03 ,7.56e-03 }; - std::vector ebarl = // E coefficient for asymmetry parameter - { 0.829, 0.794, 0.754, 0.826 }; - std::vector fbarl = // F coefficient for asymmetry parameter - { 2.482e-03, 4.226e-03, 6.560e-03, 4.353e-03 }; - - real abarli; // A coefficient for current spectral band - real bbarli; // B coefficient for current spectral band - real cbarli; // C coefficient for current spectral band - real dbarli; // D coefficient for current spectral band - real ebarli; // E coefficient for current spectral band - real fbarli; // F coefficient for current spectral band + real1d abarl("abarl", 4); // A coefficient for extinction optical depth + real1d bbarl("bbarl", 4); // B coefficient for extinction optical depth + real1d cbarl("cbarl", 4); // C coefficient for extinction optical depth + real1d dbarl("dbarl", 4); // D coefficient for extinction optical depth + real1d ebarl("ebarl", 4); // E coefficient for extinction optical depth + real1d fbarl("fbarl", 4); // F coefficient for extinction optical depth + parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int i) { + abarl(1) = 2.817e-02; + abarl(2) = 2.682e-02; + abarl(3) = 2.264e-02; + abarl(4) = 1.281e-02; + bbarl(1) = 1.305; + bbarl(2) = 1.346; + bbarl(3) = 1.454; + bbarl(4) = 1.641; + cbarl(1) = -5.62e-08; + cbarl(2) = -6.94e-06; + cbarl(3) = 4.64e-04; + abarl(4) = 0.201; + dbarl(1) = 1.63e-07; + dbarl(2) = 2.35e-05; + dbarl(3) = 1.24e-03; + dbarl(4) = 7.56e-03; + ebarl(1) = 0.829; + ebarl(2) = 0.794; + ebarl(3) = 0.754; + ebarl(4) = 0.826; + fbarl(1) = 2.482e-03; + fbarl(2) = 4.226e-03; + fbarl(3) = 6.560e-03; + fbarl(4) = 4.353e-03; + }); // Caution... A. Slingo recommends no less than 4.0 micro-meters nor // greater than 20 micro-meters. Here we set effective radius limits @@ -52,7 +65,6 @@ class Slingo { const real rel_max = 16.; int indxsl; - real tmp1l, tmp2l, tmp3l, g; RadConstants::get_sw_spectral_boundaries(wavmin,wavmax,RadConstants::micrometer); @@ -75,41 +87,39 @@ class Slingo { // Set cloud extinction optical depth, single scatter albedo, // asymmetry parameter, and forward scattered fraction: - abarli = abarl[indxsl]; - bbarli = bbarl[indxsl]; - cbarli = cbarl[indxsl]; - dbarli = dbarl[indxsl]; - ebarli = ebarl[indxsl]; - fbarli = fbarl[indxsl]; - - for (auto k=0; k rel > 16 micron (Slingo 89) - if (cldn(i,k) >= cldmin && cldn(i,k) >= cldeps) { - tmp1l = abarli + bbarli/std::min(std::max(rel_min,rel(i,k)),rel_max); - liq_tau(ns,i,k) = 1000.*cliqwp(i,k)*tmp1l; - } else { - liq_tau(ns,i,k) = 0.0; - } - - tmp2l = 1. - cbarli - dbarli*std::min(std::max(rel_min,rel(i,k)),rel_max); - tmp3l = fbarli*std::min(std::max(rel_min,rel(i,k)),rel_max); - // Do not let single scatter albedo be 1. Delta-eddington solution - // for non-conservative case has different analytic form from solution - // for conservative case, and raddedmx is written for non-conservative case. - liq_tau_w(ns,i,k) = liq_tau(ns,i,k) * std::min(tmp2l,.999999); - g = ebarli + tmp3l; - liq_tau_w_g(ns,i,k) = liq_tau_w(ns,i,k) * g; - liq_tau_w_f(ns,i,k) = liq_tau_w(ns,i,k) * g * g; - - } // End do i=1,ncol - } // End do k=1,nlev + parallel_for(SimpleBounds<2>(nlev, ncol), YAKL_LAMBDA (int k, int i) { + auto abarli = abarl(indxsl); + auto bbarli = bbarl(indxsl); + auto cbarli = cbarl(indxsl); + auto dbarli = dbarl(indxsl); + auto ebarli = ebarl(indxsl); + auto fbarli = fbarl(indxsl); + real tmp1l, tmp2l, tmp3l, g; + + // note that optical properties for liquid valid only + // in range of 4.2 > rel > 16 micron (Slingo 89) + if (cldn(i,k) >= cldmin && cldn(i,k) >= cldeps) { + tmp1l = abarli + bbarli/std::min(std::max(rel_min,rel(i,k)),rel_max); + liq_tau(ns,i,k) = 1000.*cliqwp(i,k)*tmp1l; + } else { + liq_tau(ns,i,k) = 0.0; + } + + tmp2l = 1. - cbarli - dbarli*std::min(std::max(rel_min,rel(i,k)),rel_max); + tmp3l = fbarli*std::min(std::max(rel_min,rel(i,k)),rel_max); + // Do not let single scatter albedo be 1. Delta-eddington solution + // for non-conservative case has different analytic form from solution + // for conservative case, and raddedmx is written for non-conservative case. + liq_tau_w(ns,i,k) = liq_tau(ns,i,k) * std::min(tmp2l,.999999); + g = ebarli + tmp3l; + liq_tau_w_g(ns,i,k) = liq_tau_w(ns,i,k) * g; + liq_tau_w_f(ns,i,k) = liq_tau_w(ns,i,k) * g * g; + }); } // nswbands } - static void slingo_liq_optics_lw(int ncol, int nlev, int nlwbands, real2d& cldn, real2d& iclwpth, real2d& iciwpth, real3d& abs_od) { + static void slingo_liq_optics_lw(int ncol, int nlev, int nlwbands, const real2d& cldn, + const real2d& iclwpth, const real2d& iciwpth, const real3d& abs_od) { real2d ficemr("ficemr", ncol, nlev); real2d cwp("cwp", ncol, nlev); real2d cldtau("cldtau", ncol, nlev); From 6655431d089cb60aad2dbac96c9a60be6817d99b Mon Sep 17 00:00:00 2001 From: xyuan Date: Wed, 22 Nov 2023 20:16:01 -0500 Subject: [PATCH 2/2] fix the tabs etc. --- Source/Radiation/Aero_rad_props.cpp | 18 +++++++-------- Source/Radiation/Ebert_curry.H | 34 ++++++++++++++--------------- Source/Radiation/Mam4_aero.H | 4 ++-- Source/Radiation/Optics.H | 10 ++++----- Source/Radiation/Optics.cpp | 6 ++--- Source/Radiation/Radiation.cpp | 8 +++---- Source/Radiation/Slingo.H | 26 +++++++++++----------- 7 files changed, 53 insertions(+), 53 deletions(-) diff --git a/Source/Radiation/Aero_rad_props.cpp b/Source/Radiation/Aero_rad_props.cpp index 69e97ac58..9894d358a 100644 --- a/Source/Radiation/Aero_rad_props.cpp +++ b/Source/Radiation/Aero_rad_props.cpp @@ -11,25 +11,25 @@ using yakl::fortran::SimpleBounds; void AerRadProps::initialize(int ngas_, int nmodes_, int num_aeroes_, int nswbands_, int nlwbands_, - int ncol_, int nlev_, int nrh_, int top_lev_, + int ncol_, int nlev_, int nrh_, int top_lev_, const std::vector& aero_names_, - const real2d& zi_, const real2d& pmid_, const real2d& temp_, + const real2d& zi_, const real2d& pmid_, const real2d& temp_, const real2d& qt_, const real2d& geom_radius) { nmodes = nmodes_; - ngas = ngas_; + ngas = ngas_; num_aeroes = num_aeroes_; aero_names = aero_names_; - nswbands = nswbands_; + nswbands = nswbands_; nlwbands = nlwbands_; - ncol = ncol_; - nlev = nlev_; - nrh = nrh_; + ncol = ncol_; + nlev = nlev_; + nrh = nrh_; top_lev = top_lev_; pmid = pmid_; - pdeldry = pmid_; - temp = temp_; + pdeldry = pmid_; + temp = temp_; qt = qt_; // geometric radius geometric_radius = geom_radius; diff --git a/Source/Radiation/Ebert_curry.H b/Source/Radiation/Ebert_curry.H index 5124a0855..93d3ad238 100644 --- a/Source/Radiation/Ebert_curry.H +++ b/Source/Radiation/Ebert_curry.H @@ -24,28 +24,28 @@ class EbertCurry { parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int i) { abari(1) = 3.448e-03; - abari(2) = 3.448e-03; - abari(3) = 3.448e-03; + abari(2) = 3.448e-03; + abari(3) = 3.448e-03; abari(4) = 3.448e-03; - bbari(1) = 2.431; - bbari(2) = 2.431; - bbari(3) = 2.431; + bbari(1) = 2.431; + bbari(2) = 2.431; + bbari(3) = 2.431; bbari(4) = 2.431; - cbari(1) = 1.00e-05; - cbari(2) = 1.10e-04; + cbari(1) = 1.00e-05; + cbari(2) = 1.10e-04; cbari(3) = 1.861e-02; cbari(4) = 0.46658; - dbari(1) = 0.0; - dbari(2) = 1.405e-05; - dbari(3) = 8.328e-04; + dbari(1) = 0.0; + dbari(2) = 1.405e-05; + dbari(3) = 8.328e-04; dbari(4) = 2.05e-05; - ebari(1) = 0.7661; - ebari(2) = 0.7730; - ebari(3) = 0.794; + ebari(1) = 0.7661; + ebari(2) = 0.7730; + ebari(3) = 0.794; ebari(4) = 0.9595; - fbari(1) = 5.851e-04; - fbari(2) = 5.665e-04; - fbari(3) = 7.267e-04; + fbari(1) = 5.851e-04; + fbari(2) = 5.665e-04; + fbari(3) = 7.267e-04; fbari(4) = 1.076e-04; }); @@ -111,7 +111,7 @@ class EbertCurry { } - static void ec_ice_optics_lw(int ncol, int nlev, int nlwbands, const real2d& cldn, const real2d& iclwpth, + static void ec_ice_optics_lw(int ncol, int nlev, int nlwbands, const real2d& cldn, const real2d& iclwpth, const real2d& iciwpth, const real2d& rei, const real3d& abs_od) { real2d ficemr("ficemr",ncol,nlev); real2d cwp("cwp",ncol,nlev); diff --git a/Source/Radiation/Mam4_aero.H b/Source/Radiation/Mam4_aero.H index c6add80f6..b5706f24b 100644 --- a/Source/Radiation/Mam4_aero.H +++ b/Source/Radiation/Mam4_aero.H @@ -57,7 +57,7 @@ class Mam4_aer { MamConstituents mam_consti; // initialization - inline + inline void initialize(int ncol_, int nlev_, int top_lev_, int nswbands_, int nlwbands_) { ncol = ncol_; nlev = nlev_; @@ -321,7 +321,7 @@ class Mam4_aer { void modal_aero_sw(int list_idx, real dt, int nnite, const int1d& idxnite, const bool& is_cmip6_volc, const real2d& ext_cmip6_sw, const int1d& trop_level, - const real3d& tauxar, const real3d& wa, const real3d& ga, + const real3d& tauxar, const real3d& wa, const real3d& ga, const real3d& fa, const real2d& clear_rh) { int ilev_tropp; diff --git a/Source/Radiation/Optics.H b/Source/Radiation/Optics.H index 51722c433..7c884a316 100644 --- a/Source/Radiation/Optics.H +++ b/Source/Radiation/Optics.H @@ -38,8 +38,8 @@ class Optics { const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out); // get long wave cloud optics property - void get_cloud_optics_lw(int ncol, int nlev, int nbnd, bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, - const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, const real2d& dei, const real2d& des, + void get_cloud_optics_lw(int ncol, int nlev, int nbnd, bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, + const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, const real2d& dei, const real2d& des, const real2d& rei, const real3d& tau_out, const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out); // sample short wave cloud optics property @@ -55,11 +55,11 @@ class Optics { // set the short wave aerosol optics property void set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, real dt, const int1d& night_indices, - bool is_cmip6_volc, const real3d& tau_out, const real3d& ssa_out, const real3d& asm_out, + bool is_cmip6_volc, const real3d& tau_out, const real3d& ssa_out, const real3d& asm_out, const real2d& clear_rh); // set the long wave aerosol optics property - void set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, const real2d& zi, + void set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, const real2d& zi, const real3d& tau, const real2d& clear_rh); // mcica subcol mask @@ -78,7 +78,7 @@ class Optics { const std::vector& aero_names, const real2d& zi, const real2d& pmid, const real2d& temp, const real2d& qi, const real2d& geom_radius); - + // finalize/clean up void finalize(); diff --git a/Source/Radiation/Optics.cpp b/Source/Radiation/Optics.cpp index 85c861f7a..f754b58b7 100644 --- a/Source/Radiation/Optics.cpp +++ b/Source/Radiation/Optics.cpp @@ -175,8 +175,8 @@ void Optics::get_cloud_optics_sw(int ncol, int nlev, int nbnd, //---------------------------------------------------------------------------- void Optics::get_cloud_optics_lw( - int ncol, int nlev, int nbnd, bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, - const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, const real2d& dei, const real2d& des, + int ncol, int nlev, int nbnd, bool do_snow, const real2d& cld, const real2d& cldfsnow, const real2d& iclwp, + const real2d& iciwp, const real2d& icswp, const real2d& lambdac, const real2d& mu, const real2d& dei, const real2d& des, const real2d& rei, const real3d& tau_out, const real3d& liq_tau_out, const real3d& ice_tau_out, const real3d& snw_tau_out) { // Temporary variables to hold absorption optical depth real3d ice_tau("ice_tau", nbnd, ncol, nlev); @@ -393,7 +393,7 @@ void Optics::set_aerosol_optics_sw(int icall, int ncol, int nlev, int nswbands, } //---------------------------------------------------------------------------- -void Optics::set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, const real2d& zi, +void Optics::set_aerosol_optics_lw(int icall, real dt, bool is_cmip6_volc, const real2d& zi, const real3d& tau, const real2d& clear_rh) { // Get aerosol absorption optical depth from CAM routine yakl::memset(tau, 0.); diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index 55ddb4f5b..8767ea9a2 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -163,8 +163,8 @@ void Radiation::initialize(const MultiFab& cons_in, MultiFab& qmoist, auto geom_radius = real2d("geom_radius", ncol, nlev); yakl::memset(geom_radius, 0.1); - optics.initialize(ngas, nmodes, naer, nswbands, nlwbands, - ncol, nlev, nrh, top_lev, aero_names, zi, + optics.initialize(ngas, nmodes, naer, nswbands, nlwbands, + ncol, nlev, nrh, top_lev, aero_names, zi, pmid, tmid, qt, geom_radius); amrex::Print() << " LW coefficents file: \n" @@ -275,7 +275,7 @@ void Radiation::run() { real1d cld_tau_bnd_sw_o_1d("cld_tau_bnd_sw_1d", nswbands); real1d cld_ssa_bnd_sw_o_1d("cld_ssa_bnd_sw_1d", nswbands); real1d cld_asm_bnd_sw_o_1d("cld_asm_bnd_sw_1d", nswbands); - + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay) { for (auto ibnd = 1; ibnd <= nswbands; ++ibnd) { cld_tau_bnd_sw_1d(ibnd) = cld_tau_bnd_sw(icol,ilay,ibnd); @@ -336,7 +336,7 @@ void Radiation::run() { real1d aer_tau_bnd_sw_o_1d("cld_tau_bnd_sw_1d", nswbands); real1d aer_ssa_bnd_sw_o_1d("cld_ssa_bnd_sw_1d", nswbands); real1d aer_asm_bnd_sw_o_1d("cld_asm_bnd_sw_1d", nswbands); - + parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilay) { for (auto ibnd = 1; ibnd < nswbands; ++ibnd) { aer_tau_bnd_sw_1d(ibnd) = aer_tau_bnd_sw(icol,ilay,ibnd); diff --git a/Source/Radiation/Slingo.H b/Source/Radiation/Slingo.H index dc28f8e59..6b2fe8cc7 100644 --- a/Source/Radiation/Slingo.H +++ b/Source/Radiation/Slingo.H @@ -33,29 +33,29 @@ class Slingo { real1d fbarl("fbarl", 4); // F coefficient for extinction optical depth parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int i) { abarl(1) = 2.817e-02; - abarl(2) = 2.682e-02; - abarl(3) = 2.264e-02; + abarl(2) = 2.682e-02; + abarl(3) = 2.264e-02; abarl(4) = 1.281e-02; - bbarl(1) = 1.305; - bbarl(2) = 1.346; - bbarl(3) = 1.454; + bbarl(1) = 1.305; + bbarl(2) = 1.346; + bbarl(3) = 1.454; bbarl(4) = 1.641; cbarl(1) = -5.62e-08; cbarl(2) = -6.94e-06; cbarl(3) = 4.64e-04; abarl(4) = 0.201; - dbarl(1) = 1.63e-07; + dbarl(1) = 1.63e-07; dbarl(2) = 2.35e-05; dbarl(3) = 1.24e-03; dbarl(4) = 7.56e-03; - ebarl(1) = 0.829; - ebarl(2) = 0.794; - ebarl(3) = 0.754; + ebarl(1) = 0.829; + ebarl(2) = 0.794; + ebarl(3) = 0.754; ebarl(4) = 0.826; fbarl(1) = 2.482e-03; - fbarl(2) = 4.226e-03; - fbarl(3) = 6.560e-03; - fbarl(4) = 4.353e-03; + fbarl(2) = 4.226e-03; + fbarl(3) = 6.560e-03; + fbarl(4) = 4.353e-03; }); // Caution... A. Slingo recommends no less than 4.0 micro-meters nor @@ -118,7 +118,7 @@ class Slingo { } // nswbands } - static void slingo_liq_optics_lw(int ncol, int nlev, int nlwbands, const real2d& cldn, + static void slingo_liq_optics_lw(int ncol, int nlev, int nlwbands, const real2d& cldn, const real2d& iclwpth, const real2d& iciwpth, const real3d& abs_od) { real2d ficemr("ficemr", ncol, nlev); real2d cwp("cwp", ncol, nlev);