ERF_AdvanceWSM6.cpp File Reference

#include "ERF_WSM6.H"
#include "ERF_WSM6_Fortran_Interface.H"
#include "ERF_MicrophysicsConstants.H"
#include "ERF_Constants.H"
#include <AMReX_Reduce.H>
#include <algorithm>
#include <array>
#include <cctype>
#include <cstdio>
#include <cmath>
#include <cstdint>
#include <cstring>
#include <string>
#include <vector>

Include dependency graph for ERF_AdvanceWSM6.cpp:

Namespaces
	WSM6SedCellScratch
	WSM6SedNodeScratch

Enumerations
enum	{   WSM6SedCellScratch::wd = 0 , WSM6SedCellScratch::wa , WSM6SedCellScratch::wa2 , WSM6SedCellScratch::qn ,   WSM6SedCellScratch::qn2 , WSM6SedCellScratch::dz , WSM6SedCellScratch::ww , WSM6SedCellScratch::qq ,   WSM6SedCellScratch::qq2 , WSM6SedCellScratch::was , WSM6SedCellScratch::den , WSM6SedCellScratch::denfac ,   WSM6SedCellScratch::tk , WSM6SedCellScratch::qr , WSM6SedCellScratch::qr2 , WSM6SedCellScratch::tmp ,   WSM6SedCellScratch::tmp1 , WSM6SedCellScratch::tmp2 , WSM6SedCellScratch::tmp3 , WSM6SedCellScratch::workr_col ,   WSM6SedCellScratch::worka_col , WSM6SedCellScratch::denqrs1_col , WSM6SedCellScratch::denqrs2_col , WSM6SedCellScratch::denqrs3_col ,   WSM6SedCellScratch::qsum_col , WSM6SedCellScratch::work1c_col , WSM6SedCellScratch::denqci_col , WSM6SedCellScratch::NumComps }
enum	{   WSM6SedNodeScratch::wi = 0 , WSM6SedNodeScratch::zi , WSM6SedNodeScratch::za , WSM6SedNodeScratch::dza ,   WSM6SedNodeScratch::qa , WSM6SedNodeScratch::qa2 , WSM6SedNodeScratch::qmi , WSM6SedNodeScratch::qpi ,   WSM6SedNodeScratch::NumComps }

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_cpmcal (Real x, Real qmin_arg, Real cpd_arg, Real cpv_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_xlcal (Real x, Real xlv0_arg, Real xlv1_arg, Real t0c_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_diffus (Real x, Real y)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_viscos (Real x, Real y)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_xka (Real x, Real y)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_diffac (Real a, Real b, Real c, Real d, Real e, Real rv_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_venfac (Real a, Real b, Real c, Real den0_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_conden (Real a, Real b, Real c, Real d, Real e, Real qmin_arg, Real rv_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_lamdar (Real x, Real y, Real pidn0r_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_lamdas (Real x, Real y, Real z, Real pidn0s_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wsm6_lamdag (Real x, Real y, Real pidn0g_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wsm6_slope_rain_cell (Real qr, Real den, Real denfac, Real pidn0r_arg, Real qcrmin_arg, Real rslopermax_arg, Real rsloperbmax_arg, Real rsloper2max_arg, Real rsloper3max_arg, Real bvtr_arg, Real pvtr_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wsm6_slope_snow_cell (Real qs, Real den, Real denfac, Real t, Real pidn0s_arg, Real alpha_arg, Real n0smax_arg, Real n0s_arg, Real t0c_arg, Real qcrmin_arg, Real rslopesmax_arg, Real rslopesbmax_arg, Real rslopes2max_arg, Real rslopes3max_arg, Real bvts_arg, Real pvts_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt, Real &n0sfac)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wsm6_slope_graup_cell (Real qg, Real den, Real denfac, Real pidn0g_arg, Real qcrmin_arg, Real rslopegmax_arg, Real rslopegbmax_arg, Real rslopeg2max_arg, Real rslopeg3max_arg, Real bvtg_arg, Real pvtg_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wsm6_nislfv_rain_plm_scratch (int km, int ww_comp, int rq_comp, Real *precip, Real dt, int iter, Array4< Real > const &sed_cell, Array4< Real > const &sed_node, int i_s, int j_s, int klo_s)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wsm6_nislfv_rain_plm6_scratch (int km, int ww_comp, int rq_comp, int rq2_comp, Real *precip1, Real *precip2, Real dt, int iter, Array4< Real > const &sed_cell, Array4< Real > const &sed_node, int i_s, int j_s, int klo_s)

Function Documentation

wsm6_conden()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_conden	(	Real	a,
		Real	b,
		Real	c,
		Real	d,
		Real	e,
		Real	qmin_arg,
		Real	rv_arg
	)

                                               {
    return (amrex::max(b,qmin_arg)-c)
          /(Real(1.0)+d*d/(rv_arg*e)*c/(a*a));
}

Referenced by WSM6::AdvanceColumns().

Here is the caller graph for this function:

wsm6_cpmcal()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_cpmcal	(	Real	x,
		Real	qmin_arg,
		Real	cpd_arg,
		Real	cpv_arg
	)

                                                                     {
    return cpd_arg*(Real(1.0)-amrex::max(x,qmin_arg))
          +amrex::max(x,qmin_arg)*cpv_arg;
}

Referenced by WSM6::AdvanceColumns().

Here is the caller graph for this function:

wsm6_diffac()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_diffac	(	Real	a,
		Real	b,
		Real	c,
		Real	d,
		Real	e,
		Real	rv_arg
	)

                                {
    return d*a*a/(wsm6_xka(c,d)*rv_arg*c*c)
          +Real(1.0)/(e*wsm6_diffus(c,b));
}

Referenced by WSM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_diffus()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_diffus	(	Real	x,
		Real	y
	)

                                  {
    return Real(8.794e-5)*std::exp(std::log(x)*Real(1.81))/y;
}

Referenced by wsm6_diffac(), and wsm6_venfac().

Here is the caller graph for this function:

wsm6_lamdag()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_lamdag	(	Real	x,
		Real	y,
		Real	pidn0g_arg
	)

                                                   {
    return std::sqrt(std::sqrt(pidn0g_arg/(x*y)));
}

Referenced by wsm6_slope_graup_cell().

Here is the caller graph for this function:

wsm6_lamdar()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_lamdar	(	Real	x,
		Real	y,
		Real	pidn0r_arg
	)

                                                   {
    return std::sqrt(std::sqrt(pidn0r_arg/(x*y)));
}

Referenced by wsm6_slope_rain_cell().

Here is the caller graph for this function:

wsm6_lamdas()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_lamdas	(	Real	x,
		Real	y,
		Real	z,
		Real	pidn0s_arg
	)

                                                           {
    return std::sqrt(std::sqrt(pidn0s_arg*z/(x*y)));
}

Referenced by wsm6_slope_snow_cell().

Here is the caller graph for this function:

wsm6_nislfv_rain_plm6_scratch()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_nislfv_rain_plm6_scratch	(	int	km,
		int	ww_comp,
		int	rq_comp,
		int	rq2_comp,
		Real *	precip1,
		Real *	precip2,
		Real	dt,
		int	iter,
		Array4< Real > const &	sed_cell,
		Array4< Real > const &	sed_node,
		int	i_s,
		int	j_s,
		int	klo_s
	)

{
    auto DZ = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::dz);
    };
    auto WW = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::ww);
    };
    auto QQ = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qq);
    };
    auto QQ2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qq2);
    };
    auto WD = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::wd);
    };
    auto WA = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::wa);
    };
    auto WA2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::wa2);
    };
    auto WAS = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::was);
    };
    auto DEN = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::den);
    };
    auto DENFAC = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::denfac);
    };
    auto TK = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tk);
    };
    auto QN = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qn);
    };
    auto QN2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qn2);
    };
    auto QR = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qr);
    };
    auto QR2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qr2);
    };
    auto TMP = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp);
    };
    auto TMP1 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp1);
    };
    auto TMP2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp2);
    };
    auto TMP3 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp3);
    };
 
    auto WI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::wi);
    };
    auto ZI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::zi);
    };
    auto ZA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::za);
    };
    auto DZA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::dza);
    };
    auto QA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qa);
    };
    auto QA2 = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qa2);
    };
    auto QMI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qmi);
    };
    auto QPI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qpi);
    };
 
    constexpr Real pi = Real(3.141592653589793238462643383279502884);
    auto rgmma = [](Real x) -> Real {
        if (x == Real(1.0)) return Real(0.0);
        constexpr Real euler = Real(0.577215664901532);
        Real rg = x * std::exp(euler * x);
        for (int ii = 1; ii <= 10000; ++ii) {
            const Real y = static_cast<Real>(ii);
            rg = rg * (Real(1.0) + x / y) * std::exp(-x / y);
        }
        return Real(1.0) / rg;
    };
 
    const Real dens = WSM6::dens_snow;
    const Real n0g = Real(4.0e6);
    const Real deng = Real(500.0);
    const Real avtg = Real(330.0);
    const Real bvtg = Real(0.8);
    const Real lamdagmax = Real(6.0e4);
 
    const Real pidn0s = pi * dens * WSM6::n0s;
    const Real pidn0g = pi * deng * n0g;
    const Real rslopesmax = Real(1.0) / WSM6::lamdasmax;
    const Real rslopesbmax = std::pow(rslopesmax, WSM6::bvts);
    const Real rslopes2max = rslopesmax * rslopesmax;
    const Real rslopes3max = rslopes2max * rslopesmax;
    const Real rslopegmax = Real(1.0) / lamdagmax;
    const Real rslopegbmax = std::pow(rslopegmax, bvtg);
    const Real rslopeg2max = rslopegmax * rslopegmax;
    const Real rslopeg3max = rslopeg2max * rslopegmax;
    const Real pvts = WSM6::avts * rgmma(Real(4.0) + WSM6::bvts) / Real(6.0);
    const Real pvtg = avtg * rgmma(Real(4.0) + bvtg) / Real(6.0);
 
    Real allold = Real(0.0);
    for (int k = 0; k < km; ++k) {
        QQ(k) = sed_cell(i_s, j_s, klo_s + k, rq_comp);
        QQ2(k) = sed_cell(i_s, j_s, klo_s + k, rq2_comp);
        WW(k) = sed_cell(i_s, j_s, klo_s + k, ww_comp);
        WD(k) = WW(k);
        allold += QQ(k) + QQ2(k);
    }
 
    precip1[0] = Real(0.0);
    precip2[0] = Real(0.0);
    if (allold <= Real(0.0)) {
        return;
    }
 
    ZI(0) = Real(0.0);
    for (int k = 0; k < km; ++k) {
        ZI(k + 1) = ZI(k) + DZ(k);
    }
 
    auto update_wind_and_state = [&](void) {
        WI(0) = WW(0);
        WI(km) = WW(km - 1);
        for (int k = 1; k < km; ++k) {
            WI(k) = (WW(k) * DZ(k - 1) + WW(k - 1) * DZ(k)) / (DZ(k - 1) + DZ(k));
        }
 
        WI(0) = WW(0);
        WI(1) = Real(0.5) * (WW(1) + WW(0));
        for (int k = 2; k < km - 1; ++k) {
            WI(k) = Real(9.0) / Real(16.0) * (WW(k) + WW(k - 1))
                  - Real(1.0) / Real(16.0) * (WW(k + 1) + WW(k - 2));
        }
        if (km > 1) {
            WI(km - 1) = Real(0.5) * (WW(km - 1) + WW(km - 2));
            WI(km) = WW(km - 1);
        }
 
        for (int k = 1; k < km; ++k) {
            if (WW(k) == Real(0.0)) WI(k) = WW(k - 1);
        }
 
        const Real con1 = Real(0.05);
        for (int k = km - 1; k >= 0; --k) {
            const Real decfl = (WI(k + 1) - WI(k)) * dt / DZ(k);
            if (decfl > con1) {
                WI(k) = WI(k + 1) - con1 * DZ(k) / dt;
            }
        }
 
        for (int k = 0; k <= km; ++k) {
            ZA(k) = ZI(k) - WI(k) * dt;
        }
 
        for (int k = 0; k < km; ++k) {
            DZA(k) = ZA(k + 1) - ZA(k);
            if (DZA(k) <= Real(0.0)) DZA(k) = DZ(k);
        }
        DZA(km) = ZI(km) - ZA(km);
        if (DZA(km) <= Real(0.0)) DZA(km) = DZ(km > 0 ? km - 1 : 0);
        for (int k = 0; k < km; ++k) {
            QA(k) = QQ(k) * DZ(k) / DZA(k);
            QA2(k) = QQ2(k) * DZ(k) / DZA(k);
            QR(k) = QA(k) / DEN(k);
            QR2(k) = QA2(k) / DEN(k);
        }
        QA(km) = Real(0.0);
        QA2(km) = Real(0.0);
    };
 
    update_wind_and_state();
 
    if (iter > 0) {
        Real n0sfac_dummy;
        for (int k = 0; k < km; ++k) {
            wsm6_slope_snow_cell(QR(k), DEN(k), DENFAC(k), TK(k), pidn0s,
                                 Real(0.12), Real(1.0e11), Real(2.0e6),
                                 Real(273.15), Real(WSM6::qcrmin),
                                 rslopesmax, rslopesbmax, rslopes2max,
                                 rslopes3max, WSM6::bvts, pvts,
                                 TMP(k), TMP1(k), TMP2(k), TMP3(k),
                                 WA(k), n0sfac_dummy);
            wsm6_slope_graup_cell(QR2(k), DEN(k), DENFAC(k), pidn0g,
                                  Real(WSM6::qcrmin), rslopegmax,
                                  rslopegbmax, rslopeg2max, rslopeg3max,
                                  Real(0.8), pvtg,
                                  TMP(k), TMP1(k), TMP2(k), TMP3(k),
                                  WA2(k));
        }
        for (int k = 0; k < km; ++k) {
            const Real tmpq = amrex::max(QR(k) + QR2(k), Real(1.0e-15));
            if (tmpq > Real(1.0e-15)) {
                WA(k) = (WA(k) * QR(k) + WA2(k) * QR2(k)) / tmpq;
            } else {
                WA(k) = Real(0.0);
            }
        }
        for (int k = 0; k < km; ++k) {
            WW(k) = Real(0.5) * (WD(k) + WA(k));
            WAS(k) = WA(k);
        }
        update_wind_and_state();
    }
 
    for (int ist = 0; ist < 2; ++ist) {
        const int qn_comp = (ist == 0)
            ? WSM6SedCellScratch::qn
            : WSM6SedCellScratch::qn2;
        const int qa_comp = (ist == 0)
            ? WSM6SedNodeScratch::qa
            : WSM6SedNodeScratch::qa2;
        auto QN_DST = [&](int k) -> Real& {
            return sed_cell(i_s,j_s,klo_s+k,qn_comp);
        };
        auto QA_SRC = [&](int k) -> Real& {
            return sed_node(i_s,j_s,klo_s+k,qa_comp);
        };
        Real* precip_dst = (ist == 0) ? &precip1[0] : &precip2[0];
 
        for (int k = 1; k < km; ++k) {
            const Real dip = (QA_SRC(k + 1) - QA_SRC(k)) / (DZA(k + 1) + DZA(k));
            const Real dim = (QA_SRC(k) - QA_SRC(k - 1)) / (DZA(k - 1) + DZA(k));
            if (dip * dim <= Real(0.0)) {
                QMI(k) = QA_SRC(k);
                QPI(k) = QA_SRC(k);
            } else {
                QPI(k) = QA_SRC(k) + Real(0.5) * (dip + dim) * DZA(k);
                QMI(k) = Real(2.0) * QA_SRC(k) - QPI(k);
                if (QPI(k) < Real(0.0) || QMI(k) < Real(0.0)) {
                    QPI(k) = QA_SRC(k);
                    QMI(k) = QA_SRC(k);
                }
            }
        }
        QPI(0) = QA_SRC(0);
        QMI(0) = QA_SRC(0);
        QMI(km) = QA_SRC(km);
        QPI(km) = QA_SRC(km);
 
        for (int k = 0; k < km; ++k) {
            QN_DST(k) = Real(0.0);
        }
 
        int kb = 0;
        int kt = 0;
        for (int k = 0; k < km; ++k) {
            if (ZI(k) >= ZA(km)) {
                break;
            }
 
            for (int kk = kb; kk < km; ++kk) {
                if (ZI(k) <= ZA(kk + 1)) {
                    kb = kk;
                    break;
                }
            }
 
            // Unlike the WRF original, also try the departure point of the column top, and when
            // the cell top lies above even that, take the rest from the empty region above it
            // (QA(km) = 0). Otherwise the search fails for the top cells of the column, which
            // then come out empty and lose what they held: harmless at the model top, where
            // there is nothing, but not at the top of a refined patch.
            bool kt_found = false;
            for (int kk = kt; kk <= km; ++kk) {
                if (ZI(k + 1) <= ZA(kk)) {
                    kt = kk;
                    kt_found = true;
                    break;
                }
            }
            if (!kt_found) { kt = km + 1; }
            kt = amrex::max(kt - 1, 0);
 
            if (kt == kb) {
                const Real tl = (ZI(k) - ZA(kb)) / DZA(kb);
                const Real th = (ZI(k + 1) - ZA(kb)) / DZA(kb);
                const Real tl2 = tl * tl;
                const Real th2 = th * th;
                const Real qqd = Real(0.5) * (QPI(kb) - QMI(kb));
                const Real qqh = qqd * th2 + QMI(kb) * th;
                const Real qql = qqd * tl2 + QMI(kb) * tl;
                QN_DST(k) = (qqh - qql) / (th - tl);
            } else if (kt > kb) {
                const Real tl = (ZI(k) - ZA(kb)) / DZA(kb);
                const Real tl2 = tl * tl;
                const Real qqd = Real(0.5) * (QPI(kb) - QMI(kb));
                const Real qql = qqd * tl2 + QMI(kb) * tl;
                const Real dql = QA_SRC(kb) - qql;
                Real zsum = (Real(1.0) - tl) * DZA(kb);
                Real qsum = dql * DZA(kb);
                if (kt - kb > 1) {
                    for (int m = kb + 1; m < kt; ++m) {
                        zsum += DZA(m);
                        qsum += QA_SRC(m) * DZA(m);
                    }
                }
                const Real th = (ZI(k + 1) - ZA(kt)) / DZA(kt);
                const Real th2 = th * th;
                const Real dqh = Real(0.5) * (QPI(kt) - QMI(kt)) * th2 + QMI(kt) * th;
                zsum += th * DZA(kt);
                qsum += dqh * DZA(kt);
                QN_DST(k) = qsum / zsum;
            }
        }
 
        Real precip = Real(0.0);
        for (int k = 0; k < km; ++k) {
            if (ZA(k) < Real(0.0) && ZA(k + 1) < Real(0.0)) {
                precip += QA_SRC(k) * DZA(k);
            } else if (ZA(k) < Real(0.0) && ZA(k + 1) >= Real(0.0)) {
                precip += QA_SRC(k) * (Real(0.0) - ZA(k));
                break;
            } else {
                break;
            }
        }
        *precip_dst = precip;
    }
 
    for (int k = 0; k < km; ++k) {
        sed_cell(i_s, j_s, klo_s + k, rq_comp) = QN(k);
        sed_cell(i_s, j_s, klo_s + k, rq2_comp) = QN2(k);
        sed_cell(i_s, j_s, klo_s + k, ww_comp) = WW(k);
    }
}

Referenced by WSM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_nislfv_rain_plm_scratch()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_nislfv_rain_plm_scratch	(	int	km,
		int	ww_comp,
		int	rq_comp,
		Real *	precip,
		Real	dt,
		int	iter,
		Array4< Real > const &	sed_cell,
		Array4< Real > const &	sed_node,
		int	i_s,
		int	j_s,
		int	klo_s
	)

{
    auto DZ = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::dz);
    };
    auto WW = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::ww);
    };
    auto QQ = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qq);
    };
    auto WD = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::wd);
    };
    auto WA = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::wa);
    };
    auto WAS = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::was);
    };
    auto DEN = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::den);
    };
    auto DENFAC = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::denfac);
    };
    auto QN = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qn);
    };
    auto QR = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::qr);
    };
    auto TMP = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp);
    };
    auto TMP1 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp1);
    };
    auto TMP2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp2);
    };
    auto TMP3 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WSM6SedCellScratch::tmp3);
    };
 
    auto WI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::wi);
    };
    auto ZI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::zi);
    };
    auto ZA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::za);
    };
    auto DZA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::dza);
    };
    auto QA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qa);
    };
    auto QMI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qmi);
    };
    auto QPI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WSM6SedNodeScratch::qpi);
    };
 
    constexpr Real pi = Real(3.141592653589793238462643383279502884);
    auto rgmma = [](Real x) -> Real {
        if (x == Real(1.0)) return Real(0.0);
        constexpr Real euler = Real(0.577215664901532);
        Real rg = x * std::exp(euler * x);
        for (int ii = 1; ii <= 10000; ++ii) {
            const Real y = static_cast<Real>(ii);
            rg = rg * (Real(1.0) + x / y) * std::exp(-x / y);
        }
        return Real(1.0) / rg;
    };
 
    const Real pidn0r = pi * Real(rhoh2o) * WSM6::n0r;
    const Real rslopermax = Real(1.0) / WSM6::lamdarmax;
    const Real rsloperbmax = std::pow(rslopermax, WSM6::bvtr);
    const Real rsloper2max = rslopermax * rslopermax;
    const Real rsloper3max = rsloper2max * rslopermax;
    const Real pvtr = WSM6::avtr * rgmma(Real(4.0) + WSM6::bvtr) / Real(6.0);
 
    Real allold = Real(0.0);
    for (int k = 0; k < km; ++k) {
        QQ(k) = sed_cell(i_s, j_s, klo_s + k, rq_comp);
        WW(k) = sed_cell(i_s, j_s, klo_s + k, ww_comp);
        WD(k) = WW(k);
        allold += QQ(k);
    }
 
    precip[0] = Real(0.0);
    if (allold <= Real(0.0)) {
        return;
    }
 
    ZI(0) = Real(0.0);
    for (int k = 0; k < km; ++k) {
        ZI(k + 1) = ZI(k) + DZ(k);
    }
 
    auto update_wind_and_state = [&](void) {
        WI(0) = WW(0);
        WI(km) = WW(km - 1);
        for (int k = 1; k < km; ++k) {
            WI(k) = (WW(k) * DZ(k - 1) + WW(k - 1) * DZ(k)) / (DZ(k - 1) + DZ(k));
        }
 
        WI(0) = WW(0);
        WI(1) = Real(0.5) * (WW(1) + WW(0));
        for (int k = 2; k < km - 1; ++k) {
            WI(k) = Real(9.0) / Real(16.0) * (WW(k) + WW(k - 1))
                  - Real(1.0) / Real(16.0) * (WW(k + 1) + WW(k - 2));
        }
        if (km > 1) {
            WI(km - 1) = Real(0.5) * (WW(km - 1) + WW(km - 2));
            WI(km) = WW(km - 1);
        }
 
        for (int k = 1; k < km; ++k) {
            if (WW(k) == Real(0.0)) WI(k) = WW(k - 1);
        }
 
        const Real con1 = Real(0.05);
        for (int k = km - 1; k >= 0; --k) {
            const Real decfl = (WI(k + 1) - WI(k)) * dt / DZ(k);
            if (decfl > con1) {
                WI(k) = WI(k + 1) - con1 * DZ(k) / dt;
            }
        }
 
        for (int k = 0; k <= km; ++k) {
            ZA(k) = ZI(k) - WI(k) * dt;
        }
 
        for (int k = 0; k < km; ++k) {
            DZA(k) = ZA(k + 1) - ZA(k);
            if (DZA(k) <= Real(0.0)) DZA(k) = DZ(k);
        }
        DZA(km) = ZI(km) - ZA(km);
        if (DZA(km) <= Real(0.0)) DZA(km) = DZ(km > 0 ? km - 1 : 0);
        for (int k = 0; k < km; ++k) {
            QA(k) = QQ(k) * DZ(k) / DZA(k);
            QR(k) = QA(k) / DEN(k);
        }
        QA(km) = Real(0.0);
    };
 
    update_wind_and_state();
 
    if (iter > 0) {
        for (int k = 0; k < km; ++k) {
            wsm6_slope_rain_cell(QR(k), DEN(k), DENFAC(k), pidn0r,
                                 WSM6::qcrmin,
                                 rslopermax, rsloperbmax, rsloper2max,
                                 rsloper3max, WSM6::bvtr, pvtr,
                                 TMP(k), TMP1(k), TMP2(k), TMP3(k), WA(k));
        }
        for (int k = 0; k < km; ++k) {
            WW(k) = Real(0.5) * (WD(k) + WA(k));
            WAS(k) = WA(k);
        }
        update_wind_and_state();
    }
 
    for (int k = 1; k < km; ++k) {
        const Real dip = (QA(k + 1) - QA(k)) / (DZA(k + 1) + DZA(k));
        const Real dim = (QA(k) - QA(k - 1)) / (DZA(k - 1) + DZA(k));
        if (dip * dim <= Real(0.0)) {
            QMI(k) = QA(k);
            QPI(k) = QA(k);
        } else {
            QPI(k) = QA(k) + Real(0.5) * (dip + dim) * DZA(k);
            QMI(k) = Real(2.0) * QA(k) - QPI(k);
            if (QPI(k) < Real(0.0) || QMI(k) < Real(0.0)) {
                QPI(k) = QA(k);
                QMI(k) = QA(k);
            }
        }
    }
    QPI(0) = QA(0);
    QMI(0) = QA(0);
    QMI(km) = QA(km);
    QPI(km) = QA(km);
 
    for (int k = 0; k < km; ++k) {
        QN(k) = Real(0.0);
    }
 
    int kb = 0;
    int kt = 0;
    for (int k = 0; k < km; ++k) {
        kb = amrex::max(kb - 1, 0);
        kt = amrex::max(kt - 1, 0);
 
        if (ZI(k) >= ZA(km)) {
            break;
        }
 
        for (int kk = kb; kk < km; ++kk) {
            if (ZI(k) <= ZA(kk + 1)) {
                kb = kk;
                break;
            }
        }
 
        // Unlike the WRF original, also try the departure point of the column top, and when
        // the cell top lies above even that, take the rest from the empty region above it
        // (QA(km) = 0). Otherwise the search fails for the top cells of the column, which
        // then come out empty and lose what they held: harmless at the model top, where
        // there is nothing, but not at the top of a refined patch.
        bool kt_found = false;
        for (int kk = kt; kk <= km; ++kk) {
            if (ZI(k + 1) <= ZA(kk)) {
                kt = kk;
                kt_found = true;
                break;
            }
        }
        if (!kt_found) { kt = km + 1; }
        kt = amrex::max(kt - 1, 0);
 
        if (kt == kb) {
            const Real tl = (ZI(k) - ZA(kb)) / DZA(kb);
            const Real th = (ZI(k + 1) - ZA(kb)) / DZA(kb);
            const Real tl2 = tl * tl;
            const Real th2 = th * th;
            const Real qqd = Real(0.5) * (QPI(kb) - QMI(kb));
            const Real qqh = qqd * th2 + QMI(kb) * th;
            const Real qql = qqd * tl2 + QMI(kb) * tl;
            QN(k) = (qqh - qql) / (th - tl);
        } else if (kt > kb) {
            const Real tl = (ZI(k) - ZA(kb)) / DZA(kb);
            const Real tl2 = tl * tl;
            const Real qqd = Real(0.5) * (QPI(kb) - QMI(kb));
            const Real qql = qqd * tl2 + QMI(kb) * tl;
            const Real dql = QA(kb) - qql;
            Real zsum = (Real(1.0) - tl) * DZA(kb);
            Real qsum = dql * DZA(kb);
            if (kt - kb > 1) {
                for (int m = kb + 1; m < kt; ++m) {
                    zsum += DZA(m);
                    qsum += QA(m) * DZA(m);
                }
            }
            const Real th = (ZI(k + 1) - ZA(kt)) / DZA(kt);
            const Real th2 = th * th;
            const Real dqh = Real(0.5) * (QPI(kt) - QMI(kt)) * th2 + QMI(kt) * th;
            zsum += th * DZA(kt);
            qsum += dqh * DZA(kt);
            QN(k) = qsum / zsum;
        }
    }
 
    for (int k = 0; k < km; ++k) {
        if (ZA(k) < Real(0.0) && ZA(k + 1) < Real(0.0)) {
            precip[0] += QA(k) * DZA(k);
        } else if (ZA(k) < Real(0.0) && ZA(k + 1) >= Real(0.0)) {
            precip[0] += QA(k) * (Real(0.0) - ZA(k));
            break;
        } else {
            break;
        }
    }
 
    for (int k = 0; k < km; ++k) {
        sed_cell(i_s, j_s, klo_s + k, rq_comp) = QN(k);
        sed_cell(i_s, j_s, klo_s + k, ww_comp) = WW(k);
    }
}

