Use Native SHOC cloud fraction with RRTMGP

Docs/sphinx_doc/Inputs.rst

@@ -2036,6 +2036,9 @@ List of Parameters
 +-------------------------------------+----------------------------------------+-------------------+-----------------------------------+
 | **erf.rad_do_subcol_sampling**      | Enable MCICA subcolumn sampling        | true / false      | true                              |
 +-------------------------------------+----------------------------------------+-------------------+-----------------------------------+
+| **erf.rad_use_shoc_cldfrac**        | Use Native SHOC diagnosed liquid-cloud | true / false      | true                              |
+|                                     | fraction in RRTMGP when available      |                   |                                   |
++-------------------------------------+----------------------------------------+-------------------+-----------------------------------+
 | **erf.rad_orbital_year**            | Fixed orbital year for zenith calcs    | Integer           | < 0 uses timestamp year           |
 +-------------------------------------+----------------------------------------+-------------------+-----------------------------------+
 | **erf.rad_orbital_eccentricity**    | Override orbital eccentricity          | Real              | < 0 uses computed value           |

Source/ERF.H

@@ -462,7 +462,8 @@ public:
 
     void advance_radiation (int lev,
                             amrex::MultiFab& cons_in,
-                            const double& dt_advance);
+                            const double& dt_advance,
+                            const amrex::MultiFab* liquid_cloud_fraction = nullptr);
 
 #ifdef ERF_USE_EAMXX_SHOC
     void compute_shoc_tendencies (int lev,

Source/PBL/Shoc/README.md

@@ -456,6 +456,16 @@ Also run representative non-SHOC regression tests after changing shared coupling
 microphysics, plotfile, or time-integration paths. Native SHOC must not alter
 non-SHOC results when it is not selected at runtime.
 
+## RRTMGP radiation coupling
+
+When Native SHOC is coupled to RRTMGP, ERF passes SHOC's diagnosed liquid-cloud
+fraction to radiation when `erf.rad_use_shoc_cldfrac = true` (the default).
+Radiation is called after the current Native SHOC update so the diagnostic and
+the host thermodynamic state are from the same step. The diagnostic describes
+liquid-cloud macrophysics only; cloud ice retains a binary fraction. ERF's host
+`qc` field remains the liquid input to RRTMGP, and `shoc_ql` is not substituted
+for it by this coupling.
+
 ## See also
 
 User documentation:

Source/PhysicsInterfaces/Radiation/ERF_Radiation.H

@@ -39,6 +39,7 @@
 #include <ERF_IndexDefines.H>
 #include <ERF_DataStruct.H>
 #include <ERF_EOS.H>
+#include <ERF_RadiationCloudFraction.H>
 
 #include <ERF_LandSurface.H>
 
@@ -100,6 +101,7 @@ public:
          const amrex::BoxArray& ba,
          amrex::Geometry& geom,
          amrex::MultiFab* cons_in,
+         const amrex::MultiFab* liquid_cloud_fraction,
          amrex::iMultiFab* lmask,
          amrex::MultiFab*  t_surf,
          amrex::Vector<amrex::MultiFab*>& lsm_input_ptrs,
@@ -112,7 +114,7 @@ public:
          const bool updated_lsm) override
     {
         set_grids(level, step, time, dt, ba, geom,
-                  cons_in, lmask, t_surf,
+                  cons_in, liquid_cloud_fraction, lmask, t_surf,
                   lsm_input_ptrs, qheating_rates,
                   rad_fluxes, z_phys, lat_ptr, lon_ptr,
                   updated_lsm);
@@ -128,6 +130,7 @@ public:
                const amrex::BoxArray& ba,
                amrex::Geometry& geom,
                amrex::MultiFab* cons_in,
+               const amrex::MultiFab* liquid_cloud_fraction,
                amrex::iMultiFab* lmask,
                amrex::MultiFab*  t_surf,
                amrex::Vector<amrex::MultiFab*>& lsm_input_ptrs,
@@ -245,6 +248,7 @@ private:
     // Do we have moisture and cold comps?
     bool m_moist = false;
     bool m_ice   = false;
+    bool m_use_shoc_cldfrac = true;
 
     // Do we have a land surface model?
     bool m_lsm   = false;
@@ -266,6 +270,9 @@ private:
     // Pointer to the CC conserved vars
     amrex::MultiFab* m_cons_in = nullptr;
 
+    // Non-owning Native SHOC liquid-cloud fraction input.
+    const amrex::MultiFab* m_liquid_cloud_fraction = nullptr;
+
     // Pointer to the radiation source terms
     amrex::MultiFab* m_qheating_rates = nullptr;
 
@@ -411,6 +418,8 @@ private:
     real2d_k qv_lay;
     real2d_k qc_lay;
     real2d_k qi_lay;
+    real2d_k cldfrac_liq;
+    real2d_k cldfrac_ice;
     real2d_k cldfrac_tot;
     real2d_k eff_radius_qc;
     real2d_k eff_radius_qi;

Source/PhysicsInterfaces/Radiation/ERF_Radiation.cpp

@@ -40,6 +40,9 @@ Radiation::Radiation (const int& lev,
     // Radiation timestep, as a number of atm steps
     pp.query("rad_freq_in_steps", m_rad_freq_in_steps);
 
+    // Use Native SHOC's diagnosed liquid-cloud fraction when supplied.
+    pp.query("rad_use_shoc_cldfrac", m_use_shoc_cldfrac);
+
     // Get nvar if specified
     pp.query("rad_nvar", m_rad_nvar);
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_rad_nvar >= 0,
@@ -150,6 +153,9 @@ Radiation::Radiation (const int& lev,
                 << m_nswbands << " " << m_nlwbands << "\n";
         Print() << "Number of short/longwave gauss points: "
                 << m_nswgpts  << " " << m_nlwgpts  << "\n";
+        Print() << "Native SHOC cloud fraction for RRTMGP: "
+                << ((m_use_shoc_cldfrac && sc.turbChoice[lev].uses_native_shoc()) ? "enabled" : "disabled")
+                << "\n";
         Print() << "========================================================\n";
     }
 }
@@ -162,6 +168,7 @@ Radiation::set_grids (int& level,
                       const BoxArray& ba,
                       Geometry& geom,
                       MultiFab* cons_in,
+                      const MultiFab* liquid_cloud_fraction,
                       iMultiFab* lmask,
                       MultiFab*  t_surf,
                       Vector<MultiFab*>& lsm_input_ptrs,
@@ -180,12 +187,24 @@ Radiation::set_grids (int& level,
     m_dt             = dt;
     m_geom           = geom;
     m_cons_in        = cons_in;
+    m_liquid_cloud_fraction = liquid_cloud_fraction;
     m_qheating_rates = qheating_rates;
     m_rad_fluxes     = rad_fluxes;
     m_z_phys         = z_phys;
     m_lat            = lat;
     m_lon            = lon;
 
+    if (m_liquid_cloud_fraction != nullptr) {
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->boxArray().ixType().cellCentered(),
+                                         "Radiation liquid cloud fraction must be cell-centered.");
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->nComp() > 0,
+                                         "Radiation liquid cloud fraction must have at least one component.");
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->boxArray() == cons_in->boxArray(),
+                                         "Radiation liquid cloud fraction must match the conserved state BoxArray.");
+        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_liquid_cloud_fraction->DistributionMap() == cons_in->DistributionMap(),
+                                         "Radiation liquid cloud fraction must match the conserved state DistributionMap.");
+    }
+
     // Update the day and month
     time_t timestamp = time_t(time);
     struct tm *timeinfo = gmtime(&timestamp);
@@ -280,6 +299,8 @@ Radiation::alloc_buffers ()
     qv_lay        = real2d_k("qv"           , m_ncol, m_nlay);
     qc_lay        = real2d_k("qc"           , m_ncol, m_nlay);
     qi_lay        = real2d_k("qi"           , m_ncol, m_nlay);
+    cldfrac_liq   = real2d_k("cldfrac_liq"  , m_ncol, m_nlay);
+    cldfrac_ice   = real2d_k("cldfrac_ice"  , m_ncol, m_nlay);
     cldfrac_tot   = real2d_k("cldfrac_tot"  , m_ncol, m_nlay);
     eff_radius_qc = real2d_k("eff_radius_qc", m_ncol, m_nlay);
     eff_radius_qi = real2d_k("eff_radius_qi", m_ncol, m_nlay);
@@ -409,6 +430,8 @@ Radiation::dealloc_buffers ()
     qv_lay            = real2d_k();
     qc_lay            = real2d_k();
     qi_lay            = real2d_k();
+    cldfrac_liq       = real2d_k();
+    cldfrac_ice       = real2d_k();
     cldfrac_tot       = real2d_k();
     eff_radius_qc     = real2d_k();
     eff_radius_qi     = real2d_k();
@@ -478,6 +501,8 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask,
     Table2D<Real,Order::C> qv_lay_tab(qv_lay.data(), {0,0}, {static_cast<int>(qv_lay.extent(0)),static_cast<int>(qv_lay.extent(1))});
     Table2D<Real,Order::C> qc_lay_tab(qc_lay.data(), {0,0}, {static_cast<int>(qc_lay.extent(0)),static_cast<int>(qc_lay.extent(1))});
     Table2D<Real,Order::C> qi_lay_tab(qi_lay.data(), {0,0}, {static_cast<int>(qi_lay.extent(0)),static_cast<int>(qi_lay.extent(1))});
+    Table2D<Real,Order::C> cldfrac_liq_tab(cldfrac_liq.data(), {0,0}, {static_cast<int>(cldfrac_liq.extent(0)),static_cast<int>(cldfrac_liq.extent(1))});
+    Table2D<Real,Order::C> cldfrac_ice_tab(cldfrac_ice.data(), {0,0}, {static_cast<int>(cldfrac_ice.extent(0)),static_cast<int>(cldfrac_ice.extent(1))});
     Table2D<Real,Order::C> cldfrac_tot_tab(cldfrac_tot.data(), {0,0}, {static_cast<int>(cldfrac_tot.extent(0)),static_cast<int>(cldfrac_tot.extent(1))});
 
     Table2D<Real,Order::C> lwp_tab(lwp.data(), {0,0}, {static_cast<int>(lwp.extent(0)),static_cast<int>(lwp.extent(1))});
@@ -497,6 +522,7 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask,
     const bool has_lsm = m_lsm;
     const bool has_lat = m_lat;
     const bool has_lon = m_lon;
+    const bool has_shoc_cldfrac = m_use_shoc_cldfrac && (m_liquid_cloud_fraction != nullptr);
     const bool has_surflayer = (t_surf);
     int  ncol  = m_ncol;
     int  nlay  = m_nlay;
@@ -518,6 +544,8 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask,
                                                           Array4<const Real>{};
         const Array4<const Real>& lon_arr  = (m_lon)    ? m_lon->const_array(mfi) :
                                                           Array4<const Real>{};
+        const Array4<const Real>& shoc_cf_arr = has_shoc_cldfrac
+            ? m_liquid_cloud_fraction->const_array(mfi) : Array4<const Real>{};
         ParallelFor(vbx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             // map [i,j,k] 0-based to [icol, ilay] 0-based
@@ -530,6 +558,10 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask,
             Real qv = (moist) ? std::max(cons_arr(i,j,k,RhoQ1_comp)/r,Real(0.)) : Real(0.);
             Real qc = (moist) ? std::max(cons_arr(i,j,k,RhoQ2_comp)/r,Real(0.)) : Real(0.);
             Real qi = (ice)   ? std::max(cons_arr(i,j,k,RhoQ3_comp)/r,Real(0.)) : Real(0.);
+            const RadiationCloudFractions cloud_fractions =
+                radiation_cloud_fractions(qc, qi,
+                                           has_shoc_cldfrac ? shoc_cf_arr(i,j,k) : Real(0.),
+                                           has_shoc_cldfrac);
 
             // EOS avg to z-face
             Real r_lo   = cons_arr(i,j,k-1,Rho_comp);
@@ -559,7 +591,9 @@ Radiation::mf_to_kokkos_buffers (iMultiFab* lmask,
             qv_lay_tab(icol,ilay) = qv;
             qc_lay_tab(icol,ilay) = qc;
             qi_lay_tab(icol,ilay) = qi;
-            cldfrac_tot_tab(icol,ilay) = ((qc+qi)>Real(0.)) ? Real(1.) : Real(0.);
+            cldfrac_liq_tab(icol,ilay) = cloud_fractions.liquid;
+            cldfrac_ice_tab(icol,ilay) = cloud_fractions.ice;
+            cldfrac_tot_tab(icol,ilay) = cloud_fractions.total;
 
             // NOTE: These are populated in 'mixing_ratio_to_cloud_mass'
             lwp_tab(icol,ilay) = Real(0.);
@@ -1177,8 +1211,8 @@ Radiation::run_impl ()
     Kokkos::deep_copy(mu0, h_mu0);
 
     // Compute layer cloud mass per unit area (populates lwp/iwp)
-    rrtmgp::mixing_ratio_to_cloud_mass(qc_lay, cldfrac_tot, r_lay, z_del, lwp);
-    rrtmgp::mixing_ratio_to_cloud_mass(qi_lay, cldfrac_tot, r_lay, z_del, iwp);
+    rrtmgp::mixing_ratio_to_cloud_mass(qc_lay, cldfrac_liq, r_lay, z_del, lwp);
+    rrtmgp::mixing_ratio_to_cloud_mass(qi_lay, cldfrac_ice, r_lay, z_del, iwp);
 
     // Convert to g/m2 (needed by RRTMGP)
     Table2D<Real,Order::C> lwp_tab(lwp.data(), {0,0}, {static_cast<int>(lwp.extent(0)),static_cast<int>(lwp.extent(1))});

Source/PhysicsInterfaces/Radiation/ERF_RadiationCloudFraction.H

@@ -0,0 +1,40 @@
+#ifndef ERF_RADIATION_CLOUD_FRACTION_H
+#define ERF_RADIATION_CLOUD_FRACTION_H
+
+#include <AMReX_GpuQualifiers.H>
+#include <AMReX_REAL.H>
+
+struct RadiationCloudFractions
+{
+    amrex::Real liquid;
+    amrex::Real ice;
+    amrex::Real total;
+};
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+RadiationCloudFractions
+radiation_cloud_fractions (const amrex::Real qc,
+                           const amrex::Real qi,
+                           const amrex::Real supplied_liquid_fraction,
+                           const bool use_supplied_liquid_fraction) noexcept
+{
+    constexpr amrex::Real min_cloud_fraction = amrex::Real(1.0e-4);
+
+    amrex::Real liquid_cf = amrex::Real(0.0);
+    if (qc > amrex::Real(0.0)) {
+        if (use_supplied_liquid_fraction) {
+            liquid_cf = supplied_liquid_fraction;
+            liquid_cf = (liquid_cf < amrex::Real(0.0)) ? amrex::Real(0.0) : liquid_cf;
+            liquid_cf = (liquid_cf > amrex::Real(1.0)) ? amrex::Real(1.0) : liquid_cf;
+            liquid_cf = (liquid_cf < min_cloud_fraction) ? min_cloud_fraction : liquid_cf;
+        } else {
+            liquid_cf = amrex::Real(1.0);
+        }
+    }
+
+    const amrex::Real ice_cf = (qi > amrex::Real(0.0)) ? amrex::Real(1.0) : amrex::Real(0.0);
+    return {liquid_cf, ice_cf,
+            (liquid_cf > ice_cf) ? liquid_cf : ice_cf};
+}
+
+#endif

Source/PhysicsInterfaces/Radiation/ERF_RadiationInterface.H

@@ -30,6 +30,7 @@ public:
          const amrex::BoxArray& ba,
          amrex::Geometry& geom,
          amrex::MultiFab* cons_in,
+         const amrex::MultiFab* liquid_cloud_fraction,
          amrex::iMultiFab* lmask,
          amrex::MultiFab*  t_surf,
          amrex::Vector<amrex::MultiFab*>& lsm_input_ptrs,

Source/TimeIntegration/ERF_Advance.cpp

@@ -147,7 +147,13 @@ ERF::Advance (int lev, double time, double dt_lev, int iteration, int /*ncycle*/
     // **************************************************************************************
     // Update the radiation sources with the "old" state
     // **************************************************************************************
-    advance_radiation(lev, S_old, dt_lev);
+    const bool use_native_shoc_rad_coupling =
+        solverChoice.rad_type == RadiationType::RRTMGP &&
+        solverChoice.turbChoice[lev].uses_native_shoc();
+
+    if (!use_native_shoc_rad_coupling) {
+        advance_radiation(lev, S_old, dt_lev);
+    }
 
     // **************************************************************************************
     // Update the "old" state using SHOC
@@ -202,6 +208,13 @@ ERF::Advance (int lev, double time, double dt_lev, int iteration, int /*ncycle*/
         }
     }
 
+    if (use_native_shoc_rad_coupling) {
+        AMREX_ALWAYS_ASSERT(native_shoc_driver[lev]);
+        AMREX_ALWAYS_ASSERT(native_shoc_driver[lev]->has_native_diagnostics());
+        advance_radiation(lev, S_old, dt_lev,
+                          &native_shoc_driver[lev]->shoc_cldfrac_diagnostics());
+    }
+
     const BoxArray&            ba = S_old.boxArray();
     const DistributionMapping& dm = S_old.DistributionMap();
 

Source/TimeIntegration/ERF_AdvanceRadiation.cpp

@@ -4,7 +4,8 @@ using namespace amrex;
 
 void ERF::advance_radiation (int lev,
                              MultiFab& cons,
-                             const double& dt_advance)
+                             const double& dt_advance,
+                             const MultiFab* liquid_cloud_fraction)
 {
     if (solverChoice.rad_type != RadiationType::None) {
 #ifdef ERF_USE_NETCDF
@@ -40,6 +41,7 @@ void ERF::advance_radiation (int lev,
         double time_for_rad = t_old[lev] + start_time;
         rad[lev]->Run(lev, istep[lev], time_for_rad, dt_advance,
                       cons.boxArray(), geom[lev], &(cons),
+                      liquid_cloud_fraction,
                       lmask_lev[lev][0].get(), t_surf,
                       lsm_input_ptrs, lsm_output_ptrs,
                       qheating_rates[lev].get(), rad_fluxes[lev].get(),

Tests/Unit/Shoc/CMakeLists.txt

@@ -20,6 +20,7 @@ add_executable(erf_shoc_unit_tests
   ERF_ShocDriverTests.cpp
   ERF_ShocPhysicalPropertyTests.cpp
   ERF_ShocThermoTests.cpp
+  ERF_ShocRadiationCouplingTests.cpp
 )
 
 build_erf_exe(erf_shoc_unit_tests NO_ERF_MAIN)
@@ -30,6 +31,7 @@ target_include_directories(erf_shoc_unit_tests PRIVATE
   ${CMAKE_SOURCE_DIR}/Source/PBL
   ${CMAKE_SOURCE_DIR}/Source/PBL/Shoc
   ${CMAKE_SOURCE_DIR}/Source/Utils
+  ${CMAKE_SOURCE_DIR}/Source/PhysicsInterfaces/Radiation
 )
 
 target_compile_definitions(erf_shoc_unit_tests PRIVATE

Tests/Unit/Shoc/ERF_ShocRadiationCouplingTests.cpp

@@ -0,0 +1,48 @@
+#include "ERF_RadiationCloudFraction.H"
+
+#include <gtest/gtest.h>
+
+TEST(ShocRadiationCloudFraction, BinaryClearAndLiquidFallback)
+{
+    const auto clear = radiation_cloud_fractions(0.0, 0.0, 0.0, false);
+    EXPECT_EQ(clear.liquid, 0.0);
+    EXPECT_EQ(clear.ice, 0.0);
+    EXPECT_EQ(clear.total, 0.0);
+
+    const auto liquid = radiation_cloud_fractions(1.0e-4, 0.0, 0.0, false);
+    EXPECT_EQ(liquid.liquid, 1.0);
+    EXPECT_EQ(liquid.ice, 0.0);
+    EXPECT_EQ(liquid.total, 1.0);
+}
+
+TEST(ShocRadiationCloudFraction, SuppliedLiquidFractionIsBounded)
+{
+    const auto fractional = radiation_cloud_fractions(1.0e-4, 0.0, 0.25, true);
+    EXPECT_DOUBLE_EQ(fractional.liquid, 0.25);
+    EXPECT_EQ(fractional.ice, 0.0);
+    EXPECT_DOUBLE_EQ(fractional.total, 0.25);
+
+    const auto below_zero = radiation_cloud_fractions(1.0e-4, 0.0, -0.5, true);
+    EXPECT_DOUBLE_EQ(below_zero.liquid, 1.0e-4);
+
+    const auto above_one = radiation_cloud_fractions(1.0e-4, 0.0, 1.5, true);
+    EXPECT_DOUBLE_EQ(above_one.liquid, 1.0);
+}
+
+TEST(ShocRadiationCloudFraction, IceRemainsBinaryAndPreservesIceOnlyClouds)
+{
+    const auto ice_only = radiation_cloud_fractions(0.0, 1.0e-4, 0.0, true);
+    EXPECT_EQ(ice_only.liquid, 0.0);
+    EXPECT_EQ(ice_only.ice, 1.0);
+    EXPECT_EQ(ice_only.total, 1.0);
+
+    const auto mixed = radiation_cloud_fractions(1.0e-4, 1.0e-4, 0.3, true);
+    EXPECT_DOUBLE_EQ(mixed.liquid, 0.3);
+    EXPECT_EQ(mixed.ice, 1.0);
+    EXPECT_EQ(mixed.total, 1.0);
+
+    const auto diagnosed_clear = radiation_cloud_fractions(0.0, 0.0, 0.7, true);
+    EXPECT_EQ(diagnosed_clear.liquid, 0.0);
+    EXPECT_EQ(diagnosed_clear.ice, 0.0);
+    EXPECT_EQ(diagnosed_clear.total, 0.0);
+}