RadiationSimple Class Reference

#include <ERF_RadiationSimple.H>

Inheritance diagram for RadiationSimple:

Collaboration diagram for RadiationSimple:

Public Member Functions
	RadiationSimple (const int &, const SolverChoice &sc)
	~RadiationSimple ()=default
virtual void	Init (const amrex::Geometry &geom, const amrex::BoxArray &ba, amrex::MultiFab *cons_in) override
virtual void	Run (int &level, int &step, double &time, const double &dt, const amrex::BoxArray &ba, amrex::Geometry &geom, amrex::MultiFab *cons_in, amrex::iMultiFab *lmask, amrex::MultiFab *t_surf, const amrex::Vector< const amrex::MultiFab * > &lsm_input_ptrs, const amrex::Vector< amrex::MultiFab * > &lsm_output_ptrs, amrex::MultiFab *qheating_rates, amrex::MultiFab *rad_fluxes, amrex::MultiFab *z_phys, amrex::MultiFab *lat, amrex::MultiFab *lon, const bool force_update, amrex::Real &declin, amrex::Real &calday) override
virtual void	WriteDataLog (const double &time) override
bool	is_nested_patch () const override
const amrex::MultiFab *	get_radlwdn ()
const amrex::MultiFab *	get_radqrlw ()
Public Member Functions inherited from IRadiation
virtual	~IRadiation ()=default
virtual bool	radiation_updated () const
	Reports whether the last radiation call updated its outputs. More...
virtual amrex::Vector< std::string >	get_lsm_input_varnames ()
virtual amrex::Vector< std::string >	get_lsm_output_varnames ()
void	setupDataLog ()
void	setDataLogFrequency (const int nstep)
bool	hasDatalog ()

Private Attributes
bool	m_have_qv = false
bool	m_have_qc = false
bool	m_ice = false
int	m_qv_comp = -1
int	m_qc_comp = -1
int	m_qi_comp = -1
amrex::Geometry	m_geom
amrex::BoxArray	m_ba
bool	m_partial_columns = false
amrex::MultiFab	deltaq
amrex::MultiFab	flux
std::unique_ptr< amrex::MultiFab >	radlwdn
std::unique_ptr< amrex::MultiFab >	radqrlw

Additional Inherited Members
Protected Attributes inherited from IRadiation
std::unique_ptr< std::fstream >	datalog = nullptr
std::string	datalogname
int	datalog_int = -1

Constructor & Destructor Documentation

RadiationSimple()

RadiationSimple::RadiationSimple ( const int &  , const SolverChoice &  sc  )

inline

                                             {
        // Every moist species we read is taken from the configured moisture-component
        // mapping, so each scheme is indexed by its own layout rather than by an
        // enumerated model list: every scheme carrying cloud ice (SAM, Morrison, WSM6,
        // WDM6, ...) is picked up, and the warm-rain schemes are not.
        m_qv_comp = sc.moisture_indices.qv;
        m_qc_comp = sc.moisture_indices.qc;
        m_qi_comp = sc.moisture_indices.qi;
 
        m_have_qv = (m_qv_comp >= 0);
        m_have_qc = (m_qc_comp >= 0);
        m_ice     = (m_qi_comp >= 0);
    }

~RadiationSimple()

RadiationSimple::~RadiationSimple ( )

default

Member Function Documentation

get_radlwdn()

const amrex::MultiFab* RadiationSimple::get_radlwdn ( )

inline

{ return radlwdn.get(); }

get_radqrlw()

const amrex::MultiFab* RadiationSimple::get_radqrlw ( )

inline

{ return radqrlw.get(); }

Init()

void RadiationSimple::Init ( const amrex::Geometry &  geom, const amrex::BoxArray &  ba, amrex::MultiFab *  cons_in  )

overridevirtual

Implements IRadiation.

{
    m_geom = geom;
    m_ba = ba;
 
    // Set (or cleared) on every Init, so a level that spans the column again after a
    // regrid resumes the sweep.
    m_partial_columns = false;
    for (int ib = 0; ib < ba.size(); ++ib) {
        if (ba[ib].smallEnd(2) != geom.Domain().smallEnd(2) || ba[ib].bigEnd(2) != geom.Domain().bigEnd(2)) {
            m_partial_columns = true;
            break;
        }
    }
 
    DistributionMapping dm = cons_in->DistributionMap();
    deltaq.define(ba, dm, 1, 0);
    flux.define(ba, dm, 1, IntVect(0, 0, 1));
    radlwdn = std::make_unique<amrex::MultiFab>(ba, dm, 1, 0);
    radqrlw = std::make_unique<amrex::MultiFab>(ba, dm, 1, 0);
 
    deltaq.setVal(zero);
    flux.setVal(zero);
    radlwdn->setVal(zero);
    radqrlw->setVal(zero);
}

is_nested_patch()

bool RadiationSimple::is_nested_patch ( ) const

inlineoverridevirtual

Reimplemented from IRadiation.

{ return m_partial_columns; }

Run()

void RadiationSimple::Run ( int &  level, int &  step, double &  time, const double &  dt, const amrex::BoxArray &  ba, amrex::Geometry &  geom, amrex::MultiFab *  cons_in, amrex::iMultiFab *  lmask, amrex::MultiFab *  t_surf, const amrex::Vector< const amrex::MultiFab * > &  lsm_input_ptrs, const amrex::Vector< amrex::MultiFab * > &  lsm_output_ptrs, amrex::MultiFab *  qheating_rates, amrex::MultiFab *  rad_fluxes, amrex::MultiFab *  z_phys, amrex::MultiFab *  lat, amrex::MultiFab *  lon, const bool  force_update, amrex::Real &  declin, amrex::Real &  calday  )

overridevirtual

Implements IRadiation.

