ERF_WSM6.H

Go to the documentation of this file.

#ifndef ERF_WSM6_H
#define ERF_WSM6_H
 
#include <memory>
#include <string>
#include <vector>
 
#include <AMReX_FArrayBox.H>
#include <AMReX_Geometry.H>
#include <AMReX_MultiFabUtil.H>
 
#include "ERF_Constants.H"
#include "ERF_MicrophysicsConstants.H"
#include "ERF_DataStruct.H"
#include "ERF_IndexDefines.H"
#include "ERF_MicrophysicsUtils.H"
#include "ERF_NullMoist.H"
#include "ERF_SurfacePrecipitation.H"
 
namespace MicVar_WSM6 {
   enum {
      rho = 0,
      theta,
      tabs,
      pres,
      qv,
      qc,
      qi,
      qr,
      qs,
      qg,
      rain_accum,
      snow_accum,
      graup_accum,
      NumVars
   };
}
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void wsm6_copy_state_to_micro_cell (
    const amrex::Array4<const amrex::Real>& states,
    const amrex::Array4<const amrex::Real>& base,
    const amrex::Array4<amrex::Real>& rho,
    const amrex::Array4<amrex::Real>& theta,
    const amrex::Array4<amrex::Real>& tabs,
    const amrex::Array4<amrex::Real>& pres,
    const amrex::Array4<amrex::Real>& qv,
    const amrex::Array4<amrex::Real>& qc,
    const amrex::Array4<amrex::Real>& qi,
    const amrex::Array4<amrex::Real>& qr,
    const amrex::Array4<amrex::Real>& qs,
    const amrex::Array4<amrex::Real>& qg,
    const amrex::Real rdOcp,
    const bool use_anelastic_reference_pressure,
    const int i, const int j, const int k) noexcept
{
    const amrex::Real rho_value = states(i,j,k,Rho_comp);
    const amrex::Real rho_theta = states(i,j,k,RhoTheta_comp);
    rho(i,j,k) = rho_value;
    theta(i,j,k) = rho_theta / rho_value;
    qv(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ1_comp) / rho_value);
    qc(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ2_comp) / rho_value);
    qi(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ3_comp) / rho_value);
    qr(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ4_comp) / rho_value);
    qs(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ5_comp) / rho_value);
    qg(i,j,k) = amrex::max(amrex::Real(0.0), states(i,j,k,RhoQ6_comp) / rho_value);
 
    const amrex::Real p0 = use_anelastic_reference_pressure
        ? base(i,j,k,BaseState::p0_comp) : amrex::Real(0.0);
    const MicrophysicsThermoState thermo = diagnose_microphysics_thermo_state(
        rho_value, rho_theta, qv(i,j,k), rdOcp,
        use_anelastic_reference_pressure, p0);
    tabs(i,j,k) = thermo.temperature;
    pres(i,j,k) = thermo.pressure_pa;
}
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
amrex::Real wsm6_theta_from_temperature_and_pressure (
    const amrex::Real rho,
    const amrex::Real temperature,
    const amrex::Real pressure,
    const amrex::Real qv,
    const amrex::Real configured_rdOcp,
    const bool use_anelastic_reference_pressure) noexcept
{
    return use_anelastic_reference_pressure
        ? getThgivenTandP(temperature, pressure, configured_rdOcp)
        : getThgivenRandT(rho, temperature, RdoCp, qv);
}
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void wsm6_copy_micro_to_state_cell (
    const amrex::Array4<amrex::Real>& states,
    const amrex::Array4<amrex::Real>& theta,
    const amrex::Array4<const amrex::Real>& rho,
    const amrex::Array4<const amrex::Real>& tabs,
    const amrex::Array4<const amrex::Real>& pres,
    const amrex::Array4<const amrex::Real>& qv,
    const amrex::Array4<const amrex::Real>& qc,
    const amrex::Array4<const amrex::Real>& qi,
    const amrex::Array4<const amrex::Real>& qr,
    const amrex::Array4<const amrex::Real>& qs,
    const amrex::Array4<const amrex::Real>& qg,
    const bool use_anelastic_reference_pressure,
    const amrex::Real rdOcp,
    const int i, const int j, const int k) noexcept
{
    const amrex::Real theta_value = wsm6_theta_from_temperature_and_pressure(
        rho(i,j,k), tabs(i,j,k), pres(i,j,k), qv(i,j,k), rdOcp,
        use_anelastic_reference_pressure);
    theta(i,j,k) = theta_value;
    states(i,j,k,RhoTheta_comp) = rho(i,j,k) * theta_value;
    states(i,j,k,RhoQ1_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qv(i,j,k));
    states(i,j,k,RhoQ2_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qc(i,j,k));
    states(i,j,k,RhoQ3_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qi(i,j,k));
    states(i,j,k,RhoQ4_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qr(i,j,k));
    states(i,j,k,RhoQ5_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qs(i,j,k));
    states(i,j,k,RhoQ6_comp) = rho(i,j,k) * amrex::max(amrex::Real(0), qg(i,j,k));
}
 
