Morrison Class Reference

#include <ERF_Morrison.H>

Inheritance diagram for Morrison:

Collaboration diagram for Morrison:

Public Member Functions
	Morrison ()
virtual	~Morrison ()=default
void	Cloud (const SolverChoice &sc)
void	IceFall (const SolverChoice &sc)
void	Precip (const SolverChoice &sc)
void	PrecipFall (const SolverChoice &sc)
void	Define (SolverChoice &sc) override
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	Copy_State_to_Micro (const amrex::MultiFab &cons_in) override
void	Copy_Micro_to_State (amrex::MultiFab &cons_in) override
void	Update_Micro_Vars (amrex::MultiFab &cons_in) override
void	Update_State_Vars (amrex::MultiFab &cons_in) override
void	Advance (const amrex::Real &dt_advance, const SolverChoice &sc) override
amrex::MultiFab *	Qmoist_Ptr (const int &varIdx) override
void	Compute_Coefficients ()
int	Qmoist_Size () override
int	Qstate_Size () override
void	Qmoist_Restart_Vars (const SolverChoice &, std::vector< int > &a_idx, std::vector< std::string > &a_names) const override
Public Member Functions inherited from NullMoist
	NullMoist ()
virtual	~NullMoist ()=default

Static Public Member Functions
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real	NewtonIterSat (int &i, int &j, int &k, const amrex::Real &fac_cond, const amrex::Real &fac_fus, const amrex::Real &fac_sub, const amrex::Real &an, const amrex::Real &bn, const amrex::Array4< amrex::Real > &tabs_array, const amrex::Array4< amrex::Real > &pres_array, const amrex::Array4< amrex::Real > &qv_array, const amrex::Array4< amrex::Real > &qc_array, const amrex::Array4< amrex::Real > &qi_array, const amrex::Array4< amrex::Real > &qn_array, const amrex::Array4< amrex::Real > &qt_array)

Private Types
using	FabPtr = std::shared_ptr< amrex::MultiFab >

Private Attributes
int	m_qmoist_size = 3
int	m_qstate_size = 6
amrex::Vector< int >	MicVarMap
amrex::Geometry	m_geom
amrex::BoxArray	m_gtoe
amrex::Real	dt
int	nlev
int	zlo
int	zhi
int	m_axis
amrex::Real	m_fac_cond
amrex::Real	m_fac_fus
amrex::Real	m_fac_sub
amrex::Real	m_gOcp
amrex::Real	m_rdOcp
amrex::MultiFab *	m_z_phys_nd
amrex::MultiFab *	m_detJ_cc
amrex::Array< FabPtr, MicVar_Morr::NumVars >	mic_fab_vars
amrex::TableData< amrex::Real, 1 >	accrrc
amrex::TableData< amrex::Real, 1 >	accrsi
amrex::TableData< amrex::Real, 1 >	accrsc
amrex::TableData< amrex::Real, 1 >	coefice
amrex::TableData< amrex::Real, 1 >	evaps1
amrex::TableData< amrex::Real, 1 >	evaps2
amrex::TableData< amrex::Real, 1 >	accrgi
amrex::TableData< amrex::Real, 1 >	accrgc
amrex::TableData< amrex::Real, 1 >	evapg1
amrex::TableData< amrex::Real, 1 >	evapg2
amrex::TableData< amrex::Real, 1 >	evapr1
amrex::TableData< amrex::Real, 1 >	evapr2
amrex::TableData< amrex::Real, 1 >	rho1d
amrex::TableData< amrex::Real, 1 >	pres1d
amrex::TableData< amrex::Real, 1 >	tabs1d
amrex::TableData< amrex::Real, 1 >	qt1d
amrex::TableData< amrex::Real, 1 >	qv1d
amrex::TableData< amrex::Real, 1 >	qn1d
amrex::TableData< amrex::Real, 1 >	gamaz
amrex::TableData< amrex::Real, 1 >	zmid

Static Private Attributes
static constexpr amrex::Real	CFL_MAX = 0.5

Member Typedef Documentation

FabPtr

using Morrison::FabPtr = std::shared_ptr<amrex::MultiFab>

private

Constructor & Destructor Documentation

Morrison()

Morrison::Morrison ( )

inline

{}

~Morrison()

virtual Morrison::~Morrison ( )

virtualdefault

Member Function Documentation

Advance()

void Morrison::Advance ( const amrex::Real &  dt_advance, const SolverChoice &  sc  )

inlineoverridevirtual

Reimplemented from NullMoist.

    {
        dt = dt_advance;
 
        this->Cloud(sc);
        this->IceFall(sc);
        this->Precip(sc);
        this->PrecipFall(sc);
    }

Here is the call graph for this function:

Cloud()

void Morrison::Cloud

(

const SolverChoice &

sc

)

Split cloud components according to saturation pressures; source theta from latent heat.

