#include <ERF_SAM.H>

Inheritance diagram for SAM:

Collaboration diagram for SAM:

[legend]

Public Member Functions
	SAM ()

virtual	~SAM ()=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	Set_dzmin (const amrex::Real dz_min) override

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

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, const amrex::MultiFab &) override

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

const SedimentationFlux *	Sed_Flux () const override
	The sedimentation face fluxes of the most recent Advance, or null for a scheme that does not sediment (see ERF_SedimentationFlux.H) More...

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

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

void	Compute_Coefficients (bool use_anelastic_reference_pressure=false)

SAMCoefficientRow	CoefficientRowAt (int k) const

int	Qmoist_Size () override

int	Qstate_Moist_Size () override

int	Qstate_Moist_NumConc_Size () override

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

SurfacePrecipAccumulationSources	Get_Surface_Precip_Accumulation_Ptrs (const int &) const override

Public Member Functions inherited from NullMoist
	NullMoist ()

virtual	~NullMoist ()=default

virtual int	Qstate_NonMoist_Size ()

virtual void	GetPlotVarNames (amrex::Vector< std::string > &a_vec) const

virtual void	GetPlotVar (const std::string &, amrex::MultiFab &) const

virtual void	GetPlotVar (const std::string &a_name, amrex::MultiFab &a_mf, const int) const

virtual void	SetCurrentLevel (const int &lev)

virtual void	InitLevel (const int, const amrex::MultiFab &)

virtual int	getDiagnosticsInterval () const

virtual void	Set_Lmask (amrex::iMultiFab *)
	Import ERF's land/water mask. Only schemes whose physics branches on land vs water need to override this. More...

virtual void	Set_Sed_Inflow (const amrex::MultiFab *inflow)
	The inflow through the top of a fine patch for the next Advance, or null (see ERF_SedimentationFlux.H); the caller keeps it alive for that Advance. More...

virtual void	Set_RealWidth (const int)

