ERF_AdvanceWDM6.cpp File Reference

#include "ERF_WDM6.H"
#include "ERF_MicrophysicsConstants.H"
#include "ERF_Constants.H"
#include <AMReX_IArrayBox.H>
#include <AMReX_Reduce.H>
#include <algorithm>
#include <array>
#include <cctype>
#include <cstdio>
#include <cmath>
#include <limits>
#include <cstdint>
#include <cstring>
#include <string>
#include <vector>

Include dependency graph for ERF_AdvanceWDM6.cpp:

Namespaces
	WDM6SedCellScratch
	WDM6SedNodeScratch

Enumerations
enum	{   WDM6SedCellScratch::wd , WDM6SedCellScratch::wa , WDM6SedCellScratch::wa2 , WDM6SedCellScratch::was ,   WDM6SedCellScratch::qn , WDM6SedCellScratch::qn2 , WDM6SedCellScratch::qq , WDM6SedCellScratch::qq2 ,   WDM6SedCellScratch::qr , WDM6SedCellScratch::qr2 , WDM6SedCellScratch::dz , WDM6SedCellScratch::den ,   WDM6SedCellScratch::denfac , WDM6SedCellScratch::tk , WDM6SedCellScratch::work_col , WDM6SedCellScratch::rq_col ,   WDM6SedCellScratch::rq2_col , WDM6SedCellScratch::NumComps }
enum	{   WDM6SedNodeScratch::wi , WDM6SedNodeScratch::zi , WDM6SedNodeScratch::za , WDM6SedNodeScratch::dza ,   WDM6SedNodeScratch::qa , WDM6SedNodeScratch::qa2 , WDM6SedNodeScratch::qmi , WDM6SedNodeScratch::qpi ,   WDM6SedNodeScratch::NumComps }

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_cpmcal (Real x, Real qmin_arg, Real cpd_arg, Real cpv_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_xlcal (Real x, Real xlv0_arg, Real xlv1_arg, Real t0c_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_default_real_pow (double base, double exponent)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_diffus (Real x, Real y)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_viscos (Real x, Real y)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_xka (Real x, Real y)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_diffac (Real a, Real b, Real c, Real d, Real e, Real rv_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_venfac (Real a, Real b, Real c, Real den0_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_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	wdm6_lamdar (Real qr, Real den, Real nr, Real pidnr_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_lamdac_exact (Real qc, Real den, Real nc, Real pidnc_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_rslopec_exact (Real qc, Real den, Real nc, Real pidnc_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_xni_exact (Real qi, Real den, Real qmin_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wdm6_slope_rain_cell (Real qr, Real nr, Real den, Real denfac, Real qcrmin_arg, Real nrmin_arg, Real rslopermax_arg, Real rsloperbmax_arg, Real rsloper2max_arg, Real rsloper3max_arg, Real bvtr_arg, Real pvtr_arg, Real pvtrn_arg, Real pidnr_arg, Real &rslope, Real &rslopeb, Real &rslope2, Real &rslope3, Real &vt, Real &vtn)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_lamdas (Real x, Real y, Real z, Real pidn0s_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	wdm6_lamdag (Real x, Real y, Real pidn0g_arg)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wdm6_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	wdm6_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	wdm6_nislfv_rain_plm6_column (int km, Real *precip1, Real *precip2, Real dt, int iter, Real pidn0s, Real pidn0g, Real qcrmin, Real alpha, Real n0smax, Real n0s, Real t0c, Real rslopesmax, Real rslopesbmax, Real rslopes2max, Real rslopes3max, Real bvts, Real pvts, Real rslopegmax, Real rslopegbmax, Real rslopeg2max, Real rslopeg3max, Real bvtg, Real pvtg, 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 Real	wdm6_mean_droplet_diameter (Real qc, Real nc, Real den, Real)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void	wdm6_ccn_activation (Real &nc, Real &nn, const Real qv, const Real qc, const Real qvs, const Real, const Real w, const Real, const Real, const Real)

Variables
constexpr amrex::Real	wdm6_slope_t0c = wdm6_literal(273.15)

Function Documentation

wdm6_ccn_activation()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_ccn_activation	(	Real &	nc,
		Real &	nn,
		const Real	qv,
		const Real	qc,
		const Real	qvs,
		const	Real,
		const Real	w,
		const	Real,
		const	Real,
		const	Real
	)

  {
    // Only activate if supersaturated and cloud exists
    if (qv <= qvs || qc < Real(1.e-8)) return;
 
    // Supersaturation
    Real supersaturation = qv / qvs - Real(1.0);
    supersaturation = amrex::min(supersaturation, Real(0.0048));  // satmax - 1
 
    // Critical supersaturation based on updraft strength
    // From WRF: stronger updrafts -> more activation
    Real s_crit = Real(0.6) * std::pow(amrex::max(w, Real(0.01)), Real(0.5));
 
    if (supersaturation > s_crit && nn > Real(0.0)) {
        // Activation rate: convert aerosols to droplets
        Real nc_activate = amrex::min(nn, Real(0.1) * (supersaturation - s_crit) / s_crit * nn);
 
        // Update concentrations
        nc += nc_activate;
        nn -= nc_activate;
 
        // Enforce physical bounds
        nc = amrex::max(nc, Real(1.e1));
        nn = amrex::max(nn, Real(0.0));
    }
}

wdm6_conden()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_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 WDM6::AdvanceColumns().

Here is the caller graph for this function:

wdm6_cpmcal()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_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 WDM6::AdvanceColumns().

Here is the caller graph for this function:

wdm6_default_real_pow()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_default_real_pow	(	double	base,
		double	exponent
	)

                                                          {
#ifdef ERF_WDM6_F32_LITERALS
    return Real(std::pow(static_cast<float>(base), static_cast<float>(exponent)));
#else
    return Real(std::pow(base, exponent));
#endif
}

Referenced by WDM6::AdvanceColumns().

Here is the caller graph for this function:

wdm6_diffac()

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

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

Referenced by WDM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_diffus()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_diffus	(	Real	x,
		Real	y
	)

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

Referenced by wdm6_diffac(), and wdm6_venfac().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_lamdac_exact()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_lamdac_exact	(	Real	qc,
		Real	den,
		Real	nc,
		Real	pidnc_arg
	)

                                                                    {
    // Fortran statement function, ERF_module_mp_wdm6.F90:562:
    //     lamdac(x,y,z) = exp(log(((pidnc*z)/(x*y)))*((.33333333)))
    // called as lamdac(qci(i,k,1), den(i,k), ncr(i,k,2)), so x=qc, y=den, z=nc.
    // The exponent literal is unsuffixed in the Fortran and therefore obeys the
    // LITERAL PRECISION CONTRACT; it is not 1/3.
    return std::exp(std::log((pidnc_arg * nc) / (qc * den)) * wdm6_literal(0.33333333));
}

Referenced by wdm6_rslopec_exact().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_lamdag()

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

                                                   {
    // Graupel slope parameter (single-moment)
    return std::sqrt(std::sqrt(pidn0g_arg/(x*y)));
}

Referenced by wdm6_slope_graup_cell().

Here is the caller graph for this function:

wdm6_lamdar()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_lamdar	(	Real	qr,
		Real	den,
		Real	nr,
		Real	pidnr_arg
	)

                                                              {
    // Fortran statement function, ERF_module_mp_wdm6.F90:3277 and :3357:
    //     lamdar(x,y,z) = exp(log(((pidnr*z)/(x*y)))*((.33333333)))
    // called as lamdar(qrs(i,k,1), den(i,k), ncr(i,k)), so x=qr, y=den, z=nr:
    // the denominator is qr*den, NOT qr alone.
    //
    // This carried the identical defect already fixed in wdm6_lamdac: it wrote
    // (pidnr*nr*den)/(den*qr), in which den cancels, dropping the /den entirely
    // and inflating rslope by den^(-1/3). Measured at (199,3,43): the rain slope
    // was 3.098e-05 native against 2.833e-05 on the bridge, a ratio of 1.0936
    // against den = 0.765, matching den^(-1/3) exactly. That pushed avedia
    // across the di82 threshold in G16a, 8.9365e-05 native against 8.1717e-05,
    // so the bridge collapsed rain back into cloud (qr and nr to zero, qc and nc
    // incremented) while the native path left the cell untouched.
    //
    // The old qr/nr guard is dropped: every caller already gates on
    // qr <= qcrmin || nr <= nrmin with the same values, and the branch returned
    // 1/5.0e4, a slope, where the caller expects a lambda and inverts it again.
    return std::exp(std::log((pidnr_arg * nr) / (qr * den)) * wdm6_literal(0.33333333));
}

Referenced by wdm6_slope_rain_cell().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_lamdas()

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

                                                           {
    // Snow slope parameter (single-moment)
    return std::sqrt(std::sqrt(pidn0s_arg*z/(x*y)));
}

Referenced by wdm6_slope_snow_cell().

Here is the caller graph for this function:

wdm6_mean_droplet_diameter()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_mean_droplet_diameter	(	Real	qc,
		Real	nc,
		Real	den,
		Real
	)

                                                                     {
    // Volume-weighted mean diameter of cloud droplets
    // Returns diameter in meters
    if (nc < Real(1.e1) || qc < Real(1.e-9)) return Real(0.0);
 
    // Mean mass = rho * qc / nc
    // Mean volume = mean mass / rho_water
    // Mean diameter = (6 * volume / pi)^(1/3)
    Real mean_mass = den * qc / (nc * den);  // kg per droplet
    Real mean_volume = mean_mass / Real(1000.0);  // m^3 per droplet
    Real diameter = std::pow(Real(6.0) * mean_volume / Real(3.14159265359), Real(1.0)/Real(3.0));
 
    return diameter;
}

wdm6_nislfv_rain_plm6_column()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_nislfv_rain_plm6_column	(	int	km,
		Real *	precip1,
		Real *	precip2,
		Real	dt,
		int	iter,
		Real	pidn0s,
		Real	pidn0g,
		Real	qcrmin,
		Real	alpha,
		Real	n0smax,
		Real	n0s,
		Real	t0c,
		Real	rslopesmax,
		Real	rslopesbmax,
		Real	rslopes2max,
		Real	rslopes3max,
		Real	bvts,
		Real	pvts,
		Real	rslopegmax,
		Real	rslopegbmax,
		Real	rslopeg2max,
		Real	rslopeg3max,
		Real	bvtg,
		Real	pvtg,
		Array4< Real > const &	sed_cell,
		Array4< Real > const &	sed_node,
		int	i_s,
		int	j_s,
		int	klo_s
	)

{
    auto WD = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::wd);
    };
    auto WA = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::wa);
    };
    auto WA2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::wa2);
    };
    auto WAS = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::was);
    };
    auto QN = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::qn);
    };
    auto QN2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::qn2);
    };
    auto QQ = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::qq);
    };
    auto QQ2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::qq2);
    };
    auto QR = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::qr);
    };
    auto QR2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::qr2);
    };
    auto DZ = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::dz);
    };
    auto DEN = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::den);
    };
    auto DENFAC = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::denfac);
    };
    auto TK = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::tk);
    };
    auto WW = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::work_col);
    };
    auto RQL = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::rq_col);
    };
    auto RQL2 = [&](int k) -> amrex::Real& {
        return sed_cell(i_s, j_s, klo_s + k, WDM6SedCellScratch::rq2_col);
    };
 
    auto WI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::wi);
    };
    auto ZI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::zi);
    };
    auto ZA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::za);
    };
    auto DZA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::dza);
    };
    auto QA = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::qa);
    };
    auto QA2 = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::qa2);
    };
    auto QMI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::qmi);
    };
    auto QPI = [&](int k) -> amrex::Real& {
        return sed_node(i_s, j_s, klo_s + k, WDM6SedNodeScratch::qpi);
    };
 
    Real allold = Real(0.0);
    for (int k = 0; k < km; ++k) {
        QQ(k) = RQL(k);
        QQ2(k) = RQL2(k);
        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 rebuild = [&]() {
        WI(0) = WW(0);
        if (km > 1) {
            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));
            }
            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);
        }
 
        constexpr Real con1 = wdm6_literal(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);
            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);
        }
        DZA(km) = ZI(km) - ZA(km);
        QA(km) = Real(0.0);
        QA2(km) = Real(0.0);
    };
 
    rebuild();
 
    if (iter > 0) {
        for (int k = 0; k < km; ++k) {
            Real rslope, rslopeb, rslope2, rslope3, vt, n0sfac_dummy;
            wdm6_slope_snow_cell(QR(k), DEN(k), DENFAC(k), TK(k),
                                 pidn0s, alpha, n0smax, n0s, t0c, qcrmin,
                                 rslopesmax, rslopesbmax, rslopes2max, rslopes3max,
                                 bvts, pvts, rslope, rslopeb, rslope2, rslope3, vt,
                                 n0sfac_dummy);
            WA(k) = vt;
            wdm6_slope_graup_cell(QR2(k), DEN(k), DENFAC(k),
                                  pidn0g, qcrmin, rslopegmax, rslopegbmax,
                                  rslopeg2max, rslopeg3max, bvtg, pvtg,
                                  rslope, rslopeb, rslope2, rslope3, vt);
            WA2(k) = vt;
        }
        for (int k = 0; k < km; ++k) {
            const Real tmpq = amrex::max(QR(k) + QR2(k), wdm6_literal(1.0e-15));
            WA(k) = (tmpq > wdm6_literal(1.0e-15))
                ? (WA(k) * QR(k) + WA2(k) * QR2(k)) / tmpq
                : Real(0.0);
            WW(k) = Real(0.5) * (WD(k) + WA(k));
            WAS(k) = WA(k);
        }
        rebuild();
    }
 
    for (int ist = 0; ist < 2; ++ist) {
        const int qn_comp = (ist == 0)
            ? WDM6SedCellScratch::qn
            : WDM6SedCellScratch::qn2;
        const int qa_comp = (ist == 0)
            ? WDM6SedNodeScratch::qa
            : WDM6SedNodeScratch::qa2;
        auto qn_dst = [&](int k) -> amrex::Real& {
            return sed_cell(i_s, j_s, klo_s + k, qn_comp);
        };
        auto qa_src = [&](int k) -> amrex::Real& {
            return sed_node(i_s, j_s, klo_s + k, qa_comp);
        };
        Real* precip_dst = (ist == 0) ? precip1 : precip2;
 
        QPI(0) = qa_src(0);
        QMI(0) = qa_src(0);
        QMI(km) = qa_src(km);
        QPI(km) = qa_src(km);
        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);
                }
            }
        }
 
        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) {
            // Fortran backs both brackets off by one at the TOP of every k
            // before searching forward:
            //     kb=max(kb-1,1) ; kt=max(kt-1,1)
            // Omitting this let the forward search start above the correct
            // bracket, so cells whose departure interval reaches back below the
            // previous k's bracket were remapped from the wrong donor.
            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; }
            // Fortran writes a bare kt = kt - 1 with no floor. Clamping it to 0
            // is not equivalent: when the decrement should drop kt BELOW kb,
            // the Fortran falls through both branches and leaves qn(k) at zero,
            // whereas the clamp can make kt equal kb and take the single-cell
            // branch against the wrong donor. kt is only used as an index
            // inside the two branches, both of which require kt >= kb >= 0, so
            // a negative value here is safe.
            kt = kt - 1;
 
            if (kt == kb) {
                const Real tl = (ZI(k) - ZA(kb)) / DZA(kb);
                const Real th = (ZI(k + 1) - ZA(kb)) / DZA(kb);
                const Real qqd = Real(0.5) * (QPI(kb) - QMI(kb));
                qn_dst(k) = (qqd * th * th + QMI(kb) * th
                           - qqd * tl * tl - QMI(kb) * tl) / (th - tl);
            } else if (kt > kb) {
                const Real tl = (ZI(k) - ZA(kb)) / DZA(kb);
                const Real qqd = Real(0.5) * (QPI(kb) - QMI(kb));
                const Real qql = qqd * tl * tl + QMI(kb) * tl;
                const Real dql = qa_src(kb) - qql;
                Real zsum = (Real(1.0) - tl) * DZA(kb);
                Real qsum = dql * DZA(kb);
                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 dqh = Real(0.5) * (QPI(kt) - QMI(kt)) * th * th + QMI(kt) * th;
                zsum += th * DZA(kt);
                qsum += dqh * DZA(kt);
                qn_dst(k) = qsum / zsum;
            }
        }
 
        precip_dst[0] = Real(0.0);
        for (int k = 0; k < km; ++k) {
            if (ZA(k) < Real(0.0) && ZA(k + 1) < Real(0.0)) {
                precip_dst[0] += qa_src(k) * DZA(k);
            } else if (ZA(k) < Real(0.0) && ZA(k + 1) >= Real(0.0)) {
                precip_dst[0] += qa_src(k) * (Real(0.0) - ZA(k));
                break;
            } else {
                break;
            }
        }
    }
 
    for (int k = 0; k < km; ++k) {
        RQL(k) = QN(k);
        RQL2(k) = QN2(k);
    }
}