Referenced by WSM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_slope_graup_cell()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_slope_graup_cell	(	Real	qg,
		Real	den,
		Real	denfac,
		Real	pidn0g_arg,
		Real	qcrmin_arg,
		Real	rslopegmax_arg,
		Real	rslopegbmax_arg,
		Real	rslopeg2max_arg,
		Real	rslopeg3max_arg,
		Real	bvtg_arg,
		Real	pvtg_arg,
		Real &	rslope,
		Real &	rslopeb,
		Real &	rslope2,
		Real &	rslope3,
		Real &	vt
	)

{
    if (qg <= qcrmin_arg) {
        rslope  = rslopegmax_arg;
        rslopeb = rslopegbmax_arg;
        rslope2 = rslopeg2max_arg;
        rslope3 = rslopeg3max_arg;
    } else {
        rslope  = Real(1.0)/wsm6_lamdag(qg,den,pidn0g_arg);
        rslopeb = std::pow(rslope,bvtg_arg);
        rslope2 = rslope*rslope;
        rslope3 = rslope2*rslope;
    }
    vt = pvtg_arg*rslopeb*denfac;
    if (qg <= Real(0.0)) vt = Real(0.0);
}

Referenced by WSM6::AdvanceColumns(), and wsm6_nislfv_rain_plm6_scratch().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_slope_rain_cell()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_slope_rain_cell	(	Real	qr,
		Real	den,
		Real	denfac,
		Real	pidn0r_arg,
		Real	qcrmin_arg,
		Real	rslopermax_arg,
		Real	rsloperbmax_arg,
		Real	rsloper2max_arg,
		Real	rsloper3max_arg,
		Real	bvtr_arg,
		Real	pvtr_arg,
		Real &	rslope,
		Real &	rslopeb,
		Real &	rslope2,
		Real &	rslope3,
		Real &	vt
	)

{
    if (qr <= qcrmin_arg) {
        rslope  = rslopermax_arg;
        rslopeb = rsloperbmax_arg;
        rslope2 = rsloper2max_arg;
        rslope3 = rsloper3max_arg;
    } else {
        rslope  = Real(1.0)/wsm6_lamdar(qr,den,pidn0r_arg);
        rslopeb = std::pow(rslope,bvtr_arg);
        rslope2 = rslope*rslope;
        rslope3 = rslope2*rslope;
    }
    vt = pvtr_arg*rslopeb*denfac;
    if (qr <= Real(0.0)) vt = Real(0.0);
}

Referenced by WSM6::AdvanceColumns(), and wsm6_nislfv_rain_plm_scratch().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_slope_snow_cell()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wsm6_slope_snow_cell	(	Real	qs,
		Real	den,
		Real	denfac,
		Real	t,
		Real	pidn0s_arg,
		Real	alpha_arg,
		Real	n0smax_arg,
		Real	n0s_arg,
		Real	t0c_arg,
		Real	qcrmin_arg,
		Real	rslopesmax_arg,
		Real	rslopesbmax_arg,
		Real	rslopes2max_arg,
		Real	rslopes3max_arg,
		Real	bvts_arg,
		Real	pvts_arg,
		Real &	rslope,
		Real &	rslopeb,
		Real &	rslope2,
		Real &	rslope3,
		Real &	vt,
		Real &	n0sfac
	)