Static Public Member Functions
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real	NewtonIterSat (int &i, int &j, int &k, const int &SAM_moisture_type, const amrex::Real &fac_cond, const amrex::Real &fac_fus, const amrex::Real &, 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 = 4

int	n_qstate_moist_size = 6

int	n_qstate_moist_numconc_size = 0

amrex::Vector< int >	MicVarMap

amrex::Geometry	m_geom

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_rdOcp {RdoCp}

bool	m_do_cond

MoistureType	m_moisture_type = MoistureType::None

amrex::Real	m_dzmin {0.0}

amrex::MultiFab *	m_z_phys_nd

amrex::MultiFab *	m_detJ_cc

amrex::Array< FabPtr, MicVar::NumVars >	mic_fab_vars

SedimentationFlux	m_sed

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

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

Additional Inherited Members
Protected Member Functions inherited from NullMoist
void	set_anelastic_reference_pressure_mode (const SolverChoice &sc)

void	assert_base_state_available (const amrex::MultiFab *base_state) const

Protected Attributes inherited from NullMoist
const amrex::MultiFab *	m_sed_inflow = nullptr

int	m_level {0}

bool	m_use_anelastic_reference_pressure {false}

Member Typedef Documentation

◆ FabPtr

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

private

Constructor & Destructor Documentation

◆ SAM()

SAM::SAM ( )

inline

63 {}

◆ ~SAM()

virtual SAM::~SAM ( )

virtualdefault

Member Function Documentation

◆ Advance()

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

inlineoverridevirtual

Reimplemented from NullMoist.

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

Here is the call graph for this function:

◆ Cloud()

void SAM::Cloud ( const SolverChoice & sc )

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

 {
     if (!m_do_cond) { return; }
  
     constexpr Real an = one/(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;
  
     int SAM_moisture_type = 1;
     if (sam_is_no_ice(sc.moisture_type)) {
         SAM_moisture_type = 2;
     }
  
     for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
         auto  qt_array = mic_fab_vars[MicVar::qt]->array(mfi);
         auto  qn_array = mic_fab_vars[MicVar::qn]->array(mfi);
         auto  qv_array = mic_fab_vars[MicVar::qv]->array(mfi);
         auto qcl_array = mic_fab_vars[MicVar::qcl]->array(mfi);
         auto qci_array = mic_fab_vars[MicVar::qci]->array(mfi);
  
         auto  tabs_array = mic_fab_vars[MicVar::tabs]->array(mfi);
         auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi);
         auto  pres_array = mic_fab_vars[MicVar::pres]->array(mfi);
  
         auto tbx = mfi.tilebox();
  
         ParallelFor(tbx, [=] 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.
  
             const SAMCloudPhaseChange phase_change =
                 sam_partition_cloud_phase(SAM_moisture_type,
                                           tabs_array(i,j,k),
                                           qn_array(i,j,k),
                                           qcl_array(i,j,k),
                                           qci_array(i,j,k),
                                           fac_cond, fac_fus,
                                           an, bn);
             omn = phase_change.omn;
             delta_qc = phase_change.delta_qc;
             delta_qi = phase_change.delta_qi;
             qcl_array(i,j,k) = phase_change.qcl;
             qci_array(i,j,k) = phase_change.qci;
             tabs_array(i,j,k) = phase_change.tabs;
             // Cloud adjustment updates tabs under the held-pressure SAM source
             // convention, then refreshes theta from the same stored pressure.
             theta_array(i,j,k) = sam_theta_from_stored_mbar_converted_to_pa(tabs_array(i,j,k),
                                                                              pres_array(i,j,k), rdOcp);
  
             // 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 = sam_mixed_qsat(omn, qsatw, 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   , SAM_moisture_type   ,
                                                   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) = sam_theta_from_stored_mbar_converted_to_pa(tabs_array(i,j,k),
                                                                                  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)  = zero;
                 qci_array(i,j,k)  = zero;
                  qn_array(i,j,k)  = zero;
                  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) = sam_theta_from_stored_mbar_converted_to_pa(tabs_array(i,j,k),
                                                                                  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 = sam_mixed_qsat(omn, qsatw, qsati);
                 if (qt_array(i,j,k) > qsat) {
  
                     // Update temperature
                     tabs_array(i,j,k) = NewtonIterSat(i, j, k   , SAM_moisture_type   ,
                                                       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) = sam_theta_from_stored_mbar_converted_to_pa(tabs_array(i,j,k),
                                                                                      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:

◆ CoefficientRowAt()

SAMCoefficientRow SAM::CoefficientRowAt ( int k ) const

 {
     AMREX_ALWAYS_ASSERT(k >= zlo && k <= zhi);
  
     SAMCoefficientRow row{};
     row.accrrc = copy_table_value_to_host(accrrc, zlo, zhi, k);
     row.accrsi = copy_table_value_to_host(accrsi, zlo, zhi, k);
     row.accrsc = copy_table_value_to_host(accrsc, zlo, zhi, k);
     row.coefice = copy_table_value_to_host(coefice, zlo, zhi, k);
     row.evaps1 = copy_table_value_to_host(evaps1, zlo, zhi, k);
     row.evaps2 = copy_table_value_to_host(evaps2, zlo, zhi, k);
     row.accrgi = copy_table_value_to_host(accrgi, zlo, zhi, k);
     row.accrgc = copy_table_value_to_host(accrgc, zlo, zhi, k);
     row.evapg1 = copy_table_value_to_host(evapg1, zlo, zhi, k);
     row.evapg2 = copy_table_value_to_host(evapg2, zlo, zhi, k);
     row.evapr1 = copy_table_value_to_host(evapr1, zlo, zhi, k);
     row.evapr2 = copy_table_value_to_host(evapr2, zlo, zhi, k);
     return row;
 }

Here is the call graph for this function:

◆ Compute_Coefficients()

void SAM::Compute_Coefficients ( bool use_anelastic_reference_pressure = false )

 {
     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();
  
     Real gam3  = erf_gammafff(three             );
     Real gamr1 = erf_gammafff(three+b_rain      );
     Real gamr2 = erf_gammafff((Real(5.0)+b_rain)/two);
     Real gams1 = erf_gammafff(three+b_snow      );
     Real gams2 = erf_gammafff((Real(5.0)+b_snow)/two);
     Real gamg1 = erf_gammafff(three+b_grau      );
     Real gamg2 = erf_gammafff((Real(5.0)+b_grau)/two);
  
     // calculate the plane average variables
     PlaneAverage rho_ave(mic_fab_vars[MicVar::rho].get(), m_geom, m_axis);
     PlaneAverage theta_ave(mic_fab_vars[MicVar::theta].get(), m_geom, m_axis);
     PlaneAverage qv_ave(mic_fab_vars[MicVar::qv].get(), m_geom, m_axis);
     PlaneAverage pres_ave(mic_fab_vars[MicVar::pres].get(), m_geom, m_axis);
     PlaneAverage tabs_ave(mic_fab_vars[MicVar::tabs].get(), m_geom, m_axis);
     rho_ave.compute_averages(ZDir(), rho_ave.field());
     if (use_anelastic_reference_pressure) {
         pres_ave.compute_averages(ZDir(), pres_ave.field());
         tabs_ave.compute_averages(ZDir(), tabs_ave.field());
     } else {
         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);
     Gpu::HostVector<Real> pres_h(ncell), tabs_h(ncell);
     rho_ave.line_average(0, rho_h);
     if (use_anelastic_reference_pressure) {
         pres_ave.line_average(0, pres_h);
         tabs_ave.line_average(0, tabs_h);
     } else {
         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::DeviceVector<Real> pres_d(ncell), tabs_d(ncell);
     Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin());
     if (use_anelastic_reference_pressure) {
         Gpu::copyAsync(Gpu::hostToDevice, pres_h.begin(), pres_h.end(), pres_d.begin());
         Gpu::copyAsync(Gpu::hostToDevice, tabs_h.begin(), tabs_h.end(), tabs_d.begin());
     } else {
         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* pres_dptr  = pres_d.data();
     Real* tabs_dptr  = tabs_d.data();
  
     ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
     {
         rho1d_t(k)    = rho_dptr[k];
         if (use_anelastic_reference_pressure) {
             pres1d_t(k) = pres_dptr[k];
             tabs1d_t(k) = amrex::min(tabs_dptr[k], Real(273.16));
         } else {
             Real RhoTheta = rho_dptr[k]*theta_dptr[k];
             Real pressure = getPgivenRTh(RhoTheta, qv_dptr[k]);
             pres1d_t(k) = sam_pa_to_mbar(pressure);
             // NOTE: Limit the temperature to the melting point of ice to avoid
             // a divide by 0 condition when computing cold coefficients.
             tabs1d_t(k) = amrex::min(getTgivenRandRTh(rho_dptr[k], RhoTheta, qv_dptr[k]),
                                      Real(273.16));
         }
     });
  
     if(round(gam3) != 2) {
         amrex::Abort("cannot compute gamma-function in Microphysics::Init");
     }
  
     // Populate all the coefficients
     ParallelFor(nlev, [=] AMREX_GPU_DEVICE (int k) noexcept
     {
         const SAMCoefficientRow row =
             sam_compute_coefficient_row(rho1d_t(k), tabs1d_t(k),
                                         gamr1, gamr2, gams1, gams2, gamg1, gamg2);
         accrrc_t(k) = row.accrrc;
         accrsi_t(k) = row.accrsi;
         accrsc_t(k) = row.accrsc;
         coefice_t(k) = row.coefice;
         evaps1_t(k) = row.evaps1;
         evaps2_t(k) = row.evaps2;
         accrgi_t(k) = row.accrgi;
         accrgc_t(k) = row.accrgc;
         evapg1_t(k) = row.evapg1;
         evapg2_t(k) = row.evapg2;
         evapr1_t(k) = row.evapr1;
         evapr2_t(k) = row.evapr2;
     });
 }

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 SAM::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::rho]->array(mfi);
         auto theta_arr  = mic_fab_vars[MicVar::theta]->array(mfi);
  
         auto qv_arr     = mic_fab_vars[MicVar::qv]->array(mfi);
         auto qc_arr     = mic_fab_vars[MicVar::qcl]->array(mfi);
         auto qi_arr     = mic_fab_vars[MicVar::qci]->array(mfi);
  
         auto qpr_arr     = mic_fab_vars[MicVar::qpr]->array(mfi);
         auto qps_arr     = mic_fab_vars[MicVar::qps]->array(mfi);
         auto qpg_arr     = mic_fab_vars[MicVar::qpg]->array(mfi);
  
         // get potential total density, temperature, qt, qp
         ParallelFor( box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
         {
             SAMPrimitiveCell primitive{};
             primitive.rho = rho_arr(i,j,k);
             primitive.theta = theta_arr(i,j,k);
             primitive.qv = qv_arr(i,j,k);
             primitive.qcl = qc_arr(i,j,k);
             primitive.qci = qi_arr(i,j,k);
             primitive.qpr = qpr_arr(i,j,k);
             primitive.qps = qps_arr(i,j,k);
             primitive.qpg = qpg_arr(i,j,k);
             sam_primitive_to_cons(primitive, states_arr, i, j, k);
         });
     }
  
     // Fill interior ghost cells and periodic boundaries
     cons.FillBoundary(m_geom.periodicity());
 }

Referenced by Update_State_Vars().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Copy_State_to_Micro() [1/2]

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

overridevirtual

Reimplemented from NullMoist.

Referenced by Update_Micro_Vars().

Here is the caller graph for this function:

◆ Copy_State_to_Micro() [2/2]

void SAM::Copy_State_to_Micro	(	const amrex::MultiFab &	cons_in,
		const amrex::MultiFab *	base_state
	)

◆ Define()

void SAM::Define ( SolverChoice & sc )

inlineoverridevirtual

Reimplemented from NullMoist.

     {
         m_moisture_type = sc.moisture_type;
         m_fac_cond = lcond / sc.c_p;
         m_fac_fus  = lfus / sc.c_p;
         m_fac_sub  = lsub / sc.c_p;
         m_axis     = sc.ave_plane;
         m_rdOcp    = sc.rdOcp;
         m_do_cond  = (!sc.uses_shoc_family());
         set_anelastic_reference_pressure_mode(sc);
     }

Here is the call graph for this function:

◆ Get_Surface_Precip_Accumulation_Ptrs()

SurfacePrecipAccumulationSources SAM::Get_Surface_Precip_Accumulation_Ptrs ( const int & ) const

inlineoverridevirtual

Reimplemented from NullMoist.

     {
         if (m_moisture_type == MoistureType::SAM_NoPrecip_NoIce) {
             return {};
         }
  
         SurfacePrecipAccumulationSources sources;
         // No total is offered: the land and fire couplings take the total where a
         // scheme offers one and the rain otherwise, and SAM has always given them rain
         sources.rain = {mic_fab_vars[MicVar::rain_accum].get(), rhor / amrex::Real(1000.0)};
  
         if (m_moisture_type == MoistureType::SAM) {
             sources.snow    = {mic_fab_vars[MicVar::snow_accum].get(), rhos / amrex::Real(1000.0)};
             sources.graupel = {mic_fab_vars[MicVar::graupel_accum].get(), rhog / amrex::Real(1000.0)};
         }
  
         return sources;
     }

◆ IceFall()

void SAM::IceFall ( const SolverChoice & sc )

Sedimentation of cloud ice (A32)

                                          {
  
     if (sam_is_no_ice(sc.moisture_type))
       return;
  
     Real dtn  = dt;
     // NOTE: coef carries only the reference vertical spacing. The physical
     //       spacing is supplied by the inverse Jacobian in the tendency below,
     //       since detJ = dz_phys/CellSize(2). Dividing by m_dzmin here as well
     //       would apply the vertical metric twice. m_dzmin is still the correct
     //       length scale for the substep (CFL) count, which must bound the
     //       thinnest cell in the domain.
     Real coef = dtn * m_geom.InvCellSize(2);
  
     auto domain = m_geom.Domain();
     int k_lo = domain.smallEnd(2);
     int k_hi = domain.bigEnd(2);
  
     auto qcl   = mic_fab_vars[MicVar::qcl];
     auto qci   = mic_fab_vars[MicVar::qci];
     auto qn    = mic_fab_vars[MicVar::qn];
     auto qt    = mic_fab_vars[MicVar::qt];
     auto rho   = mic_fab_vars[MicVar::rho];
     auto tabs  = mic_fab_vars[MicVar::tabs];
  
     // Cloud ice that falls out of the bottom reaches the ground, so it counts in
     // the precipitation total (though in none of the rain, snow or graupel ones)
     auto precip_accum = mic_fab_vars[MicVar::precip_accum];
  
     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.);
  
     // The face state reads qci at k-1, a ghost cell whenever a box boundary is interior in
     // z, and Cloud has just updated qci in valid cells only: refresh before the fluxes that
     // size the substeps are computed (the substep loop refreshes it again as it falls).
     qci->FillBoundary(m_geom.periodicity());
  
     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
         {
             const Real rho_km1 = (k == k_lo) ? rho_array(i,j,k) : rho_array(i,j,k-1);
             const Real rho_k = (k == k_hi + 1) ? rho_array(i,j,k-1) : rho_array(i,j,k);
             const Real qci_km1 = (k == k_lo) ? qci_array(i,j,k) : qci_array(i,j,k-1);
             const Real qci_k = (k == k_hi + 1) ? qci_array(i,j,k-1) : qci_array(i,j,k);
             const SAMFaceState face_state = sam_face_average_state(k, k_lo, k_hi,
                                                                    rho_km1, rho_k,
                                                                    zero, zero,
                                                                    qci_km1, qci_k,
                                                                    zero, zero);
             const Real vt_ice = sam_cloud_ice_terminal_velocity(face_state.qci_avg);
  
             // 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) = face_state.rho_avg * vt_ice * face_state.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);
     // ParReduce over a MultiFab gives the local rank maximum here. IceFall
     // needs one global substep count across ranks, matching PrecipFall.
     ParallelDescriptor::ReduceRealMax(wt_max);
     wt_max += std::numeric_limits<Real>::epsilon();
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_dzmin > Real(0),
         "SAM::IceFall: the minimum cell height was never set (Set_dzmin), so the fall substeps cannot be sized");
     n_substep = sam_substep_count_from_reduced_flux(wt_max, dtn, m_dzmin);
     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) {
         // The face state and the donor-cap limiter read qci at k-1, which is a ghost cell
         // whenever a box boundary is interior in z. The previous substep updated qci in
         // valid cells only, so refresh ghosts here, as PrecipFall does; otherwise the two
         // sides of a shared z face see different donor values.
         qci->FillBoundary(m_geom.periodicity());
  
         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);
             auto sed_array   = m_sed.down->array(mfi);
             auto precip_accum_array = precip_accum->array(mfi);
             const auto inflow_array = (m_sed_inflow) ? m_sed_inflow->const_array(mfi)
                                                      : Array4<const Real>{};
             const int inflow_flag = m_sed.inflow_flag_comp();
             const Real inv_dt = one / dt;
  
             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
             {
                 const Real rho_km1 = (k == k_lo) ? rho_array(i,j,k) : rho_array(i,j,k-1);
                 const Real rho_k = (k == k_hi + 1) ? rho_array(i,j,k-1) : rho_array(i,j,k);
                 const Real qci_km1 = (k == k_lo) ? qci_array(i,j,k) : qci_array(i,j,k-1);
                 const Real qci_k = (k == k_hi + 1) ? qci_array(i,j,k-1) : qci_array(i,j,k);
                 const SAMFaceState face_state = sam_face_average_state(k, k_lo, k_hi,
                                                                        rho_km1, rho_k,
                                                                        zero, zero,
                                                                        qci_km1, qci_k,
                                                                        zero, zero);
                 const Real vt_ice = sam_cloud_ice_terminal_velocity(face_state.qci_avg);
                 const int donor_k = sam_precip_face_donor_k(k, k_lo, k_hi);
                 const Real detJ_donor = (dJ_array) ? dJ_array(i,j,donor_k) : one;
  
                 // 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.
                 //
                 // Cap the flux by the ice the donor cell holds, as PrecipFall does, so a
                 // substep cannot export more than the cell contains. The clip on dqi
                 // below then binds only at roundoff, so the column budget closes.
                 fz_array(i,j,k) = sam_limit_precip_component_flux(face_state.rho_avg * vt_ice * face_state.qci_avg,
                                                                   rho_array(i,j,donor_k),
                                                                   qci_array(i,j,donor_k),
                                                                   detJ_donor,
                                                                   coef);
                 // At the top of a fine patch the ice coming in is the coarse level's flux
                 // through that face, not one built from interpolated ghost cells
                 if (inflow_array && inflow_array(i,j,k,inflow_flag) > Real(0.5)) {
                     fz_array(i,j,k) = inflow_array(i,j,k,3) * inv_dt;
                 }
                 sed_array(i,j,k,3) += fz_array(i,j,k) * dtn;
  
                 // In mm of liquid water, as the other accumulations
                 if (k == k_lo) {
                     precip_accum_array(i,j,k) += fz_array(i,j,k) * dtn / rhor * Real(1000.0);
                 }
             });
  
             // Update precip every substep
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 // Jacobian determinant
                 Real dJinv = (dJ_array) ? one/dJ_array(i,j,k) : one;
  
                 //==================================================
                 // Cloud ice sedimentation (A32)
                 //==================================================
                 Real dqi  = sam_sedimentation_tendency(fz_array(i,j,k+1), fz_array(i,j,k),
                                                        rho_array(i,j,k), dJinv, 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 SAM::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.

 {
     dt = dt_advance;
     m_geom = geom;
  
     m_z_phys_nd = z_phys_nd.get();
     m_detJ_cc   = detJ_cc.get();
  
     MicVarMap.resize(m_qmoist_size);
     MicVarMap = {MicVar::rain_accum, MicVar::snow_accum, MicVar::graupel_accum, MicVar::precip_accum};
  
     // initialize microphysics variables
     for (auto ivar = 0; ivar < MicVar::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.);
     }
  
     m_sed.define(cons_in.boxArray(), cons_in.DistributionMap(),
                  {RhoQ4_comp, RhoQ5_comp, RhoQ6_comp, RhoQ3_comp});
  
     // NOTE: For multi-level not all ranks will own a box.
     //       Furthermore, the plane average allocates space
     //       for the entire domain. We make this consistent.
     nlev = m_geom.Domain().length(2);
     zlo  = m_geom.Domain().smallEnd(2);
     zhi  = m_geom.Domain().bigEnd(2);
  
     // 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});
 }

◆ NewtonIterSat()

AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE amrex::Real SAM::NewtonIterSat	(	int &	i,
		int &	j,
		int &	k,
		const int &	SAM_moisture_type,
		const amrex::Real &	fac_cond,
		const amrex::Real &	fac_fus,
		const amrex::Real &	,
		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
 #ifdef AMREX_USE_FLOAT
         constexpr amrex::Real tol = amrex::Real(1.e-4);
 #else
         constexpr amrex::Real tol = amrex::Real(1.e-8);
 #endif
         constexpr int max_iter = 20;
  
         // Saturation moisture fractions
         amrex::Real omn, domn;
         amrex::Real qsatm, dqsatm;
         amrex::Real qsatw, dqsatw;
         amrex::Real qsati, dqsati;
         amrex::Real qn, qif, dqif;
  
         // Newton iteration vars
         int niter    = 0;
         amrex::Real F, dFdT;
  
         //=======================================================================
         // Mixed-phase Newton iteration to qv=qsat (cloud phase only)
         //=======================================================================
  
         // Initial guess for temperature & pressure
         amrex::Real tabs = tabs_array(i,j,k);
         amrex::Real pres = pres_array(i,j,k);
  
         // Mixed saturation moisture fractions
         erf_qsatw(tabs, pres, qsatw);
         erf_qsati(tabs, pres, qsati);
         omn   = sam_cloud_liquid_fraction(SAM_moisture_type, tabs, an, bn);
         qsatm = sam_mixed_qsat(omn, qsatw, qsati);
  
         // Saturation derivative WRT temperature
         domn = zero;
         if ( (SAM_moisture_type == 1) &&
              (tabs < tbgmax) && (tabs > tbgmin) ) {
             domn = an;
         }
         erf_dtqsatw(tabs, pres, dqsatw);
         erf_dtqsati(tabs, pres, dqsati);
         dqsatm = sam_mixed_dqsat_dT(omn, domn, qsatw, qsati, dqsatw, dqsati);
  
         // Enforce omega splitting
         qn   = (qt_array(i,j,k) - qsatm);
         qif  = qn * (one - omn);
         dqif = sam_mixed_dqif_dT(omn, domn, qn, dqsatm);
  
         // Function for root finding (conserve latent proxy):
         // 0 = -T_new + T_old + L_v/C_p * (qv,0 - qs,m) - L_f/C_p * (qi,0 - qi,f)
         F = sam_newton_residual(tabs, tabs_array(i,j,k),
                                 fac_cond, qv_array(i,j,k), qsatm,
                                 fac_fus , qi_array(i,j,k), qif  );
         while (std::abs(F) > tol && niter < max_iter) {
             // Derivative of function (T_new iterated on)
             dFdT  = sam_newton_residual_derivative(fac_cond, fac_fus, dqsatm, dqif);
  
             // Update the temperature
             tabs -= F/dFdT;
  
             // Mixed saturation moisture fractions
             erf_qsatw(tabs, pres, qsatw);
             erf_qsati(tabs, pres, qsati);
             omn   = sam_cloud_liquid_fraction(SAM_moisture_type, tabs, an, bn);
             qsatm = sam_mixed_qsat(omn, qsatw, qsati);
  
             // Saturation derivative WRT temperature
             domn = zero;
             if ( (SAM_moisture_type == 1) &&
                  (tabs < tbgmax) && (tabs > tbgmin) ) {
                 domn = an;
             }
             erf_dtqsatw(tabs, pres, dqsatw);
             erf_dtqsati(tabs, pres, dqsati);
             dqsatm = sam_mixed_dqsat_dT(omn, domn, qsatw, qsati, dqsatw, dqsati);
  
             // Enforce omega splitting
             qn   = (qt_array(i,j,k) - qsatm);
             qif  = qn * (one - omn);
             dqif = sam_mixed_dqif_dT(omn, domn, qn, dqsatm);
  
             // Function for root finding (conserve latent proxy):
             // 0 = -T_new + T_old + L_v/C_p * (qv,0 - qs,m) - L_f/C_p * (qi,0 - qi,f)
             F = sam_newton_residual(tabs, tabs_array(i,j,k),
                                     fac_cond, qv_array(i,j,k), qsatm,
                                     fac_fus , qi_array(i,j,k), qif  );
  
             // Update iteration
             ++niter;
         }
  
         // NOTE: the while loop already enforces the iteration limit, so the
         //       residual alone determines whether we found a valid root; testing
         //       niter < max_iter here would spuriously fail a solve that
         //       converged on exactly the last allowed iteration.
         bool valid_solution = (std::abs(F) <= tol);
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(valid_solution, "SAM mixed phase saturation adjustment failed to converge.");
  
         // Partition the change in non-precipitating q
         qv_array(i,j,k) = qsatm;
         qc_array(i,j,k) = qn * omn;
         qi_array(i,j,k) = qn * (one - omn);
         qn_array(i,j,k) = qn;
         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 SAM::Precip ( const SolverChoice & sc )

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

 {
     if (sam_is_no_precip(sc.moisture_type)) return;
  
     Real powr1 = (three + b_rain) / Real(4.0);
     Real powr2 = (Real(5.0) + b_rain) / Real(8.0);
     Real pows1 = (three + b_snow) / Real(4.0);
     Real pows2 = (Real(5.0) + b_snow) / Real(8.0);
     Real powg1 = (three + b_grau) / Real(4.0);
     Real powg2 = (Real(5.0) + b_grau) / Real(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;
  
     int SAM_moisture_type = 1;
     if (sc.moisture_type == MoistureType::SAM_NoIce) {
         SAM_moisture_type = 2;
     }
  
     const SAMPrecipConfig precip_config{
         SAM_moisture_type,
         true,
         dtn,
         rdOcp,
         fac_cond,
         fac_fus,
         fac_sub,
         eps,
         powr1,
         pows1,
         powg1,
         powr2,
         pows2,
         powg2};
     // get the temperature, density, theta, qt and qp from input
     for ( MFIter mfi(*(mic_fab_vars[MicVar::tabs]),TilingIfNotGPU()); mfi.isValid(); ++mfi) {
         auto rho_array   = mic_fab_vars[MicVar::rho]->array(mfi);
         auto theta_array = mic_fab_vars[MicVar::theta]->array(mfi);
         auto tabs_array  = mic_fab_vars[MicVar::tabs]->array(mfi);
         auto pres_array  = mic_fab_vars[MicVar::pres]->array(mfi);
  
         // Non-precipitating
         auto qv_array    = mic_fab_vars[MicVar::qv]->array(mfi);
         auto qcl_array   = mic_fab_vars[MicVar::qcl]->array(mfi);
         auto qci_array   = mic_fab_vars[MicVar::qci]->array(mfi);
         auto qn_array    = mic_fab_vars[MicVar::qn]->array(mfi);
         auto qt_array    = mic_fab_vars[MicVar::qt]->array(mfi);
  
         // Precipitating
         auto qpr_array   = mic_fab_vars[MicVar::qpr]->array(mfi);
         auto qps_array   = mic_fab_vars[MicVar::qps]->array(mfi);
         auto qpg_array   = mic_fab_vars[MicVar::qpg]->array(mfi);
         auto qp_array    = mic_fab_vars[MicVar::qp]->array(mfi);
  
         auto tbx = mfi.tilebox();
  
         ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
         {
             SAMCellState state{};
             state.rho = rho_array(i,j,k);
             state.theta = theta_array(i,j,k);
             state.tabs = tabs_array(i,j,k);
             state.pres_mbar = pres_array(i,j,k);
             state.qv = qv_array(i,j,k);
             state.qcl = qcl_array(i,j,k);
             state.qci = qci_array(i,j,k);
             state.qn = qn_array(i,j,k);
             state.qt = qt_array(i,j,k);
             state.qpr = qpr_array(i,j,k);
             state.qps = qps_array(i,j,k);
             state.qpg = qpg_array(i,j,k);
             state.qp = qp_array(i,j,k);
  
             const SAMCoefficientRow coeffs{
                 accrrc_t(k),
                 accrsi_t(k),
                 accrsc_t(k),
                 coefice_t(k),
                 evaps1_t(k),
                 evaps2_t(k),
                 accrgi_t(k),
                 accrgc_t(k),
                 evapg1_t(k),
                 evapg2_t(k),
                 evapr1_t(k),
                 evapr2_t(k)};
  
             // The scalar helper owns the local held-pressure source update;
             // this public kernel handles state staging and coefficient lookup.
             state = sam_precip_cell_update(state, coeffs, precip_config);
  
             theta_array(i,j,k) = state.theta;
             tabs_array(i,j,k) = state.tabs;
             qv_array(i,j,k) = state.qv;
             qcl_array(i,j,k) = state.qcl;
             qci_array(i,j,k) = state.qci;
             qn_array(i,j,k) = state.qn;
             qt_array(i,j,k) = state.qt;
             qpr_array(i,j,k) = state.qpr;
             qps_array(i,j,k) = state.qps;
             qpg_array(i,j,k) = state.qpg;
             qp_array(i,j,k) = state.qp;
         });
     }
 }

Referenced by Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ PrecipFall()

void SAM::PrecipFall ( const SolverChoice & sc )

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

Code modified from SAMXX, the C++ version of the SAM code.

Parameters

[in] hydro_type Type selection for the precipitation advection hydrodynamics scheme (0-3)

 {
     if (sam_is_no_precip(sc.moisture_type)) return;
  
     // Conservative component-sedimentation contract: rain, snow, and graupel
     // sediment with separate component fluxes. Each face flux is limited by
     // the donor cell's available component mass, including detJ when present.
     // The same limited bottom fluxes update surface accumulation and the
     // in-column qpr/qps/qpg flux-divergence update. This closes component and
     // total precipitation column budgets up to roundoff. Sedimentation does
     // not update theta.
  
     Real rho_0 = Real(1.29);
  
     Real gamr3 = erf_gammafff(Real(4.0)+b_rain);
     Real gams3 = erf_gammafff(Real(4.0)+b_snow);
     Real gamg3 = erf_gammafff(Real(4.0)+b_grau);
  
     Real vrain = (a_rain*gamr3/Real(6.0))*std::pow((PI*rhor*nzeror),-crain);
     Real vsnow = (a_snow*gams3/Real(6.0))*std::pow((PI*rhos*nzeros),-csnow);
     Real vgrau = (a_grau*gamg3/Real(6.0))*std::pow((PI*rhog*nzerog),-cgrau);
  
     Real dtn  = dt;
     // NOTE: coef carries only the reference vertical spacing. The physical
     //       spacing is supplied by the inverse Jacobian in the tendency below,
     //       since detJ = dz_phys/CellSize(2). Dividing by m_dzmin here as well
     //       would apply the vertical metric twice. m_dzmin is still the correct
     //       length scale for the substep (CFL) count, which must bound the
     //       thinnest cell in the domain.
     Real coef = dtn * m_geom.InvCellSize(2);
  
     auto domain = m_geom.Domain();
     int k_lo = domain.smallEnd(2);
     int k_hi = domain.bigEnd(2);
  
     auto qpr   = mic_fab_vars[MicVar::qpr];
     auto qps   = mic_fab_vars[MicVar::qps];
     auto qpg   = mic_fab_vars[MicVar::qpg];
     auto qp    = mic_fab_vars[MicVar::qp];
     auto rho   = mic_fab_vars[MicVar::rho];
     auto tabs  = mic_fab_vars[MicVar::tabs];
     auto rain_accum = mic_fab_vars[MicVar::rain_accum];
     auto snow_accum = mic_fab_vars[MicVar::snow_accum];
     auto graupel_accum = mic_fab_vars[MicVar::graupel_accum];
     auto precip_accum = mic_fab_vars[MicVar::precip_accum];
  
     auto ba    = tabs->boxArray();
     auto dm    = tabs->DistributionMap();
     auto ngrow = tabs->nGrowVect();
  
     MultiFab fz;
     fz.define(convert(ba, IntVect(0,0,1)), dm, 3, ngrow);
  
     // The face states read tabs, rho and the precipitating species at k-1, a ghost cell
     // whenever a box boundary is interior in z. Cloud and Precip have just updated them in
     // valid cells only, so refresh the ghosts before the fluxes that size the substeps are
     // computed; the substep loop below refreshes the species again as they fall. tabs does
     // not change here, so once is enough for it.
     tabs->FillBoundary(m_geom.periodicity());
     rho->FillBoundary(m_geom.periodicity());
     qpr->FillBoundary(m_geom.periodicity());
     qps->FillBoundary(m_geom.periodicity());
     qpg->FillBoundary(m_geom.periodicity());
  
     //  Precompute the vertical component fluxes for the legacy reduced-flux substep rule.
     for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
         auto qpr_array  = qpr->const_array(mfi);
         auto qps_array  = qps->const_array(mfi);
         auto qpg_array  = qpg->const_array(mfi);
         auto rho_array  = rho->const_array(mfi);
         auto tabs_array = tabs->const_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
         {
             const SAMPrecipComponentFaceState face_state =
                 sam_precip_component_face_state(rho_array, tabs_array,
                                                 qpr_array, qps_array, qpg_array,
                                                 i, j, k, k_lo, k_hi);
             const SAMPrecipFluxComponents raw_fluxes =
                 sam_precip_component_fluxes_from_face_state(face_state,
                                                             vrain, vsnow, vgrau);
             const SAMPrecipFluxComponents corrected_fluxes =
                 sam_precip_flux_components_density_corrected(raw_fluxes,
                                                              rho_0,
                                                              face_state.rho_avg);
  
             fz_array(i,j,k,0) = corrected_fluxes.rain;
             fz_array(i,j,k,1) = corrected_fluxes.snow;
             fz_array(i,j,k,2) = corrected_fluxes.graupel;
         });
     }
  
     // Compute the legacy reduced-flux substep count from the maximum
     // density-corrected sedimentation flux rather than a direct fall speed.
     Real max_reduced_flux;
     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,0)
                                  + ma_fz_arr[box_no](i,j,k,1)
                                  + ma_fz_arr[box_no](i,j,k,2) };
                          });
     max_reduced_flux = get<0>(max);
     // ParReduce over a MultiFab gives the local rank maximum here. PrecipFall
     // needs one global substep count across ranks, including ranks that own no
     // active tiles for this reduction, so keep the explicit MPI max reduction.
     ParallelDescriptor::ReduceRealMax(max_reduced_flux);
     max_reduced_flux += std::numeric_limits<Real>::epsilon();
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_dzmin > Real(0),
         "SAM::PrecipFall: the minimum cell height was never set (Set_dzmin), so the fall substeps cannot be sized");
     n_substep = sam_substep_count_from_reduced_flux(max_reduced_flux, dtn, m_dzmin);
     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)
     {
         // Refresh ghosts so both sides of a shared z face use the same donor values
         qpr->FillBoundary(m_geom.periodicity());
         qps->FillBoundary(m_geom.periodicity());
         qpg->FillBoundary(m_geom.periodicity());
         mic_fab_vars[MicVar::rho]->FillBoundary(m_geom.periodicity());
  
         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->const_array(mfi);
             auto tabs_array   = tabs->const_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 graupel_accum_array = graupel_accum->array(mfi);
             auto precip_accum_array = precip_accum->array(mfi);
             auto sed_array = m_sed.down->array(mfi);
             const auto inflow_array = (m_sed_inflow) ? m_sed_inflow->const_array(mfi)
                                                      : Array4<const Real>{};
             const int inflow_flag = m_sed.inflow_flag_comp();
             const Real inv_dt = one / dt;
  
             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
             {
                 const SAMPrecipComponentFaceState face_state =
                     sam_precip_component_face_state(rho_array, tabs_array,
                                                     qpr_array, qps_array, qpg_array,
                                                     i, j, k, k_lo, k_hi);
                 const SAMPrecipFluxComponents raw_fluxes =
                     sam_precip_component_fluxes_from_face_state(face_state,
                                                                 vrain, vsnow, vgrau);
                 const SAMPrecipFluxComponents corrected_fluxes =
                     sam_precip_flux_components_density_corrected(raw_fluxes,
                                                                  rho_0,
                                                                  face_state.rho_avg);
                 const int donor_k = sam_precip_face_donor_k(k, k_lo, k_hi);
                 const Real detJ_donor = (dJ_array) ? dJ_array(i,j,donor_k) : one;
  
                 // 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,0) = sam_limit_precip_component_flux(corrected_fluxes.rain,
                                                                     rho_array(i,j,donor_k),
                                                                     qpr_array(i,j,donor_k),
                                                                     detJ_donor,
                                                                     coef);
                 fz_array(i,j,k,1) = sam_limit_precip_component_flux(corrected_fluxes.snow,
                                                                     rho_array(i,j,donor_k),
                                                                     qps_array(i,j,donor_k),
                                                                     detJ_donor,
                                                                     coef);
                 fz_array(i,j,k,2) = sam_limit_precip_component_flux(corrected_fluxes.graupel,
                                                                     rho_array(i,j,donor_k),
                                                                     qpg_array(i,j,donor_k),
                                                                     detJ_donor,
                                                                     coef);
  
                 // At the top of a fine patch the precipitation coming in is the coarse
                 // level's flux through that face, not one built from interpolated ghost cells
                 if (inflow_array && inflow_array(i,j,k,inflow_flag) > Real(0.5)) {
                     fz_array(i,j,k,0) = inflow_array(i,j,k,0) * inv_dt;
                     fz_array(i,j,k,1) = inflow_array(i,j,k,1) * inv_dt;
                     fz_array(i,j,k,2) = inflow_array(i,j,k,2) * inv_dt;
                 }
  
                 sed_array(i,j,k,0) += fz_array(i,j,k,0) * dtn;
                 sed_array(i,j,k,1) += fz_array(i,j,k,1) * dtn;
                 sed_array(i,j,k,2) += fz_array(i,j,k,2) * dtn;
  
                 if (k == k_lo) {
                     const SAMSurfaceAccumulation surface_accum =
                         sam_surface_accumulation_from_component_fluxes(
                             {fz_array(i,j,k,0), fz_array(i,j,k,1), fz_array(i,j,k,2)},
                             dtn);
                     rain_accum_array(i,j,k)  = rain_accum_array(i,j,k)  + surface_accum.rain;
                     snow_accum_array(i,j,k)  = snow_accum_array(i,j,k)  + surface_accum.snow;
                     graupel_accum_array(i,j,k) = graupel_accum_array(i,j,k) + surface_accum.graupel;
                     // The total is in mm of liquid water, whatever the species' own densities
                     precip_accum_array(i,j,k) += (fz_array(i,j,k,0) + fz_array(i,j,k,1) + fz_array(i,j,k,2))
                                                * dtn / rhor * Real(1000.0);
                 }
             });
  
             // Update precip every substep
             ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
             {
                 // Jacobian determinant
                 Real dJinv = (dJ_array) ? one/dJ_array(i,j,k) : one;
  
                 //==================================================
                 // Precipitating sedimentation (A19)
                 //==================================================
                 const Real dqpr = sam_sedimentation_tendency(fz_array(i,j,k+1,0), fz_array(i,j,k,0),
                                                              rho_array(i,j,k), dJinv, coef);
                 const Real dqps = sam_sedimentation_tendency(fz_array(i,j,k+1,1), fz_array(i,j,k,1),
                                                              rho_array(i,j,k), dJinv, coef);
                 const Real dqpg = sam_sedimentation_tendency(fz_array(i,j,k+1,2), fz_array(i,j,k,2),
                                                              rho_array(i,j,k), dJinv, coef);
  
                 qpr_array(i,j,k) = std::max(Real(0.0), qpr_array(i,j,k) + dqpr);
                 qps_array(i,j,k) = std::max(Real(0.0), qps_array(i,j,k) + dqps);
                 qpg_array(i,j,k) = std::max(Real(0.0), qpg_array(i,j,k) + dqpg);
                  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() [1/2]

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

inline

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

Here is the call graph for this function:

◆ Qmoist_Ptr() [2/2]

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

inlineoverridevirtual

Reimplemented from NullMoist.

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

Here is the call graph for this function:

◆ Qmoist_Restart_Vars()

void SAM::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");
         a_idx.push_back(3); a_names.push_back("PrecipAccum");
     }