Referenced by WDM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_rslopec_exact()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_rslopec_exact	(	Real	qc,
		Real	den,
		Real	nc,
		Real	pidnc_arg
	)

                                                                     {
    // rslopec is the RECIPROCAL of lamdac: the Fortran writes
    //     rslopec(i,k) = 1./lamdac(qci(i,k,1),den(i,k),ncr(i,k,2))
    // at both of its two sites (G3 at :915 and G11 at :1544). Returning lamdac
    // itself here put rslopec off by 1/lamdac^2, about nine orders of magnitude
    // once cloud water exists. It was invisible until cloud water first formed,
    // because with qc <= qmin both legs take the rslopecmax branch instead.
    return Real(1.0) / wdm6_lamdac_exact(qc, den, nc, pidnc_arg);
}

Referenced by WDM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_slope_graup_cell()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_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)/wdm6_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 WDM6::AdvanceColumns(), and wdm6_nislfv_rain_plm6_column().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_slope_rain_cell()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_slope_rain_cell	(	Real	qr,
		Real	nr,
		Real	den,
		Real	denfac,
		Real	qcrmin_arg,
		Real	nrmin_arg,
		Real	rslopermax_arg,
		Real	rsloperbmax_arg,
		Real	rsloper2max_arg,
		Real	rsloper3max_arg,
		Real	bvtr_arg,
		Real	pvtr_arg,
		Real	pvtrn_arg,
		Real	pidnr_arg,
		Real &	rslope,
		Real &	rslopeb,
		Real &	rslope2,
		Real &	rslope3,
		Real &	vt,
		Real &	vtn
	)

