Radiation Debug (#1552)

* fix cmake compile for radiation.

* Clean up for debug compile and bump YAKL hash to get NETCDF read fix.

---------

Co-authored-by: Aaron Lattanzi <[email protected]>

Author: AMLattanzi

Source/Radiation/Aero_rad_props.cpp

@@ -309,8 +309,8 @@ void AerRadProps::aer_rad_props_lw (const bool& is_cmip6_volc,
     real2d r_lw_abs("r_lw_abs", ncol, nlev);  // radius dependent mass-specific absorption coefficient
     real1d r_mu("r_mu", ncol);        // log(geometric_mean_radius) domain samples of r_lw_abs(:,:)
     real2d mu("mu", ncol, nlev);          // log(geometric_radius)
-    real r_mu_min, r_mu_max, wmu, mutrunc;
-    int nmu, kmu;
+    //real r_mu_min, r_mu_max, wmu, mutrunc;
+    //int nmu, kmu;
 
     // for table lookup into rh grid
     real2d es("es", ncol, nlev);     // saturation vapor pressure
@@ -327,7 +327,7 @@ void AerRadProps::aer_rad_props_lw (const bool& is_cmip6_volc,
 
     //For cmip6 volcanic file
     int1d trop_level("trop_level", ncol);
-    real lyr_thk;
+    //real lyr_thk;
     real3d ext_cmip6_lw("ext_cmip6_lw", ncol, nlev, nlwbands);
     real3d ext_cmip6_lw_inv_m("ext_cmip6_lw_inv_m",ncol, nlev, nlwbands); //long wave extinction in the units of 1/m
 
@@ -421,7 +421,7 @@ void AerRadProps::aer_rad_props_lw (const bool& is_cmip6_volc,
         {
             int ilev_tropp = trop_level(icol); //tropopause level
             if (ilev < ilev_tropp) {
-                auto lyr_thk = zi(icol,ilev) - zi(icol,ilev+1);
+                //auto lyr_thk = zi(icol,ilev) - zi(icol,ilev+1);
                 // odap_aer(icol,ilev,ilw) = lyr_thk * ext_cmip6_lw_inv_m(icol,ilev,ilw);
             }
         });
@@ -675,7 +675,7 @@ void AerRadProps::get_volcanic_rad_props(const int& ncol,
                                          const real3d& tau_w,
                                          const real3d& tau_w_g,
                                          const real3d& tau_w_f) {
-   int nswbands;
+   //int nswbands;
    parallel_for(SimpleBounds<3>(nswbands, ncol,nlev), YAKL_LAMBDA (int iswband, int icol, int ilev) {
      real g = 0;
      if (scat(iswband) > 0.)

Source/Radiation/Ebert_curry.H

@@ -122,7 +122,7 @@ class EbertCurry {
         real2d cldtau("cldtau",ncol,nlev);
 
         // longwave liquid absorption coeff (m**2/g)
-        const real kabsl = 0.090361;
+        //const real kabsl = 0.090361;
 
         // Optical properties for ice are valid only in the range of
         // 13 < rei < 130 micron (Ebert and Curry 92)

Source/Radiation/Mam4_aero.H

@@ -326,7 +326,6 @@ class Mam4_aer {
                          const real3d& fa, const real2d& clear_rh)
     {
         // local variables
-        int ilev_tropp;
         real2d mass("mass",ncol,nlev);               // layer mass
         real2d air_density("air_density",ncol,nlev);  // (kg/m3)
 
@@ -656,7 +655,7 @@ class Mam4_aer {
                 }
             }
 
-            int nlwbands; // add by xyuan to be compiled
+            //int nlwbands; // add by xyuan to be compiled
             for(auto ilw = 1; ilw <= nlwbands; ++ilw) {
                 for(auto k = top_lev; k > nlev; --k) {
                     // form bulk refractive index. Use volume mixing for infrared

Source/Radiation/Mam4_constituents.H

@@ -279,8 +279,8 @@ class MamConstituents {
     inline
     void rad_cnst_init ()
     {
-        int num_aerosols;
-        constexpr bool stricttest = true;
+        //int num_aerosols;
+        //constexpr bool stricttest = true;
 
         // memory to point to if zero value requested
         //allocate(zero_cols(pcols,pver))
@@ -567,7 +567,7 @@ class MamConstituents {
     const
     void rad_cnst_out (int list_idx)
     {
-        int ncol;
+        //int ncol;
         int idx;
         std::string name, cbname, source;
 
@@ -645,7 +645,7 @@ class MamConstituents {
 
         // Get data source
         auto source = aerlist.aer[aer_idx].source;
-        auto idx    = aerlist.aer[aer_idx].idx;
+        //auto idx    = aerlist.aer[aer_idx].idx;
     }
 
     // Return pointer to mass mixing ratio for the modal aerosol specie from the specified
@@ -654,7 +654,7 @@ class MamConstituents {
     void rad_cnst_get_mam_mmr_by_idx (int list_idx, int mode_idx, int spec_idx,
                                       const std::string& phase, real2d& mmr)
     {
-        int idx;
+        //int idx;
         std::string source;
         modelist_t mlist;
 
@@ -729,7 +729,7 @@ class MamConstituents {
     {
         modelist_t mlist;
         std::string source;
-        int idx;
+        //int idx;
 
         if (list_idx >= 0 && list_idx <= N_DIAG) {
             mlist = ma_list[list_idx];

Source/Radiation/Modal_aero_wateruptake.H

@@ -47,10 +47,7 @@ class ModalAeroWateruptake {
 
         real1d es("es",ncol);             // saturation vapor pressure
         real1d qs("qs",ncol);             // saturation specific humidity
-        real cldn_thresh;
         real2d aerosol_water("aerosol_water",ncol,nlev); //sum of aerosol water (wat_a1 + wat_a2 + wat_a3 + wat_a4)
-        bool history_aerosol;      // Output the MAM aerosol variables and tendencies
-        bool history_verbose;      // produce verbose history output
         bool compute_wetdens;
 
         std::string trnum;       // used to hold mode number (as characters)
@@ -276,7 +273,6 @@ class ModalAeroWateruptake {
     {
         const real eps = 1.e-4;
         const real mw = 18.;
-        const real pi = 3.14159;
         const real rhow = 1.;
         const real surften = 76.;
         const real tair = 273.;

Source/Radiation/Phys_prop.H

@@ -485,7 +485,7 @@ class PhysProp {
         real2d flw_abs;
 
         auto nrh       = RadConstants::nrh;
-        auto nbnd      = prop.getDimSize( "lw_band" );
+        //auto nbnd      = prop.getDimSize( "lw_band" );
         auto nswbands  = prop.getDimSize( "sw_band" );
 
         prop.read( fsw_ext, "ext_sw");
@@ -578,9 +578,9 @@ class PhysProp {
     // Read optics data of type 'volcanic_radius'
     void volcanic_radius_optics_init (physprop_t& phys_prop, yakl::SimpleNetCDF& prop)
     {
-        auto n_mu_samples = prop.getDimSize( "mu_samples" );
-        auto nbnd         = prop.getDimSize( "lw_band" );
-        auto nswbands     = prop.getDimSize( "sw_band" );
+        //auto n_mu_samples = prop.getDimSize( "mu_samples" );
+        //auto nbnd         = prop.getDimSize( "lw_band" );
+        //auto nswbands     = prop.getDimSize( "sw_band" );
 
         prop.read( phys_prop.r_sw_ext,   "bext_sw");
         prop.read( phys_prop.r_sw_scat,  "bsca_sw");
@@ -595,8 +595,8 @@ class PhysProp {
     // Read optics data of type 'volcanic'
     void volcanic_optics_init (physprop_t& phys_prop, yakl::SimpleNetCDF& prop)
     {
-        auto nbnd     = prop.getDimSize( "lw_band" );
-        auto nswbands = prop.getDimSize( "sw_band" );
+        //auto nbnd     = prop.getDimSize( "lw_band" );
+        //auto nswbands = prop.getDimSize( "sw_band" );
 
         prop.read( phys_prop.sw_nonhygro_ext,   "bext_sw");
         prop.read( phys_prop.sw_nonhygro_scat,  "bsca_sw");
@@ -611,7 +611,7 @@ class PhysProp {
     void hygroscopic_optics_init (physprop_t& phys_prop, yakl::SimpleNetCDF& prop)
     {
         // temp data from hygroscopic file before interpolation onto cam-rh-mesh
-        int nfilerh; // number of rh values in file
+        //int nfilerh; // number of rh values in file
         real1d frh;
         real2d fsw_ext, fsw_ssa, fsw_asm, flw_abs;
 
@@ -631,7 +631,7 @@ class PhysProp {
         prop.read( frh,     "rh");
 
         real1d fswe("",nrh), fsws("",nrh),
-            fswa("",nrh), flwa("",nrh);
+               fswa("",nrh), flwa("",nrh);
 
         // interpolate onto cam's rh mesh
         parallel_for (SimpleBounds<2> (nswbands, nrh), YAKL_LAMBDA (int kbnd, int krh)
@@ -671,8 +671,8 @@ class PhysProp {
     // Read optics data of type 'nonhygro'
     void nonhygro_optics_init (physprop_t& phys_prop, yakl::SimpleNetCDF& prop)
     {
-        auto nlwbands = prop.getDimSize( "lw_band" );
-        auto nswbands = prop.getDimSize( "sw_band" );
+        //auto nlwbands = prop.getDimSize( "lw_band" );
+        //auto nswbands = prop.getDimSize( "sw_band" );
 
         prop.read( phys_prop.sw_nonhygro_ext, "ext_sw");
         prop.read( phys_prop.sw_nonhygro_ssa, "ssa_sw");
@@ -833,9 +833,9 @@ class PhysProp {
     void aer_optics_log_rh (std::string name, const real1d& ext, const real1d& ssa, const real1d& asmin)
     {
         const int nrh_test = 36;
-        int krh;
+        //int krh;
         real1d rh_test("rh_test", nrh_test);
-        auto nrh = ext.extent(0);
+        //auto nrh = ext.extent(0);
 
         parallel_for (SimpleBounds<1> (nrh_test), YAKL_LAMBDA (int krh_test)
         {
@@ -845,13 +845,15 @@ class PhysProp {
         // loop through test rh values
         parallel_for (SimpleBounds<1> (nrh_test), YAKL_LAMBDA (int krh_test)
         {
+            /*
             // find corresponding rh index
             auto rh = rh_test(krh_test);
             auto krh = std::min(floor( (rh) * nrh ) + 1, static_cast<real>(nrh - 1));
             auto wrh = (rh) *nrh - krh;
             auto exti = ext(krh + 1) * (wrh + 1) - ext(krh) * wrh;
             auto ssai = ssa(krh + 1) * (wrh + 1) - ssa(krh) * wrh;
             auto asmi = asmin(krh + 1) * (wrh + 1) - asmin(krh) * wrh;
+            */
         });
     }
 };

Source/Radiation/Rad_constants.H

@@ -188,7 +188,7 @@ class RadConstants {
     {
         real1d lower_boundaries("lower_boundaries", nlwbands);
         real1d upper_boundaries("upper_boundaries", nlwbands);
-        int iband;
+        //int iband;
 
         // Get band boundaries
         get_lw_spectral_boundaries(lower_boundaries, upper_boundaries, units);

Source/Radiation/Radiation.cpp

@@ -172,7 +172,7 @@ void Radiation::initialize (const MultiFab& cons_in,
 
         const auto& box3d = mfi.tilebox();
         auto nx = box3d.length(0);
-        auto ny = box3d.length(1);
+        //auto ny = box3d.length(1);
         // Get pressure, theta, temperature, density, and qt, qp
         amrex::ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
@@ -283,7 +283,7 @@ void Radiation::run ()
     real3d gas_vmr("gas_vmr", ngas, ncol, nlev);
 
     // Needed for shortwave aerosol;
-    int nday, nnight;     // Number of daylight columns
+    //int nday, nnight;     // Number of daylight columns
     int1d day_indices("day_indices", ncol), night_indices("night_indices", ncol);   // Indicies of daylight coumns
 
     // Flag to carry (QRS,QRL)*dp across time steps.
@@ -783,7 +783,7 @@ void Radiation::get_gas_vmr (const std::vector<std::string>& gas_names, const re
     const std::vector<real> mol_weight_gas = {18.01528, 44.0095, 47.9982, 44.0128,
                                               28.0101, 16.04246, 31.998, 28.0134}; // g/mol
     // Molar weight of air
-    const real mol_weight_air = 28.97; // g/mol
+    //const real mol_weight_air = 28.97; // g/mol
     // Defaults for gases that are not available (TODO: is this still accurate?)
     const real co_vol_mix_ratio = 1.0e-7;
     const real n2_vol_mix_ratio = 0.7906;

Source/Radiation/Run_shortwave_rrtmgp.cpp

@@ -35,7 +35,7 @@ void Rrtmgp::run_shortwave_rrtmgp (int ngas, int ncol, int nlay,
     boolHost1d top_at_1_h("top_at_1_h",1);
     bool top_at_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);
+    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
 
     // Apply TOA flux scaling
     parallel_for(SimpleBounds<2>(nswgpts,ncol), YAKL_LAMBDA (int igpt, int icol) {

Source/Radiation/m2005_effradius.H

@@ -98,7 +98,7 @@ void m2005_effradius (const real2d& ql, const real2d& nl, const real2d& qi,
 
             // set the minimum droplet number as 20/cm3.
             // nlic   = max(nlic,20.e6_r8/rho) ! sghan minimum in #/cm3
-            auto tempnc = nlic/rho/1.0e6;    // #/kg --> #/cm3
+            //auto tempnc = nlic/rho/1.0e6;    // #/kg --> #/cm3
 
             // Should be the in-cloud dropelt number calculated as nlic*rho/1.0e6_r8 ????!!!! +++mhwang
             auto pgam   = 0.0005714*(nlic*rho/1.e6) + 0.2714;

Source/TimeIntegration/ERF_slow_rhs_post.cpp

@@ -117,12 +117,6 @@ void erf_slow_rhs_post (int level, int finest_level,
                                     tc.pbl_type == PBLType::MYNN25      ||
                                     tc.pbl_type == PBLType::YSU );
 
-    // Open bc will be imposed upon all vars (we only access cons here for simplicity)
-    const bool xlo_open = (bc_ptr_h[BCVars::cons_bc].lo(0) == ERFBCType::open);
-    const bool xhi_open = (bc_ptr_h[BCVars::cons_bc].hi(0) == ERFBCType::open);
-    const bool ylo_open = (bc_ptr_h[BCVars::cons_bc].lo(1) == ERFBCType::open);
-    const bool yhi_open = (bc_ptr_h[BCVars::cons_bc].hi(1) == ERFBCType::open);
-
     const Box& domain = geom.Domain();
 
     const GpuArray<Real, AMREX_SPACEDIM> dxInv = geom.InvCellSizeArray();