{
 
    constexpr Real an = 1.0/(tbgmax-tbgmin);
    constexpr Real bn = tbgmin*an;
 
    Real fac_cond = m_fac_cond;
    Real fac_sub  = m_fac_sub;
    Real fac_fus  = m_fac_fus;
    Real rdOcp    = m_rdOcp;
 
    for ( MFIter mfi(*(mic_fab_vars[MicVar_Morr::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        auto  qt_array = mic_fab_vars[MicVar_Morr::qt]->array(mfi);
        auto  qn_array = mic_fab_vars[MicVar_Morr::qn]->array(mfi);
        auto  qv_array = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
        auto qcl_array = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
        auto qci_array = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
 
        auto   rho_array = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
        auto  tabs_array = mic_fab_vars[MicVar_Morr::tabs]->array(mfi);
        auto theta_array = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
        auto  pres_array = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
 
        const auto& box3d = mfi.tilebox();
 
        ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            // Saturation moisture fractions
            Real omn;
            Real qsat;
            Real qsatw;
            Real qsati;
 
            // Newton iteration vars
            Real delta_qv, delta_qc, delta_qi;
 
            // NOTE: Conversion before iterations is necessary to
            //       convert cloud water to ice or vice versa.
            //       This ensures the omn splitting is enforced
            //       before the Newton iteration, which assumes it is.
 
            omn = 1.0;
            // Cloud ice not permitted (melt to form water)
            if (tabs_array(i,j,k) >= tbgmax) {
                omn = 1.0;
                delta_qi = qci_array(i,j,k);
                qci_array(i,j,k)   = 0.0;
                qcl_array(i,j,k)  += delta_qi;
                tabs_array(i,j,k) -= fac_fus * delta_qi;
                pres_array(i,j,k)  = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
                                     * (1.0 + R_v/R_d * qv_array(i,j,k));
                theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
                pres_array(i,j,k) *= 0.01;
            }
            // Cloud water not permitted (freeze to form ice)
            else if (tabs_array(i,j,k) <= tbgmin) {
                omn = 0.0;
                delta_qc = qcl_array(i,j,k);
                qcl_array(i,j,k)   = 0.0;
                qci_array(i,j,k)  += delta_qc;
                tabs_array(i,j,k) += fac_fus * delta_qc;
                pres_array(i,j,k)  = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
                                     * (1.0 + R_v/R_d * qv_array(i,j,k));
                theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
                pres_array(i,j,k) *= 0.01;
            }
            // Mixed cloud phase (split according to omn)
            else {
                omn = an*tabs_array(i,j,k)-bn;
                delta_qc = qcl_array(i,j,k) - qn_array(i,j,k) * omn;
                delta_qi = qci_array(i,j,k) - qn_array(i,j,k) * (1.0 - omn);
                qcl_array(i,j,k)   = qn_array(i,j,k) * omn;
                qci_array(i,j,k)   = qn_array(i,j,k) * (1.0 - omn);
                tabs_array(i,j,k) += fac_fus * delta_qc;
                pres_array(i,j,k)  = rho_array(i,j,k) * R_d * tabs_array(i,j,k)
                                     * (1.0 + R_v/R_d * qv_array(i,j,k));
                theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), pres_array(i,j,k), rdOcp);
                pres_array(i,j,k) *= 0.01;
            }
 
            // Saturation moisture fractions
            erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
            erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
            qsat = omn * qsatw  + (1.0-omn) * qsati;
 
            // We have enough total moisture to relax to equilibrium
            if (qt_array(i,j,k) > qsat) {
 
                // Update temperature
                tabs_array(i,j,k) = NewtonIterSat(i, j, k   ,
                                                  fac_cond  , fac_fus   , fac_sub ,
                                                  an        , bn        ,
                                                  tabs_array, pres_array,
                                                  qv_array  , qcl_array  , qci_array,
                                                  qn_array  , qt_array);
 
                // Update theta
                theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
 
            //
            // We cannot blindly relax to qsat, but we can convert qc/qi -> qv.
            // The concept here is that if we put all the moisture into qv and modify
            // the temperature, we can then check if qv > qsat occurs (for final T/P/qv).
            // If the reduction in T/qsat and increase in qv does trigger the
            // aforementioned condition, we can do Newton iteration to drive qv = qsat.
            //
            } else {
                // Changes in each component
                delta_qv = qcl_array(i,j,k) + qci_array(i,j,k);
                delta_qc = qcl_array(i,j,k);
                delta_qi = qci_array(i,j,k);
 
                // Partition the change in non-precipitating q
                 qv_array(i,j,k) += delta_qv;
                qcl_array(i,j,k)  = 0.0;
                qci_array(i,j,k)  = 0.0;
                 qn_array(i,j,k)  = 0.0;
                 qt_array(i,j,k)  = qv_array(i,j,k);
 
                // Update temperature (endothermic since we evap/sublime)
                tabs_array(i,j,k) -= fac_cond * delta_qc + fac_sub * delta_qi;
 
                // Update theta
                theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
 
                // Verify assumption that qv > qsat does not occur
                erf_qsatw(tabs_array(i,j,k), pres_array(i,j,k), qsatw);
                erf_qsati(tabs_array(i,j,k), pres_array(i,j,k), qsati);
                qsat = omn * qsatw  + (1.0-omn) * qsati;
                if (qt_array(i,j,k) > qsat) {
 
                    // Update temperature
                    tabs_array(i,j,k) = NewtonIterSat(i, j, k   ,
                                                      fac_cond  , fac_fus   , fac_sub ,
                                                      an        , bn        ,
                                                      tabs_array, pres_array,
                                                      qv_array  , qcl_array  , qci_array,
                                                      qn_array  , qt_array);
 
                    // Update theta
                    theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
 
                }
            }
        });
    } // mfi
}

Referenced by Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

Compute_Coefficients()

void Morrison::Compute_Coefficients

(

)

{
    auto dz   = m_geom.CellSize(2);
    auto lowz = m_geom.ProbLo(2);
 
    auto accrrc_t  = accrrc.table();
    auto accrsi_t  = accrsi.table();
    auto accrsc_t  = accrsc.table();
    auto coefice_t = coefice.table();
    auto evaps1_t  = evaps1.table();
    auto evaps2_t  = evaps2.table();
    auto accrgi_t  = accrgi.table();
    auto accrgc_t  = accrgc.table();
    auto evapg1_t  = evapg1.table();
    auto evapg2_t  = evapg2.table();
    auto evapr1_t  = evapr1.table();
    auto evapr2_t  = evapr2.table();
 
    auto rho1d_t  = rho1d.table();
    auto pres1d_t = pres1d.table();
    auto tabs1d_t = tabs1d.table();
 
    auto gamaz_t  = gamaz.table();
    auto zmid_t   = zmid.table();
 
    Real gam3  = erf_gammafff(3.0             );
    Real gamr1 = erf_gammafff(3.0+b_rain      );
    Real gamr2 = erf_gammafff((5.0+b_rain)/2.0);
    Real gams1 = erf_gammafff(3.0+b_snow      );
    Real gams2 = erf_gammafff((5.0+b_snow)/2.0);
    Real gamg1 = erf_gammafff(3.0+b_grau      );
    Real gamg2 = erf_gammafff((5.0+b_grau)/2.0);
 
    // calculate the plane average variables
    PlaneAverage rho_ave(mic_fab_vars[MicVar_Morr::rho].get(), m_geom, m_axis);
    PlaneAverage theta_ave(mic_fab_vars[MicVar_Morr::theta].get(), m_geom, m_axis);
    PlaneAverage qv_ave(mic_fab_vars[MicVar_Morr::qv].get(), m_geom, m_axis);
    rho_ave.compute_averages(ZDir(), rho_ave.field());
    theta_ave.compute_averages(ZDir(), theta_ave.field());
    qv_ave.compute_averages(ZDir(), qv_ave.field());
 
    // get host variable rho, and rhotheta
    int ncell = rho_ave.ncell_line();
 
    Gpu::HostVector<Real> rho_h(ncell), theta_h(ncell), qv_h(ncell);
    rho_ave.line_average(0, rho_h);
    theta_ave.line_average(0, theta_h);
    qv_ave.line_average(0, qv_h);
 
    // copy data to device
    Gpu::DeviceVector<Real> rho_d(ncell), theta_d(ncell), qv_d(ncell);
    Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
    Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin());
    Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin());
    Gpu::streamSynchronize();
 
    Real* rho_dptr   = rho_d.data();
    Real* theta_dptr = theta_d.data();
    Real* qv_dptr    = qv_d.data();
 
    Real gOcp = m_gOcp;
 
    ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
    {
        Real RhoTheta = rho_dptr[k]*theta_dptr[k];
        Real pressure = getPgivenRTh(RhoTheta, qv_dptr[k]);
        rho1d_t(k)    = rho_dptr[k];
        pres1d_t(k)   = pressure*0.01;
        // NOTE: Limit the temperature to the melting point of ice to avoid a divide by
        //       0 condition when computing the cold evaporation coefficients. This should
        //       not affect results since evporation requires snow/graupel to be present
        //       and thus T<273.16
        tabs1d_t(k)   = std::min(getTgivenRandRTh(rho_dptr[k], RhoTheta, qv_dptr[k]),273.16);
        zmid_t(k)     = lowz + (k+0.5)*dz;
        gamaz_t(k)    = gOcp*zmid_t(k);
    });
 
    if(round(gam3) != 2) {
        std::cout << "cannot compute gamma-function in Microphysics::Init" << std::endl;
        std::exit(-1);
    }
 
    // Populate all the coefficients
    ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
    {
        Real Prefactor;
        Real pratio = sqrt(1.29 / rho1d_t(k));
        //Real rrr1   = 393.0/(tabs1d_t(k)+120.0)*std::pow((tabs1d_t(k)/273.0),1.5);
        //Real rrr2   = std::pow((tabs1d_t(k)/273.0),1.94)*(1000.0/pres1d_t(k));
        Real estw   = 100.0*erf_esatw(tabs1d_t(k));
        Real esti   = 100.0*erf_esati(tabs1d_t(k));
 
        // accretion by snow:
        Real coef1   = 0.25 * PI * nzeros * a_snow * gams1 * pratio/pow((PI * rhos * nzeros/rho1d_t(k) ) , ((3.0+b_snow)/4.0));
        Real coef2   = exp(0.025*(tabs1d_t(k) - 273.15));
        accrsi_t(k)  =  coef1 * coef2 * esicoef;
        accrsc_t(k)  =  coef1 * esccoef;
        coefice_t(k) =  coef2;
 
        // evaporation of snow:
        coef1 = (lsub/(tabs1d_t(k)*R_v)-1.0)*lsub/(therco*tabs1d_t(k));
        coef2 = R_v * R_d / (diffelq * esti);
        Prefactor = 2.0 * PI * nzeros / (rho1d_t(k) * (coef1 + coef2));
        Prefactor *= (2.0/PI); // Shape factor snow
        evaps1_t(k) = Prefactor * 0.65 * sqrt(rho1d_t(k) / (PI * rhos * nzeros));
        evaps2_t(k) = Prefactor * 0.44 * sqrt(a_snow * rho1d_t(k) / muelq) * gams2
                    * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhos * nzeros) , ((5.0+b_snow)/8.0));
 
        // accretion by graupel:
        coef1 = 0.25*PI*nzerog*a_grau*gamg1*pratio/pow((PI*rhog*nzerog/rho1d_t(k)) , ((3.0+b_grau)/4.0));
        coef2 = exp(0.025*(tabs1d_t(k) - 273.15));
        accrgi_t(k) = coef1 * coef2 * egicoef;
        accrgc_t(k) = coef1 * egccoef;
 
        // evaporation of graupel:
        coef1 = (lsub/(tabs1d_t(k)*R_v)-1.0)*lsub/(therco*tabs1d_t(k));
        coef2 = R_v * R_d / (diffelq * esti);
        Prefactor = 2.0 * PI * nzerog / (rho1d_t(k) * (coef1 + coef2)); // Shape factor for graupel is 1
        evapg1_t(k) = Prefactor * 0.78 * sqrt(rho1d_t(k) / (PI * rhog * nzerog));
        evapg2_t(k) = Prefactor * 0.31 * sqrt(a_grau * rho1d_t(k) / muelq) * gamg2
                    * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhog * nzerog) , ((5.0+b_grau)/8.0));
 
        // accretion by rain:
        accrrc_t(k) = 0.25 * PI * nzeror * a_rain * gamr1 * pratio/pow((PI * rhor * nzeror / rho1d_t(k)) , ((3.0+b_rain)/4.))* erccoef;
 
        // evaporation of rain:
        coef1 = (lcond/(tabs1d_t(k)*R_v)-1.0)*lcond/(therco*tabs1d_t(k));
        coef2 = R_v * R_d / (diffelq * estw);
        Prefactor = 2.0 * PI * nzeror / (rho1d_t(k) * (coef1 + coef2)); // Shape factor for rain is 1
        evapr1_t(k) = Prefactor * 0.78 * sqrt(rho1d_t(k) / (PI * rhor * nzeror));
        evapr2_t(k) = Prefactor * 0.31 * sqrt(a_rain * rho1d_t(k) / muelq) * gamr2
                    * sqrt(pratio) * pow(rho1d_t(k) / (PI * rhor * nzeror) , ((5.0+b_rain)/8.0));
    });
}

Referenced by Update_Micro_Vars().

Here is the call graph for this function:

Here is the caller graph for this function:

Copy_Micro_to_State()

void Morrison::Copy_Micro_to_State ( amrex::MultiFab &  cons_in )

overridevirtual

Updates conserved and microphysics variables in the provided MultiFabs from the internal MultiFabs that store Microphysics module data.

Parameters

[out]	cons	Conserved variables
[out]	qmoist	qv, qc, qi, qr, qs, qg

Reimplemented from NullMoist.

{
    // Get the temperature, density, theta, qt and qp from input
    for ( MFIter mfi(cons,TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        const auto& box3d = mfi.tilebox();
 
        auto states_arr = cons.array(mfi);
 
        auto rho_arr    = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
        auto theta_arr  = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
 
        auto qv_arr     = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
        auto qc_arr     = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
        auto qi_arr     = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
 
        auto qpr_arr     = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
        auto qps_arr     = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
        auto qpg_arr     = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
 
        // get potential total density, temperature, qt, qp
        ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            states_arr(i,j,k,RhoTheta_comp) = rho_arr(i,j,k)*theta_arr(i,j,k);
 
            states_arr(i,j,k,RhoQ1_comp)    = rho_arr(i,j,k)*std::max(0.0,qv_arr(i,j,k));
            states_arr(i,j,k,RhoQ2_comp)    = rho_arr(i,j,k)*std::max(0.0,qc_arr(i,j,k));
            states_arr(i,j,k,RhoQ3_comp)    = rho_arr(i,j,k)*std::max(0.0,qi_arr(i,j,k));
 
            states_arr(i,j,k,RhoQ4_comp)    = rho_arr(i,j,k)*std::max(0.0,qpr_arr(i,j,k));
            states_arr(i,j,k,RhoQ5_comp)    = rho_arr(i,j,k)*std::max(0.0,qps_arr(i,j,k));
            states_arr(i,j,k,RhoQ6_comp)    = rho_arr(i,j,k)*std::max(0.0,qpg_arr(i,j,k));
        });
    }
 
    // Fill interior ghost cells and periodic boundaries
    cons.FillBoundary(m_geom.periodicity());
}

Referenced by Update_State_Vars().

Here is the caller graph for this function:

Copy_State_to_Micro()

void Morrison::Copy_State_to_Micro ( const amrex::MultiFab &  cons_in )

overridevirtual

Initializes the Microphysics module.

Parameters

[in]

cons_in

Conserved variables input

Reimplemented from NullMoist.

{
    // Get the temperature, density, theta, qt and qp from input
    for ( MFIter mfi(cons_in); mfi.isValid(); ++mfi) {
        const auto& box3d = mfi.growntilebox();
 
        auto states_array = cons_in.array(mfi);
 
        // Non-precipitating
        auto qv_array    = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
        auto qc_array    = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
        auto qi_array    = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
        auto qn_array    = mic_fab_vars[MicVar_Morr::qn]->array(mfi);
        auto qt_array    = mic_fab_vars[MicVar_Morr::qt]->array(mfi);
 
        // Precipitating
        auto qpr_array   = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
        auto qps_array   = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
        auto qpg_array   = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
        auto qp_array    = mic_fab_vars[MicVar_Morr::qp]->array(mfi);
 
        auto rho_array   = mic_fab_vars[MicVar_Morr::rho]->array(mfi);
        auto theta_array = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
        auto tabs_array  = mic_fab_vars[MicVar_Morr::tabs]->array(mfi);
        auto pres_array  = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
 
        // Get pressure, theta, temperature, density, and qt, qp
        ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
        {
            rho_array(i,j,k)   = states_array(i,j,k,Rho_comp);
            theta_array(i,j,k) = states_array(i,j,k,RhoTheta_comp)/states_array(i,j,k,Rho_comp);
 
            qv_array(i,j,k)    = std::max(0.0,states_array(i,j,k,RhoQ1_comp)/states_array(i,j,k,Rho_comp));
            qc_array(i,j,k)    = std::max(0.0,states_array(i,j,k,RhoQ2_comp)/states_array(i,j,k,Rho_comp));
            qi_array(i,j,k)    = std::max(0.0,states_array(i,j,k,RhoQ3_comp)/states_array(i,j,k,Rho_comp));
            qn_array(i,j,k)    = qc_array(i,j,k) + qi_array(i,j,k);
            qt_array(i,j,k)    = qv_array(i,j,k) + qn_array(i,j,k);
 
            qpr_array(i,j,k)   = std::max(0.0,states_array(i,j,k,RhoQ4_comp)/states_array(i,j,k,Rho_comp));
            qps_array(i,j,k)   = std::max(0.0,states_array(i,j,k,RhoQ5_comp)/states_array(i,j,k,Rho_comp));
            qpg_array(i,j,k)   = std::max(0.0,states_array(i,j,k,RhoQ6_comp)/states_array(i,j,k,Rho_comp));
             qp_array(i,j,k)   = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
 
            tabs_array(i,j,k)  = getTgivenRandRTh(states_array(i,j,k,Rho_comp),
                                                  states_array(i,j,k,RhoTheta_comp),
                                                  qv_array(i,j,k));
            pres_array(i,j,k)  = getPgivenRTh(states_array(i,j,k,RhoTheta_comp), qv_array(i,j,k)) * 0.01;
        });
    }
}

Referenced by Update_Micro_Vars().

Here is the call graph for this function:

Here is the caller graph for this function:

Define()

void Morrison::Define ( SolverChoice &  sc )

inlineoverridevirtual

Reimplemented from NullMoist.

    {
        m_fac_cond = lcond / sc.c_p;
        m_fac_fus  = lfus / sc.c_p;
        m_fac_sub  = lsub / sc.c_p;
        m_gOcp     = CONST_GRAV / sc.c_p;
        m_axis     = sc.ave_plane;
        m_rdOcp    = sc.rdOcp;
    }

IceFall()

void Morrison::IceFall

(

const SolverChoice &

sc

)

Sedimentation of cloud ice (A32)

{
    Real dz   = m_geom.CellSize(2);
    Real dtn  = dt;
    Real coef = dtn/dz;
 
    auto domain = m_geom.Domain();
    int k_lo = domain.smallEnd(2);
    int k_hi = domain.bigEnd(2);
 
    auto qcl   = mic_fab_vars[MicVar_Morr::qcl];
    auto qci   = mic_fab_vars[MicVar_Morr::qci];
    auto qn    = mic_fab_vars[MicVar_Morr::qn];
    auto qt    = mic_fab_vars[MicVar_Morr::qt];
    auto rho   = mic_fab_vars[MicVar_Morr::rho];
    auto tabs  = mic_fab_vars[MicVar_Morr::tabs];
 
    MultiFab fz;
    IntVect  ng = qcl->nGrowVect();
    BoxArray ba = qcl->boxArray();
    DistributionMapping dm = qcl->DistributionMap();
    fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ng);
    fz.setVal(0.);
 
    for (MFIter mfi(fz, TileNoZ()); mfi.isValid(); ++mfi) {
        auto qci_array = qci->array(mfi);
        auto rho_array = rho->array(mfi);
        auto fz_array  = fz.array(mfi);
 
        const auto& box3d  = mfi.tilebox();
 
        ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            Real rho_avg, qci_avg;
            if (k==k_lo) {
                rho_avg = rho_array(i,j,k);
                qci_avg = qci_array(i,j,k);
            } else if (k==k_hi+1) {
                rho_avg = rho_array(i,j,k-1);
                qci_avg = qci_array(i,j,k-1);
            } else {
                rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
                qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k));
            }
            Real vt_ice = min( 0.4 , 8.66 * pow( (std::max(0.,qci_avg)+1.e-10) , 0.24) );
 
            // NOTE: Fz is the sedimentation flux from the advective operator.
            //       In the terrain-following coordinate system, the z-deriv in
            //       the divergence uses the normal velocity (Omega). However,
            //       there are no u/v components to the sedimentation velocity.
            //       Therefore, we simply end up with a division by detJ when
            //       evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
            fz_array(i,j,k) = rho_avg*vt_ice*qci_avg;
        });
    }
 
    // Compute number of substeps from maximum terminal velocity
    Real wt_max;
    int n_substep;
    auto const& ma_fz_arr = fz.const_arrays();
    GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
                                   TypeList<Real>{},
                                   fz, IntVect(0),
                         [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
                         -> GpuTuple<Real>
                         {
                             return { ma_fz_arr[box_no](i,j,k) };
                         });
    wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
    n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
    AMREX_ALWAYS_ASSERT(n_substep >= 1);
    coef /= Real(n_substep);
    dtn  /= Real(n_substep);
 
    // Substep the vertical advection
    for (int nsub(0); nsub<n_substep; ++nsub) {
        for (MFIter mfi(*qci, TileNoZ()); mfi.isValid(); ++mfi) {
            auto qci_array   = qci->array(mfi);
            auto qn_array    = qn->array(mfi);
            auto qt_array    = qt->array(mfi);
            auto rho_array   = rho->array(mfi);
            auto fz_array    = fz.array(mfi);
 
            const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};
 
            const auto& tbx  = mfi.tilebox();
            const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));
 
            // Update vertical flux every substep
            ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                Real rho_avg, qci_avg;
                if (k==k_lo) {
                    rho_avg = rho_array(i,j,k);
                    qci_avg = qci_array(i,j,k);
                } else if (k==k_hi+1) {
                    rho_avg = rho_array(i,j,k-1);
                    qci_avg = qci_array(i,j,k-1);
                } else {
                    rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k));
                    qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k));
                }
                Real vt_ice = min( 0.4 , 8.66 * pow( (std::max(0.,qci_avg)+1.e-10) , 0.24) );
 
                // NOTE: Fz is the sedimentation flux from the advective operator.
                //       In the terrain-following coordinate system, the z-deriv in
                //       the divergence uses the normal velocity (Omega). However,
                //       there are no u/v components to the sedimentation velocity.
                //       Therefore, we simply end up with a division by detJ when
                //       evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
                fz_array(i,j,k) = rho_avg*vt_ice*qci_avg;
            });
 
            // Update precip every substep
            ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                // Jacobian determinant
                Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;
 
                //==================================================
                // Cloud ice sedimentation (A32)
                //==================================================
                Real dqi  = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
                dqi = std::max(-qci_array(i,j,k), dqi);
 
                // Add this increment to both non-precipitating and total water.
                qci_array(i,j,k) += dqi;
                 qn_array(i,j,k) += dqi;
                 qt_array(i,j,k) += dqi;
 
                // NOTE: Sedimentation does not affect the potential temperature,
                //       but it does affect the liquid/ice static energy.
                //       No source to Theta occurs here.
            });
        } // mfi
    } // nsub
}

Referenced by Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

Init()

void Morrison::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  )

overridevirtual

Initializes the Microphysics module.

Parameters

[in]	cons_in	Conserved variables input
[in]	qc_in	Cloud variables input
[in,out]	qv_in	Vapor variables input
[in]	qi_in	Ice variables input
[in]	grids	The boxes on which we will evolve the solution
[in]	geom	Geometry associated with these MultiFabs and grids
[in]	dt_advance	Timestep for the advance

Reimplemented from NullMoist.