{
    if (qr <= qcrmin_arg || nr <= nrmin_arg) {
        rslope  = rslopermax_arg;
        rslopeb = rsloperbmax_arg;
        rslope2 = rsloper2max_arg;
        rslope3 = rsloper3max_arg;
    } else {
        // The rain-slope cap is unsuffixed in the Fortran at BOTH of its sites,
        // slope_wdm6 :3420 and slope_rain :3493:
        //     rslope(i,k,1) = min(1./lamdar(...),1.e-3)
        // so it obeys the LITERAL PRECISION CONTRACT: float32(1.e-3) widens to
        // 0.0010000000474974513, not the exact double 1e-3. The cap therefore
        // sits +4.749745e-08 relative ABOVE the port's, and the gap propagates
        // as a power: rslopeb = rslope**bvtr, rslope2, rslope3, so rslope3
        // carries 3x it at 1.424924e-07.
        //
        // This is a min(), so it binds only where 1./lamdar exceeds the cap --
        // rare, and identically zero for the first fifteen steps of the Bubble
        // case, which is why a bitwise-clean 10-step run did not expose it.
        rslope  = amrex::min(Real(1.0) / wdm6_lamdar(qr, den, nr, pidnr_arg),
                             wdm6_literal(1.e-3));
        rslopeb = std::pow(rslope, bvtr_arg);
        rslope2 = rslope  * rslope;
        rslope3 = rslope2 * rslope;
    }
    vt = pvtr_arg * rslopeb * denfac;
    vtn = pvtrn_arg * rslopeb * denfac;
    if (qr <= Real(0.0)) { vt  = Real(0.0); }
    if (nr <= Real(0.0)) { vtn = Real(0.0); }
}

