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:

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 = 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 = 0.0;
    amrex::Real del_T = T - 273.15;  // Convert to Celsius
 
    if (type == 1) {  // Ice (lines ~5631-5644)
        if (T >= 195.8) {
            // Flatau et al. formula for ice
            const amrex::Real a0i = 6.11147274;
            const amrex::Real a1i = 0.503160820;
            const amrex::Real a2i = 0.188439774e-1;
            const amrex::Real a3i = 0.420895665e-3;
            const amrex::Real a4i = 0.615021634e-5;
            const amrex::Real a5i = 0.602588177e-7;
            const amrex::Real a6i = 0.385852041e-9;
            const amrex::Real a7i = 0.146898966e-11;
            const amrex::Real a8i = 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 *= 100.0;  // Convert from hPa to Pa
        } else {
            // Goff-Gratch formula for ice at cold temperatures
            polysvp = std::pow(10.0, (-9.09718*(273.16/T-1.0) - 3.56654*std::log10(273.16/T) +
                                      0.876793*(1.0-T/273.16) + std::log10(6.1071))) * 100.0;
            polysvp = 0.0;
        } // T
    } else {  // Water (lines ~5648-5665)
      if (T >= 202.0) {
        // Flatau et al. formula for liquid water
        const amrex::Real a0 = 6.11239921;
        const amrex::Real a1 = 0.443987641;
        const amrex::Real a2 = 0.142986287e-1;
        const amrex::Real a3 = 0.264847430e-3;
        const amrex::Real a4 = 0.302950461e-5;
        const amrex::Real a5 = 0.206739458e-7;
        const amrex::Real a6 = 0.640689451e-10;
        const amrex::Real a7 = -0.952447341e-13;
        const amrex::Real a8 = -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 *= 100.0;  // Convert from hPa to Pa
      } else {
        // Goff-Gratch formula for water at cold temperatures
        polysvp = std::pow(10.0, (-7.90298*(373.16/T-1.0) + 5.02808*std::log10(373.16/T) -
                                  1.3816e-7*(std::pow(10.0, (11.344*(1.0-T/373.16)))-1.0) +
                                  8.1328e-3*(std::pow(10.0, (-3.49149*(373.16/T-1.0)))-1.0) +
                                  std::log10(1013.246))) * 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, half, one, res, sum, twelve, two, xbig, xden, xinf, xminin;
    amrex::Real xnum, y, y1, ysq, z, zero;
    amrex::Real c[7];
    amrex::Real p[8];
    amrex::Real q[8];
 
    // Mathematical constants
    one = 1.0;
    half = 0.5;
    twelve = 12.0;
    two = 2.0;
    zero = 0.0;
 
    // Machine dependent parameters
    xbig = 35.040;
    xminin = 1.18e-38;
    amrex::Real eps = 1.19e-7;
    xinf = 3.4e38;
 
    // Numerator and denominator coefficients for rational minimax approximation over (1,2)
    p[0] = -1.71618513886549492533811e+0;
    p[1] = 2.47656508055759199108314e+1;
    p[2] = -3.79804256470945635097577e+2;
    p[3] = 6.29331155312818442661052e+2;
    p[4] = 8.66966202790413211295064e+2;
    p[5] = -3.14512729688483675254357e+4;
    p[6] = -3.61444134186911729807069e+4;
    p[7] = 6.64561438202405440627855e+4;
 
    q[0] = -3.08402300119738975254353e+1;
    q[1] = 3.15350626979604161529144e+2;
    q[2] = -1.01515636749021914166146e+3;
    q[3] = -3.10777167157231109440444e+3;
    q[4] = 2.25381184209801510330112e+4;
    q[5] = 4.75584627752788110767815e+3;
    q[6] = -1.34659959864969306392456e+5;
    q[7] = -1.15132259675553483497211e+5;
 
    // Coefficients for minimax approximation over (12, inf)
    c[0] = -1.910444077728e-03;
    c[1] = 8.4171387781295e-04;
    c[2] = -5.952379913043012e-04;
    c[3] = 7.93650793500350248e-04;
    c[4] = -2.777777777777681622553e-03;
    c[5] = 8.333333333333333331554247e-02;
    c[6] = 5.7083835261e-03;
 
    // Initialize variables
    parity = false;
    fact = one;
    n = 0;
    y = x;
 
//    printf("wrf_gamma: Initial y = %g\n", y);
 
    if (y <= zero) {
        // Argument is negative
//        printf("wrf_gamma: Handling negative argument\n");
        y = -x;
        y1 = std::floor(y);
        res = y - y1;
        if (res != zero) {
            if (y1 != std::floor(y1 * half) * two)
                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) {
            // 0.0 < argument < 1.0
//            printf("wrf_gamma: Sub-branch: 0 < y < 1\n");
            z = y;
            y = y + one;
        }
        else {
            // 1.0 < argument < 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 = zero;
        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 0.0 < argument < 1.0
            res = res / y1;
//            printf("wrf_gamma: Adjusted for y < 1: res = %g\n", res);
        }
        else if (y1 > y) {
            // Adjust for 2.0 < argument < 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 - half) * 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 = 0.9189385332046727417803297

constexpr

Referenced by wrf_gamma().