ERF_AdvanceMorrison.cpp File Reference

#include <string>
#include <vector>
#include <memory>
#include <complex>
#include <cmath>
#include <AMReX_Math.H>
#include <AMReX_FArrayBox.H>
#include <AMReX_Geometry.H>
#include <AMReX_TableData.H>
#include <AMReX_MultiFabUtil.H>
#include "ERF.H"
#include "ERF_Constants.H"
#include "ERF_MicrophysicsUtils.H"
#include "ERF_IndexDefines.H"
#include "ERF_DataStruct.H"
#include "ERF_NullMoist.H"
#include "ERF_Morrison.H"

Include dependency graph for ERF_AdvanceMorrison.cpp:

Namespaces
	MORRInd

Enumerations
enum	{   MORRInd::lamc = 0 , MORRInd::lami , MORRInd::lams , MORRInd::lamr ,   MORRInd::lamg , MORRInd::cdist1 , MORRInd::n0i , MORRInd::n0s ,   MORRInd::n0r , MORRInd::n0g , MORRInd::pgam , MORRInd::qc3dten ,   MORRInd::qi3dten , MORRInd::qni3dten , MORRInd::qr3dten , MORRInd::ni3dten ,   MORRInd::ns3dten , MORRInd::nr3dten , MORRInd::qc3d , MORRInd::qi3d ,   MORRInd::qni3d , MORRInd::qr3d , MORRInd::ni3d , MORRInd::ns3d ,   MORRInd::nr3d , MORRInd::t3dten , MORRInd::qv3dten , MORRInd::t3d ,   MORRInd::qv3d , MORRInd::pres , MORRInd::dzq , MORRInd::w3d ,   MORRInd::nc3d , MORRInd::nc3dten , MORRInd::qg3dten , MORRInd::ng3dten ,   MORRInd::qg3d , MORRInd::ng3d , MORRInd::qgsten , MORRInd::qrsten ,   MORRInd::qisten , MORRInd::qnisten , MORRInd::qcsten , MORRInd::qrcu1d ,   MORRInd::qscu1d , MORRInd::qicu1d , MORRInd::precrt , MORRInd::snowrt ,   MORRInd::snowprt , MORRInd::grplprt , MORRInd::effc , MORRInd::effi ,   MORRInd::effs , MORRInd::effr , MORRInd::effg , MORRInd::rho ,   MORRInd::mu , MORRInd::ain , MORRInd::arn , MORRInd::asn ,   MORRInd::acn , MORRInd::agn , MORRInd::dumi , MORRInd::dumr ,   MORRInd::dumfni , MORRInd::dumg , MORRInd::dumfng , MORRInd::uni ,   MORRInd::umi , MORRInd::umr , MORRInd::fr , MORRInd::fi ,   MORRInd::fni , MORRInd::fg , MORRInd::fng , MORRInd::rgvm ,   MORRInd::faloutr , MORRInd::falouti , MORRInd::faloutni , MORRInd::faltndr ,   MORRInd::faltndi , MORRInd::faltndni , MORRInd::dumqs , MORRInd::dumfns ,   MORRInd::ums , MORRInd::uns , MORRInd::fs , MORRInd::fns ,   MORRInd::falouts , MORRInd::faloutns , MORRInd::faloutg , MORRInd::faloutng ,   MORRInd::faltnds , MORRInd::faltndns , MORRInd::unr , MORRInd::faltndg ,   MORRInd::faltndng , MORRInd::dumc , MORRInd::dumfnc , MORRInd::unc ,   MORRInd::umc , MORRInd::ung , MORRInd::umg , MORRInd::fc ,   MORRInd::faloutc , MORRInd::faloutnc , MORRInd::faltndc , MORRInd::faltndnc ,   MORRInd::fnc , MORRInd::dumfnr , MORRInd::faloutnr , MORRInd::faltndnr ,   MORRInd::fnr , MORRInd::dlams , MORRInd::dlamr , MORRInd::dlami ,   MORRInd::dlamc , MORRInd::dlamg , MORRInd::xxls , MORRInd::xxlv ,   MORRInd::cpm , MORRInd::xlf , MORRInd::NumInds }

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	wrf_gamma (amrex::Real x)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real	gamma_function (Real x)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	calc_saturation_vapor_pressure (const amrex::Real T, const int type)

Variables
constexpr Real	xxx = Real(0.9189385332046727417803297)

Function Documentation

calc_saturation_vapor_pressure()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real calc_saturation_vapor_pressure	(	const amrex::Real	T,
		const int	type
	)

Helper function to calculate saturation vapor pressure for water or ice. This corresponds to the POLYSVP function in the Fortran code (line ~5580).

Parameters

[in]	T	Temperature in Kelvin
[in]	type	0 for liquid water, 1 for ice

Returns