{
    amrex::Abort("Morrison not actually implemented yet; this is SAM");
 
    dt     = dt_advance;
    m_geom = geom;
    m_gtoe = grids;
 
    m_z_phys_nd = z_phys_nd.get();
    m_detJ_cc   = detJ_cc.get();
 
    MicVarMap.resize(m_qmoist_size);
    MicVarMap = {MicVar_Morr::rain_accum, MicVar_Morr::snow_accum, MicVar_Morr::graup_accum};
 
    // initialize microphysics variables
    for (auto ivar = 0; ivar < MicVar_Morr::NumVars; ++ivar) {
        mic_fab_vars[ivar] = std::make_shared<MultiFab>(cons_in.boxArray(), cons_in.DistributionMap(),
                                                        1, cons_in.nGrowVect());
        mic_fab_vars[ivar]->setVal(0.);
    }
 
    // Set class data members
    for ( MFIter mfi(cons_in, TileNoZ()); mfi.isValid(); ++mfi) {
        const auto& box3d = mfi.tilebox();
 
        const auto& lo = lbound(box3d);
        const auto& hi = ubound(box3d);
 
        nlev = box3d.length(2);
        zlo  = lo.z;
        zhi  = hi.z;
 
        // parameters
        accrrc.resize({zlo},  {zhi});
        accrsi.resize({zlo},  {zhi});
        accrsc.resize({zlo},  {zhi});
        coefice.resize({zlo}, {zhi});
        evaps1.resize({zlo},  {zhi});
        evaps2.resize({zlo},  {zhi});
        accrgi.resize({zlo},  {zhi});
        accrgc.resize({zlo},  {zhi});
        evapg1.resize({zlo},  {zhi});
        evapg2.resize({zlo},  {zhi});
        evapr1.resize({zlo},  {zhi});
        evapr2.resize({zlo},  {zhi});
 
        // data (input)
        rho1d.resize({zlo}, {zhi});
        pres1d.resize({zlo}, {zhi});
        tabs1d.resize({zlo}, {zhi});
        gamaz.resize({zlo}, {zhi});
        zmid.resize({zlo}, {zhi});
    }
}

Here is the call graph for this function:

NewtonIterSat()

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real Morrison::NewtonIterSat ( int &  i, int &  j, int &  k, const amrex::Real &  fac_cond, const amrex::Real &  fac_fus, const amrex::Real &  fac_sub, const amrex::Real &  an, const amrex::Real &  bn, const amrex::Array4< amrex::Real > &  tabs_array, const amrex::Array4< amrex::Real > &  pres_array, const amrex::Array4< amrex::Real > &  qv_array, const amrex::Array4< amrex::Real > &  qc_array, const amrex::Array4< amrex::Real > &  qi_array, const amrex::Array4< amrex::Real > &  qn_array, const amrex::Array4< amrex::Real > &  qt_array  )

inlinestatic

    {
        // Solution tolerance
        amrex::Real tol = 1.0e-4;
 
        // Saturation moisture fractions
        amrex::Real omn, domn;
        amrex::Real qsat, dqsat;
        amrex::Real qsatw, dqsatw;
        amrex::Real qsati, dqsati;
 
        // Newton iteration vars
        int  niter;
        amrex::Real fff, dfff, dtabs;
        amrex::Real lstar, dlstar;
        amrex::Real lstarw, lstari;
        amrex::Real delta_qv, delta_qc, delta_qi;
 
        // Initial guess for temperature & pressure
        amrex::Real tabs = tabs_array(i,j,k);
        amrex::Real pres = pres_array(i,j,k);
 
        niter = 0;
        dtabs = 1;
        //==================================================
        // Newton iteration to qv=qsat (cloud phase only)
        //==================================================
        do {
            // Latent heats and their derivatives wrt to T
            lstarw  = fac_cond;
            lstari  = fac_fus;
            domn    = 0.0;
 
            // Saturation moisture fractions
            erf_qsatw(tabs, pres, qsatw);
            erf_qsati(tabs, pres, qsati);
            erf_dtqsatw(tabs, pres, dqsatw);
            erf_dtqsati(tabs, pres, dqsati);
 
            // Cloud ice not permitted (condensation & fusion)
            if(tabs >= tbgmax) {
                omn   = 1.0;
            }
            // Cloud water not permitted (sublimation & fusion)
            else if(tabs <= tbgmin) {
                omn    = 0.0;
                lstarw = fac_sub;
            }
            // Mixed cloud phase (condensation & fusion)
            else {
                omn   = an*tabs-bn;
                domn  = an;
            }
 
            // Linear combination of each component
            qsat   =  omn * qsatw  + (1.0-omn) * qsati;
            dqsat  =  omn * dqsatw + (1.0-omn) * dqsati
                   + domn *  qsatw -     domn  *  qsati;
            lstar  =  omn * lstarw + (1.0-omn) * lstari;
            dlstar = domn * lstarw -     domn  * lstari;
 
            // Function for root finding:
            // 0 = -T_new + T_old + L_eff/C_p * (qv - qsat)
            fff   = -tabs + tabs_array(i,j,k) +  lstar*(qv_array(i,j,k) - qsat);
 
            // Derivative of function (T_new iterated on)
            dfff  = -1.0 + dlstar*(qv_array(i,j,k) - qsat) - lstar*dqsat;
 
            // Update the temperature
            dtabs = -fff/dfff;
            tabs += dtabs;
 
            // Update iteration
            niter = niter+1;
        } while(std::abs(dtabs) > tol && niter < 20);
 
        // Update qsat from last iteration (dq = dq/dt * dt)
        qsat += dqsat*dtabs;
 
        // Changes in each component
        delta_qv = qv_array(i,j,k) - qsat;
        delta_qc = std::max(-qc_array(i,j,k), delta_qv * omn);
        delta_qi = std::max(-qi_array(i,j,k), delta_qv * (1.0-omn));
 
        // Partition the change in non-precipitating q
        qv_array(i,j,k)  = qsat;
        qc_array(i,j,k) += delta_qc;
        qi_array(i,j,k) += delta_qi;
        qn_array(i,j,k)  = qc_array(i,j,k) + qi_array(i,j,k);
        qt_array(i,j,k)  = qv_array(i,j,k) + qn_array(i,j,k);
 
        // Return to temperature
        return tabs;
    }

Here is the call graph for this function:

Precip()

void Morrison::Precip

(

const SolverChoice &

sc

)

Autoconversion (A30), Accretion (A28), Evaporation (A24)

{
    Real powr1 = (3.0 + b_rain) / 4.0;
    Real powr2 = (5.0 + b_rain) / 8.0;
    Real pows1 = (3.0 + b_snow) / 4.0;
    Real pows2 = (5.0 + b_snow) / 8.0;
    Real powg1 = (3.0 + b_grau) / 4.0;
    Real powg2 = (5.0 + b_grau) / 8.0;
 
    auto accrrc_t  = accrrc.table();
    auto accrsc_t  = accrsc.table();
    auto accrsi_t  = accrsi.table();
    auto accrgc_t  = accrgc.table();
    auto accrgi_t  = accrgi.table();
    auto coefice_t = coefice.table();
    auto evapr1_t  = evapr1.table();
    auto evapr2_t  = evapr2.table();
    auto evaps1_t  = evaps1.table();
    auto evaps2_t  = evaps2.table();
    auto evapg1_t  = evapg1.table();
    auto evapg2_t  = evapg2.table();
 
    Real fac_cond = m_fac_cond;
    Real fac_sub  = m_fac_sub;
    Real fac_fus  = m_fac_fus;
    Real rdOcp    = m_rdOcp;
 
    Real eps = std::numeric_limits<Real>::epsilon();
 
    Real dtn = dt;
 
    // get the temperature, density, theta, qt and qp from input
    for ( MFIter mfi(*(mic_fab_vars[MicVar_Morr::tabs]),TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        auto theta_array = mic_fab_vars[MicVar_Morr::theta]->array(mfi);
        auto tabs_array  = mic_fab_vars[MicVar_Morr::tabs]->array(mfi);
        auto pres_array  = mic_fab_vars[MicVar_Morr::pres]->array(mfi);
 
        // Non-precipitating
        auto qv_array    = mic_fab_vars[MicVar_Morr::qv]->array(mfi);
        auto qcl_array   = mic_fab_vars[MicVar_Morr::qcl]->array(mfi);
        auto qci_array   = mic_fab_vars[MicVar_Morr::qci]->array(mfi);
        auto qn_array    = mic_fab_vars[MicVar_Morr::qn]->array(mfi);
        auto qt_array    = mic_fab_vars[MicVar_Morr::qt]->array(mfi);
 
        // Precipitating
        auto qpr_array   = mic_fab_vars[MicVar_Morr::qpr]->array(mfi);
        auto qps_array   = mic_fab_vars[MicVar_Morr::qps]->array(mfi);
        auto qpg_array   = mic_fab_vars[MicVar_Morr::qpg]->array(mfi);
        auto qp_array    = mic_fab_vars[MicVar_Morr::qp]->array(mfi);
 
        const auto& box3d = mfi.tilebox();
 
        ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            //------- Autoconversion/accretion
            Real omn, omp, omg;
            Real qsat, qsatw, qsati;
 
            Real qcc, qii, qpr, qps, qpg;
            Real dprc, dpsc, dpgc;
            Real dpsi, dpgi;
 
            Real dqc, dqca, dqi, dqia, dqp;
            Real dqpr, dqps, dqpg;
 
            Real auto_r, autos;
            Real accrcr, accrcs, accris, accrcg, accrig;
 
            // Work to be done for autoc/accr or evap
            if (qn_array(i,j,k)+qp_array(i,j,k) > 0.0) {
                omn = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tbgmin)*a_bg));
                omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
                omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr));
 
                qcc = qcl_array(i,j,k);
                qii = qci_array(i,j,k);
 
                qpr = qpr_array(i,j,k);
                qps = qps_array(i,j,k);
                qpg = qpg_array(i,j,k);
 
                //==================================================
                // Autoconversion (A30/A31) and accretion (A27)
                //==================================================
                if (qn_array(i,j,k) > 0.0) {
                    accrcr = 0.0;
                    accrcs = 0.0;
                    accris = 0.0;
                    accrcg = 0.0;
                    accrig = 0.0;
 
                    if (qcc > qcw0) {
                        auto_r = alphaelq;
                    } else {
                        auto_r = 0.0;
                    }
 
                    if (qii > qci0) {
                        autos = betaelq*coefice_t(k);
                    } else {
                        autos = 0.0;
                    }
 
                    if (omp > 0.001) {
                        accrcr = accrrc_t(k);
                    }
 
                    if (omp < 0.999 && omg < 0.999) {
                        accrcs = accrsc_t(k);
                        accris = accrsi_t(k);
                    }
 
                    if (omp < 0.999 && omg > 0.001) {
                        accrcg = accrgc_t(k);
                        accrig = accrgi_t(k);
                    }
 
                    // Autoconversion & accretion (sink for cloud comps)
                    dqca = dtn * auto_r  * (qcc-qcw0);
                    dprc = dtn * accrcr * qcc * std::pow(qpr, powr1);
                    dpsc = dtn * accrcs * qcc * std::pow(qps, pows1);
                    dpgc = dtn * accrcg * qcc * std::pow(qpg, powg1);
 
                    dqia = dtn * autos  * (qii-qci0);
                    dpsi = dtn * accris * qii * std::pow(qps, pows1);
                    dpgi = dtn * accrig * qii * std::pow(qpg, powg1);
 
                    // Rescale sinks to avoid negative cloud fractions
                    dqc  = dqca + dprc + dpsc + dpgc;
                    dqi  = dqia + dpsi + dpgi;
                    Real scalec = std::min(qcl_array(i,j,k),dqc) / (dqc + eps);
                    Real scalei = std::min(qci_array(i,j,k),dqi) / (dqi + eps);
                    dqca *= scalec; dprc *= scalec; dpsc *= scalec; dpgc *= scalec;
                    dqia *= scalei; dpsi *= scalei; dpgi *= scalei;
                    dqc   = dqca + dprc + dpsc + dpgc;
                    dqi   = dqia + dpsi + dpgi;
 
                    // NOTE: Autoconversion of cloud water and ice are sources
                    //       to qp, while accretion is a source to an individual
                    //       precipitating component (e.g., qpr/qps/qpg). So we
                    //       only split autoconversion with omega. The omega
                    //       splitting does imply a latent heat source.
 
                    // Partition formed precip componentss
                    dqpr = (dqca + dqia) * omp + dprc;
                    dqps = (dqca + dqia) * (1.0 - omp) * (1.0 - omg) + dpsc + dpsi;
                    dqpg = (dqca + dqia) * (1.0 - omp) * omg         + dpgc + dpgi;
 
                    // Update the primitive state variables
                    qcl_array(i,j,k) -= dqc;
                    qci_array(i,j,k) -= dqi;
                    qpr_array(i,j,k) += dqpr;
                    qps_array(i,j,k) += dqps;
                    qpg_array(i,j,k) += dqpg;
 
                    // Update the primitive derived vars
                    qn_array(i,j,k) = qcl_array(i,j,k) + qci_array(i,j,k);
                    qt_array(i,j,k) =  qv_array(i,j,k) +  qn_array(i,j,k);
                    qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
 
                    // Update temperature
                    tabs_array(i,j,k) += fac_fus * ( dqca * (1.0 - omp) - dqia * omp );
 
                    // Update theta
                    theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
                }
 
                //==================================================
                // Evaporation (A24)
                //==================================================
                erf_qsatw(tabs_array(i,j,k),pres_array(i,j,k),qsatw);
                erf_qsati(tabs_array(i,j,k),pres_array(i,j,k),qsati);
                qsat = qsatw * omn + qsati * (1.0-omn);
                if((qp_array(i,j,k) > 0.0) && (qv_array(i,j,k) < qsat)) {
 
                    dqpr = evapr1_t(k)*sqrt(qpr) + evapr2_t(k)*pow(qpr,powr2);
                    dqps = evaps1_t(k)*sqrt(qps) + evaps2_t(k)*pow(qps,pows2);
                    dqpg = evapg1_t(k)*sqrt(qpg) + evapg2_t(k)*pow(qpg,powg2);
 
                    // NOTE: This is always a sink for precipitating comps
                    //       since qv<qsat and thus (1 - qv/qsat)>0. If we are
                    //       in a super-saturated state (qv>qsat) the Newton
                    //       iterations in Cloud() will have handled condensation.
                    dqpr *= dtn * (1.0 - qv_array(i,j,k)/qsat);
                    dqps *= dtn * (1.0 - qv_array(i,j,k)/qsat);
                    dqpg *= dtn * (1.0 - qv_array(i,j,k)/qsat);
 
                    // Limit to avoid negative moisture fractions
                    dqpr = std::min(qpr_array(i,j,k),dqpr);
                    dqps = std::min(qps_array(i,j,k),dqps);
                    dqpg = std::min(qpg_array(i,j,k),dqpg);
                    dqp  = dqpr + dqps + dqpg;
 
                    // Update the primitive state variables
                     qv_array(i,j,k) += dqp;
                    qpr_array(i,j,k) -= dqpr;
                    qps_array(i,j,k) -= dqps;
                    qpg_array(i,j,k) -= dqpg;
 
                    // Update the primitive derived vars
                    qt_array(i,j,k) =  qv_array(i,j,k) +  qn_array(i,j,k);
                    qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
 
                    // Update temperature
                    tabs_array(i,j,k) -= fac_cond * dqpr + fac_sub * (dqps + dqpg);
 
                    // Update theta
                    theta_array(i,j,k) = getThgivenPandT(tabs_array(i,j,k), 100.0*pres_array(i,j,k), rdOcp);
                }
            }
        });
    }
}