◆ Qmoist_Size()

int SAM::Qmoist_Size ( )

inlineoverridevirtual

Reimplemented from NullMoist.

177 { return SAM::m_qmoist_size; }

◆ Qstate_Moist_NumConc_Size()

int SAM::Qstate_Moist_NumConc_Size ( )

inlineoverridevirtual

Reimplemented from NullMoist.

183 { return SAM::n_qstate_moist_numconc_size; }

SAM::n_qstate_moist_numconc_size

int n_qstate_moist_numconc_size

Definition: ERF_SAM.H:355

◆ Qstate_Moist_Size()

int SAM::Qstate_Moist_Size ( )

inlineoverridevirtual

Reimplemented from NullMoist.

180 { return SAM::n_qstate_moist_size; }

SAM::n_qstate_moist_size

int n_qstate_moist_size

Definition: ERF_SAM.H:352

◆ Sed_Flux()

const SedimentationFlux* SAM::Sed_Flux ( ) const

inlineoverridevirtual

The sedimentation face fluxes of the most recent Advance, or null for a scheme that does not sediment (see ERF_SedimentationFlux.H)

Reimplemented from NullMoist.

154 { return &m_sed; }

◆ Set_dzmin()

void SAM::Set_dzmin ( const amrex::Real dz_min )

inlineoverridevirtual

Reimplemented from NullMoist.

     {
         m_dzmin = dz_min;
     }

◆ Update_Micro_Vars() [1/2]

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

inlineoverridevirtual

Reimplemented from NullMoist.

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

Here is the call graph for this function:

◆ Update_Micro_Vars() [2/2]

void SAM::Update_Micro_Vars	(	amrex::MultiFab &	cons_in,
		const amrex::MultiFab *	base_state
	)

overridevirtual

Reimplemented from NullMoist.

◆ Update_State_Vars()

void SAM::Update_State_Vars	(	amrex::MultiFab &	cons_in,
		const amrex::MultiFab &
	)

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> SAM::accrgc

private

◆ accrgi

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

private

◆ accrrc

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

private

◆ accrsc

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

private

◆ accrsi

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

private

◆ CFL_MAX

constexpr amrex::Real SAM::CFL_MAX = myhalf

staticconstexprprivate

◆ coefice

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

private

◆ dt

amrex::Real SAM::dt

private

Referenced by Advance().

◆ evapg1

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

private

◆ evapg2

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

private

◆ evapr1

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

private

◆ evapr2

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

private

◆ evaps1

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

private

◆ evaps2

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

private

◆ m_axis

int SAM::m_axis

private

Referenced by Define().

