ERF_Kessler.H

Go to the documentation of this file.

#ifndef ERF_Kessler_H
#define ERF_Kessler_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_Constants.H"
#include "ERF_MicrophysicsConstants.H"
#include "ERF_KesslerUtils.H"
#include "ERF_MicrophysicsUtils.H"
#include "ERF_IndexDefines.H"
#include "ERF_DataStruct.H"
#include "ERF_NullMoist.H"
#include "ERF_SurfacePrecipitation.H"
 
namespace MicVar_Kess {
   enum {
      // independent variables
      rho=0, // density
      theta, // liquid/ice water potential temperature
      tabs,  // temperature
      pres,  // pressure
      // non-precipitating vars
      qt,    // total cloud
      qv,    // cloud vapor
      qcl,   // cloud water
      // precipitating vars
      qp,    // total precip
      // derived vars
      rain_accum,
      NumVars
   };
}
 
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
void kessler_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>& qv,
    const amrex::Array4<amrex::Real>& qc,
    const amrex::Array4<amrex::Real>& qp,
    const amrex::Array4<amrex::Real>& qt,
    const amrex::Array4<amrex::Real>& tabs,
    const amrex::Array4<amrex::Real>& pres,
    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);
    const amrex::Real qv_value = states(i,j,k,RhoQ1_comp) / rho_value;
 
    rho(i,j,k) = rho_value;
    theta(i,j,k) = rho_theta / rho_value;
    qv(i,j,k) = qv_value;
    qc(i,j,k) = states(i,j,k,RhoQ2_comp) / rho_value;
    qp(i,j,k) = states(i,j,k,RhoQ3_comp) / rho_value;
    qt(i,j,k) = qv(i,j,k) + qc(i,j,k);
 
    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_value, rdOcp,
        use_anelastic_reference_pressure, p0);
    tabs(i,j,k) = thermo.temperature;
    // Kessler's saturation helpers use mbar / hPa.
    pres(i,j,k) = amrex::Real(0.01) * thermo.pressure_pa;
}
 
class Kessler : public NullMoist {
 
    using FabPtr = std::shared_ptr<amrex::MultiFab>;
 
public:
    // constructor
    Kessler () {}
 
    // destructor
    virtual ~Kessler () = default;
 
    // cloud physics
    void AdvanceKessler (const SolverChoice &solverChoice);
 
    // Set up for first time
    void
    Define (SolverChoice& sc) override
    {
        m_moisture_type = sc.moisture_type;
        // Use one latent-over-cp factor for both the saturation solve and theta
        // update. L_v comes from ERF_Constants.H; c_p comes from SolverChoice.
        m_fac_cond = L_v / sc.c_p;
        m_rdOcp    = sc.rdOcp;
        m_do_cond  = (!sc.uses_shoc_family());
        set_anelastic_reference_pressure_mode(sc);
    }
 
    // init
    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;
 
    // Import minimum dz at this level
    void
    Set_dzmin (const amrex::Real dz_min) override
    {
        m_dzmin = dz_min;
    }
 
    // Copy state into micro vars
    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
    Update_Micro_Vars (amrex::MultiFab& cons_in,
                       const amrex::MultiFab* base_state) override;
 
    // Copy state into micro vars
    void
    Copy_Micro_to_State (amrex::MultiFab& cons_in) override;
 
    // update micro vars
    void
    Update_Micro_Vars (amrex::MultiFab& cons_in) override
    {
        this->Copy_State_to_Micro(cons_in);
    }
 
    // update state vars
    void
    Update_State_Vars (amrex::MultiFab& cons_in,
                       const amrex::MultiFab& /*z_phys_nd*/) override
    {
        this->Copy_Micro_to_State(cons_in);
    }
 
    // wrapper to advance micro vars
    void
    Advance (const amrex::Real& dt_advance,
             const SolverChoice& solverChoice) override
    {
        dt = dt_advance;
 
        this->AdvanceKessler(solverChoice);
    }
 
    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 Kessler::m_qmoist_size; }
 
    int
    Qstate_Moist_Size () override { return Kessler::n_qstate_moist_size; }
 
    int
    Qstate_Moist_NumConc_Size () override { return Kessler::n_qstate_moist_numconc_size; }
 
    void
    Qmoist_Restart_Vars (const SolverChoice& /*a_sc*/,
                         std::vector<int>& a_idx,
                         std::vector<std::string>& a_names) const override
    {
        a_idx.clear();
        a_names.clear();
 
        // NOTE: These are the indices to access into qmoist (not mic_fab_vars)
        a_idx.push_back(0); a_names.push_back("RainAccum");
    }
 
    SurfacePrecipAccumulationSources
    Get_Surface_Precip_Accumulation_Ptrs (const int&) const override
    {
        if (m_moisture_type == MoistureType::Kessler_NoRain) {
            return {};
        }
 
        SurfacePrecipAccumulationSources sources;
        sources.rain = {mic_fab_vars[MicVar_Kess::rain_accum].get(), rhor / amrex::Real(1000.0)};
        return sources;
    }
 
private:
    // Number of qmoist variables (rain_accum)
    int m_qmoist_size = 1;
 
    // Number of qstate variables
    int n_qstate_moist_size = 3;
 
    // Number of qstate variables that are number concentrations
    int n_qstate_moist_numconc_size = 0;
 
    // CFL MAX for vertical advection
    static constexpr amrex::Real CFL_MAX = kessler_sedimentation_cfl_max;
 
    // MicVar map (Qmoist indices -> MicVar enum)
    amrex::Vector<int> MicVarMap;
 
    // geometry
    amrex::Geometry m_geom;
 
    // timestep
    amrex::Real dt;
 
    // Minimum dz at this level
    amrex::Real m_dzmin{0.0};   // set by Set_dzmin; the advance asserts it is positive
 
    // number of vertical levels
    int nlev, zlo, zhi;
 
    // constants
    amrex::Real m_fac_cond;
    amrex::Real m_rdOcp{RdoCp};
    bool m_do_cond;
    MoistureType m_moisture_type = MoistureType::None;
 
    // Pointer to terrain data
    amrex::MultiFab* m_z_phys_nd;
    amrex::MultiFab* m_detJ_cc;
 
    // independent variables
    amrex::Array<FabPtr, MicVar_Kess::NumVars> mic_fab_vars;
};
#endif