Referenced by Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

PrecipFall()

void Morrison::PrecipFall

(

const SolverChoice &

sc

)

Precipitation fluxes P_{r/s/g} (A19)

Code modified from MorrisonXX, the C++ version of the Morrison code.

Parameters

[in]

hydro_type

Type selection for the precipitation advection hydrodynamics scheme (0-3)

{
    Real rho_0 = 1.29;
 
    Real gamr3 = erf_gammafff(4.0+b_rain);
    Real gams3 = erf_gammafff(4.0+b_snow);
    Real gamg3 = erf_gammafff(4.0+b_grau);
 
    Real vrain = (a_rain*gamr3/6.0)*pow((PI*rhor*nzeror),-crain);
    Real vsnow = (a_snow*gams3/6.0)*pow((PI*rhos*nzeros),-csnow);
    Real vgrau = (a_grau*gamg3/6.0)*pow((PI*rhog*nzerog),-cgrau);
 
    auto dz   = m_geom.CellSize(2);
    Real dtn  = dt;
    Real coef = dtn/dz;
 
    auto domain = m_geom.Domain();
    int k_lo = domain.smallEnd(2);
    int k_hi = domain.bigEnd(2);
 
    auto qpr   = mic_fab_vars[MicVar_Morr::qpr];
    auto qps   = mic_fab_vars[MicVar_Morr::qps];
    auto qpg   = mic_fab_vars[MicVar_Morr::qpg];
    auto qp    = mic_fab_vars[MicVar_Morr::qp];
    auto rho   = mic_fab_vars[MicVar_Morr::rho];
    auto tabs  = mic_fab_vars[MicVar_Morr::tabs];
    auto theta = mic_fab_vars[MicVar_Morr::theta];
    auto rain_accum = mic_fab_vars[MicVar_Morr::rain_accum];
    auto snow_accum = mic_fab_vars[MicVar_Morr::snow_accum];
    auto graup_accum = mic_fab_vars[MicVar_Morr::graup_accum];
 
    auto ba    = tabs->boxArray();
    auto dm    = tabs->DistributionMap();
    auto ngrow = tabs->nGrowVect();
 
    MultiFab fz;
    fz.define(convert(ba, IntVect(0,0,1)), dm, 1, ngrow);
 
    //  Precompute the vertical fluxes for CFL constraint
    for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
        auto qp_array   = qp->array(mfi);
        auto rho_array  = rho->array(mfi);
        auto tabs_array = tabs->array(mfi);
        auto fz_array   = fz.array(mfi);
 
        const auto& box3d = mfi.tilebox();
 
        ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
        {
            Real rho_avg, tab_avg, qp_avg;
            if (k==k_lo) {
                rho_avg =  rho_array(i,j,k);
                tab_avg = tabs_array(i,j,k);
                 qp_avg =   qp_array(i,j,k);
            } else if (k==k_hi+1) {
                rho_avg =  rho_array(i,j,k-1);
                tab_avg = tabs_array(i,j,k-1);
                 qp_avg =   qp_array(i,j,k-1);
            } else {
                rho_avg = 0.5*( rho_array(i,j,k-1) +  rho_array(i,j,k));
                tab_avg = 0.5*(tabs_array(i,j,k-1) + tabs_array(i,j,k));
                 qp_avg = 0.5*(  qp_array(i,j,k-1) +   qp_array(i,j,k));
            }
 
            Real Pprecip = 0.0;
            if(qp_avg > qp_threshold) {
                Real omp, omg;
                omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr));
                omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr));
                Real qrr = omp*qp_avg;
                Real qss = (1.0-omp)*(1.0-omg)*qp_avg;
                Real qgg = (1.0-omp)*(omg)*qp_avg;
                Pprecip = omp*vrain*std::pow(rho_avg*qrr,1.0+crain)
                        + (1.0-omp)*( (1.0-omg)*vsnow*std::pow(rho_avg*qss,1.0+csnow)
                                    +      omg *vgrau*std::pow(rho_avg*qgg,1.0+cgrau) );
            }
 
            // NOTE: Fz is the sedimentation flux from the advective operator.
            //       In the terrain-following coordinate system, the z-deriv in
            //       the divergence uses the normal velocity (Omega). However,
            //       there are no u/v components to the sedimentation velocity.
            //       Therefore, we simply end up with a division by detJ when
            //       evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
            fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg);
        });
    }
 
    // Compute number of substeps from maximum terminal velocity
    Real wt_max;
    int n_substep;
    auto const& ma_fz_arr = fz.const_arrays();
    GpuTuple<Real> max = ParReduce(TypeList<ReduceOpMax>{},
                                   TypeList<Real>{},
                                   fz, IntVect(0),
                         [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
                         -> GpuTuple<Real>
                         {
                             return { ma_fz_arr[box_no](i,j,k) };
                         });
    wt_max = get<0>(max) + std::numeric_limits<Real>::epsilon();
    n_substep = int( std::ceil(wt_max * coef / CFL_MAX) );
    AMREX_ALWAYS_ASSERT(n_substep >= 1);
    coef /= Real(n_substep);
    dtn  /= Real(n_substep);
 
    // Substep the vertical advection
    for (int nsub(0); nsub<n_substep; ++nsub) {
        for (MFIter mfi(*qp, TileNoZ()); mfi.isValid(); ++mfi) {
            auto qpr_array    = qpr->array(mfi);
            auto qps_array    = qps->array(mfi);
            auto qpg_array    = qpg->array(mfi);
            auto qp_array     = qp->array(mfi);
            auto rho_array    = rho->array(mfi);
            auto tabs_array   = tabs->array(mfi);
            auto fz_array     = fz.array(mfi);
 
            auto rain_accum_array = rain_accum->array(mfi);
            auto snow_accum_array = snow_accum->array(mfi);
            auto graup_accum_array = graup_accum->array(mfi);
 
            const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4<const Real>{};
 
            const auto& tbx = mfi.tilebox();
            const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0));
 
            // Update vertical flux every substep
            ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                Real rho_avg, tab_avg, qp_avg;
                if (k==k_lo) {
                    rho_avg =  rho_array(i,j,k);
                    tab_avg = tabs_array(i,j,k);
                     qp_avg =   qp_array(i,j,k);
                } else if (k==k_hi+1) {
                    rho_avg =  rho_array(i,j,k-1);
                    tab_avg = tabs_array(i,j,k-1);
                     qp_avg =   qp_array(i,j,k-1);
                } else {
                    rho_avg = 0.5*( rho_array(i,j,k-1) +  rho_array(i,j,k));
                    tab_avg = 0.5*(tabs_array(i,j,k-1) + tabs_array(i,j,k));
                     qp_avg = 0.5*(  qp_array(i,j,k-1) +   qp_array(i,j,k));
                }
 
                Real Pprecip = 0.0;
                if(qp_avg > qp_threshold) {
                    Real omp, omg;
                    omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr));
                    omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr));
                    Real qrr = omp*qp_avg;
                    Real qss = (1.0-omp)*(1.0-omg)*qp_avg;
                    Real qgg = (1.0-omp)*(omg)*qp_avg;
                    Pprecip = omp*vrain*std::pow(rho_avg*qrr,1.0+crain)
                            + (1.0-omp)*( (1.0-omg)*vsnow*std::pow(rho_avg*qss,1.0+csnow)
                                        +      omg *vgrau*std::pow(rho_avg*qgg,1.0+cgrau) );
                }
 
                // NOTE: Fz is the sedimentation flux from the advective operator.
                //       In the terrain-following coordinate system, the z-deriv in
                //       the divergence uses the normal velocity (Omega). However,
                //       there are no u/v components to the sedimentation velocity.
                //       Therefore, we simply end up with a division by detJ when
                //       evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv.
                fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg);
 
                if(k==k_lo){
                    Real omp, omg;
                    omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr));
                    omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr));
                    rain_accum_array(i,j,k)  = rain_accum_array(i,j,k)  + rho_avg*(omp*qp_avg)*vrain*dtn/rhor*1000.0; // Divide by rho_water and convert to mm
                    snow_accum_array(i,j,k)  = snow_accum_array(i,j,k)  + rho_avg*(1.0-omp)*(1.0-omg)*qp_avg*vrain*dtn/rhos*1000.0; // Divide by rho_snow and convert to mm
                    graup_accum_array(i,j,k) = graup_accum_array(i,j,k) + rho_avg*(1.0-omp)*(omg)*qp_avg*vrain*dtn/rhog*1000.0; // Divide by rho_graupel and convert to mm
                }
            });
 
            // Update precip every substep
            ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
            {
                // Jacobian determinant
                Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0;
 
                //==================================================
                // Precipitating sedimentation (A19)
                //==================================================
                Real dqp = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef;
                Real omp, omg;
                omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr));
                omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr));
 
                qpr_array(i,j,k) = std::max(0.0, qpr_array(i,j,k) + dqp*omp);
                qps_array(i,j,k) = std::max(0.0, qps_array(i,j,k) + dqp*(1.0-omp)*(1.0-omg));
                qpg_array(i,j,k) = std::max(0.0, qpg_array(i,j,k) + dqp*(1.0-omp)*omg);
                 qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k);
 
                 // NOTE: Sedimentation does not affect the potential temperature,
                 //       but it does affect the liquid/ice static energy.
                 //       No source to Theta occurs here.
            });
        } // mfi
    } // nsub
}