◆ m_detJ_cc

amrex::MultiFab* SAM::m_detJ_cc

private

◆ m_do_cond

bool SAM::m_do_cond

private

Referenced by Define().

◆ m_dzmin

amrex::Real SAM::m_dzmin {0.0}

private

Referenced by Set_dzmin().

◆ m_fac_cond

amrex::Real SAM::m_fac_cond

private

Referenced by Define().

◆ m_fac_fus

amrex::Real SAM::m_fac_fus

private

Referenced by Define().

◆ m_fac_sub

amrex::Real SAM::m_fac_sub

private

Referenced by Define().

◆ m_geom

amrex::Geometry SAM::m_geom

private

◆ m_moisture_type

MoistureType SAM::m_moisture_type = MoistureType::None

private

Referenced by Define(), and Get_Surface_Precip_Accumulation_Ptrs().

◆ m_qmoist_size

int SAM::m_qmoist_size = 4

private

Referenced by Qmoist_Ptr(), and Qmoist_Size().

◆ m_rdOcp

amrex::Real SAM::m_rdOcp {RdoCp}

private

Referenced by Define().

◆ m_sed

SedimentationFlux SAM::m_sed

private

Referenced by Advance(), and Sed_Flux().

◆ m_z_phys_nd

amrex::MultiFab* SAM::m_z_phys_nd

private

◆ mic_fab_vars

amrex::Array<FabPtr, MicVar::NumVars> SAM::mic_fab_vars

private

Referenced by Get_Surface_Precip_Accumulation_Ptrs(), and Qmoist_Ptr().

◆ MicVarMap

amrex::Vector<int> SAM::MicVarMap

private

Referenced by Qmoist_Ptr().

◆ n_qstate_moist_numconc_size

int SAM::n_qstate_moist_numconc_size = 0

private

Referenced by Qstate_Moist_NumConc_Size().

◆ n_qstate_moist_size

int SAM::n_qstate_moist_size = 6

private

Referenced by Qstate_Moist_Size().

◆ nlev

int SAM::nlev

private

◆ pres1d

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

private

◆ qn1d

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

private

◆ qt1d

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

private

◆ qv1d

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

private

◆ rho1d

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

private

◆ tabs1d

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

private

◆ zhi

int SAM::zhi

private

◆ zlo

int SAM::zlo

private

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

Source/Microphysics/SAM/ERF_SAM.H
Source/Microphysics/SAM/ERF_CloudSAM.cpp
Source/Microphysics/SAM/ERF_IceFall.cpp
Source/Microphysics/SAM/ERF_InitSAM.cpp
Source/Microphysics/SAM/ERF_Precip.cpp
Source/Microphysics/SAM/ERF_PrecipFall.cpp
Source/Microphysics/SAM/ERF_UpdateSAM.cpp

Public Member Functions

Static Public Member Functions

Private Types

Private Attributes

Static Private Attributes

Additional Inherited Members

Member Typedef Documentation

◆ FabPtr

Constructor & Destructor Documentation

◆ SAM()

◆ ~SAM()

Member Function Documentation

◆ Advance()

◆ Cloud()

◆ CoefficientRowAt()

◆ Compute_Coefficients()

◆ Copy_Micro_to_State()

◆ Copy_State_to_Micro() [1/2]

◆ Copy_State_to_Micro() [2/2]

◆ Define()

◆ Get_Surface_Precip_Accumulation_Ptrs()

◆ IceFall()

◆ Init()

◆ NewtonIterSat()

◆ Precip()

◆ PrecipFall()

◆ Qmoist_Ptr() [1/2]

◆ Qmoist_Ptr() [2/2]

◆ Qmoist_Restart_Vars()

◆ Qmoist_Size()

◆ Qstate_Moist_NumConc_Size()

◆ Qstate_Moist_Size()

◆ Sed_Flux()

◆ Set_dzmin()

◆ Update_Micro_Vars() [1/2]

◆ Update_Micro_Vars() [2/2]

◆ Update_State_Vars()

Member Data Documentation

◆ accrgc

◆ accrgi

◆ accrrc

◆ accrsc

◆ accrsi

◆ CFL_MAX

◆ coefice

◆ dt

◆ evapg1

◆ evapg2

◆ evapr1

◆ evapr2

◆ evaps1

◆ evaps2

◆ m_axis

◆ m_detJ_cc

◆ m_do_cond

◆ m_dzmin

◆ m_fac_cond

◆ m_fac_fus

◆ m_fac_sub

◆ m_geom

◆ m_moisture_type

◆ m_qmoist_size

◆ m_rdOcp

◆ m_sed

◆ m_z_phys_nd

◆ mic_fab_vars

◆ MicVarMap

◆ n_qstate_moist_numconc_size

◆ n_qstate_moist_size

◆ nlev

◆ pres1d

◆ qn1d

◆ qt1d

◆ qv1d

◆ rho1d

◆ tabs1d

◆ zhi

◆ zlo