class WSM6 : public NullMoist {
    using FabPtr = std::shared_ptr<amrex::MultiFab>;
 
public:
    WSM6() {}
    virtual ~WSM6() = default;
 
    void Define(SolverChoice& sc) override
    {
        m_moisture_type = sc.moisture_type;
        m_rdOcp = sc.rdOcp;
        set_anelastic_reference_pressure_mode(sc);
        m_axis = sc.ave_plane;
        m_do_cond = (!sc.uses_shoc_family());
    }
 
    // WSM6 module parameter constants (Fortran: parameter, private)
    static constexpr amrex::Real dtcldcr    = amrex::Real(120.0);
    static constexpr amrex::Real n0r        = amrex::Real(8.0e6);
    static constexpr amrex::Real avtr       = amrex::Real(841.9);
    static constexpr amrex::Real bvtr       = amrex::Real(0.8);
    static constexpr amrex::Real r0         = amrex::Real(0.8e-5);
    static constexpr amrex::Real peaut      = amrex::Real(0.55);
    static constexpr amrex::Real xncr       = amrex::Real(3.0e8);
    static constexpr amrex::Real xmyu       = amrex::Real(1.718e-5);
    static constexpr amrex::Real avts       = amrex::Real(11.72);
    static constexpr amrex::Real bvts       = amrex::Real(0.41);
    static constexpr amrex::Real lamdarmax  = amrex::Real(8.0e4);
    static constexpr amrex::Real lamdasmax  = amrex::Real(1.0e5);
    static constexpr amrex::Real dicon      = amrex::Real(11.9);
    static constexpr amrex::Real dimax      = amrex::Real(500.0e-6);
    static constexpr amrex::Real pfrz1      = amrex::Real(100.0);
    static constexpr amrex::Real pfrz2      = amrex::Real(0.66);
    static constexpr amrex::Real qcrmin     = amrex::Real(1.0e-9);
    static constexpr amrex::Real eacrc      = amrex::Real(1.0);
    static constexpr amrex::Real dens_snow  = amrex::Real(100.0);
    static constexpr amrex::Real qs0        = amrex::Real(6.0e-4);
    static constexpr amrex::Real n0smax     = amrex::Real(1.0e11);
    static constexpr amrex::Real n0s        = amrex::Real(2.0e6);
    static constexpr amrex::Real alpha_wsm6 = amrex::Real(0.12);
 
    void Init(const amrex::MultiFab& cons_in,
              const amrex::BoxArray& grids,
              const amrex::Geometry& geom,
              const amrex::Real& dt_advance,
              std::unique_ptr<amrex::MultiFab>& z_phys_nd,
              std::unique_ptr<amrex::MultiFab>& detJ_cc) override;
 
    void Set_dzmin(const amrex::Real dz_min) override { m_dzmin = dz_min; }
 