Saturation vapor pressure in Pascals

  {
    amrex::Real polysvp = Real(0);
    amrex::Real del_T = T - Real(273.15);  // Convert to Celsius
 
    if (type == 1) {  // Ice (lines ~5631-5644)
        if (T >= Real(195.8)) {
            // Flatau et al. formula for ice
            const amrex::Real a0i = Real(6.11147274);
            const amrex::Real a1i = Real(0.503160820);
            const amrex::Real a2i = Real(0.188439774e-1);
            const amrex::Real a3i = Real(0.420895665e-3);
            const amrex::Real a4i = Real(0.615021634e-5);
            const amrex::Real a5i = Real(0.602588177e-7);
            const amrex::Real a6i = Real(0.385852041e-9);
            const amrex::Real a7i = Real(0.146898966e-11);
            const amrex::Real a8i = Real(0.252751365e-14);
 
            polysvp = a0i + del_T*(a1i + del_T*(a2i + del_T*(a3i + del_T*(a4i + del_T*(a5i + del_T*(a6i + del_T*(a7i + a8i*del_T)))))));
            polysvp *= Real(100.0);  // Convert from hPa to Pa
        } else {
            // Goff-Gratch formula for ice at cold temperatures
            polysvp = std::pow(Real(10.0), (-Real(9.09718)*(Real(273.16)/T-one) - Real(3.56654)*std::log10(Real(273.16)/T) +
                                             Real(0.876793)*(one-T/Real(273.16)) + std::log10(Real(6.1071)))) * Real(100.0);
        } // T
    } else {  // Water (lines ~5648-5665)
      if (T >= Real(202.0)) {
        // Flatau et al. formula for liquid water
        const amrex::Real a0 = Real(6.11239921);
        const amrex::Real a1 = Real(0.443987641);
        const amrex::Real a2 = Real(0.142986287e-1);
        const amrex::Real a3 = Real(0.264847430e-3);
        const amrex::Real a4 = Real(0.302950461e-5);
        const amrex::Real a5 = Real(0.206739458e-7);
        const amrex::Real a6 = Real(0.640689451e-10);
        const amrex::Real a7 = -Real(0.952447341e-13);
        const amrex::Real a8 = -Real(0.976195544e-15);
 
        polysvp = a0 + del_T*(a1 + del_T*(a2 + del_T*(a3 + del_T*(a4 + del_T*(a5 + del_T*(a6 + del_T*(a7 + a8*del_T)))))));
        polysvp *= Real(100.0);  // Convert from hPa to Pa
      } else {
        // Goff-Gratch formula for water at cold temperatures
        polysvp = std::pow(Real(10.0), (-Real(7.90298)*(Real(373.16)/T-one) + Real(5.02808)*std::log10(Real(373.16)/T) -
                                  Real(1.3816e-7)*(std::pow(Real(10.0), (Real(11.344)*(one-T/Real(373.16))))-one) +
                                  Real(8.1328e-3)*(std::pow(Real(10.0), (-Real(3.49149)*(Real(373.16)/T-one)))-one) +
                                  std::log10(Real(1013.246)))) * Real(100.0);
      }
    }
 
    return polysvp;
  }

Referenced by Morrison::Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

gamma_function()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Real gamma_function

(

Real

x

)

                            {
  return wrf_gamma(x);
}

Referenced by Morrison::Advance().

Here is the call graph for this function:

Here is the caller graph for this function:

wrf_gamma()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real wrf_gamma

(

amrex::Real

x

)