{
    if (m_partial_columns) {
        // Level 0 has no parent to interpolate from: grids split in z are unsupported here.
        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(level > 0,
            "RadiationSimple: every level-0 box must span the domain in z (do not split grids in z)");
        return;
    }
 
    constexpr amrex::Real cp_spec = 1015.0;
    constexpr amrex::Real coef = 70.0;
    constexpr amrex::Real coef1 = 22.0;
    constexpr amrex::Real f0=3.75e-6;
    constexpr amrex::Real xk = 85.0;
 
    const Real fixed_dz = geom.CellSize(2);
    const int nz = geom.Domain().length(2);
 
    flux.setVal(zero);
    deltaq.setVal(zero);
 
    const bool have_qv = m_have_qv;
    const bool have_qc = m_have_qc;
    const bool ice     = m_ice;
    const int  qv_comp = m_qv_comp;
    const int  qc_comp = m_qc_comp;
    const int  qi_comp = m_qi_comp;
 
    for (MFIter mfi(*cons_in, TileNoZ()); mfi.isValid(); ++mfi)
    {
        Box box = mfi.tilebox();
        box.makeSlab(2, 0);
 
        const Array4<const Real>& cons_arr = cons_in->const_array(mfi);
        const Array4<const Real>& z_nd_arr = (z_phys) ? z_phys->const_array(mfi) : Array4<const Real>{};
        const Array4<Real>& deltaq_arr = deltaq.array(mfi);
        const Array4<Real>& flux_arr = flux.array(mfi);
        const Array4<Real>& radlwdn_arr = radlwdn->array(mfi);
        const Array4<Real>& radqrlw_arr = radqrlw->array(mfi);
        const Array4<Real>& qheating_arr = qheating_rates->array(mfi);
 
        ParallelFor(box, [=] AMREX_GPU_DEVICE (int i, int j, int)
        {
            int itop = 0;
            Real qzinf = 0.0;  // holds accumulated optical depth between z and domain top
            Real qzeroz = 0.0; // holds accumulated optical depth between surface and z
 
            for (int k = 0; k < nz; k++)
            {
                Real rho = cons_arr(i, j, k, Rho_comp);
                Real qv = (have_qv) ? cons_arr(i, j, k, qv_comp) / rho : zero;
                Real qc = (have_qc) ? cons_arr(i, j, k, qc_comp) / rho : zero;
                Real qi = (ice)     ? cons_arr(i, j, k, qi_comp) / rho : zero;
 
                Real dz = (z_nd_arr) ? Real(0.25) * ( (z_nd_arr(i  ,j  ,k+1) - z_nd_arr(i  ,j  ,k))
                                            + (z_nd_arr(i+1,j  ,k+1) - z_nd_arr(i+1,j  ,k))
                                            + (z_nd_arr(i  ,j+1,k+1) - z_nd_arr(i  ,j+1,k))
                                            + (z_nd_arr(i+1,j+1,k+1) - z_nd_arr(i+1,j+1,k)) ) : fixed_dz;
 
                // optical depth only includes that due to liquid water
                if (qc + qi > 0.0)
                {
                    deltaq_arr(i, j, k) = xk * rho * (qc + qi) * dz;
                }
 
                // inversion height is top of highest layer w/q>8g/kg
                if (qv + qc + qi > 0.008)
                {
                    itop = max(itop, k+1); // note zi(k+1) is inversion height
                }
 
                // accumulate optical depth in qzinf
                qzinf = qzinf + deltaq_arr(i, j, k);
            }
 
            // note: qzinf now holds total optical depth of column (due to cloud)
            //       qzeroz is initialized to zero since first level is at surface
 
            // compute net upward longwave flux up to inversion height.  When the
            // moist layer reaches the top cell, itop = nz: the face at nz is then
            // written below, where its above-inversion term vanishes (zi = z_itop),
            // and deltaq has no cell at nz to read.
            for (int k = 0; k <= amrex::min(itop, nz - 1); k++)
            {
                flux_arr(i, j, k) = coef*exp(-qzinf) + coef1*exp(-qzeroz);
                qzinf = qzinf - deltaq_arr(i, j, k);
                qzeroz = qzeroz + deltaq_arr(i, j, k);
            }
 
            // compute net upward longwave flux above inversion height
            // this includes correction for clearsky fluxes which balances
            // the prescribed subsidence heating above the inversion
            Real z_itop = (z_nd_arr) ? z_nd_arr(i, j, itop) : itop * fixed_dz;
            for (int k = itop + 1; k < nz; k++) {
                Real rho_nd = 0.5 * (cons_arr(i, j, k - 1, Rho_comp) + cons_arr(i, j, k, Rho_comp)); // rho at interface
                Real zi = (z_nd_arr) ? z_nd_arr(i, j, k) : k * fixed_dz;
                flux_arr(i, j, k) = coef * exp(-qzinf) + coef1 * exp(-qzeroz)
                                    + cp_spec * rho_nd * f0 * (0.25 * zi + 0.75 * z_itop)
                                    * pow(zi - z_itop, 1.0 / 3.0);
                qzinf = qzinf - deltaq_arr(i, j, k);
                qzeroz = qzeroz + deltaq_arr(i, j, k);
            }
            Real rho_nd = 0.5 * (cons_arr(i, j, nz - 1, Rho_comp) + cons_arr(i, j, nz, Rho_comp));
            Real zi = (z_nd_arr) ? z_nd_arr(i, j, nz) : nz * fixed_dz;
            flux_arr(i, j, nz) = coef * exp(-qzinf) + coef1 * exp(-qzeroz)
                                 + cp_spec * rho_nd * f0 * (0.25 * zi + 0.75 * z_itop)
                                 * pow(zi - z_itop, 1.0 / 3.0);
 
            // note that our flux differs from the specification in that it
            // uses the local density (rather than the interface density) in
            // computing the correction to the flux above the inversion.
            // Formulating things in this way ensures a rough balance between
            // the subsidence heating and radiative cooling above the inversion.
 
            // compute radiative heating as divergence of net upward lw flux
            for (int k = 0; k < nz; k++)
            {
                Real dz = (z_nd_arr) ? 0.25 * ( (z_nd_arr(i  ,j  ,k+1) - z_nd_arr(i  ,j  ,k))
                                            + (z_nd_arr(i+1,j  ,k+1) - z_nd_arr(i+1,j  ,k))
                                            + (z_nd_arr(i  ,j+1,k+1) - z_nd_arr(i  ,j+1,k))
                                            + (z_nd_arr(i+1,j+1,k+1) - z_nd_arr(i+1,j+1,k)) ) : fixed_dz;
 
                Real cpmassl = cp_spec * cons_arr(i, j, k, Rho_comp) * dz; // thermal mass
                Real FTHRL = -(flux_arr(i, j, k+1) - flux_arr(i, j, k)) / cpmassl;
                // convert heating rate FTHRL for theta_d
                Real qv_k = (have_qv) ? cons_arr(i,j,k,qv_comp)/cons_arr(i,j,k,Rho_comp) : zero;
                Real pres = getPgivenRTh(cons_arr(i,j,k,RhoTheta_comp), qv_k);
                Real exner = getExnergivenP(pres, R_d/Cp_d);
                FTHRL /= exner;
                qheating_arr(i, j, k, 1) = FTHRL;          // radiative heating for source term
                radlwdn_arr(i, j, k) = flux_arr(i, j, k);  // net lw flux
                radqrlw_arr(i, j, k) = FTHRL;              // net lw heating
 
                // The net flux is not written to rad_fluxes: that array holds
                // (SW up, SW down, LW up, LW down), and this model computes
                // none of those. It stays zero, and HasRadiation() is false.
            }
        });
    }
}