    void Copy_State_to_Micro(const amrex::MultiFab& cons_in) override;
    void Copy_State_to_Micro(const amrex::MultiFab& cons_in,
                             const amrex::MultiFab* base_state);
    void Copy_Micro_to_State(amrex::MultiFab& cons_in) override;
 
    void Update_Micro_Vars(amrex::MultiFab& cons_in,
                           const amrex::MultiFab* base_state) override;
 
    void Update_Micro_Vars(amrex::MultiFab& cons_in) override
    {
        Copy_State_to_Micro(cons_in);
    }
 
    void Update_State_Vars(amrex::MultiFab& cons_in,
                           const amrex::MultiFab& /*z_phys_nd*/) override
    {
        Copy_Micro_to_State(cons_in);
    }
 
    void Advance(const amrex::Real& dt_advance,
                 const SolverChoice& solverChoice) override;
 
    amrex::MultiFab* Qmoist_Ptr(const int& varIdx) override
    {
        AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
        return mic_fab_vars[MicVarMap[varIdx]].get();
    }
 
    int Qmoist_Size() override { return m_qmoist_size; }
    int Qstate_Moist_Size() override { return n_qstate_moist_size; }
    int Qstate_Moist_NumConc_Size () override { return n_qstate_moist_numconc_size; }
 
    void Qmoist_Restart_Vars(const SolverChoice&,
                             std::vector<int>& a_idx,
                             std::vector<std::string>& a_names) const override
    {
        a_idx = {0, 1, 2};
        a_names = {"RainAccum", "SnowAccum", "GraupAccum"};
    }
 
    SurfacePrecipAccumulationSources
    Get_Surface_Precip_Accumulation_Ptrs (const int&) const override
    {
        SurfacePrecipAccumulationSources sources;
        sources.total   = {mic_fab_vars[MicVar_WSM6::rain_accum].get(), rhoh2o / amrex::Real(1000.0)};
        sources.snow    = {mic_fab_vars[MicVar_WSM6::snow_accum].get(), rhoh2o / amrex::Real(1000.0)};
        sources.graupel = {mic_fab_vars[MicVar_WSM6::graup_accum].get(), rhoh2o / amrex::Real(1000.0)};
        return sources;
    }
 
private:
    int m_qmoist_size = 3;
    int n_qstate_moist_size = 6;
    int n_qstate_moist_numconc_size = 0;
 
    amrex::Vector<int> MicVarMap;
 
    amrex::Geometry m_geom;
    amrex::Real dt{0.0};
    amrex::Real m_dzmin{0.0};
    int nlev{0}, zlo{0}, zhi{0};
    int m_axis{2};
    bool m_do_cond{true};
    amrex::Real m_rdOcp{RdoCp};
    MoistureType m_moisture_type{MoistureType::None};
 
    amrex::MultiFab* m_z_phys_nd{nullptr};
    amrex::MultiFab* m_detJ_cc{nullptr};
 
    amrex::Array<FabPtr, MicVar_WSM6::NumVars> mic_fab_vars;
 
    // ---- WSM6 init coefficients (from mp_wsm6_init save variables) ----
    // hail_opt-dependent (set in Initialize() based on m_hail_opt)
    bool m_hail_opt{false};
    amrex::Real m_n0g{0}, m_deng{0}, m_avtg{0}, m_bvtg{0}, m_lamdagmax{0};
 
    void initialize_coeffs();
 
    // Derived coefficients (computed in initialize_coeffs(), called from Init())
    amrex::Real m_pi_wsm6{0}, m_xlv1{0};
    amrex::Real m_qc0{0}, m_qck1{0}, m_pidnc{0};
    amrex::Real m_bvtr1{0}, m_bvtr2{0}, m_bvtr3{0}, m_bvtr4{0}, m_bvtr6{0};
    amrex::Real m_g1pbr{0}, m_g3pbr{0}, m_g4pbr{0}, m_g5pbro2{0}, m_g6pbr{0};
    amrex::Real m_pvtr{0}, m_eacrr{0}, m_pacrr{0};
    amrex::Real m_precr1{0}, m_precr2{0}, m_roqimax{0};
    amrex::Real m_bvts1{0}, m_bvts2{0}, m_bvts3{0}, m_bvts4{0};
    amrex::Real m_g1pbs{0}, m_g3pbs{0}, m_g4pbs{0}, m_g5pbso2{0};
    amrex::Real m_pvts{0}, m_pacrs{0}, m_precs1{0}, m_precs2{0};
    amrex::Real m_pidn0r{0}, m_pidn0s{0}, m_pacrc{0};
    amrex::Real m_bvtg1{0}, m_bvtg2{0}, m_bvtg3{0}, m_bvtg4{0};
    amrex::Real m_g1pbg{0}, m_g3pbg{0}, m_g4pbg{0}, m_g5pbgo2{0};
    amrex::Real m_pvtg{0}, m_pacrg{0}, m_precg1{0}, m_precg2{0}, m_pidn0g{0};
    amrex::Real m_rslopermax{0}, m_rslopesmax{0}, m_rslopegmax{0};
    amrex::Real m_rsloperbmax{0}, m_rslopesbmax{0}, m_rslopegbmax{0};
    amrex::Real m_rsloper2max{0}, m_rslopes2max{0}, m_rslopeg2max{0};
    amrex::Real m_rsloper3max{0}, m_rslopes3max{0}, m_rslopeg3max{0};
};
 
namespace WSM6Ind {
    enum {
        // --- inout state (mirrors mp_wsm6_run arguments) ---
        t = 0,
        q,
        qc,
        qi,
        qr,
        qs,
        qg,
        den,
        p,
        delz,
 
        // --- sedimentation working arrays ---
        fallc,
        falkc,
        work1c,
        work2c,
        workr,
        worka,
        den_tmp,
        delz_tmp,
 
        // --- microphysics process rates ---
        pigen,
        pidep,
        pcond,
        prevp,
        psevp,
        pgevp,
        psdep,
        pgdep,
        praut,
        psaut,
        pgaut,
        piacr,
        pracw,
        praci,
        pracs,
        psacw,
        psaci,
        psacr,
        pgacw,
        pgaci,
        pgacr,
        pgacs,
        paacw,
        psmlt,
        pgmlt,
        pseml,
        pgeml,
 
        // --- derived thermodynamic fields ---
        qsum,
        xl,
        cpm,
        work2,
        denfac,
        xni,
        denqrs1,
        denqrs2,
        denqrs3,
        denqci,
        n0sfac,
 
        // --- multi-component arrays expanded from Fortran dimension(:,:,N) ---
        // rh/qsat: N=2, w=over liquid(1), i=over ice(2) — matches Morrison qsatw/qsati
        // slope/fall/qrs/falk/work1: N=3, _r=rain(1) _s=snow(2) _g=graupel(3)
        // workdiff: w=liquid diffusion work1(1), i=ice diffusion work1(2) — separate from work1_r/s/g
        rhw,  rhi,       // relative humidity: w=over liquid, i=over ice (Fortran rh(:,:,1/2))
        qsatw, qsati,    // saturation mixing ratio: w=over liquid, i=over ice (Fortran qsat(:,:,1/2))
        rslope_r,  rslope_s,  rslope_g,
        rslope2_r, rslope2_s, rslope2_g,
        rslope3_r, rslope3_s, rslope3_g,
        rslopeb_r, rslopeb_s, rslopeb_g,
        qrs_tmp_r, qrs_tmp_s, qrs_tmp_g,
        falk_r,  falk_s,  falk_g,
        fall_r,  fall_s,  fall_g,
        work1_r, work1_s, work1_g,  // fall speeds post slope_wsm6: 1=rain 2=snow 3=graupel
        workdiffw, workdiffi,        // diffusion denominators post slope_wsm6 call 2: w=liq i=ice
 
        NumInds
    };
}
 
#endif