diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/.gitignore b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/.gitignore index b39e14de8..5dc9fae88 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/.gitignore +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/.gitignore @@ -7,4 +7,4 @@ test_data/ .gt_cache_* .translate-errors/ .vscode -test_data/ \ No newline at end of file +test_data/ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/aer_activation.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/aer_activation.py new file mode 100644 index 000000000..f61b6b226 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/aer_activation.py @@ -0,0 +1,443 @@ +from gt4py.cartesian.gtscript import PARALLEL, computation, exp, interval, log, sqrt + +import pyMoist.aer_activation_constants as constants +from ndsl import QuantityFactory, StencilFactory, orchestrate +from ndsl.dsl.typing import Float, FloatField, FloatFieldIJ, Int +from pyMoist.numerical_recipes import Erf +from pyMoist.types import FloatField_NModes + + +def aer_activation_stencil( + aero_dgn: FloatField_NModes, + aero_num: FloatField_NModes, + nacti: FloatField, + t: FloatField, + plo: FloatField, + qicn: FloatField, + qils: FloatField, + qlcn: FloatField, + qlls: FloatField, + nn_land: Float, + frland: FloatFieldIJ, + nn_ocean: Float, + aero_hygroscopicity: FloatField_NModes, + nwfa: FloatField, + nactl: FloatField, + vvel: FloatField, + tke: FloatField, + aero_sigma: FloatField_NModes, + nact: FloatField_NModes, + ni: FloatField_NModes, + rg: FloatField_NModes, + sig0: FloatField_NModes, + bibar: FloatField_NModes, +): + """ + Compute the aerosol activation stencil. + + Parameters: + aero_dgn (4D in): AeroProps aerosol geometric mean diameter. + aero_num (4D in): AeroProps aerosol number concentration. + nacti (3D out): Activated aerosol number concentration. + t (3D in): Temperature field. + plo (3D in): Pressure field. + qicn (3D in): Ice cloud number concentration. + qils (3D in): Ice liquid water content. + qlcn (3D in): Cloud number concentration. + qlls (3D in): Liquid water content. + nn_land (1D in): Land-based aerosol activation number. + frland (2D in): Fraction of land. + nn_ocean (1D in): Ocean-based aerosol activation number. + aero_hygroscopicity (4D in): Aerosol hygroscopicity parameter. + nwfa (3D out): Number of activated aerosols. + nactl (3D out): Number of activated aerosols in liquid clouds. + vvel (3D in): Vertical velocity field. + tke (3D in): Turbulent kinetic energy field. + aero_sigma (4D in): AeroProps aerosol geometric standard deviation. + nact (4D temporary in): Activated aerosol number concentration field. + ni (4D temporary in): AeroProp ice crystal number concentration field. + rg (4D temporary in): AeroProp geometric mean radius of aerosols. + sig0 (4D temporary in): AeroProp aerosol geometric standard deviation field. + bibar (4D temporary in): AeroProp Hygroscopicity parameter field. + + Returns: + None + """ + with computation(PARALLEL), interval(...): + + # Compute nwfa + # Fortran AeroProps aero_kap is aero_hygroscopicity + nfaux = 0.0 + n = 0 + while n < constants.n_modes: + if aero_hygroscopicity[0, 0, 0][n] > 0.4: + nfaux += aero_num[0, 0, 0][n] + n += 1 + nwfa = nfaux + + # Determine aerosol number concentration at cloud base + nactl = nn_land * frland + nn_ocean * (1.0 - frland) + nacti = nn_land * frland + nn_ocean * (1.0 - frland) + + tk = t # [K] + press = plo # [Pa] + air_den = press * 28.8e-3 / 8.31 / tk # [kg/m3] + qi = (qicn + qils) * 1.0e3 # [g/kg] + ql = (qlcn + qlls) * 1.0e3 # [g/kg] + wupdraft = vvel + sqrt(tke) + + # Liquid Clouds Calculation + if ( + (tk >= (constants.MAPL_TICE - 40.0)) + and (plo > 10000.0) + and (0.1 < wupdraft) + and (wupdraft < 100.0) + ): + n = 0 + while n < constants.n_modes: + ni[0, 0, 0][n] = max( + aero_num[0, 0, 0][n] * air_den, constants.ZERO_PAR + ) # unit: [m-3] + rg[0, 0, 0][n] = max( + aero_dgn[0, 0, 0][n] * 0.5 * 1.0e6, constants.ZERO_PAR + ) # unit: [um] + sig0[0, 0, 0][n] = aero_sigma[0, 0, 0][n] # unit: [um] + bibar[0, 0, 0][n] = max( + aero_hygroscopicity[0, 0, 0][n], constants.ZERO_PAR + ) + n += 1 + + """ + 12-12-06, DLW: Routine to calculate the activated fraction of the number + and mass concentrations, as well as the number and mass + concentrations activated for each of nmodes modes. The + minimum dry radius for activation for each mode is also returned. + + Each mode is assumed to potentially contains 5 chemical species: + (1) sulfate + (2) BC + (3) OC + (4) mineral dust + (5) sea salt + + The aerosol activation parameterizations are described in + + 1. Abdul-Razzak et al. 1998, JGR, vol.103, p.6123-6131. + 2. Abdul-Razzak and Ghan 2000, JGR, vol.105, p.6837-6844. + + and values for many of the required parameters were taken from + + 3. Ghan et al. 2001, JGR vol 106, p.5295-5316. + + With the density of sea salt set to the value used in ref. 3 (1900 kg/m^3), + this routine yields values for the hygroscopicity parameters Bi in + agreement with ref. 3. + + This routine is for the multiple-aerosol type parameterization. + """ + # rdrp is the radius value used in eqs.(17) & (18) and was adjusted to + # yield eta and zeta values close to those given in + # a-z et al. 1998 figure 5. + + # tuned to approximate the results in figures 1-5 in a-z et al. 1998. + rdrp = 0.105e-06 # [m] + + # These variables are common to all modes and need only be computed once. + dv = ( + constants.DIJH2O0 + * (constants.P0DIJ / plo) + * (tk / constants.T0DIJ) ** 1.94e00 + ) # [m^2/s] (p&k,2nd ed., p.503) + surten = 76.10e-3 - 0.155e-3 * (tk - 273.15e00) # [j/m^2] + wpe = exp( + 77.34491296 - 7235.424651 / tk - 8.2 * log(tk) + tk * 5.7113e-3 + ) # [pa] + dumw = sqrt( + constants.TWOPI * constants.WMOLMASS / constants.RGASJMOL / tk + ) # [s/m] + dvprime = dv / ( + (rdrp / (rdrp + constants.DELTAV)) + + (dv * dumw / (rdrp * constants.ALPHAC)) + ) # [m^2/s] - eq. (17) + xka = ( + 5.69 + 0.017 * (tk - 273.15) + ) * 418.4e-5 # [j/m/s/k] (0.0238 j/m/s/k at 273.15 k) + duma = sqrt( + constants.TWOPI * constants.AMOLMASS / constants.RGASJMOL / tk + ) # [s/m] + xkaprime = xka / ( + (rdrp / (rdrp + constants.DELTAT)) + + ( + xka + * duma + / (rdrp * constants.ALPHAT * constants.DENH2O * constants.CPAIR) + ) + ) # [j/m/s/k] + g = 1.0 / ( + (constants.DENH2O * constants.RGASJMOL * tk) + / (wpe * dvprime * constants.WMOLMASS) + + ((constants.HEATVAP * constants.DENH2O) / (xkaprime * tk)) + * ( + (constants.HEATVAP * constants.WMOLMASS) / (constants.RGASJMOL * tk) + - 1.0 + ) + ) # [m^2/s] + a = (2.0 * surten * constants.WMOLMASS) / ( + constants.DENH2O * constants.RGASJMOL * tk + ) # [m] + alpha = (constants.GRAVITY / (constants.RGASJMOL * tk)) * ( + (constants.WMOLMASS * constants.HEATVAP) / (constants.CPAIR * tk) + - constants.AMOLMASS + ) # [1/m] + gamma = (constants.RGASJMOL * tk) / (wpe * constants.WMOLMASS) + ( + constants.WMOLMASS * constants.HEATVAP * constants.HEATVAP + ) / ( + constants.CPAIR * plo * constants.AMOLMASS * tk + ) # [m^3/kg] + dum = sqrt(alpha * wupdraft / g) # [1/m] + zeta = 2.0 * a * dum / 3.0 # [1] + + # These variables must be computed for each mode + n = 0 + while n < constants.n_modes: + xlogsigm = log(sig0[0, 0, 0][n]) # [1] + smax = 0.0 # [1] + sm = (2.0 / sqrt(bibar[0, 0, 0][n])) * ( + a / (3.0 * rg[0, 0, 0][n]) + ) ** 1.5 # [1] + eta = dum ** 3 / ( + constants.TWOPI * constants.DENH2O * gamma * ni[0, 0, 0][n] + ) # [1] + f1 = 0.5 * exp(2.50 * xlogsigm ** 2) # [1] + f2 = 1.0 + 0.25 * xlogsigm # [1] + smax = ( + smax + + ( + f1 * (zeta / eta) ** 1.5 + + f2 * (sm ** 2 / (eta + 3.0 * zeta)) ** 0.75 + ) + / sm ** 2 + ) # [1] - eq. (6) + n += 1 + + smax = 1.0e00 / sqrt(smax) # [1] + n = 0 + u = 0.0 + while n < constants.n_modes: + sm = (2.0 / sqrt(bibar[0, 0, 0][n])) * ( + a / (3.0 * rg[0, 0, 0][n]) + ) ** 1.5 # [1] + xlogsigm = log(sig0[0, 0, 0][n]) # [1] + ac = rg[0, 0, 0][n] * (sm / smax) ** 0.66666666666666667 # [um] + u = log(ac / rg[0, 0, 0][n]) / (constants.SQRT2 * xlogsigm) # [1] + fracactn = 0.5 * (1.0 - Erf(u)) # [1] + nact[0, 0, 0][n] = fracactn * ni[0, 0, 0][n] # [#/m^3] + n += 1 + + numbinit = 0.0 + nactl = 0.0 + n = 0 + while n < constants.n_modes: + numbinit += aero_num[0, 0, 0][n] * air_den + nactl += nact[0, 0, 0][n] + n += 1 + nactl = min(nactl, 0.99 * numbinit) + + # Ice Clouds Calculation + if (tk <= constants.MAPL_TICE) and ( + qi > constants.FLOAT_TINY or ql > constants.FLOAT_TINY + ): + numbinit = 0.0 + n = 0 + while n < constants.n_modes: + # diameters > 0.5 microns + if aero_dgn[0, 0, 0][n] >= 0.5e-6: + numbinit += aero_num[0, 0, 0][n] + n += 1 + numbinit *= air_den # [#/m3] + # Number of activated IN following deMott (2010) [#/m3] + nacti = ( + constants.AI + * ((constants.MAPL_TICE - tk) ** constants.BI) + * ( + numbinit + ** (constants.CI * (constants.MAPL_TICE - tk) + constants.DI) + ) + ) + + # apply limits for NACTL/NACTI + if nactl < constants.NN_MIN: + nactl = constants.NN_MIN + if nactl > constants.NN_MAX: + nactl = constants.NN_MAX + if nacti < constants.NN_MIN: + nacti = constants.NN_MIN + if nacti > constants.NN_MAX: + nacti = constants.NN_MAX + + +class AerActivation: + """ + Class for aerosol activation computation. + + Attributes: + stencil_factory (StencilFactory): Factory for creating stencil computations. + quantity_factory (QuantityFactory): Factory for creating quantities. + n_modes (Int): Number of aerosol modes. + USE_AERSOL_NN (bool): Flag indicating whether to use neural network for aerosol. + """ + + def __init__( + self, + stencil_factory: StencilFactory, + quantity_factory: QuantityFactory, + n_modes: Int, + USE_AERSOL_NN: bool, + ) -> None: + """ + Initialize the AerActivation class. + + Parameters: + stencil_factory (StencilFactory): Factory for creating stencil computations. + quantity_factory (QuantityFactory): Factory for creating quantities. + n_modes (Int): Number of aerosol modes. + USE_AERSOL_NN (bool): Flag indicating whether to use neural network for aerosol. + + Raises: + NotImplementedError: If the number of modes is not equal to the expected number. + NotImplementedError: If the neural network for aerosol is not implemented. + """ + self._nact = quantity_factory._numpy.zeros( + ( + stencil_factory.grid_indexing.domain[0], + stencil_factory.grid_indexing.domain[1], + stencil_factory.grid_indexing.domain[2], + constants.n_modes, + ), + dtype=Float, + ) + self._ni = quantity_factory._numpy.zeros( + ( + stencil_factory.grid_indexing.domain[0], + stencil_factory.grid_indexing.domain[1], + stencil_factory.grid_indexing.domain[2], + constants.n_modes, + ), + dtype=Float, + ) + self._rg = quantity_factory._numpy.zeros( + ( + stencil_factory.grid_indexing.domain[0], + stencil_factory.grid_indexing.domain[1], + stencil_factory.grid_indexing.domain[2], + constants.n_modes, + ), + dtype=Float, + ) + self._sig0 = quantity_factory._numpy.zeros( + ( + stencil_factory.grid_indexing.domain[0], + stencil_factory.grid_indexing.domain[1], + stencil_factory.grid_indexing.domain[2], + constants.n_modes, + ), + dtype=Float, + ) + self._bibar = quantity_factory._numpy.zeros( + ( + stencil_factory.grid_indexing.domain[0], + stencil_factory.grid_indexing.domain[1], + stencil_factory.grid_indexing.domain[2], + constants.n_modes, + ), + dtype=Float, + ) + + if constants.n_modes != n_modes: + raise NotImplementedError( + f"Coding limitation: 14 modes are expected, getting {n_modes}" + ) + + if not USE_AERSOL_NN: + raise NotImplementedError("Non NN Aerosol not implemented") + + orchestrate(obj=self, config=stencil_factory.config.dace_config) + + self.aer_activation = stencil_factory.from_origin_domain( + func=aer_activation_stencil, + origin=(0, 0, 0), + domain=stencil_factory.grid_indexing.domain, + ) + + def __call__( + self, + aero_dgn: FloatField_NModes, + aero_num: FloatField_NModes, + nacti: FloatField, + t: FloatField, + plo: FloatField, + qicn: FloatField, + qils: FloatField, + qlcn: FloatField, + qlls: FloatField, + nn_land: Float, + frland: FloatFieldIJ, + nn_ocean: Float, + aero_hygroscopicity: FloatField_NModes, + nwfa: FloatField, + nactl: FloatField, + vvel: FloatField, + tke: FloatField, + aero_sigma: FloatField_NModes, + ) -> None: + """ + Compute aerosol activation by calling the stencil function. + + Parameters: + aero_dgn (4D in): AeroProps aerosol geometric mean diameter. + aero_num (4D in): AeroProps aerosol number concentration. + nacti (3D out): Activated aerosol number concentration. + t (3D in): Temperature field. + plo (3D in): Pressure field. + qicn (3D in): Ice cloud number concentration. + qils (3D in): Ice liquid water content. + qlcn (3D in): Cloud number concentration. + qlls (3D in): Liquid water content. + nn_land (1D in): Land-based aerosol activation number. + frland (2D in): Fraction of land. + nn_ocean (1D in): Ocean-based aerosol activation number. + aero_hygroscopicity (4D in): Aerosol hygroscopicity parameter. + nwfa (3D out): Number of activated aerosols. + nactl (3D out): Number of activated aerosols in liquid clouds. + vvel (3D in): Vertical velocity field. + tke (3D in): Turbulent kinetic energy field. + aero_sigma (4D in): AeroProps aerosol geometric standard deviation. + + Returns: + None + """ + self.aer_activation( + aero_dgn, + aero_num, + nacti, + t, + plo, + qicn, + qils, + qlcn, + qlls, + nn_land, + frland, + nn_ocean, + aero_hygroscopicity, + nwfa, + nactl, + vvel, + tke, + aero_sigma, + self._nact, + self._ni, + self._rg, + self._sig0, + self._bibar, + ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/aer_activation_constants.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/aer_activation_constants.py new file mode 100644 index 000000000..0b386982b --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/aer_activation_constants.py @@ -0,0 +1,62 @@ +import numpy as np + +from ndsl.dsl.typing import Float + + +""" +All global constants used for aer_actv_single_moment.F90 +""" +# MAPL Constants +MAPL_TICE = 273.16 # K + +# Aerosol Activation Constants +# 32 bit +R_AIR = 3.47e-3 # m3 Pa kg-1K-1 + +# 64 bit +ZERO_PAR = 1.0e-6 # small non-zero value +AI = 0.0000594 +BI = 3.33 +CI = 0.0264 +DI = 0.0033 + +BETAAI = -2.262e3 +GAMAI = 5.113e6 +DELTAAI = 2.809e3 +DENSIC = 917.0 # Ice crystal density in kgm-3 + +# 32 bit +NN_MIN = 100.0e6 +NN_MAX = 1000.0e6 + +# ACTFRAC_Mat constants - all 64 bit +PI = 3.141592653589793e00 +TWOPI = 2.0e00 * PI +SQRT2 = 1.414213562e00 +THREESQRT2BY2 = 1.5e00 * SQRT2 + +AVGNUM = 6.0221367e23 # [1/mol] +RGASJMOL = 8.31451e00 # [j/mol/k] +WMOLMASS = 18.01528e-03 # molar mass of h2o [kg/mol] +AMOLMASS = 28.966e-03 # molar mass of air [kg/mol] +ASMOLMSS = 132.1406e-03 # molar mass of nh42so4 [kg/mol] +DENH2O = 1.00e03 # density of water [kg/m^3] +DENAMSUL = 1.77e03 # density of pure ammonium sulfate [kg/m^3] +XNUAMSUL = 3.00e00 # number of ions formed when the salt is dissolved in water [1] +PHIAMSUL = 1.000e00 # osmotic coefficient value in a-r 1998. [1] +GRAVITY = 9.81e00 # grav. accel. at the earth's surface [m/s/s] +HEATVAP = 40.66e03 / WMOLMASS # latent heat of vap. for water and tnbp [j/kg] +CPAIR = 1006.0e00 # heat capacity of air [j/kg/k] +T0DIJ = 273.15e00 # reference temp. for dv [k] +P0DIJ = 101325.0e00 # reference pressure for dv [pa] +DIJH2O0 = 0.211e-04 # reference value of dv [m^2/s] (p&k,2nd ed., p.503) +DELTAV = 1.096e-07 # vapor jump length [m] +DELTAT = 2.160e-07 # thermal jump length [m] +ALPHAC = 1.000e00 # condensation mass accommodation coefficient [1] +ALPHAT = 0.960e00 # thermal accommodation coefficient [1] + +# Define how many modes in an aerosol +n_modes = 14 + +# Python euqivalent of Fortran's tiny(X) +FLOAT_TINY = np.finfo(Float).tiny diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/numerical_recipes.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/numerical_recipes.py new file mode 100644 index 000000000..1b14136e7 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/numerical_recipes.py @@ -0,0 +1,168 @@ +import gt4py.cartesian.gtscript as gtscript +from gt4py.cartesian.gtscript import exp, log + +from ndsl.dsl.typing import Float + + +@gtscript.function +def GammLn(xx: Float) -> Float: + """ + See numerical recipes, w. press et al., 2nd edition. + + Compute the natural logarithm of the gamma function for a given value xx. + + Parameters: + xx (Float in): Input value for which the natural logarithm of the + gamma function is to be computed. + + Returns: + Float: The natural logarithm of the gamma function value for the input xx. + """ + stp = 2.5066282746310005 + + x = xx + y = x + tmp = x + 5.5 + tmp = (x + 0.5) * log(tmp) - tmp + ser = 1.000000000190015 + + ser += 76.18009172947146 / (y + 1) + ser += -86.50532032941677 / (y + 1) + ser += 24.01409824083091 / (y + 1) + ser += -1.231739572450155 / (y + 1) + ser += 0.001208650973866179 / (y + 1) + ser += -0.000005395239384953 / (y + 1) + + gammln = tmp + log(stp * ser / x) + return gammln + + +@gtscript.function +def gser(a: Float, x: Float, gln: Float) -> Float: + """ + See numerical recipes, w. press et al., 2nd edition. + + Compute the series representation of the incomplete gamma function. + + Parameters: + a (Float in): Parameter a for the incomplete gamma function. + x (Float in): Parameter x for the incomplete gamma function. + gln (Float in): Natural logarithm of the gamma function. + + Returns: + Float: The series representation of the incomplete gamma function. + """ + eps = 3.0e-9 # was eps=3.0d-07 in press et al. + itmax = 10000 # was itmax=100 in press et al. + gln = GammLn(a) + if x <= 0: + # Fortran messages here x < 0 in gser + # TODO: Allow print in GT4Py + # if x < 0: + # raise ValueError('aero_actv: subroutine gser: x < 0 in gser') + gamser = 0.0 + else: + ap = a + sum_ = 1.0 / a + del_ = sum_ + n = 0 + while n < itmax: + ap += 1.0 + del_ *= x / ap + sum_ += del_ + if abs(del_) < abs(sum_) * eps: + gamser = sum_ * exp(-x + a * log(x) - gln) + n = itmax + n += 1 + gamser = sum_ * exp(-x + a * log(x) - gln) + return gamser + + +@gtscript.function +def gcf_matrix(a: Float, x: Float, gln: Float) -> Float: + """ + See numerical recipes, w. press et al., 2nd edition. + + Compute the continued fraction representation of the incomplete gamma function. + + Parameters: + a (Float in): Parameter a for the incomplete gamma function. + x (Float in): Parameter x for the incomplete gamma function. + gln (Float in): Natural logarithm of the gamma function. + + Returns: + Float: The continued fraction representation of the incomplete gamma function. + """ + itmax = 10000 + eps = 3.0e-7 + fpmin = 1.0e-30 + gln = GammLn(a) + b = x + 1.0 - a + c = 1.0 / fpmin + d = 1.0 / b + h = d + + i = 1 + while i <= itmax: + an = -i * (i - a) + b += 2.0 + d = an * d + b + if abs(d) < fpmin: + d = fpmin + c = b + an / c + if abs(c) < fpmin: + c = fpmin + d = 1.0 / d + del_ = d * c + h *= del_ + if abs(del_ - 1.0) < eps: + i = itmax + 1 + i += 1 + return exp(-x + a * log(x) - gln) * h + + +@gtscript.function +def GammP(a: Float, x: Float) -> Float: + """ + See numerical recipes, w. press et al., 2nd edition. + + Compute the incomplete gamma function for given values a and x. + + Parameters: + a (Float in): Parameter a for the incomplete gamma function. + x (Float in): Parameter x for the incomplete gamma function. + + Returns: + Float: The incomplete gamma function value for the input parameters a and x. + """ + # Fortran messages here potential bad arguments + # TODO: Allow print in GT4Py + # if (x < 0.0) or (a <= 0.0): + # raise ValueError("aero_actv: function gammp: bad arguments") + gln = GammLn(a) + if x < a + 1.0: + gammp = gser(a, x, gln) + else: + gammp = 1.0 - gcf_matrix(a, x, gln) + return gammp + + +@gtscript.function +def Erf(x: Float) -> Float: + """ + See numerical recipes, w. press et al., 2nd edition. + + Compute the error function for a given value x. + + Parameters: + x (Float in): Input value for which the error function is to be computed. + + Returns: + Float: The error function value for the input x. + """ + erf = 0.0 + if x < 0.0e00: + erf = -1.0 * GammP(0.5, x ** 2) + else: + erf = GammP(0.5, x ** 2) + return erf diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling.py index a31640328..ce0f0449c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling.py @@ -1,16 +1,15 @@ -from ndsl.constants import X_DIM, Y_DIM, Z_DIM -from ndsl.dsl.typing import FloatField, Int, Float -from ndsl import QuantityFactory, StencilFactory import gt4py.cartesian.gtscript as gtscript -from gt4py.cartesian.gtscript import computation, interval, PARALLEL, log10 +from gt4py.cartesian.gtscript import PARALLEL, computation, interval, log10 + import pyMoist.radiation_coupling_constants as radconstants +from ndsl import QuantityFactory, StencilFactory +from ndsl.constants import X_DIM, Y_DIM, Z_DIM +from ndsl.dsl.typing import Float, FloatField, Int + @gtscript.function -def air_density( - PL: Float, - TE: Float - )-> Float: - ''' +def air_density(PL: Float, TE: Float) -> Float: + """ Calculate air density [kg/m^3] Parameters: @@ -19,19 +18,20 @@ def air_density( Returns: Float: Calculated air density. - ''' + """ air_density = (100.0 * PL) / (radconstants.MAPL_RGAS * TE) return air_density + @gtscript.function def cloud_effective_radius_ice( PL: Float, - TE: Float, - QC: Float, + TE: Float, + QC: Float, NNL: Float, NNI: Float, -)-> Float: - ''' +) -> Float: + """ Calculate the effective radius of ice clouds [m] Implementation of LDRADIUS4 for Ice clouds @@ -44,38 +44,43 @@ def cloud_effective_radius_ice( Returns: Float: Effective radius of ice clouds. - ''' - #Calculate ice water content - WC = 1.0e3 * air_density(PL, TE) * QC #air density [g/m3] * ice cloud mixing ratio [kg/kg] - #Calculate radius in meters [m] + """ + # Calculate ice water content + WC = ( + 1.0e3 * air_density(PL, TE) * QC + ) # air density [g/m3] * ice cloud mixing ratio [kg/kg] + # Calculate radius in meters [m] if radconstants.ICE_RADII_PARAM == 1: - #Ice cloud effective radius -- [klaus wyser, 1998] + # Ice cloud effective radius -- [klaus wyser, 1998] if TE > radconstants.MAPL_TICE or QC <= 0.0: BB = -2.0 else: - BB = -2.0 + log10(WC / 50.0) * (1.0e-3 * (radconstants.MAPL_TICE - TE)**1.5) + BB = -2.0 + log10(WC / 50.0) * ( + 1.0e-3 * (radconstants.MAPL_TICE - TE) ** 1.5 + ) BB = min(max(BB, -6.0), -2.0) - RADIUS = 377.4 + 203.3 * BB + 37.91 * BB**2 + 2.3696 * BB**3 + RADIUS = 377.4 + 203.3 * BB + 37.91 * BB ** 2 + 2.3696 * BB ** 3 RADIUS = min(150.0e-6, max(5.0e-6, 1.0e-6 * RADIUS)) - else: - #Ice cloud effective radius ----- [Sun, 2001] + else: + # Ice cloud effective radius ----- [Sun, 2001] TC = TE - radconstants.MAPL_TICE ZFSR = 1.2351 + 0.0105 * TC - AA = 45.8966 * (WC**0.2214) - BB = 0.79570 * (WC**0.2535) + AA = 45.8966 * (WC ** 0.2214) + BB = 0.79570 * (WC ** 0.2535) RADIUS = ZFSR * (AA + BB * (TE - 83.15)) RADIUS = min(150.0e-6, max(5.0e-6, 1.0e-6 * RADIUS * 0.64952)) return RADIUS + @gtscript.function def cloud_effective_radius_liquid( PL: Float, - TE: Float, - QC: Float, + TE: Float, + QC: Float, NNL: Float, NNI: Float, -)-> Float: - ''' +) -> Float: + """ Calculate the effective radius of liquid clouds [m] Implementation of LDRADIUS4 for liquid clouds @@ -88,31 +93,47 @@ def cloud_effective_radius_liquid( Returns: Float: Effective radius of liquid clouds. - ''' - #Calculate liquid water content - WC = 1.0e3 * air_density(PL,TE) * QC #air density [g/m3] * liquid cloud mixing ratio [kg/kg] - #Calculate cloud drop number concentration from the aerosol model + .... + """ + # Calculate liquid water content + WC = ( + 1.0e3 * air_density(PL, TE) * QC + ) # air density [g/m3] * liquid cloud mixing ratio [kg/kg] + # Calculate cloud drop number concentration from the aerosol model + .... NNX = max(NNL * 1.0e-6, 10.0) - #Calculate Radius in meters [m] + # Calculate Radius in meters [m] if radconstants.LIQ_RADII_PARAM == 1: - #Jason Version - RADIUS = min(60.0e-6, max(2.5e-6, 1.0e-6 * radconstants.BX * (WC / NNX)**radconstants.R13BBETA * radconstants.ABETA * 6.92)) + # Jason Version + RADIUS = min( + 60.0e-6, + max( + 2.5e-6, + 1.0e-6 + * radconstants.BX + * (WC / NNX) ** radconstants.R13BBETA + * radconstants.ABETA + * 6.92, + ), + ) else: - #[liu&daum, 2000 and 2005. liu et al 2008] - RADIUS = min(60.0e-6, max(2.5e-6, 1.0e-6 * radconstants.LBX * (WC / NNX)**radconstants.LBE)) + # [liu&daum, 2000 and 2005. liu et al 2008] + RADIUS = min( + 60.0e-6, + max(2.5e-6, 1.0e-6 * radconstants.LBX * (WC / NNX) ** radconstants.LBE), + ) return RADIUS + def _fix_up_clouds_stencil( QV: FloatField, TE: FloatField, - QLC: FloatField, + QLC: FloatField, QIC: FloatField, - CF: FloatField, - QLA: FloatField, + CF: FloatField, + QLA: FloatField, QIA: FloatField, - AF: FloatField, -)->None: - ''' + AF: FloatField, +) -> None: + """ Fix cloud variables to ensure physical consistency. Parameters: @@ -124,89 +145,90 @@ def _fix_up_clouds_stencil( QLA (FloatField): Anvil liquid cloud mixing ratio. QIA (FloatField): Anvil ice cloud mixing ratio. AF (FloatField): Anvil cloud fraction. - ''' + """ with computation(PARALLEL), interval(...): - #Fix if Anvil cloud fraction too small - if AF < 1.E-5: + # Fix if Anvil cloud fraction too small + if AF < 1.0e-5: QV = QV + QLA + QIA TE = TE - (radconstants.ALHLBCP) * QLA - (radconstants.ALHSBCP) * QIA AF = 0.0 QLA = 0.0 QIA = 0.0 - #Fix if LS cloud fraction too small - if CF < 1.E-5: + # Fix if LS cloud fraction too small + if CF < 1.0e-5: QV = QV + QLC + QIC TE = TE - (radconstants.ALHLBCP) * QLC - (radconstants.ALHSBCP) * QIC CF = 0.0 QLC = 0.0 QIC = 0.0 - #LS LIQUID too small - if QLC < 1.E-8: + # LS LIQUID too small + if QLC < 1.0e-8: QV = QV + QLC TE = TE - (radconstants.ALHLBCP) * QLC QLC = 0.0 - #LS ICE too small - if QIC < 1.E-8: + # LS ICE too small + if QIC < 1.0e-8: QV = QV + QIC TE = TE - (radconstants.ALHSBCP) * QIC QIC = 0.0 - #Anvil LIQUID too small - if QLA < 1.E-8: + # Anvil LIQUID too small + if QLA < 1.0e-8: QV = QV + QLA TE = TE - (radconstants.ALHLBCP) * QLA QLA = 0.0 - #Anvil ICE too small - if QIA < 1.E-8: + # Anvil ICE too small + if QIA < 1.0e-8: QV = QV + QIA TE = TE - (radconstants.ALHSBCP) * QIA QIA = 0.0 - #Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small - if (QLA + QIA) < 1.E-8: + # Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small + if (QLA + QIA) < 1.0e-8: QV = QV + QLA + QIA TE = TE - (radconstants.ALHLBCP) * QLA - (radconstants.ALHSBCP) * QIA AF = 0.0 QLA = 0.0 QIA = 0.0 - #Fix ALL cloud quants if LS cloud LIQUID+ICE too small - if (QLC + QIC) < 1.E-8: + # Fix ALL cloud quants if LS cloud LIQUID+ICE too small + if (QLC + QIC) < 1.0e-8: QV = QV + QLC + QIC TE = TE - (radconstants.ALHLBCP) * QLC - (radconstants.ALHSBCP) * QIC CF = 0.0 QLC = 0.0 QIC = 0.0 + def _radcouple_stencil( - TE: FloatField, - PL: FloatField, - CF: FloatField, - AF: FloatField, - QV: FloatField, - QClLS: FloatField, - QCiLS: FloatField, - QClAN: FloatField, - QCiAN: FloatField, - QRN_ALL: FloatField, - QSN_ALL: FloatField, - QGR_ALL: FloatField, - NL: FloatField, - NI: FloatField, - RAD_QV: FloatField, - RAD_QL: FloatField, - RAD_QI: FloatField, - RAD_QR: FloatField, - RAD_QS: FloatField, - RAD_QG: FloatField, - RAD_CF: FloatField, - RAD_RL: FloatField, - RAD_RI: FloatField, - FAC_RL: Float, - MIN_RL: Float, - MAX_RL: Float, - FAC_RI: Float, - MIN_RI: Float, - MAX_RI: Float, -)-> None: - ''' + TE: FloatField, + PL: FloatField, + CF: FloatField, + AF: FloatField, + QV: FloatField, + QClLS: FloatField, + QCiLS: FloatField, + QClAN: FloatField, + QCiAN: FloatField, + QRN_ALL: FloatField, + QSN_ALL: FloatField, + QGR_ALL: FloatField, + NL: FloatField, + NI: FloatField, + RAD_QV: FloatField, + RAD_QL: FloatField, + RAD_QI: FloatField, + RAD_QR: FloatField, + RAD_QS: FloatField, + RAD_QG: FloatField, + RAD_CF: FloatField, + RAD_RL: FloatField, + RAD_RI: FloatField, + FAC_RL: Float, + MIN_RL: Float, + MAX_RL: Float, + FAC_RI: Float, + MIN_RI: Float, + MAX_RI: Float, +) -> None: + """ Couple radiation with cloud variables to ensure physical consistency. Parameters: @@ -239,19 +261,19 @@ def _radcouple_stencil( FAC_RI (Float): Factor for ice effective radius. MIN_RI (Float): Minimum ice effective radius. MAX_RI (Float): Maximum ice effective radius. - ''' + """ with computation(PARALLEL), interval(...): - #water vapor + # water vapor RAD_QV = QV - #total cloud fraction + # total cloud fraction RAD_CF = max(min(CF + AF, 1.0), 0.0) - if RAD_CF >= 1.e-5: - RAD_QL = (QClLS + QClAN) / RAD_CF if (QClLS + QClAN) >= 1.e-8 else 0.0 - RAD_QI = (QCiLS + QCiAN) / RAD_CF if (QCiLS + QCiAN) >= 1.e-8 else 0.0 - RAD_QR = QRN_ALL / RAD_CF if QRN_ALL >= 1.e-8 else 0.0 - RAD_QS = QSN_ALL / RAD_CF if QSN_ALL >= 1.e-8 else 0.0 - RAD_QG = QGR_ALL / RAD_CF if QGR_ALL >= 1.e-8 else 0.0 + if RAD_CF >= 1.0e-5: + RAD_QL = (QClLS + QClAN) / RAD_CF if (QClLS + QClAN) >= 1.0e-8 else 0.0 + RAD_QI = (QCiLS + QCiAN) / RAD_CF if (QCiLS + QCiAN) >= 1.0e-8 else 0.0 + RAD_QR = QRN_ALL / RAD_CF if QRN_ALL >= 1.0e-8 else 0.0 + RAD_QS = QSN_ALL / RAD_CF if QSN_ALL >= 1.0e-8 else 0.0 + RAD_QG = QGR_ALL / RAD_CF if QGR_ALL >= 1.0e-8 else 0.0 else: RAD_CF = 0.0 RAD_QL = 0.0 @@ -260,17 +282,24 @@ def _radcouple_stencil( RAD_QS = 0.0 RAD_QG = 0.0 - #Cap the high end of condensates + # Cap the high end of condensates RAD_QL = min(RAD_QL, 0.01) RAD_QI = min(RAD_QI, 0.01) RAD_QR = min(RAD_QR, 0.01) RAD_QS = min(RAD_QS, 0.01) RAD_QG = min(RAD_QG, 0.01) - #Liquid radii - Brams formulation with limits - RAD_RL = max(MIN_RL, min(cloud_effective_radius_liquid(PL, TE, RAD_QL, NL, NI) * FAC_RL, MAX_RL)) - #Ice radii - Brams formulation with limits - RAD_RI = max(MIN_RI, min(cloud_effective_radius_ice(PL, TE, RAD_QI, NL, NI) * FAC_RI, MAX_RI)) + # Liquid radii - Brams formulation with limits + RAD_RL = max( + MIN_RL, + min(cloud_effective_radius_liquid(PL, TE, RAD_QL, NL, NI) * FAC_RL, MAX_RL), + ) + # Ice radii - Brams formulation with limits + RAD_RI = max( + MIN_RI, + min(cloud_effective_radius_ice(PL, TE, RAD_QI, NL, NI) * FAC_RI, MAX_RI), + ) + class RadiationCoupling: def __init__( @@ -278,28 +307,31 @@ def __init__( stencil_factory: StencilFactory, quantity_factory: QuantityFactory, do_qa: bool, - )-> None: - ''' + ) -> None: + """ Initialize the RadiationCoupling class. Parameters: stencil_factory (StencilFactory): Factory to create stencils. quantity_factory (QuantityFactory): Factory to create quantities. do_qa (bool): Flag to indicate if QA should be performed. - ''' + """ self._fix_up_clouds = stencil_factory.from_dims_halo( - func=_fix_up_clouds_stencil, + func=_fix_up_clouds_stencil, compute_dims=[X_DIM, Y_DIM, Z_DIM], ) self._radcouple = stencil_factory.from_dims_halo( - func=_radcouple_stencil, + func=_radcouple_stencil, compute_dims=[X_DIM, Y_DIM, Z_DIM], ) self.do_qa = do_qa - #GEOS_GFDL_1M_InterfaceMod.F90:866 not implemented. Implement QSAT0 and QSAT3 logic if diagnostics are needed, found in GEOS/src/Shared/@GMAO_Shared/GEOS_Shared/GEOS_Utilities.F90. + # GEOS_GFDL_1M_InterfaceMod.F90:866 not implemented. Implement QSAT0 and QSAT3 logic if diagnostics are needed, found in GEOS/src/Shared/@GMAO_Shared/GEOS_Shared/GEOS_Utilities.F90. if self.do_qa: - #RHX = Q/GEOS_QSAT( T, PLmb) - raise NotImplementedError("[Radiation Coupling] Diagnostic (do_qa) not implemented. (GEOS_QSAT missing)") + # RHX = Q/GEOS_QSAT( T, PLmb) + raise NotImplementedError( + "[Radiation Coupling] Diagnostic (do_qa) not implemented. (GEOS_QSAT missing)" + ) + def __call__( self, Q: FloatField, @@ -332,7 +364,7 @@ def __call__( MIN_RI: Float, MAX_RI: Float, ): - ''' + """ Perform the radiation coupling process. Parameters: @@ -365,43 +397,45 @@ def __call__( FAC_RI (Float): Factor for ice effective radius. MIN_RI (Float): Minimum ice effective radius. MAX_RI (Float): Maximum ice effective radius. - ''' - self._fix_up_clouds(QV = Q, - TE = T, - QLC = QLLS, - QIC = QILS, - CF = CLLS, - QLA = QLCN, - QIA = QICN, - AF = CLCN, - ) - self._radcouple(TE = T, - PL = PLmb, - CF = CLLS, - AF = CLCN, - QV = Q, - QClLS = QLLS, - QCiLS = QILS, - QClAN = QLCN, - QCiAN = QICN, - QRN_ALL = QRAIN, - QSN_ALL = QSNOW, - QGR_ALL = QGRAUPEL, - NL = NACTL, - NI = NACTI, - RAD_QV = RAD_QV, - RAD_QL = RAD_QL, - RAD_QI = RAD_QI, - RAD_QR = RAD_QR, - RAD_QS = RAD_QS, - RAD_QG = RAD_QG, - RAD_CF = RAD_CF, - RAD_RL = CLDREFFL, - RAD_RI = CLDREFFI, - FAC_RL = FAC_RL, - MIN_RL = MIN_RL, - MAX_RL = MAX_RL, - FAC_RI = FAC_RI, - MIN_RI = MIN_RI, - MAX_RI = MAX_RI, - ) \ No newline at end of file + """ + self._fix_up_clouds( + QV=Q, + TE=T, + QLC=QLLS, + QIC=QILS, + CF=CLLS, + QLA=QLCN, + QIA=QICN, + AF=CLCN, + ) + self._radcouple( + TE=T, + PL=PLmb, + CF=CLLS, + AF=CLCN, + QV=Q, + QClLS=QLLS, + QCiLS=QILS, + QClAN=QLCN, + QCiAN=QICN, + QRN_ALL=QRAIN, + QSN_ALL=QSNOW, + QGR_ALL=QGRAUPEL, + NL=NACTL, + NI=NACTI, + RAD_QV=RAD_QV, + RAD_QL=RAD_QL, + RAD_QI=RAD_QI, + RAD_QR=RAD_QR, + RAD_QS=RAD_QS, + RAD_QG=RAD_QG, + RAD_CF=RAD_CF, + RAD_RL=CLDREFFL, + RAD_RI=CLDREFFI, + FAC_RL=FAC_RL, + MIN_RL=MIN_RL, + MAX_RL=MAX_RL, + FAC_RI=FAC_RI, + MIN_RI=MIN_RI, + MAX_RI=MAX_RI, + ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling_constants.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling_constants.py index d7b3e7e21..172ba8231 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling_constants.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/radiation_coupling_constants.py @@ -1,24 +1,24 @@ -'''All constants used for Radiation Coupling Port. Global constants taken from GEOS/src/Shared/@MAPL/shared/Constants''' +"""All constants used for Radiation Coupling Port. Global constants taken from GEOS/src/Shared/@MAPL/shared/Constants""" -#Define whether or not to use CODATA 2018 constants +# Define whether or not to use CODATA 2018 constants CODATA_2018_CONSTANTS = True -#Math Constants -MAPL_PI_R8 = 3.14159265358979323846e0 -MAPL_PI = MAPL_PI_R8 -MAPL_DEGREES_TO_RADIANS_R8 = MAPL_PI_R8 / 180.e0 +# Math Constants +MAPL_PI_R8 = 3.14159265358979323846e0 +MAPL_PI = MAPL_PI_R8 +MAPL_DEGREES_TO_RADIANS_R8 = MAPL_PI_R8 / 180.0e0 MAPL_DEGREES_TO_RADIANS = MAPL_PI / 180.0 -MAPL_RADIANS_TO_DEGREES = 180.e0 / MAPL_PI_R8 +MAPL_RADIANS_TO_DEGREES = 180.0e0 / MAPL_PI_R8 # Universal Constants if CODATA_2018_CONSTANTS: - MAPL_RUNIV = 8314.462618 # J/(Kmole K) -else: - MAPL_RUNIV = 8314.47 # J/(Kmole K) + MAPL_RUNIV = 8314.462618 # J/(Kmole K) +else: + MAPL_RUNIV = 8314.47 # J/(Kmole K) -MAPL_LATENT_HEAT_VAPORIZATION = 2.4665E6 # J/kg @15C @1atm +MAPL_LATENT_HEAT_VAPORIZATION = 2.4665e6 # J/kg @15C @1atm MAPL_ALHL = MAPL_LATENT_HEAT_VAPORIZATION # J/kg -MAPL_LATENT_HEAT_FUSION = 3.3370E5 # J/kg @1atm +MAPL_LATENT_HEAT_FUSION = 3.3370e5 # J/kg @1atm MAPL_ALHF = MAPL_LATENT_HEAT_FUSION # J/kg MAPL_LATENT_HEAT_SUBLIMATION = MAPL_ALHL + MAPL_ALHF # J/kg MAPL_ALHS = MAPL_LATENT_HEAT_SUBLIMATION # J/kg @@ -31,13 +31,13 @@ MAPL_RGAS = MAPL_RDRY # J/(kg K) (DEPRECATED) MAPL_CP = MAPL_RGAS / MAPL_KAPPA # J/(kg K) (DEPRECATED) -#Constants for _fix_up_clouds_stencil +# Constants for _fix_up_clouds_stencil ALHLBCP = MAPL_ALHL / MAPL_CP ALHSBCP = MAPL_ALHS / MAPL_CP -#Constants for cloud_effective_radius_liquid and cloud_effective_radius_ice +# Constants for cloud_effective_radius_liquid and cloud_effective_radius_ice MAPL_RGAS = 287.0 -MAPL_TICE = 273.16 #K +MAPL_TICE = 273.16 # K LIQ_RADII_PARAM = 1 ICE_RADII_PARAM = 1 BX = 100.0 * (3.0 / (4.0 * MAPL_PI)) ** (1.0 / 3.0) @@ -49,10 +49,10 @@ LBX = LDISS * 1.0e3 * (3.0 / (4.0 * MAPL_PI * LK * RHO_W * 1.0e-3)) ** (1.0 / 3.0) LBE = 1.0 / 3.0 - 0.14 -#Taken from GEOS_GFDL_1M_InterfaceMod.F90:269-274 - uses MAPL_GetResource to get variables. These do not change in the radiation coupling loop +# Taken from GEOS_GFDL_1M_InterfaceMod.F90:269-274 - uses MAPL_GetResource to get variables. These do not change in the radiation coupling loop MIN_RL = 2.5e-6 MAX_RL = 60.0e-6 FAC_RL = 1.0 MIN_RI = 5.0e-6 MAX_RI = 100.0e-6 -FAC_RI = 1.0 \ No newline at end of file +FAC_RI = 1.0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/types.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/types.py new file mode 100644 index 000000000..aa8f1f7d3 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/pyMoist/types.py @@ -0,0 +1,8 @@ +import gt4py.cartesian.gtscript as gtscript + +import pyMoist.aer_activation_constants as constants +from ndsl.dsl.typing import Float + + +# Global space +FloatField_NModes = gtscript.Field[gtscript.IJK, (Float, (constants.n_modes))] diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py index 8e899652b..e4055491c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/__init__.py @@ -1 +1,2 @@ from .translate_radiation_coupling import TranslateRadCouple +from .translate_aer_activation import TranslateAerActivation diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_aer_activation.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_aer_activation.py new file mode 100644 index 000000000..2ba353fcf --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_aer_activation.py @@ -0,0 +1,128 @@ +from ndsl import Namelist, StencilFactory +from ndsl.dsl.typing import Float +from ndsl.stencils.testing.translate import TranslateFortranData2Py +from pyMoist.aer_activation import AerActivation + + +class TranslateAerActivation(TranslateFortranData2Py): + def __init__(self, grid, namelist: Namelist, stencil_factory: StencilFactory): + super().__init__(grid, stencil_factory) + self.stencil_factory = stencil_factory + self.quantity_factory = grid.quantity_factory + self._grid = grid + self.max_error = 1e-9 + + # FloatField Inputs + self.in_vars["data_vars"] = { + "AERO_F_DUST": {}, + "AERO_F_ORGANIC": {}, + "TMP3D": {}, + "NACTL": {}, + "PLE": {}, + "AERO_F_SOOT": {}, + "T": {}, + "AERO_HYGROSCOPICITY": {}, + "ZLE0": {}, + "QLLS": {}, + "AERO_NUM": {}, + "ZL0": {}, + "PLmb": {}, + "QILS": {}, + "AERO_DGN": {}, + "TKE": {}, + "NACTI": {}, + "NWFA": {}, + "QLCN": {}, + "AERO_SIGMA": {}, + "AERO_DENSITY": {}, + "QICN": {}, + "FRLAND": {}, + } + + # Float/Int Inputs + self.in_vars["parameters"] = [ + "n_modes", + "SH", + "EVAP", + "KPBL", + "CCN_LND", + "CCN_OCN", + ] + + # FloatField Outputs + self.out_vars = { + "NACTL": self.grid.compute_dict(), + "NACTI": self.grid.compute_dict(), + "NWFA": self.grid.compute_dict(), + } + + def compute(self, inputs): + aer_activation = AerActivation( + self.stencil_factory, + self.quantity_factory, + int(inputs["n_modes"]), + USE_AERSOL_NN=True, + ) + + # Outputs + nactl = inputs["NACTL"].astype(Float) + nacti = inputs["NACTI"].astype(Float) + nwfa = inputs["NWFA"].astype(Float) + + # Inputs + aero_f_dust = inputs["AERO_F_DUST"].astype(Float) + aero_f_organic = inputs["AERO_F_ORGANIC"].astype(Float) + tmp3d = inputs["TMP3D"].astype(Float) + nactl = inputs["NACTL"].astype(Float) + ple = inputs["PLE"].astype(Float) + aero_f_soot = inputs["AERO_F_SOOT"].astype(Float) + t = inputs["T"].astype(Float) + aero_hygroscopicity = inputs["AERO_HYGROSCOPICITY"].astype(Float) + zle0 = inputs["ZLE0"].astype(Float) + qlls = inputs["QLLS"].astype(Float) + aero_num = inputs["AERO_NUM"].astype(Float) + zl0 = inputs["ZL0"].astype(Float) + plmb = inputs["PLmb"].astype(Float) + qils = inputs["QILS"].astype(Float) + aero_dgn = inputs["AERO_DGN"].astype(Float) + tke = inputs["TKE"].astype(Float) + nacti = inputs["NACTI"].astype(Float) + nwfa = inputs["NWFA"].astype(Float) + qlcn = inputs["QLCN"].astype(Float) + aero_sigma = inputs["AERO_SIGMA"].astype(Float) + aero_density = inputs["AERO_DENSITY"].astype(Float) + qicn = inputs["QICN"].astype(Float) + n_modes = inputs["n_modes"] + sh = inputs["SH"].astype(Float) + evap = inputs["EVAP"].astype(Float) + kpbl = inputs["KPBL"].astype(Float) + frland = inputs["FRLAND"].astype(Float) + ccn_lnd = Float(inputs["CCN_LND"]) + ccn_ocn = Float(inputs["CCN_OCN"]) + + aer_activation( + aero_dgn=aero_dgn, + aero_num=aero_num, + nacti=nacti, + t=t, + plo=plmb * 100.0, + qicn=qicn, + qils=qils, + qlcn=qlcn, + qlls=qlls, + nn_land=Float(ccn_lnd * 1.0e6), + frland=frland, + nn_ocean=Float(ccn_ocn * 1.0e6), + aero_hygroscopicity=aero_hygroscopicity, + nwfa=nwfa, + nactl=nactl, + vvel=tmp3d, + tke=tke, + aero_sigma=aero_sigma, + ) + + return { + "NACTL": nactl, + "NACTI": nacti, + "NWFA": nwfa, + } diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_radiation_coupling.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_radiation_coupling.py index 23053bbab..7aee1f047 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_radiation_coupling.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/savepoint/translate_radiation_coupling.py @@ -1,25 +1,26 @@ from ndsl import Namelist, StencilFactory +from ndsl.constants import X_DIM, Y_DIM, Z_DIM from ndsl.stencils.testing.translate import TranslateFortranData2Py from pyMoist.radiation_coupling import RadiationCoupling -from ndsl.constants import X_DIM, Y_DIM, Z_DIM + class TranslateRadCouple(TranslateFortranData2Py): def __init__( self, grid, - namelist: Namelist, + namelist: Namelist, stencil_factory: StencilFactory, ): super().__init__(grid, stencil_factory) self.compute_func = RadiationCoupling( # type: ignore self.stencil_factory, self.grid.quantity_factory, - do_qa = False, #Change to do_qa=namelist.do_qa, if QSAT module procedures QSAT0 and QSAT3 are implemented + do_qa=False, # Change to do_qa=namelist.do_qa, if QSAT module procedures QSAT0 and QSAT3 are implemented ) self._grid = grid self.max_error = 1e-9 - #FloatField Inputs + # FloatField Inputs self.in_vars["data_vars"] = { "Q": self.grid.compute_dict(), "T": self.grid.compute_dict(), @@ -35,21 +36,29 @@ def __init__( "QGRAUPEL": self.grid.compute_dict(), "NACTL": self.grid.compute_dict(), "NACTI": self.grid.compute_dict(), + "RAD_QV": self.grid.compute_dict(), + "RAD_QL": self.grid.compute_dict(), + "RAD_QI": self.grid.compute_dict(), + "RAD_QR": self.grid.compute_dict(), + "RAD_QS": self.grid.compute_dict(), + "RAD_QG": self.grid.compute_dict(), + "RAD_CF": self.grid.compute_dict(), + "CLDREFFI": self.grid.compute_dict(), + "CLDREFFL": self.grid.compute_dict(), } - #Float Inputs - self.in_vars["parameters"] = ['FAC_RL', 'MIN_RL', 'MAX_RL', 'FAC_RI', 'MIN_RI', 'MAX_RI'] + # Float Inputs + self.in_vars["parameters"] = [ + "FAC_RL", + "MIN_RL", + "MAX_RL", + "FAC_RI", + "MIN_RI", + "MAX_RI", + ] - #FloatField Outputs + # FloatField Outputs self.out_vars = { - "Q": self.grid.compute_dict(), - "T": self.grid.compute_dict(), - "QLLS": self.grid.compute_dict(), - "QILS": self.grid.compute_dict(), - "CLLS": self.grid.compute_dict(), - "QLCN": self.grid.compute_dict(), - "QICN": self.grid.compute_dict(), - "CLCN": self.grid.compute_dict(), "RAD_QV": self.grid.compute_dict(), "RAD_QL": self.grid.compute_dict(), "RAD_QI": self.grid.compute_dict(), @@ -61,19 +70,7 @@ def __init__( "CLDREFFL": self.grid.compute_dict(), } - #Calculated Outputs - def compute_from_storage(self, inputs): - #Original test data was missing these allocations. New test data will be generated that contain these allocations. - outputs = { - "RAD_QV": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "RAD_QL": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "RAD_QI": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "RAD_QR": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "RAD_QS": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "RAD_QG": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "RAD_CF": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "CLDREFFI": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown"), - "CLDREFFL": self._grid.quantity_factory.zeros(dims=[X_DIM, Y_DIM, Z_DIM], units="unknown") - } - self.compute_func(**inputs, **outputs) - return {**outputs, "T": inputs["T"], "Q": inputs["Q"], "QLLS": inputs["QLLS"], "QILS": inputs["QILS"], "CLLS": inputs["CLLS"], "QLCN": inputs["QLCN"], "QICN": inputs["QICN"], "CLCN": inputs["CLCN"]} + def compute(self, inputs): + self.make_storage_data_input_vars(inputs) + self.compute_func(**inputs) + return self.slice_output(inputs) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/get_data.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/get_data.sh index 9baac7540..f869bf494 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/get_data.sh +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/get_data.sh @@ -4,8 +4,12 @@ set -e -x TEST_DATA_PATH="../../test_data" mkdir -p $TEST_DATA_PATH +wget https://portal.nccs.nasa.gov/datashare/astg/smt/geos-fp/translate/11.5.2/Moist/AerActivation_nc.tar.gz +mv AerActivation_nc.tar.gz $TEST_DATA_PATH wget https://portal.nccs.nasa.gov/datashare/astg/smt/geos-fp/translate/11.5.2/Moist/rad_coup_data_nc.tar.gz mv rad_coup_data_nc.tar.gz $TEST_DATA_PATH cd $TEST_DATA_PATH +tar -xzvf AerActivation_nc.tar.gz +rm AerActivation_nc.tar.gz tar -xzvf rad_coup_data_nc.tar.gz rm rad_coup_data_nc.tar.gz diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh index d4b6c83dc..949b178a2 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSmoist_GridComp/pyMoist/tests/scripts/run_tests.sh @@ -1,9 +1,10 @@ #!/bin/bash export PACE_FLOAT_PRECISION=32 +export FV3_DACEMODE=Python python -m pytest -v -s -x\ - --data_path=../../test_data/rad_coup_data_nc \ - --backend=numpy \ - --which_rank=0 \ - --which_modules=RadCouple \ + --pdb \ + --data_path=../../test_data/geos_11.5.2/moist \ + --backend=dace:cpu \ + --which_modules=AerActivation \ --grid=default \ - .. \ No newline at end of file + ..