Here is the call graph for this function:

WriteDataLog()

void RadiationSimple::WriteDataLog ( const double &  time )

overridevirtual

Implements IRadiation.

{
    constexpr int datwidth = 14;
    constexpr int datprecision = 9;
    constexpr int timeprecision = 13;
 
    Gpu::HostVector<Real> h_avg_radlwdn, h_avg_radqrlw;
 
    auto domain = m_geom.Domain();
    h_avg_radlwdn = sumToLine(*radlwdn, 0, 1, domain, 2);
    h_avg_radqrlw = sumToLine(*radqrlw, 0, 1, domain, 2);
 
    Real area_z = static_cast<Real>(domain.length(0)*domain.length(1));
    int nz = domain.length(2);
    for (int k = 0; k < nz; k++) {
        h_avg_radlwdn[k] /= area_z;
        h_avg_radqrlw[k] /= area_z;
    }
 
    if (ParallelDescriptor::IOProcessor()) {
        std::ostream& log = *datalog;
 
        if (log.good()) {
            for (int k = 0; k < nz; k++)
            {
                Real z = k * m_geom.CellSize(2);
                log << std::setw(datwidth) << std::setprecision(timeprecision) << time << " "
                    << std::setw(datwidth) << std::setprecision(datprecision) << z << " "
                    << h_avg_radlwdn[k] << " "
                    << h_avg_radqrlw[k] << std::endl;
            }
            // Write top face values
            Real z = nz * m_geom.CellSize(2);
            log << std::setw(datwidth) << std::setprecision(timeprecision) << time << " "
                << std::setw(datwidth) << std::setprecision(datprecision) << z << " "
                << 0.0 << " "
                << 0.0 << std::endl;
        }
    }
}

Member Data Documentation

deltaq

amrex::MultiFab RadiationSimple::deltaq

private

flux

amrex::MultiFab RadiationSimple::flux

private

m_ba

amrex::BoxArray RadiationSimple::m_ba

private

m_geom

amrex::Geometry RadiationSimple::m_geom

private

m_have_qc

bool RadiationSimple::m_have_qc = false

private

Referenced by RadiationSimple().

m_have_qv

bool RadiationSimple::m_have_qv = false

private

Referenced by RadiationSimple().

m_ice

bool RadiationSimple::m_ice = false

private

Referenced by RadiationSimple().

m_partial_columns

bool RadiationSimple::m_partial_columns = false

private

Referenced by is_nested_patch().

m_qc_comp

int RadiationSimple::m_qc_comp = -1

private

Referenced by RadiationSimple().

m_qi_comp

int RadiationSimple::m_qi_comp = -1

private

Referenced by RadiationSimple().

m_qv_comp

int RadiationSimple::m_qv_comp = -1

private

Referenced by RadiationSimple().

radlwdn

std::unique_ptr<amrex::MultiFab> RadiationSimple::radlwdn

private

Referenced by get_radlwdn().

radqrlw

std::unique_ptr<amrex::MultiFab> RadiationSimple::radqrlw

private

Referenced by get_radqrlw().

---

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

• Source/Radiation/Simple/ERF_RadiationSimple.H
• Source/Radiation/Simple/ERF_RadiationSimple.cpp