From 94815bba22f635aaa8c74adc2c0af08ac900543c Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 23 Apr 2024 16:28:28 -0700 Subject: [PATCH] More Radiation (#1587) * Fix heating rate calculation and add qv to the EOS calls. * Fix units and pressure/temp extrap. LW fluxes look weird still. --- Source/ERF_make_new_arrays.cpp | 2 +- Source/Radiation/Radiation.cpp | 68 +++++++++++++---------- Source/Radiation/Run_shortwave_rrtmgp.cpp | 60 ++++++++++---------- Source/SourceTerms/ERF_make_sources.cpp | 3 +- 4 files changed, 74 insertions(+), 59 deletions(-) diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 1f2036bd7..3051e3d48 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -246,7 +246,7 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm, //********************************************************* // Radiation heating source terms //********************************************************* - qheating_rates[lev] = std::make_unique(ba, dm, 1, ngrow_state); + qheating_rates[lev] = std::make_unique(ba, dm, 2, ngrow_state); qheating_rates[lev]->setVal(0.); #endif } diff --git a/Source/Radiation/Radiation.cpp b/Source/Radiation/Radiation.cpp index a80718cb7..61fba0ed7 100644 --- a/Source/Radiation/Radiation.cpp +++ b/Source/Radiation/Radiation.cpp @@ -177,6 +177,7 @@ void Radiation::initialize (const MultiFab& cons_in, for ( MFIter mfi(cons_in, false); mfi.isValid(); ++mfi) { auto states_array = cons_in.array(mfi); auto qt_array = (qmoist[0]) ? qmoist[0]->array(mfi) : Array4 {}; + auto qv_array = (qmoist[1]) ? qmoist[1]->array(mfi) : Array4 {}; auto qc_array = (qmoist[2]) ? qmoist[2]->array(mfi) : Array4 {}; auto qi_array = (qmoist.size()>=8) ? qmoist[3]->array(mfi) : Array4 {}; @@ -188,26 +189,32 @@ void Radiation::initialize (const MultiFab& cons_in, { auto icol = j*nx+i+1; auto ilev = k+1; + Real qv = (qv_array) ? qv_array(i,j,k): 0.0; qt(icol,ilev) = (qt_array) ? qt_array(i,j,k): 0.0; qc(icol,ilev) = (qc_array) ? qc_array(i,j,k): 0.0; qi(icol,ilev) = (qi_array) ? qi_array(i,j,k): 0.0; qn(icol,ilev) = qc(icol,ilev) + qi(icol,ilev); - tmid(icol,ilev) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp)); - pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp))/1.0e3; + tmid(icol,ilev) = getTgivenRandRTh(states_array(i,j,k,Rho_comp),states_array(i,j,k,RhoTheta_comp),qv); + // NOTE: RRTMGP code expects pressure in mb so we convert it here + pmid(icol,ilev) = getPgivenRTh(states_array(i,j,k,RhoTheta_comp),qv) * 1.0e-2; }); } parallel_for(SimpleBounds<2>(ncol, nlev+1), YAKL_LAMBDA (int icol, int ilev) { if (ilev == 1) { - pint(icol, 1) = 2.*pmid(icol, 2) - pmid(icol, 1); - tint(icol, 1) = 2.*tmid(icol, 2) - tmid(icol, 1); + //pint(icol, 1) = 2.*pmid(icol, 2) - pmid(icol, 1); + //tint(icol, 1) = 2.*tmid(icol, 2) - tmid(icol, 1); + pint(icol, 1) = -0.5*pmid(icol, 2) + 1.5*pmid(icol, 1); + tint(icol, 1) = -0.5*tmid(icol, 2) + 1.5*tmid(icol, 1); } else if (ilev <= nlev) { pint(icol, ilev) = 0.5*(pmid(icol, ilev-1) + pmid(icol, ilev)); tint(icol, ilev) = 0.5*(tmid(icol, ilev-1) + tmid(icol, ilev)); } else { - pint(icol, nlev+1) = 2.*pmid(icol, nlev-1) - pmid(icol, nlev); - tint(icol, nlev+1) = 2.*tmid(icol, nlev-1) - tmid(icol, nlev); + //pint(icol, nlev+1) = 2.*pmid(icol, nlev-1) - pmid(icol, nlev); + //tint(icol, nlev+1) = 2.*tmid(icol, nlev-1) - tmid(icol, nlev); + pint(icol, nlev+1) = -0.5*pmid(icol, nlev-1) + 1.5*pmid(icol, nlev); + tint(icol, nlev+1) = -0.5*tmid(icol, nlev-1) + 1.5*tmid(icol, nlev); } }); @@ -252,9 +259,6 @@ void Radiation::initialize (const MultiFab& cons_in, // run radiation model void Radiation::run () { - // Temporary variable for heating rate output - real2d hr("hr", ncol, nlev); - // Cosine solar zenith angle for all columns in chunk real1d coszrs("coszrs", ncol); @@ -512,32 +516,20 @@ void Radiation::run () } } // dolw - // Compute heating rate for dtheta/dt - parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) - { - // TODO: CHECK THIS UNIT CONVERSION!! - hr(icol,ilev) = (qrs(icol,ilev) + qrl(icol,ilev)) / Cp_d * (1.e5 / std::pow(pmid(icol,ilev), R_d/Cp_d)); - }); - - // Populate theta source from hr + // Populate source term for theta dycore variable for (MFIter mfi(*(qrad_src)); mfi.isValid(); ++mfi) { auto qrad_src_array = qrad_src->array(mfi); const auto& box3d = mfi.tilebox(); auto nx = box3d.length(0); amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) { + // Map (col,lev) to (i,j,k) auto icol = j*nx+i+1; auto ilev = k+1; - qrad_src_array(i,j,k) = hr(icol,ilev); - }); - } - // convert radiative heating rates to Q*dp for energy conservation - if (conserve_energy) { - parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) - { - qrs(icol,ilev) = qrs(icol,ilev)*pdel(icol,ilev); - qrl(icol,ilev) = qrl(icol,ilev)*pdel(icol,ilev); + // SW and LW sources + qrad_src_array(i,j,k,0) = qrs(icol,ilev); + qrad_src_array(i,j,k,1) = qrl(icol,ilev); }); } } @@ -894,11 +886,29 @@ void Radiation::calculate_heating_rate (const real2d& flux_up, const real2d& pint, const real2d& heating_rate) { + // NOTE: The pressure is in [mb] for RRTMGP to use. + // The fluxes are in [W/m^2] and gravity is [m/s^2]. + // We need to convert pressure from [mb] -> [Pa] and divide by Cp [J/kg*K] + // The heating rate is {dF/dP * g / Cp} with units [K/s] + real1d heatfac("heatfac",1); + yakl::memset(heatfac, 1.0e-2/Cp_d); parallel_for(SimpleBounds<2>(ncol, nlev), YAKL_LAMBDA (int icol, int ilev) { - heating_rate(icol,ilev) = ( (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) - - (flux_up(icol,ilev ) - flux_dn(icol,ilev )) ) - * CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev)); + heating_rate(icol,ilev) = heatfac(1) * ( (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) + - (flux_up(icol,ilev ) - flux_dn(icol,ilev )) ) + * CONST_GRAV/(pint(icol,ilev+1)-pint(icol,ilev)); + /* + if (icol==1) { + amrex::Print() << "HR: " << ilev << ' ' + << heating_rate(icol,ilev) << ' ' + << (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) << ' ' + << (flux_up(icol,ilev ) - flux_dn(icol,ilev )) << ' ' + << (flux_up(icol,ilev+1) - flux_dn(icol,ilev+1)) + - (flux_up(icol,ilev ) - flux_dn(icol,ilev )) << ' ' + << (pint(icol,ilev+1)-pint(icol,ilev)) << ' ' + << CONST_GRAV << "\n"; + } + */ }); } diff --git a/Source/Radiation/Run_shortwave_rrtmgp.cpp b/Source/Radiation/Run_shortwave_rrtmgp.cpp index 9d7fcaa6e..086a90688 100644 --- a/Source/Radiation/Run_shortwave_rrtmgp.cpp +++ b/Source/Radiation/Run_shortwave_rrtmgp.cpp @@ -34,16 +34,18 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay, bool1d top_at_1_g("top_at_1_g",1); boolHost1d top_at_1_h("top_at_1_h",1); bool top_at_1; + parallel_for(SimpleBounds<1>(1), YAKL_LAMBDA (int icol) { // HACK: Single loop kernel is not efficient + top_at_1_g(1) = pmid(1, 1) < pmid (1, 2); + }); + top_at_1_g.deep_copy_to(top_at_1_h); + top_at_1 = top_at_1_h(1); real2d toa_flux("toa_flux", ncol, nswgpts); - k_dist_sw.gas_optics(ncol, nlay, &top_at_1, pmid, pint, tmid, gas_concs, combined_optics, toa_flux); // TODO: top_at_1 is not a valid address and this doesn't match longwave call + k_dist_sw.gas_optics(ncol, nlay, top_at_1, pmid, pint, tmid, gas_concs, combined_optics, toa_flux); // Apply TOA flux scaling parallel_for(SimpleBounds<2>(nswgpts,ncol), YAKL_LAMBDA (int igpt, int icol) { toa_flux(icol, igpt) = tsi_scaling * toa_flux(icol,igpt); - top_at_1_g(1) = pmid(1, 1) < pmid (1, 2); }); - top_at_1_g.deep_copy_to(top_at_1_h); - top_at_1 = top_at_1_h(1); // Add in aerosol // TODO: should we avoid allocating here? @@ -74,25 +76,26 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay, // Do the clearsky calculation before adding in clouds FluxesByband fluxes_clrsky; - fluxes_clrsky.flux_up = real2d("clrsky_flux_up", ncol, nlay+1); // clrsky_flux_up; - fluxes_clrsky.flux_dn = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn; - fluxes_clrsky.flux_dn_dir = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_dn_dir; - fluxes_clrsky.flux_net = real2d("clrsky_flux_up", ncol, nlay+1); //clrsky_flux_net; - fluxes_clrsky.bnd_flux_up = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_up; - fluxes_clrsky.bnd_flux_dn = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn; - fluxes_clrsky.bnd_flux_dn_dir = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn_dir; - fluxes_clrsky.bnd_flux_net = real3d("clrsky_flux_up", ncol, nlay+1, nswbands); //clrsky_bnd_flux_net; + fluxes_clrsky.flux_up = real2d("clrsky_flux_up" , ncol, nlay+1); // clrsky_flux_up; + fluxes_clrsky.flux_dn = real2d("clrsky_flux_nd" , ncol, nlay+1); //clrsky_flux_dn; + fluxes_clrsky.flux_net = real2d("clrsky_flux_net", ncol, nlay+1); //clrsky_flux_net; + fluxes_clrsky.flux_dn_dir = real2d("clrsky_flux_dn_dir", ncol, nlay+1); //clrsky_flux_dn_dir; + fluxes_clrsky.bnd_flux_up = real3d("clrsky_bnd_flux_up" , ncol, nlay+1, nswbands); //clrsky_bnd_flux_up; + fluxes_clrsky.bnd_flux_dn = real3d("clrsky_bnd_flux_dn" , ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn; + fluxes_clrsky.bnd_flux_net = real3d("clrsky_bnd_flux_net", ncol, nlay+1, nswbands); //clrsky_bnd_flux_net; + fluxes_clrsky.bnd_flux_dn_dir = real3d("clrsky_bnd_flux_dn_dir", ncol, nlay+1, nswbands); //clrsky_bnd_flux_dn_dir; + rte_sw(combined_optics, top_at_1, coszrs, toa_flux, albedo_dir, albedo_dif, fluxes_clrsky); // Copy fluxes back out of FluxesByband object - fluxes_clrsky.flux_up.deep_copy_to(clrsky_flux_up); - fluxes_clrsky.flux_dn.deep_copy_to(clrsky_flux_dn); - fluxes_clrsky.flux_dn_dir.deep_copy_to(clrsky_flux_dn_dir); + fluxes_clrsky.flux_up.deep_copy_to (clrsky_flux_up); + fluxes_clrsky.flux_dn.deep_copy_to (clrsky_flux_dn); fluxes_clrsky.flux_net.deep_copy_to(clrsky_flux_net); - fluxes_clrsky.bnd_flux_up.deep_copy_to(clrsky_bnd_flux_up); - fluxes_clrsky.bnd_flux_dn.deep_copy_to(clrsky_bnd_flux_dn); - fluxes_clrsky.bnd_flux_dn_dir.deep_copy_to(clrsky_bnd_flux_dn_dir); + fluxes_clrsky.flux_dn_dir.deep_copy_to(clrsky_flux_dn_dir); + fluxes_clrsky.bnd_flux_up.deep_copy_to (clrsky_bnd_flux_up); + fluxes_clrsky.bnd_flux_dn.deep_copy_to (clrsky_bnd_flux_dn); fluxes_clrsky.bnd_flux_net.deep_copy_to(clrsky_bnd_flux_net); + fluxes_clrsky.bnd_flux_dn_dir.deep_copy_to(clrsky_bnd_flux_dn_dir); // Add in clouds OpticalProps2str cloud_optics; @@ -110,25 +113,26 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay, // Call SW flux driver FluxesByband fluxes_allsky; - fluxes_allsky.flux_up = real2d("allsky_flux_up", ncol, nlay+1); // allsky_flux_up; - fluxes_allsky.flux_dn = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_dn; - fluxes_allsky.flux_dn_dir = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_dn_dir; - fluxes_allsky.flux_net = real2d("allsky_flux_up", ncol, nlay+1); //allsky_flux_net; - fluxes_allsky.bnd_flux_up = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_up; - fluxes_allsky.bnd_flux_dn = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn; - fluxes_allsky.bnd_flux_dn_dir = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn_dir; - fluxes_allsky.bnd_flux_net = real3d("allsky_flux_up", ncol, nlay+1, nswbands); //allsky_bnd_flux_net; + fluxes_allsky.flux_up = real2d("allsky_flux_up" , ncol, nlay+1); //allsky_flux_up; + fluxes_allsky.flux_dn = real2d("allsky_flux_dn" , ncol, nlay+1); //allsky_flux_dn; + fluxes_allsky.flux_net = real2d("allsky_flux_net", ncol, nlay+1); //allsky_flux_net; + fluxes_allsky.flux_dn_dir = real2d("allsky_flux_dn_dir", ncol, nlay+1); //allsky_flux_dn_dir; + fluxes_allsky.bnd_flux_up = real3d("allsky_bnd_flux_up" , ncol, nlay+1, nswbands); //allsky_bnd_flux_up; + fluxes_allsky.bnd_flux_dn = real3d("allsky_bnd_flux_dn" , ncol, nlay+1, nswbands); //allsky_bnd_flux_dn; + fluxes_allsky.bnd_flux_net = real3d("allsky_bnd_flux_net", ncol, nlay+1, nswbands); //allsky_bnd_flux_net; + fluxes_allsky.bnd_flux_dn_dir = real3d("allsky_bnd_flux_dn_dir", ncol, nlay+1, nswbands); //allsky_bnd_flux_dn_dir; + rte_sw(combined_optics, top_at_1, coszrs, toa_flux, albedo_dir, albedo_dif, fluxes_allsky); // Copy fluxes back out of FluxesByband object fluxes_allsky.flux_up.deep_copy_to(allsky_flux_up); fluxes_allsky.flux_dn.deep_copy_to(allsky_flux_dn); - fluxes_allsky.flux_dn_dir.deep_copy_to(allsky_flux_dn_dir); fluxes_allsky.flux_net.deep_copy_to(allsky_flux_net); + fluxes_allsky.flux_dn_dir.deep_copy_to(allsky_flux_dn_dir); fluxes_allsky.bnd_flux_up.deep_copy_to(allsky_bnd_flux_up); fluxes_allsky.bnd_flux_dn.deep_copy_to(allsky_bnd_flux_dn); - fluxes_allsky.bnd_flux_dn_dir.deep_copy_to(allsky_bnd_flux_dn_dir); fluxes_allsky.bnd_flux_net.deep_copy_to(allsky_bnd_flux_net); + fluxes_allsky.bnd_flux_dn_dir.deep_copy_to(allsky_bnd_flux_dn_dir); yakl::fence(); } diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index 2f0887d62..32102014e 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -153,7 +153,8 @@ void make_sources (int level, auto const& qheating_arr = qheating_rates->const_array(mfi); ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cell_src(i,j,k,RhoTheta_comp) += qheating_arr(i,j,k); + // Short-wavelength and long-wavelength radiation source terms + cell_src(i,j,k,RhoTheta_comp) += qheating_arr(i,j,k,0) + qheating_arr(i,j,k,1); }); }