{
    Real supcol = t0c_arg - t;
    n0sfac = amrex::max(amrex::min(std::exp(alpha_arg*supcol),
                                    n0smax_arg/n0s_arg), Real(1.0));
    if (qs <= qcrmin_arg) {
        rslope  = rslopesmax_arg;
        rslopeb = rslopesbmax_arg;
        rslope2 = rslopes2max_arg;
        rslope3 = rslopes3max_arg;
    } else {
        rslope  = Real(1.0)/wsm6_lamdas(qs,den,n0sfac,pidn0s_arg);
        rslopeb = std::pow(rslope,bvts_arg);
        rslope2 = rslope*rslope;
        rslope3 = rslope2*rslope;
    }
    vt = pvts_arg*rslopeb*denfac;
    if (qs <= Real(0.0)) vt = Real(0.0);
}

Referenced by WSM6::AdvanceColumns(), and wsm6_nislfv_rain_plm6_scratch().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_venfac()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_venfac	(	Real	a,
		Real	b,
		Real	c,
		Real	den0_arg
	)

                                                         {
    return std::exp(std::log(wsm6_viscos(b,c)/wsm6_diffus(b,a))
                   *Real(0.3333333))
          /std::sqrt(wsm6_viscos(b,c))
          *std::sqrt(std::sqrt(den0_arg/c));
}

Referenced by WSM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_viscos()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_viscos	(	Real	x,
		Real	y
	)

                                  {
    return Real(1.496e-6)*(x*std::sqrt(x))/(x+Real(120.0))/y;
}

Referenced by wsm6_venfac(), and wsm6_xka().

Here is the caller graph for this function:

wsm6_xka()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_xka	(	Real	x,
		Real	y
	)

                               {
    return Real(1.414e3)*wsm6_viscos(x,y)*y;
}

Referenced by WSM6::AdvanceColumns(), and wsm6_diffac().

Here is the call graph for this function:

Here is the caller graph for this function:

wsm6_xlcal()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wsm6_xlcal	(	Real	x,
		Real	xlv0_arg,
		Real	xlv1_arg,
		Real	t0c_arg
	)

                                                                     {
    return xlv0_arg - xlv1_arg*(x - t0c_arg);
}

Referenced by WSM6::AdvanceColumns().

Here is the caller graph for this function: