Merge pull request erf-model#1286 from xyuan/xyuan/rad_develop

tentatively integrate radiation code to ERF
ewquon · Nov 10, 2023 · fa25dc6 · fa25dc6
2 parents ac2e922 + 6989ae2
commit fa25dc6
Show file tree

Hide file tree

Showing 12 changed files with 155 additions and 157 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -78,6 +78,7 @@ function(build_erf_lib erf_lib_name)
                    ${SRC_DIR}/Radiation/Aero_rad_props.cpp
                    ${SRC_DIR}/Radiation/Optics.cpp
                    ${SRC_DIR}/Radiation/Radiation.cpp
+                   ${SRC_DIR}/Radiation/Albedo.cpp
                    ${CMAKE_SOURCE_DIR}/Submodules/RRTMGP/cpp/examples/mo_load_coefficients.cpp
                    ${CMAKE_SOURCE_DIR}/Submodules/RRTMGP/cpp/extensions/fluxes_byband/mo_fluxes_byband_kernels.cpp
                   )
@@ -146,6 +147,7 @@ function(build_erf_lib erf_lib_name)
        ${SRC_DIR}/TimeIntegration/ERF_TimeStep.cpp
        ${SRC_DIR}/TimeIntegration/ERF_advance_dycore.cpp
        ${SRC_DIR}/TimeIntegration/ERF_advance_microphysics.cpp
+       ${SRC_DIR}/TimeIntegration/ERF_advance_radiation.cpp
        ${SRC_DIR}/TimeIntegration/ERF_make_buoyancy.cpp
        ${SRC_DIR}/TimeIntegration/ERF_make_condensation_source.cpp
        ${SRC_DIR}/TimeIntegration/ERF_make_fast_coeffs.cpp

diff --git a/Source/ERF.H b/Source/ERF.H
@@ -48,6 +48,10 @@
 #include "Microphysics.H"
 #endif
 
+#ifdef ERF_USE_RRTMGP
+#include "Radiation.H"
+#endif
+
 #ifdef ERF_USE_NETCDF
 #include "NCWpsFile.H"
 #endif
@@ -218,6 +222,12 @@ public:
                                const amrex::Real& dt_advance);
 #endif
 
+#if defined(ERF_USE_RRTMGP)
+    void advance_radiation (int lev,
+                            amrex::MultiFab& cons_in,
+                            const amrex::Real& dt_advance);
+#endif
+
     amrex::MultiFab& build_fine_mask (int lev);
 
     void MakeHorizontalAverages ();
@@ -535,6 +545,10 @@ private:
     amrex::Vector<amrex::MultiFab> qmoist; // This has 6 components: qv, qc, qi, qr, qs, qg
 #endif
 
+#if defined(ERF_USE_RRTMGP)
+    Radiation rad;
+#endif
+
     // Fillpatcher classes for coarse-fine boundaries
     int cf_width{0};
     int cf_set_width{0};

diff --git a/Source/Radiation/Aero_rad_props.H b/Source/Radiation/Aero_rad_props.H
@@ -4,7 +4,7 @@
 #ifndef ERF_AER_RAD_PROPS_H_
 #define ERF_AER_RAD_PROPS_H_
 
-#include "ERF.H"
+#include <AMReX_MultiFabUtil.H>
 #include "Mam4_aero.H"
 #include "YAKL_netcdf.h"
 #include "rrtmgp_const.h"

diff --git a/Source/Radiation/Aero_rad_props.cpp b/Source/Radiation/Aero_rad_props.cpp
@@ -3,8 +3,8 @@
 //
 #include <cmath>
 
-#include "Aero_rad_props.H"
 #include "ERF_Constants.H"
+#include "Aero_rad_props.H"
 #include "Water_vapor_saturation.H"
 
 void AerRadProps::initialize() {

diff --git a/Source/Radiation/Albedo.H b/Source/Radiation/Albedo.H
@@ -9,14 +9,5 @@
 */
 #ifndef ERF_ALBEDO_H_
 #define ERF_ALBEDO_H_
-void set_albedo(const real1d& coszrs, real2d& albedo_dir, real2d& albedo_dif)
-{
-  // Albedos for land type I (Briegleb)
-  auto nswbands = albedo_dir.extent(0);
-  auto ncol     = albedo_dif.extent(1);
-  yakl::c::parallel_for(yakl::c::Bounds<2>(nswbands,ncol), YAKL_LAMBDA (int ibnd, int icol) {
-    albedo_dir(ibnd, icol) = 1.4 * 0.24 / ( 1. + 0.8 * coszrs(icol));
-    albedo_dif(ibnd, icol) = 1.2 * 0.24;
-  });
-}
+void set_albedo(const real1d& coszrs, real2d& albedo_dir, real2d& albedo_dif);
 #endif
diff --git a/Source/Radiation/Albedo.cpp b/Source/Radiation/Albedo.cpp
@@ -0,0 +1,15 @@
+
+#include <rrtmgp_const.h>
+#include "Albedo.H"
+
+void set_albedo(const real1d& coszrs, real2d& albedo_dir, real2d& albedo_dif)
+{
+  // Albedos for land type I (Briegleb)
+  auto nswbands = albedo_dir.extent(0);
+  auto ncol     = albedo_dif.extent(1);
+  yakl::c::parallel_for(yakl::c::Bounds<2>(nswbands,ncol), YAKL_LAMBDA (int ibnd, int icol) {
+    albedo_dir(ibnd, icol) = 1.4 * 0.24 / ( 1. + 0.8 * coszrs(icol));
+    albedo_dif(ibnd, icol) = 1.2 * 0.24;
+  });
+}
+
diff --git a/Source/Radiation/Cloud_rad_props.H b/Source/Radiation/Cloud_rad_props.H
@@ -4,9 +4,11 @@
 #ifndef ERF_CLOUD_RAD_PROPS_H_
 #define ERF_CLOUD_RAD_PROPS_H_
 
-#include "ERF.H"
-#include "Linear_interpolate.H"
+#include <AMReX_MultiFabUtil.H>
+#include "ERF_Constants.H"
 #include "rrtmgp_const.h"
+#include "Water_vapor_saturation.H"
+#include "Linear_interpolate.H"
 
 class CloudRadProps {
   public:

diff --git a/Source/Radiation/Cloud_rad_props.cpp b/Source/Radiation/Cloud_rad_props.cpp
@@ -3,8 +3,6 @@
 //
 #include "YAKL_netcdf.h"
 #include "Cloud_rad_props.H"
-#include "Water_vapor_saturation.H"
-#include "Linear_interpolate.H"
 
 void CloudRadProps::initialize() {
     realHost1d g_mu_h;         // mu samples on grid

diff --git a/Source/Radiation/Radiation.H b/Source/Radiation/Radiation.H
@@ -4,7 +4,7 @@
  *                        https://github.com/earth-system-radiation/rte-rrtmgp
  * Please reference to the following paper,
  *                        https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001621
- * NOTE: we use the C++ version of RTE-RRTMGP, which is reimplemented the original Fortran
+ * NOTE: we use the C++ version of RTE-RRTMGP, which is the implemention of the original Fortran
  * code using C++ YAKL for CUDA, HiP and SYCL application by E3SM ECP team, the C++ version
  * of the rte-rrtmgp code is located at:
  *                       https://github.com/E3SM-Project/rte-rrtmgp
@@ -28,6 +28,8 @@
 #include "Rrtmgp.H"
 #include "Optics.H"
 #include "Aero_rad_props.H"
+#include "Parameterizations.H"
+#include "Albedo.H"
 
 // Radiation code interface class
 class Radiation {
@@ -36,57 +38,54 @@ class Radiation {
    ~Radiation() = default;
 
    // init
-   void initialize (const amrex::MultiFab& cons_in,
-                    const amrex::MultiFab& qc_in,
-                    const amrex::MultiFab& qv_in,
-                    const amrex::MultiFab& qi_in,
-                    const amrex::BoxArray& grids,
-                    const amrex::Geometry& geom,
-                    const amrex::Real& dt_advance,
-                    const bool& do_sw_rad,
-                    const bool& do_lw_rad,
-                    const bool& do_aero_rad,
-                    const bool& do_snow_opt,
-                    const bool& is_cmip6_volcano);
+   void initialize(const amrex::MultiFab& cons_in, amrex::MultiFab& qmoist,
+                   const amrex::BoxArray& grids,
+                   const amrex::Geometry& geom,
+                   const amrex::Real& dt_advance,
+                   const bool& do_sw_rad,
+                   const bool& do_lw_rad,
+                   const bool& do_aero_rad,
+                   const bool& do_snow_opt,
+                   const bool& is_cmip6_volcano);
 
    // run radiation model
-   void run ();
+   void run();
 
    // call back
-   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,
-                             FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
-                             real2d& qrl, 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,
-                             FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
-                             real2d& qrs, real2d& qrsc);
-
-   void set_daynight_indices (real1d& coszrs,
-                              int1d& day_indices,
-                              int1d& night_indices);
-
-   void get_gas_vmr (std::vector<std::string>& gas_names,
-                     real3d& gas_vmr);
-
-   void calculate_heating_rate (real2d& flux_up,
-                                real2d& flux_dn,
-                                real2d& pint,
-                                real2d& heating_rate);
-private:
+   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,
+                            FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
+                            real2d& qrl, 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,
+                            FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
+                            real2d& qrs, real2d& qrsc);
+
+   void set_daynight_indices(real1d& coszrs,
+                             int1d& day_indices,
+                             int1d& night_indices);
+
+   void get_gas_vmr(const std::vector<std::string>& gas_names,
+                    real3d& gas_vmr);
+
+   void calculate_heating_rate(real2d& flux_up,
+                               real2d& flux_dn,
+                               real2d& pint,
+                               real2d& heating_rate);
+  private:
    // geometry
    amrex::Geometry m_geom;
 
    // valid boxes on which to evolve the solution
-   amrex::BoxArray m_gtoe;
+   amrex::BoxArray m_box;
 
    // number of vertical levels
    int nlev, zlo, zhi;
@@ -119,19 +118,19 @@ private:
    //
    // Net flux calculated in this routine; used to check energy conservation in
    // the physics package driver?
-   real1d net_flux; //("net_flux",ncol);
+   real1d net_flux;
 
    // This should be module data or something specific to aerosol where it is used?
    bool is_cmip6_volc;    // true if cmip6 style volcanic file is read otherwise false
 
    real dt;               // time step(s) - needed for aerosol optics call
 
    // Surface and top fluxes
-   real1d  fsns;  //(pcols)      ! Surface solar absorbed flux
-   real1d  fsnt;  //(pcols)      ! Net column abs solar flux at model top
-   real1d  flns;  //(pcols)      ! Srf longwave cooling (up-down) flux
-   real1d  flnt;  //(pcols)      ! Net outgoing lw flux at model top
-   real1d  fsds;  //(pcols)      ! Surface solar down flux
+   real1d fsns;  //(ncols)      ! Surface solar absorbed flux
+   real1d fsnt;  //(ncols)      ! Net column abs solar flux at model top
+   real1d flns;  //(ncols)      ! Srf longwave cooling (up-down) flux
+   real1d flnt;  //(ncols)      ! Net outgoing lw flux at model top
+   real1d fsds;  //(ncols)      ! Surface solar down flux
 
    // radiation data
    const std::vector<std::string> active_gases = {
@@ -186,20 +185,17 @@ private:
    int1d rrtmg_to_rrtmgp;
 
    // Pointers to heating rates on physics buffer
-   real2d qrs;  //(:,:)  ! shortwave radiative heating rate
-   real2d qrl;  //(:,:)  ! longwave  radiative heating rate
+   real2d qrs; // shortwave radiative heating rate
+   real2d qrl; // longwave  radiative heating rate
 
    // Pointers to fields on the physics buffer
    real2d zi;
    real2d clear_rh;
 
    // Clear-sky heating rates are not on the physics buffer, and we have no
    // reason to put them there, so declare these are regular arrays here
-   real2d qrsc; //(pcols,pver)
-   real2d qrlc; //(pcols,pver)
-
-   // Temporary variable for heating rate output
-   real2d hr;
+   real2d qrsc;
+   real2d qrlc;
 
    real2d tmid, pmid, pdel;
    real2d pint, tint;