Referenced by WDM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_slope_snow_cell()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void wdm6_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)/wdm6_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 WDM6::AdvanceColumns(), and wdm6_nislfv_rain_plm6_column().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_venfac()

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

                                                         {
    // Fortran: exp(log((viscos(b,c)/diffus(b,a)))*((.3333333))) / sqrt(viscos) * sqrt(sqrt(den0/c))
    return std::exp(std::log(wdm6_viscos(b,c)/wdm6_diffus(b,a))
                   *wdm6_literal(0.3333333))
          /std::sqrt(wdm6_viscos(b,c))
          *std::sqrt(std::sqrt(den0_arg/c));
}

Referenced by WDM6::AdvanceColumns().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_viscos()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_viscos	(	Real	x,
		Real	y
	)

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

Referenced by wdm6_venfac(), and wdm6_xka().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_xka()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_xka	(	Real	x,
		Real	y
	)

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

Referenced by WDM6::AdvanceColumns(), and wdm6_diffac().

Here is the call graph for this function:

Here is the caller graph for this function:

wdm6_xlcal()

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

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

Referenced by WDM6::AdvanceColumns().

Here is the caller graph for this function:

wdm6_xni_exact()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real wdm6_xni_exact	(	Real	qi,
		Real	den,
		Real	qmin_arg
	)

                                                       {
    Real temp = den * amrex::max(qi, qmin_arg);
    temp = std::sqrt(std::sqrt(temp * temp * temp));
    return amrex::min(amrex::max(Real(5.38e7) * temp, Real(1.e3)), Real(1.e6));
}

Referenced by WDM6::AdvanceColumns().

Here is the caller graph for this function:

Variable Documentation

wdm6_slope_t0c

constexpr amrex::Real wdm6_slope_t0c = wdm6_literal(273.15)

constexpr

Referenced by WDM6::AdvanceColumns().