ERF_Radiation.H

Go to the documentation of this file.

/*
 * RTE-RRTMGP radiation model interface to ERF
 * The original code is developed by RobertPincus, and the code is open source available at:
 *                        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 the implementation 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
 * The RTE-RRTMGP uses BSD-3-Clause Open Source License, if you want to make changes,
 * and modifications to the code, please refer to BSD-3-Clause Open Source License.
 */
#ifndef ERF_RADIATION_H
#define ERF_RADIATION_H
 
#include <string>
#include <vector>
#include <memory>
 
#include <AMReX_FArrayBox.H>
#include <AMReX_Geometry.H>
#include <AMReX_TableData.H>
#include <AMReX_MultiFabUtil.H>
 
#include "ERF_Config.H"
#include "ERF_Constants.H"
#include "ERF_RadConstants.H"
#include "ERF_RRTMGP.H"
#include "ERF_Optics.H"
#include "ERF_AeroRadProps.H"
#include "ERF_Parameterizations.H"
#include "ERF_Albedo.H"
 
// Radiation code interface class
class Radiation {
  public:
    Radiation () {
        // First, make sure yakl has been initialized
        if (!yakl::isInitialized()) yakl::init();
    }
 
    ~Radiation () = default;
 
    // 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,
                     amrex::Vector<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 ();
 
    // call back
    void on_complete ();
 
    void radiation_driver_lw (int ncol, int nlev,
                              const real3d& gas_vmr,
                              const real2d& pmid, const real2d& pint, const real2d& tmid, const real2d& tint,
                              const real3d& cld_tau_gpt, const real3d& aer_tau_bnd,
                              FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
                              const real2d& qrl, const real2d& qrlc);
 
    void radiation_driver_sw (int ncol,
                              const real3d& gas_vmr, const real2d& pmid, const real2d& pint, const real2d& tmid,
                              const real2d& albedo_dir, const real2d& albedo_dif, const real1d& coszrs,
                              const real3d& cld_tau_gpt, const real3d& cld_ssa_gpt, const real3d& cld_asm_gpt,
                              const real3d& aer_tau_bnd, const real3d& aer_ssa_bnd, const real3d& aer_asm_bnd,
                              FluxesByband& fluxes_clrsky, FluxesByband& fluxes_allsky,
                              const real2d& qrs, const real2d& qrsc);
 
    void set_daynight_indices (const real1d& coszrs,
                               const int1d& day_indices,
                               const int1d& night_indices);
 
    void get_gas_vmr (const std::vector<std::string>& gas_names,
                      const real3d& gas_vmr);
 
    void calculate_heating_rate (const real2d& flux_up,
                                 const real2d& flux_dn,
                                 const real2d& pint,
                                 const real2d& heating_rate);
    void
    export_surface_fluxes(FluxesByband& fluxes,
                          std::string band);
 
    void yakl_to_mf(const real2d &data, amrex::MultiFab &mf);
    void expand_yakl1d_to_mf(const real1d &data, amrex::MultiFab &mf);
 
    void writePlotfile(const std::string& plot_prefix, const amrex::Real time, const int level_step);
 
  private:
    // geometry
    amrex::Geometry m_geom;
 
    // valid boxes on which to evolve the solution
    amrex::BoxArray m_box;
 
    // Pointer to the radiation source terms
    amrex::MultiFab* qrad_src = nullptr;
 
    // Pointer to latitude and longitude
    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;
 
    std::string moisture_type = "None";
    bool has_qmoist;
 
    amrex::Vector<int> rank_offsets;
 
    // Specified uniform angle for radiation
    amrex::Real uniform_angle = 78.463;
 
    // Orbital properties (constant for now)
    static constexpr amrex::Real eccen  =   0.0;                   // Earth's eccentricity factor (unitless) (typically 0 to 0.1)
    static constexpr amrex::Real obliqr =  23.0 * PI / 180.0;      // Earth's obliquity in radians
    static constexpr amrex::Real mvelpp = 103.0 * PI / 180.0 + PI; // Earth's moving vernal equinox longitude of perihelion plus pi (radians)
    static constexpr amrex::Real lambm0 = -3.2503635878519378e-2;  // Mean longitude of perihelion at the vernal equinox (radians)
 
    // number of vertical levels
    int nlev, zlo, zhi;
 
    // number of columns in horizontal plane
    int ncol;
 
    int nlwgpts, nswgpts;
    int nlwbands, nswbands;
 
    // radiation options
    bool do_short_wave_rad;
    bool do_long_wave_rad;
    bool do_snow_optics;
 
    // Flag to indicate whether to do aerosol optical calculations. This
    // zeroes out the aerosol optical properties if False
    bool do_aerosol_rad = true;
 
    // rrtmgp
    Rrtmgp radiation;
 
    // optics radiation properties
    Optics optics;
 
    // aerosol optics properties
    AerRadProps aer_rad;
 
    // input/output of radiation model
    //
    // Net flux calculated in this routine; used to check energy conservation in
    // the physics package driver?
    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;  // Surface solar absorbed flux
    real1d fsnt;  // Net column abs solar flux at model top
    real1d flns;  // Srf longwave cooling (up-down) flux
    real1d flnt;  // Net outgoing lw flux at model top
    real1d fsds;  // Surface solar down flux
 
