add new feature for getvar, pwhi and pwlow are available

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

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
+

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)
+

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))

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)
+
+
+

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"],