Add LSM radiation vars and populate them with desired fluxes. (#1604)

Source/ERF.H

@@ -633,13 +633,18 @@ private:
     amrex::Vector<amrex::Real> wind_speed, thrust_coeff, power;
 
 //#endif
+
     LandSurface lsm;
     amrex::Vector<amrex::Vector<amrex::MultiFab*>> lsm_data; // (lev,ncomp) Components: theta, q1, q2
     amrex::Vector<amrex::Vector<amrex::MultiFab*>> lsm_flux; // (lev,ncomp) Components: theta, q1, q2
 
 #if defined(ERF_USE_RRTMGP)
     Radiation rad;
     amrex::Vector<std::unique_ptr<amrex::MultiFab>> qheating_rates;  // radiation heating rate source terms
+
+    // Containers for additional SLM inputs
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> sw_lw_fluxes; // Direct SW (visible, NIR), Diffuse SW (visible, NIR), LW flux
+    amrex::Vector<std::unique_ptr<amrex::MultiFab>> solar_zenith; // Solar zenith angle
 #endif
 
     // Fillpatcher classes for coarse-fine boundaries

Source/ERF.cpp

@@ -122,6 +122,8 @@ ERF::ERF ()
 
 #if defined(ERF_USE_RRTMGP)
     qheating_rates.resize(nlevs_max);
+    sw_lw_fluxes.resize(nlevs_max);
+    solar_zenith.resize(nlevs_max);
 #endif
 
     // NOTE: size lsm before readparams (chooses the model at all levels)
@@ -1790,6 +1792,8 @@ ERF::ERF (const RealBox& rb, int max_level_in,
 
 #if defined(ERF_USE_RRTMGP)
     qheating_rates.resize(nlevs_max);
+    sw_lw_fluxes.resize(nlevs_max);
+    solar_zenith.resize(nlevs_max);
 #endif
 
     // NOTE: size micro before readparams (chooses the model at all levels)

Source/ERF_make_new_arrays.cpp

@@ -246,8 +246,33 @@ ERF::init_stuff (int lev, const BoxArray& ba, const DistributionMapping& dm,
     //*********************************************************
     // Radiation heating source terms
     //*********************************************************
-     qheating_rates[lev] = std::make_unique<MultiFab>(ba, dm, 2, ngrow_state);
-     qheating_rates[lev]->setVal(0.);
+    qheating_rates[lev] = std::make_unique<MultiFab>(ba, dm, 2, ngrow_state);
+    qheating_rates[lev]->setVal(0.);
+
+    //*********************************************************
+    // Radiation fluxes for coupling to LSM
+    //*********************************************************
+
+    // NOTE: Finer levels do not need to coincide with the bottom domain boundary
+    //       at k=0. We make slabs here with the kmin for a given box. Therefore,
+    //       care must be taken before applying these fluxes to an LSM model. For
+
+    // Radiative fluxes for LSM
+    if (solverChoice.lsm_type != LandSurfaceType::None)
+    {
+        BoxList m_bl = ba.boxList();
+        for (auto& b : m_bl) {
+            int kmin = b.smallEnd(2);
+            b.setRange(2,kmin);
+        }
+        BoxArray m_ba(std::move(m_bl));
+
+        sw_lw_fluxes[lev] = std::make_unique<MultiFab>(m_ba, dm, 5, ngrow_state); // SW direct (2), SW diffuse (2), LW
+        solar_zenith[lev] = std::make_unique<MultiFab>(m_ba, dm, 2, ngrow_state);
+
+        sw_lw_fluxes[lev]->setVal(0.);
+        solar_zenith[lev]->setVal(0.);
+    }
 #endif
 }
 

Source/Radiation/Radiation.H

@@ -44,6 +44,8 @@ class Radiation {
 
     // init
     void initialize (const amrex::MultiFab& cons_in,
+                     amrex::MultiFab* lsm_fluxes,
+                     amrex::MultiFab* lsm_zenith,
                      amrex::MultiFab* qheating_rates,
                      amrex::MultiFab* lat,
                      amrex::MultiFab* lon,
@@ -89,11 +91,9 @@ class Radiation {
                                  const real2d& flux_dn,
                                  const real2d& pint,
                                  const real2d& heating_rate);
-#if 0
     void
     export_surface_fluxes(FluxesByband& fluxes,
                           std::string band);
-#endif
 
   private:
     // geometry
@@ -109,6 +109,10 @@ class Radiation {
     amrex::MultiFab* m_lat = nullptr;
     amrex::MultiFab* m_lon = nullptr;
 
+    // Pointer to output data for LSM
+    amrex::MultiFab* m_lsm_fluxes = nullptr;
+    amrex::MultiFab* m_lsm_zenith = nullptr;
+
     // Specified uniform angle for radiation
     amrex::Real uniform_angle = 78.463;
 

Source/Radiation/Radiation.cpp

@@ -99,6 +99,8 @@ namespace internal {
 
 // init
 void Radiation::initialize (const MultiFab& cons_in,
+                            MultiFab* lsm_fluxes,
+                            MultiFab* lsm_zenith,
                             MultiFab* qheating_rates,
                             MultiFab* lat,
                             MultiFab* lon,
@@ -131,6 +133,9 @@ void Radiation::initialize (const MultiFab& cons_in,
     m_lat = lat;
     m_lon = lon;
 
+    m_lsm_fluxes = lsm_fluxes;
+    m_lsm_zenith = lsm_zenith;
+
     rrtmgp_data_path = getRadiationDataDir() + "/";
     rrtmgp_coefficients_file_sw = rrtmgp_data_path + rrtmgp_coefficients_file_name_sw;
     rrtmgp_coefficients_file_lw = rrtmgp_data_path + rrtmgp_coefficients_file_name_lw;
@@ -478,10 +483,8 @@ void Radiation::run ()
                              fluxes_allsky, fluxes_clrsky, qrs, qrsc);
         }
 
-#if 0
         // Set surface fluxes that are used by the land model
         export_surface_fluxes(fluxes_allsky, "shortwave");
-#endif
 
     } else {
         // Conserve energy
@@ -526,10 +529,8 @@ void Radiation::run ()
         radiation_driver_lw(ncol, nlev, gas_vmr, pmid, pint, tmid, tint, cld_tau_gpt_lw, aer_tau_bnd_lw,
                             fluxes_allsky, fluxes_clrsky, qrl, qrlc);
 
-#if 0
         // Set surface fluxes that are used by the land model
         export_surface_fluxes(fluxes_allsky, "longwave");
-#endif
 
     }
     else {
@@ -946,11 +947,13 @@ void Radiation::calculate_heating_rate (const real2d& flux_up,
     });
 }
 
-#if 0
 void
-export_surface_fluxes(FluxesByband& fluxes,
-                      std::string band)
+Radiation::export_surface_fluxes(FluxesByband& fluxes,
+                                 std::string band)
 {
+    // No work to be done if we don't have valid pointers
+    if (!m_lsm_fluxes) return;
+
     if (band == "shortwave") {
         real3d flux_dn_diffuse("flux_dn_diffuse", ncol, nlev+1, nswbands);
 
@@ -963,8 +966,8 @@ export_surface_fluxes(FluxesByband& fluxes,
         });
 
         // Populate the LSM data structure (this is a 2D MF)
-        for (MFIter mfi(*(m_lsm_ptr)); mfi.isValid(); ++mfi) {
-            auto lsm_array = m_lsm_ptr->array(mfi);
+        for (MFIter mfi(*(m_lsm_fluxes)); mfi.isValid(); ++mfi) {
+            auto lsm_array = m_lsm_fluxes->array(mfi);
             const auto& box3d = mfi.tilebox();
             auto nx = box3d.length(0);
             amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
@@ -1005,8 +1008,8 @@ export_surface_fluxes(FluxesByband& fluxes,
         }
     } else if (band == "longwave") {
         // Populate the LSM data structure (this is a 2D MF)
-        for (MFIter mfi(*(m_lsm_ptr)); mfi.isValid(); ++mfi) {
-            auto lsm_array = m_lsm_ptr->array(mfi);
+        for (MFIter mfi(*(m_lsm_fluxes)); mfi.isValid(); ++mfi) {
+            auto lsm_array = m_lsm_fluxes->array(mfi);
             const auto& box3d = mfi.tilebox();
             auto nx = box3d.length(0);
             amrex::ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
@@ -1023,7 +1026,6 @@ export_surface_fluxes(FluxesByband& fluxes,
          amrex::Abort("Unknown radiation band type!");
     }
 }
-#endif
 
 // call back
 void Radiation::on_complete () { }

Source/TimeIntegration/ERF_advance_radiation.cpp

@@ -14,6 +14,8 @@ void ERF::advance_radiation (int lev,
    bool is_cmip6_volcano {false};
 
     rad.initialize(cons,
+                   sw_lw_fluxes[lev].get(),
+                   solar_zenith[lev].get(),
                    qheating_rates[lev].get(),
                    lat_m[lev].get(),
                    lon_m[lev].get(),