Referenced by Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

Qmoist_Ptr()

amrex::MultiFab* Morrison::Qmoist_Ptr ( const int &  varIdx )

inlineoverridevirtual

Reimplemented from NullMoist.

    {
        AMREX_ALWAYS_ASSERT(varIdx < m_qmoist_size);
        return mic_fab_vars[MicVarMap[varIdx]].get();
    }

Qmoist_Restart_Vars()

void Morrison::Qmoist_Restart_Vars ( const SolverChoice &  , std::vector< int > &  a_idx, std::vector< std::string > &  a_names  ) const

inlineoverridevirtual

Reimplemented from NullMoist.

    {
        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");
        a_idx.push_back(1); a_names.push_back("SnowAccum");
        a_idx.push_back(2); a_names.push_back("GraupAccum");
    }

Qmoist_Size()

int Morrison::Qmoist_Size ( )

inlineoverridevirtual

Reimplemented from NullMoist.

{ return Morrison::m_qmoist_size; }

Qstate_Size()

int Morrison::Qstate_Size ( )

inlineoverridevirtual

Reimplemented from NullMoist.

{ return Morrison::m_qstate_size; }

Update_Micro_Vars()

void Morrison::Update_Micro_Vars ( amrex::MultiFab &  cons_in )

inlineoverridevirtual

Reimplemented from NullMoist.

    {
        this->Copy_State_to_Micro(cons_in);
        this->Compute_Coefficients();
    }

Here is the call graph for this function:

Update_State_Vars()

void Morrison::Update_State_Vars ( amrex::MultiFab &  cons_in )

inlineoverridevirtual

Reimplemented from NullMoist.

    {
        this->Copy_Micro_to_State(cons_in);
    }

Here is the call graph for this function:

Member Data Documentation

accrgc

amrex::TableData<amrex::Real, 1> Morrison::accrgc

private

accrgi

amrex::TableData<amrex::Real, 1> Morrison::accrgi

private

accrrc

amrex::TableData<amrex::Real, 1> Morrison::accrrc

private

accrsc

amrex::TableData<amrex::Real, 1> Morrison::accrsc

private

accrsi

amrex::TableData<amrex::Real, 1> Morrison::accrsi

private

CFL_MAX

constexpr amrex::Real Morrison::CFL_MAX = 0.5

staticconstexprprivate

coefice

amrex::TableData<amrex::Real, 1> Morrison::coefice

private

dt

amrex::Real Morrison::dt

private

Referenced by Advance().

evapg1

amrex::TableData<amrex::Real, 1> Morrison::evapg1

private

evapg2

amrex::TableData<amrex::Real, 1> Morrison::evapg2

private

evapr1

amrex::TableData<amrex::Real, 1> Morrison::evapr1

private

evapr2

amrex::TableData<amrex::Real, 1> Morrison::evapr2

private

evaps1

amrex::TableData<amrex::Real, 1> Morrison::evaps1

private

evaps2

amrex::TableData<amrex::Real, 1> Morrison::evaps2

private

gamaz

amrex::TableData<amrex::Real, 1> Morrison::gamaz

private

m_axis

int Morrison::m_axis

private

Referenced by Define().

m_detJ_cc

amrex::MultiFab* Morrison::m_detJ_cc

private

m_fac_cond

amrex::Real Morrison::m_fac_cond

private

Referenced by Define().

m_fac_fus

amrex::Real Morrison::m_fac_fus

private

Referenced by Define().

m_fac_sub

amrex::Real Morrison::m_fac_sub

private

Referenced by Define().

m_geom

amrex::Geometry Morrison::m_geom

private

m_gOcp

amrex::Real Morrison::m_gOcp

private

Referenced by Define().

m_gtoe

amrex::BoxArray Morrison::m_gtoe

private

m_qmoist_size

int Morrison::m_qmoist_size = 3

private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

m_qstate_size

int Morrison::m_qstate_size = 6

private

Referenced by Qstate_Size().

m_rdOcp

amrex::Real Morrison::m_rdOcp

private

Referenced by Define().

m_z_phys_nd

amrex::MultiFab* Morrison::m_z_phys_nd

private

mic_fab_vars

amrex::Array<FabPtr, MicVar_Morr::NumVars> Morrison::mic_fab_vars

private

Referenced by Qmoist_Ptr().

MicVarMap

amrex::Vector<int> Morrison::MicVarMap

private

Referenced by Qmoist_Ptr().

nlev

int Morrison::nlev

private

pres1d

amrex::TableData<amrex::Real, 1> Morrison::pres1d

private

qn1d

amrex::TableData<amrex::Real, 1> Morrison::qn1d

private

qt1d

amrex::TableData<amrex::Real, 1> Morrison::qt1d

private

qv1d

amrex::TableData<amrex::Real, 1> Morrison::qv1d

private

rho1d

amrex::TableData<amrex::Real, 1> Morrison::rho1d

private

tabs1d

amrex::TableData<amrex::Real, 1> Morrison::tabs1d

private

zhi

int Morrison::zhi

private

zlo

int Morrison::zlo

private

zmid

amrex::TableData<amrex::Real, 1> Morrison::zmid

private

---

The documentation for this class was generated from the following files:

• Source/Microphysics/Morrison/ERF_Morrison.H
• Source/Microphysics/Morrison/ERF_InitMorrison.cpp
• Source/Microphysics/Morrison/ERF_Morrison_Cloud.cpp
• Source/Microphysics/Morrison/ERF_Morrison_IceFall.cpp
• Source/Microphysics/Morrison/ERF_Morrison_Precip.cpp
• Source/Microphysics/Morrison/ERF_Morrison_PrecipFall.cpp
• Source/Microphysics/Morrison/ERF_UpdateMorrison.cpp