{
    // Debug: using printf since it's GPU compatible
//    printf("wrf_gamma: Input value x = %g\n", x);
 
    // Local variables
    int i, n;
    bool parity = false;
    amrex::Real fact, res, sum, twelve, xbig, xden, xinf, xminin;
    amrex::Real xnum, y, y1, ysq, z;
    amrex::Real c[7];
    amrex::Real p[8];
    amrex::Real q[8];
 
    // Mathematical constants
    twelve = Real(12.0);
 
    // Machine dependent parameters
    xbig = Real(35.040);
    xminin = Real(1.18e-38);
    amrex::Real eps = Real(1.19e-7);
    xinf = Real(3.4e38);
 
    // Numerator and denominator coefficients for rational minimax approximation over (1,2)
    p[0] = -Real(1.71618513886549492533811e+0);
    p[1] =  Real(2.47656508055759199108314e+1);
    p[2] = -Real(3.79804256470945635097577e+2);
    p[3] =  Real(6.29331155312818442661052e+2);
    p[4] =  Real(8.66966202790413211295064e+2);
    p[5] = -Real(3.14512729688483675254357e+4);
    p[6] = -Real(3.61444134186911729807069e+4);
    p[7] =  Real(6.64561438202405440627855e+4);
 
    q[0] = -Real(3.08402300119738975254353e+1);
    q[1] =  Real(3.15350626979604161529144e+2);
    q[2] = -Real(1.01515636749021914166146e+3);
    q[3] = -Real(3.10777167157231109440444e+3);
    q[4] =  Real(2.25381184209801510330112e+4);
    q[5] =  Real(4.75584627752788110767815e+3);
    q[6] = -Real(1.34659959864969306392456e+5);
    q[7] = -Real(1.15132259675553483497211e+5);
 
    // Coefficients for minimax approximation over (12, inf)
    c[0] = -Real(1.910444077728e-03);
    c[1] = Real(8.4171387781295e-04);
    c[2] = -Real(5.952379913043012e-04);
    c[3] = Real(7.93650793500350248e-04);
    c[4] = -Real(2.777777777777681622553e-03);
    c[5] = Real(8.333333333333333331554247e-02);
    c[6] = Real(5.7083835261e-03);
 
    // Initialize variables
    parity = false;
    fact = one;
    n = 0;
    y = x;
 
//    printf("wrf_gamma: Initial y = %g\n", y);
 
    if (y <= Real(0)) {
        // Argument is negative
//        printf("wrf_gamma: Handling negative argument\n");
        y = -x;
        y1 = std::floor(y);
        res = y - y1;
        if (res != Real(0)) {
            if (y1 != std::floor(y1 * myhalf) * Real(2))
                parity = true;
            Real pi=amrex::Math::pi<Real>();
            fact = -pi / std::sin(pi * res);
            y = y + one;
//            printf("wrf_gamma: After reflection formula: y = %g, fact = %g, parity = %d\n",
//                   y, fact, parity);
        }
        else {
//            printf("wrf_gamma: Singularity detected, returning xinf = %g\n", xinf);
            res = xinf;
            return res;
        }
    }
 
    // Argument is positive
    if (y < eps) {
        // Argument < eps
//        printf("wrf_gamma: Small argument branch (y < eps)\n");
        if (y >= xminin) {
            res = one / y;
//            printf("wrf_gamma: Small argument result: res = %g\n", res);
        }
        else {
//            printf("wrf_gamma: Argument too small, returning xinf = %g\n", xinf);
            res = xinf;
            return res;
        }
    }
    else if (y < twelve) {
        // Medium range argument
//        printf("wrf_gamma: Medium range branch (eps <= y < 12)\n");
        y1 = y;
        if (y < one) {
            // Real(0) < argument < one
//            printf("wrf_gamma: Sub-branch: 0 < y < 1\n");
            z = y;
            y = y + one;
        }
        else {
            // one < argument < Real(12.0), reduce argument if necessary
            n = static_cast<int>(y) - 1;
//            printf("wrf_gamma: Sub-branch: 1 <= y < 12, n = %d\n", n);
            y = y - static_cast<amrex::Real>(n);
            z = y - one;
        }
 
        // Evaluate approximation
//        printf("wrf_gamma: Before approximation: z = %g, y = %g\n", z, y);
        xnum = Real(0);
        xden = one;
        for (i = 0; i < 8; i++) {
            xnum = (xnum + p[i]) * z;
            xden = xden * z + q[i];
        }
        res = xnum / xden + one;
//        printf("wrf_gamma: After approximation: res = %g\n", res);
 
        if (y1 < y) {
            // Adjust result for case Real(0) < argument < one
            res = res / y1;
//            printf("wrf_gamma: Adjusted for y < 1: res = %g\n", res);
        }
        else if (y1 > y) {
            // Adjust for two < argument < Real(12.0)
//            printf("wrf_gamma: Adjusting for y > 2 with %d multiplications\n", n);
            for (i = 0; i < n; i++) {
                res = res * y;
                y = y + one;
//                printf("wrf_gamma: Multiplication %d: res = %g, y = %g\n", i+1, res, y);
            }
        }
    }
    else {
        // Large argument
//        printf("wrf_gamma: Large argument branch (y >= 12)\n");
        if (y <= xbig) {
            ysq = y * y;
            sum = c[6];
            for (i = 0; i < 6; i++) {
                sum = sum / ysq + c[i];
//                printf("wrf_gamma: Sum step %d: sum = %g\n", i+1, sum);
            }
            sum = sum / y - y + xxx;
            sum = sum + (y - myhalf) * std::log(y);
//            printf("wrf_gamma: Before exp: sum = %g\n", sum);
            res = std::exp(sum);
//            printf("wrf_gamma: After exp: res = %g\n", res);
        }
        else {
//            printf("wrf_gamma: Argument too large, returning xinf = %g\n", xinf);
            res = xinf;
            return res;
        }
    }
 
    // Final adjustments
    if (parity) {
        res = -res;
//        printf("wrf_gamma: Applied parity adjustment: res = %g\n", res);
    }
    if (fact != one) {
        res = fact / res;
//        printf("wrf_gamma: Applied reflection adjustment: res = %g\n", res);
    }
 
//    printf("wrf_gamma: Final result = %g\n", res);
    return res;
}

Referenced by gamma_function().

Here is the caller graph for this function:

Variable Documentation

xxx

constexpr Real xxx = Real(0.9189385332046727417803297)

constexpr

Referenced by wrf_gamma().