diff --git a/doc/source/internal_api/index.rst b/doc/source/internal_api/index.rst index 089bbb0..736442b 100644 --- a/doc/source/internal_api/index.rst +++ b/doc/source/internal_api/index.rst @@ -32,6 +32,8 @@ The routines below are called internally by :meth:`wrf.getvar`. wrf.g_pressure.get_pressure wrf.g_pressure.get_pressure_hpa wrf.g_pw.get_pw + wrf.g_pw.get_pwhi + wrf.g_pw.get_pwlow wrf.g_rh.get_rh wrf.g_rh.get_rh_2m wrf.g_slp.get_slp diff --git a/fortran/wrf_pw.f90 b/fortran/wrf_pw.f90 index bd80876..8914f43 100644 --- a/fortran/wrf_pw.f90 +++ b/fortran/wrf_pw.f90 @@ -1,3 +1,79 @@ +!NCLFORTSTART +SUBROUTINE DCOMPUTEPWHI(p, tv, qv, ht, pwhi, nx, ny, nz, nzh) + USE wrf_constants, ONLY : RD + + IMPLICIT NONE + + !f2py threadsafe + !f2py intent(in,out) :: pwhi + + INTEGER, INTENT(IN) :: nx, ny, nz, nzh + REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: p, tv, qv + REAL(KIND=8), DIMENSION(nx,ny,nzh), INTENT(IN) :: ht + REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: pwhi + +!NCLEND + + INTEGER :: i, j, k + pwhi = 0 + + !$OMP PARALLEL + DO k=1,nz + !$OMP DO COLLAPSE(2) SCHEDULE(runtime) + DO j=1,ny + DO i=1,nx + IF (tv(i,j,k) <= 273.15) THEN + pwhi(i,j) = pwhi(i,j) + ((p(i,j,k)/(RD*tv(i,j,k)))*qv(i,j,k)*(ht(i,j,k+1) - ht(i,j,k))) + ENDIF + END DO + END DO + !$OMP END DO + END DO + !$OMP END PARALLEL + + RETURN + +END SUBROUTINE DCOMPUTEPWHI + + +!NCLFORTSTART +SUBROUTINE DCOMPUTEPWLOW(p, tv, qv, ht, pwlow, nx, ny, nz, nzh) + USE wrf_constants, ONLY : RD + + IMPLICIT NONE + + !f2py threadsafe + !f2py intent(in,out) :: pwlow + + INTEGER, INTENT(IN) :: nx, ny, nz, nzh + REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: p, tv, qv + REAL(KIND=8), DIMENSION(nx,ny,nzh), INTENT(IN) :: ht + REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: pwlow + +!NCLEND + + INTEGER :: i, j, k + pwlow = 0 + + !$OMP PARALLEL + DO k=1,nz + !$OMP DO COLLAPSE(2) SCHEDULE(runtime) + DO j=1,ny + DO i=1,nx + IF (tv(i,j,k) > 273.15) THEN + pwlow(i,j) = pwlow(i,j) + ((p(i,j,k)/(RD*tv(i,j,k)))*qv(i,j,k)*(ht(i,j,k+1) - ht(i,j,k))) + ENDIF + END DO + END DO + !$OMP END DO + END DO + + !$OMP END PARALLEL + + RETURN + +END SUBROUTINE DCOMPUTEPWLOW + !NCLFORTSTART SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh) USE wrf_constants, ONLY : RD @@ -34,5 +110,5 @@ SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh) !$OMP END PARALLEL RETURN - END SUBROUTINE DCOMPUTEPW + diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 5241230..3a33fc1 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -7,7 +7,7 @@ from .constants import default_fill from .extension import (_interpz3d, _interp2dxy, _interp1d, _slp, _tk, _td, _rh, _uvmet, _smooth2d, _cape, _cloudfrac, _ctt, _dbz, _srhel, _udhel, _avo, _pvo, _eth, _wetbulb, _tv, - _omega, _pw) + _omega, _pw, _pwhi, _pwlow) from .decorators import convert_units from .metadecorators import (set_alg_metadata, set_uvmet_alg_metadata, set_interp_metadata, set_cape_alg_metadata, @@ -1955,3 +1955,18 @@ def pw(pres, tkel, qv, height, meta=True): tv = _tv(tkel, qv) return _pw(pres, tv, qv, height) + + +@set_alg_metadata(2, "pres", refvarndims=3, units="kg m-2", + description="precipitable water above ZERO isoterm") +def pwhi(pres, tkel, qv, height, meta=True): + tv = _tv(tkel, qv) + return _pwhi(pres, tv, qv, height) + + +@set_alg_metadata(2, "pres", refvarndims=3, units="kg m-2", + description="precipitable water below ZERO isoterm") +def pwlow(pres, tkel, qv, height, meta=True): + tv = _tv(tkel, qv) + return _pwlow(pres, tv, qv, height) + diff --git a/src/wrf/extension.py b/src/wrf/extension.py index d28d758..02bb9c4 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -11,6 +11,7 @@ from wrf._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d, dcalrelhl, dcalcuh, dcomputepv, dcomputeabsvort, dlltoij, dijtoll, deqthecalc, omgcalc, virtual_temp, wetbulbcalc, dcomputepw, + dcomputepwhi, dcomputepwlow, wrf_monotonic, wrf_vintrp, dcomputewspd, dcomputewdir, dinterp3dz_2dlev, fomp_set_num_threads, fomp_get_num_threads, @@ -592,6 +593,39 @@ def _pw(p, tv, qv, ht, outview=None): return result +@check_args(0, 3, (3, 3, 3, 3), stagger=(None, None, None, -3)) +@left_iteration(3, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0, 1, 2, 3)) +@extract_and_transpose() +def _pwhi(p, tv, qv, ht, outview=None): + if outview is None: + outview = np.empty(p.shape[0:2], p.dtype, order="F") + + result = dcomputepwhi(p, + tv, + qv, + ht, + outview) + + return result + + +@check_args(0, 3, (3, 3, 3, 3), stagger=(None, None, None, -3)) +@left_iteration(3, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0, 1, 2, 3)) +@extract_and_transpose() +def _pwlow(p, tv, qv, ht, outview=None): + if outview is None: + outview = np.empty(p.shape[0:2], p.dtype, order="F") + result = dcomputepwlow(p, + tv, + qv, + ht, + outview) + + return result + + @check_args(0, 3, (3, 3, 3, 3, 3, 3)) @left_iteration(3, 3, ref_var_idx=0, ignore_args=(6, 7, 8)) @cast_type(arg_idxs=(0, 1, 2, 3, 4, 5)) diff --git a/src/wrf/g_pw.py b/src/wrf/g_pw.py index b63d252..0c0b0dc 100755 --- a/src/wrf/g_pw.py +++ b/src/wrf/g_pw.py @@ -1,6 +1,6 @@ from __future__ import (absolute_import, division, print_function) -from .extension import _pw, _tv, _tk +from .extension import _pw, _tv, _tk, _pwhi, _pwlow from .constants import Constants from .util import extract_vars from .metadecorators import copy_and_set_metadata @@ -84,3 +84,59 @@ def get_pw(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, tv = _tv(tk, qv) return _pw(full_p, tv, qv, ht) + + +@copy_and_set_metadata(copy_varname="T", name="pwhi", + remove_dims=("bottom_top",), + description="precipitable water above ZERO isoterm", + MemoryOrder="XY", + units="kg m-2") +def get_pwhi(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, + meta=True, _key=None): + varnames = ("T", "P", "PB", "PH", "PHB", "QVAPOR") + ncvars = extract_vars(wrfin, timeidx, varnames, method, squeeze, cache, + meta=False, _key=_key) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + qv = ncvars["QVAPOR"] + + full_p = p + pb + ht = (ph + phb)/Constants.G + full_t = t + Constants.T_BASE + + tk = _tk(full_p, full_t) + tv = _tv(tk, qv) + return _pwhi(full_p, tv, qv, ht) + + +@copy_and_set_metadata(copy_varname="T", name="pwlow", + remove_dims=("bottom_top",), + description="precipitable water below ZERO isoterm", + MemoryOrder="XY", + units="kg m-2") +def get_pwlow(wrfin, timeidx=0, method="cat", squeeze=True, cache=None, + meta=True, _key=None): + varnames = ("T", "P", "PB", "PH", "PHB", "QVAPOR") + ncvars = extract_vars(wrfin, timeidx, varnames, method, squeeze, cache, + meta=False, _key=_key) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + qv = ncvars["QVAPOR"] + + full_p = p + pb + ht = (ph + phb)/Constants.G + full_t = t + Constants.T_BASE + + tk = _tk(full_p, full_t) + tv = _tv(tk, qv) + + return _pwlow(full_p, tv, qv, ht) + + + diff --git a/src/wrf/routines.py b/src/wrf/routines.py index 3998209..8136828 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -14,7 +14,7 @@ from .g_helicity import get_srh, get_uh from .g_latlon import get_lat, get_lon from .g_omega import get_omega from .g_pressure import get_pressure, get_pressure_hpa -from .g_pw import get_pw +from .g_pw import get_pw, get_pwhi, get_pwlow from .g_rh import get_rh, get_rh_2m from .g_slp import get_slp from .g_temp import (get_tc, get_eth, get_temp, get_theta, get_tk, get_tv, @@ -47,6 +47,8 @@ _FUNC_MAP = {"cape2d": get_2dcape, "uhel": get_uh, "omega": get_omega, "pw": get_pw, + "pwhi": get_pwhi, + "pwlow": get_pwlow, "rh": get_rh, "rh2m": get_rh_2m, "slp": get_slp, @@ -112,6 +114,8 @@ _VALID_KARGS = {"cape2d": ["missing"], "uhel": ["bottom", "top"], "omega": [], "pw": [], + "pwhi": [], + "pwlow": [], "rh": [], "rh2m": [], "slp": ["units"],