    // radiation data
    const std::vector<std::string> active_gases = {
                                                   "H2O", "CO2", "O3", "N2O",
                                                   "CO" , "CH4", "O2", "N2" };
 
    bool spectralflux  = false;  // calculate fluxes (up and down) per band.
 
    // Flag to indicate whether or not to use the radiation timestep for solar zenith
    // angle calculations. If true, use the radiation timestep for all solar zenith
    // angle (cosz) calculations.
    // TODO: How does this differ if value is .false.?
    bool use_rad_dt_cosz  = false;
 
    // Value for prescribing an invariant solar constant (i.e. total solar
    // irradiance at TOA). Used for idealized experiments such as RCE.
    // Disabled when value is less than 0.
    real fixed_total_solar_irradiance = -1.;
 
    // The RRTMGP warnings are printed when the state variables need to be limited,
    // such as when the temperature drops too low. This is not normally an issue,
    // but in aquaplanet and RCE configurations these situations occur much more
    // frequently, so this flag was added to be able to disable those messages.
    bool rrtmgp_enable_temperature_warnings = true;
 
    // Output diagnostic brightness temperatures at the top of the
    // atmosphere for 7 TOVS/HIRS channels (2,4,6,8,10,11,12) and 4 TOVS/MSU
    // channels (1,2,3,4).
    // TODO: where are these options set?
    bool dohirs = false;
    int ihirsfq = 1;       // frequency (timesteps) of brightness temperature calcs
 
    // time step to use for the shr_orb_cosz calculation, if use_rad_dt_cosz set to true
    // TODO: where is this set, and what is shr_orb_cosz? Alternative solar zenith
    // angle calculation? What is the other behavior?
    real dt_avg = 0.0;
 
    // k-distribution coefficients files to read from. These are set via namelist
    // variables.
    std::string rrtmgp_data_path;
    std::string rrtmgp_coefficients_file_sw;
    std::string rrtmgp_coefficients_file_lw;
    std::string rrtmgp_coefficients_file_name_sw = "rrtmgp_coefficients_sw_20181204.nc";
    std::string rrtmgp_coefficients_file_name_lw = "rrtmgp_coefficients_lw_20181204.nc";
 
    // Band midpoints; these need to be module variables because of how cam_history works;
    // add_hist_coord sets up pointers to these, so they need to persist.
    real1d sw_band_midpoints;
    real1d lw_band_midpoints;
 
    // rad constituents mixing ratios
    int ngas, naer;
    std::vector<std::string> gasnames;
    std::vector<std::string> aernames;
 
    int1d rrtmg_to_rrtmgp;
 
    // Pointers to heating rates on physics buffer
    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;
    real2d qrlc;
 
    // moisture inputs from ERF
    real2d qt, qi, qc, qn; // [ncol, nlev]
    real2d tmid, pmid, pdel; // [ncol, nlev]
    real2d pint, tint; // [ncol, nlev+1]
    real2d albedo_dir, albedo_dif; // [nswbands, ncol]
 
    // timestep variables -- defined during each call to run()
    real1d coszrs; // [ncol]
    real2d cld, cldfsnow, iclwp, iciwp, icswp, dei, des, lambdac, mu, rei, rel; // [ncol, nlev]
    //   cloud, snow, and aerosol optical properties
    real3d cld_tau_gpt_sw, cld_ssa_gpt_sw, cld_asm_gpt_sw; // [ncol, nlev, nswgpts]
    real3d cld_tau_bnd_sw, cld_ssa_bnd_sw, cld_asm_bnd_sw; // [ncol, nlev, nswbands]
    real3d aer_tau_bnd_sw, aer_ssa_bnd_sw, aer_asm_bnd_sw; // [ncol, nlev, nswbands]
    real3d cld_tau_bnd_lw, aer_tau_bnd_lw; // [ncol, nlev, nlwbands]
    real3d cld_tau_gpt_lw; // [ncol, nlev, nlwgpts]
    //   diagnostic only:
    real3d liq_tau_bnd_sw, ice_tau_bnd_sw, snw_tau_bnd_sw; // [ncol, nlev, nswbands]
    real3d liq_tau_bnd_lw, ice_tau_bnd_lw, snw_tau_bnd_lw; // [ncol, nlev, nlwbands]
    //   gas volume mixing ratios
    real3d gas_vmr; // [ngas, ncol, nlev]
    int1d gpoint_bands_sw, gpoint_bands_lw; // [nswgpts], [nlwgpts]
    //  shortwave radiation vars:
    FluxesByband sw_fluxes_allsky, sw_fluxes_clrsky;
    real1d cld_tau_bnd_sw_1d, cld_ssa_bnd_sw_1d, cld_asm_bnd_sw_1d; // [nswbands]
    real1d cld_tau_bnd_sw_o_1d, cld_ssa_bnd_sw_o_1d, cld_asm_bnd_sw_o_1d; // [nswbands]
    //   aerosol optics:
    real1d aer_tau_bnd_sw_1d, aer_ssa_bnd_sw_1d, aer_asm_bnd_sw_1d; // [nswbands]
    real1d aer_tau_bnd_sw_o_1d, aer_ssa_bnd_sw_o_1d, aer_asm_bnd_sw_o_1d; // [nswbands]
    //  longwave radiation vars:
    FluxesByband lw_fluxes_allsky, lw_fluxes_clrsky;
};
#endif // ERF_RADIATION_H