ERF_M2005EffRadius.H File Reference

#include "ERF_Constants.H"
#include "ERF_RadConstants.H"
#include <cmath>

Include dependency graph for ERF_M2005EffRadius.H:

Go to the source code of this file.

Functions
YAKL_INLINE void	m2005_effradius (const real2d &ql, const real2d &nl, const real2d &qi, const real2d &ni, const real2d &qs, const real2d &ns, const real2d &cld, const real2d &pres, const real2d &tk, const real2d &effl, const real2d &effi, const real2d &deffi, const real2d &lamcrad, const real2d &pgamrad, const real2d &des)

Function Documentation

m2005_effradius()

YAKL_INLINE void m2005_effradius	(	const real2d &	ql,
		const real2d &	nl,
		const real2d &	qi,
		const real2d &	ni,
		const real2d &	qs,
		const real2d &	ns,
		const real2d &	cld,
		const real2d &	pres,
		const real2d &	tk,
		const real2d &	effl,
		const real2d &	effi,
		const real2d &	deffi,
		const real2d &	lamcrad,
		const real2d &	pgamrad,
		const real2d &	des
	)

{
    // Main computation
    const real pi = 3.1415926535897932384626434;
    const real qsmall = 1.0e-14; // in the SAM source code (module_mp_graupel)
    const real rhow = 997.;      // in module_mp_graupel, SAM
    const real rhoi = 500.;      // in both CAM and SAM
    const real dcs = 125.e-6;   // in module_mp_graupel, SAM
    const real ci = rhoi*pi/6.;
    const real di = 3.;
 
    // for snow water
    const real rhos = 100.;
    const real cs = rhos*pi/6.;
    const real ds = 3.;
    const real mincld  = 0.0001;
 
    int ncol = pres.extent(0);
    int nlev = pres.extent(1);
    // Effective diameters of snow crystals
    parallel_for (SimpleBounds<2> (ncol, nlev), YAKL_LAMBDA (int i, int j)
    {
        auto rho  = pres(i,j)/(287.15*tk(i,j));   // air density [kg/m3]
        auto cldm = std::max(cld(i,j), mincld);
        auto qlic = std::min(5.e-3, std::max(0., ql(i,j)/cldm));
        auto qiic = std::min(5.e-3, std::max(0., qi(i,j)/cldm));
        auto nlic = std::max(nl(i,j), 0.)/cldm;
        auto niic = std::max(ni(i,j), 0.)/cldm;
 
        real res;
        if(qs(i,j) > 1.0e-7) {
            auto lammaxs=1./10.e-6;
            auto lammins=1./2000.e-6;
            auto lams = std::pow((std::tgamma(1.+ds)*cs*ns(i,j)/qs(i,j)), (1./ds));
            lams = std::min(lammaxs, std::max(lams,lammins));
            res = 1.5/lams*1.0e6;
        }
        else {
            res = 500.;
        }
 
        // from Hugh Morrision: rhos/917 accounts for assumptions about
        // ice density in the Mitchell optics.
        des(i,j) = res*rhos/917.*2.;
 
        // Effective radius of cloud ice droplet
        if( qiic >= qsmall ) {
            niic   = std::min(niic, qiic*1.e20);
            auto lammax = 1./1.e-6;         // in module_mp_graupel, SAM
            auto lammin = 1./(2.*dcs+100.e-6);    // in module_mp_graupel, SAM
            auto lami   = std::pow((std::tgamma(1.+di)*ci*niic/qiic), (1./di));
            lami        = std::min(lammax, std::max(lami, lammin));
            effi(i,j)   = 1.5/lami*1.e6;
        }
        else {
            effi(i,j)   = 25.;
        }
 
        // hm ice effective radius for david mitchell's optics
        // ac morrison indicates that this is effective diameter
        // ac morrison indicates 917 (for the density of pure ice..)
        deffi(i,j)  = effi(i,j)*rhoi/917.*2.;
 
        // Effective radius of cloud liquid droplet
        if( qlic >= qsmall ) {
            // Matin et al., 1994 (JAS) formula for pgam (the same is used in both CAM and SAM).
            // See also Morrison and Grabowski (2007, JAS, Eq. (2))
            nlic   = std::min(nlic, qlic*1.e20);
 
            // set the minimum droplet number as 20/cm3.
            // nlic   = max(nlic,20.e6_r8/rho) ! sghan minimum in #/cm3
            //auto tempnc = nlic/rho/1.0e6;    // #/kg --> #/cm3
 
            // Should be the in-cloud dropelt number calculated as nlic*rho/1.0e6_r8 ????!!!! +++mhwang
            auto pgam   = 0.0005714*(nlic*rho/1.e6) + 0.2714;
            pgam   = 1./std::pow(pgam, 2.0)-1.;
            pgam   = std::min(10., std::max(pgam,2.));
            auto laml = std::pow((pi/6.*rhow*nlic*std::tgamma(pgam+4.)/(qlic*std::tgamma(pgam+1.))), (1./3.));
            auto lammin = (pgam+1.)/50.e-6;    // in cldwat2m, CAM
            auto lammax = (pgam+1.)/2.e-6;     // if lammax is too large, this will lead to crash in
            // src/physics/rrtmg/cloud_rad_props.F90 because
            // klambda-1 can be zero in gam_liquid_lw and gam_liquid_sw
            //  and g_lambda(kmu,klambda-1) will not be defined.
            laml         = std::min(std::max(laml, lammin),lammax);
            effl(i,j)    = std::tgamma(pgam+4.)/std::tgamma(pgam+3.)/laml/2.*1.e6;
            lamcrad(i,j) = laml;
            pgamrad(i,j) = pgam;
        }
        else {
            // we chose 10. over 25, since 10 is a more reasonable value for liquid droplet. +++mhwang
            effl(i,j)    = 10.;     // in cldwat2m, CAM
            lamcrad(i,j) = 0.0;
            pgamrad(i,j) = 0.0;
        }
    });
 } // m2005_effradius

Referenced by Radiation::run().

Here is the caller graph for this function: