Browse Source

Fix issues related to cloud top temperature.

Fixed indexing bug.
Fixed incorrect computation of optical depth when cloud ice is not available.
Users can use fill values for cloud free areas.
Users can now specify the optical depth threshold that triggers the calculation.
Fixes #45.
lon0
Bill Ladwig 7 years ago
parent
commit
cc86d68a35
  1. 6
      doc/source/_templates/product_table.txt
  2. 20
      doc/source/new.rst
  3. 54
      fortran/wrf_fctt.f90
  4. 5
      fortran/wrffortran.pyf
  5. 34
      src/wrf/computation.py
  6. 10
      src/wrf/extension.py
  7. 38
      src/wrf/g_ctt.py
  8. 4
      src/wrf/metadecorators.py
  9. 2
      src/wrf/routines.py
  10. 2
      src/wrf/version.py
  11. 2
      test/ci_tests/make_test_file.py
  12. 4
      test/comp_utest.py
  13. 214
      test/ipynb/Test_CTT.ipynb
  14. 214
      test/ipynb/Test_CTT_Prod.ipynb

6
doc/source/_templates/product_table.txt

@ -13,11 +13,13 @@
+--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+
| cape_3d | 3D cape and cin | J kg-1 | **missing** (float): Fill value for output only | | cape_3d | 3D cape and cin | J kg-1 | **missing** (float): Fill value for output only |
+--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+
| ctt | Cloud Top Temperature | degC | **units** (str) : Set to desired units. Default is *'degC'*. | | ctt | Cloud Top Temperature | degC | **fill_nocloud** (boolean): Set to True to use fill values for cloud free regions rather than surface temperature. Default is *False*. |
| | | | | | | | | |
| | | K | | | | | K | **missing** (float): The fill value to use when *fill_nocloud* is True. |
| | | | | | | | | |
| | | | **opt_thresh** (float): The optical depth required to trigger the cloud top temperature calculation. Default is 1.0. |
| | | degF | | | | | degF | |
| | | | **units** (str) : Set to desired units. Default is *'degC'*. |
+--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+
| cloudfrac | Cloud Fraction | % | **vert_type** (str): The vertical coordinate type for the cloud thresholds. Must be 'height_agl', 'height_msl', or 'pres'. Default is 'height_agl'. | | cloudfrac | Cloud Fraction | % | **vert_type** (str): The vertical coordinate type for the cloud thresholds. Must be 'height_agl', 'height_msl', or 'pres'. Default is 'height_agl'. |
| | | | | | | | | |

20
doc/source/new.rst

@ -4,6 +4,26 @@ What's New
Releases Releases
------------- -------------
v1.1.3
^^^^^^^^^^^^^^
- Release 1.1.3
- Fixed/Enhanced the cloud top temperature diagnostic.
- Optical depth was not being calculated correctly when
cloud ice mixing ratio was not available.
- Fixed an indexing bug that caused crashes on Windows, but should have been
crashing on all platforms.
- Users can now specify if they want cloud free regions to use fill values,
rather than the default behavior of using the surface temperature.
- Users can now specify the optical depth required to trigger the cloud
top temperature calculation. However, the default value of 1.0 should be
sufficient for most users.
- Added 'th' alias for the theta product.
- Fixed a crash issue related to updraft helicity when a dictionary is
used as the input.
- The cape_2d diagnostic can now work with a single column of data, just like
cape_3d.
v1.1.2 v1.1.2
^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^

54
fortran/wrf_fctt.f90

@ -1,5 +1,6 @@
!NCLFORTSTART !NCLFORTSTART
SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns, ew) SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci,&
fill_nocloud, missing, opt_thresh, nz, ns, ew)
USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL
IMPLICIT NONE IMPLICIT NONE
@ -7,13 +8,15 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns
!f2py threadsafe !f2py threadsafe
!f2py intent(in,out) :: ctt !f2py intent(in,out) :: ctt
INTEGER, INTENT(IN) :: nz, ns, ew, haveqci INTEGER, INTENT(IN) :: nz, ns, ew, haveqci, fill_nocloud
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght, prs, tk, qci, qcw, qvp REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght, prs, tk, qci, qcw, qvp
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter
REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt
REAL(KIND=8), INTENT(IN) :: missing
REAL(KIND=8), INTENT(IN) :: opt_thresh
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf
!NCLEND !NCLEND
! REAL(KIND=8) :: znfac(nz) ! REAL(KIND=8) :: znfac(nz)
@ -60,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns
DO i=1,ew DO i=1,ew
opdepthd = 0.D0 opdepthd = 0.D0
k = 0 k = 0
prsctt = 0 prsctt = -1
! Integrate downward from model top, calculating path at full ! Integrate downward from model top, calculating path at full
! model vertical levels. ! model vertical levels.
DO k=1, nz DO k=2,nz
opdepthu = opdepthd opdepthu = opdepthd
ripk = nz - k + 1 ripk = nz - k + 1
@ -76,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns
END IF END IF
IF (haveqci .EQ. 0) then IF (haveqci .EQ. 0) then
IF (tk(i,j,k) .LT. CELKEL) then IF (tk(i,j,ripk) .LT. CELKEL) then
! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed ! Note: abscoefi is m**2/g, qcw is g/kg, so no convrsion needed
opdepthd = opdepthu + ABSCOEFI*qcw(i,j,k) * dp/G opdepthd = opdepthu + ABSCOEFI*qcw(i,j,ripk) * dp/G
ELSE ELSE
opdepthd = opdepthu + ABSCOEF*qcw(i,j,k) * dp/G opdepthd = opdepthu + ABSCOEF*qcw(i,j,ripk) * dp/G
END IF END IF
ELSE ELSE
opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G
END IF END IF
IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN IF (opdepthd .LT. opt_thresh .AND. k .LT. nz) THEN
CYCLE CYCLE
ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN ELSE IF (opdepthd .LT. opt_thresh .AND. k .EQ. nz) THEN
prsctt = prs(i,j,1) IF (fill_nocloud .EQ. 0) THEN
prsctt = prs(i,j,1)
ENDIF
EXIT EXIT
ELSE ELSE
fac = (1. - opdepthu)/(opdepthd - opdepthu) fac = (1. - opdepthu)/(opdepthd - opdepthu)
@ -100,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns
END IF END IF
END DO END DO
DO k=2,nz ! prsctt should only be 0 if fill values are used
ripk = nz - k + 1 IF (prsctt .GT. -1) THEN
p1 = prs(i,j,ripk+1) DO k=2,nz
p2 = prs(i,j,ripk) ripk = nz - k + 1
IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN p1 = prs(i,j,ripk+1)
fac = (prsctt - p1)/(p2 - p1) p2 = prs(i,j,ripk)
arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN
ctt(i,j) = tk(i,j,ripk+1) + arg1 fac = (prsctt - p1)/(p2 - p1)
EXIT arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL
END IF ctt(i,j) = tk(i,j,ripk+1) + arg1
END DO EXIT
END IF
END DO
ELSE
ctt(i,j) = missing
END IF
END DO END DO
END DO END DO
!$OMP END DO !$OMP END DO

5
fortran/wrffortran.pyf

@ -359,7 +359,7 @@ python module _wrffortran ! in
real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285714285714 real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285714285714
real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25864379523 real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25864379523
end module wrf_constants end module wrf_constants
subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,fill_nocloud,missing,opt_thresh,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90
threadsafe threadsafe
use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel
real(kind=8) dimension(ew,ns,nz),intent(in) :: prs real(kind=8) dimension(ew,ns,nz),intent(in) :: prs
@ -372,6 +372,9 @@ python module _wrffortran ! in
real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: ctt real(kind=8) dimension(ew,ns),intent(out,in),depend(ew,ns) :: ctt
real(kind=8) dimension(ew,ns,nz),intent(inout),depend(ew,ns,nz) :: pf real(kind=8) dimension(ew,ns,nz),intent(inout),depend(ew,ns,nz) :: pf
integer intent(in) :: haveqci integer intent(in) :: haveqci
integer intent(in) :: fill_nocloud
real(kind=8) intent(in) :: missing
real(kind=8) intent(in) :: opt_thresh
integer, optional,intent(in),check(shape(prs,2)==nz),depend(prs) :: nz=shape(prs,2) integer, optional,intent(in),check(shape(prs,2)==nz),depend(prs) :: nz=shape(prs,2)
integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1) integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1)
integer, optional,intent(in),check(shape(prs,0)==ew),depend(prs) :: ew=shape(prs,0) integer, optional,intent(in),check(shape(prs,0)==ew),depend(prs) :: ew=shape(prs,0)

34
src/wrf/computation.py

@ -1015,6 +1015,11 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh,
high_thresh (:obj:`float`): The bottom vertical threshold for what is high_thresh (:obj:`float`): The bottom vertical threshold for what is
considered a high cloud. considered a high cloud.
missing (:obj:`float:`, optional): The fill value to use for areas
where the surface is higher than the cloud threshold level
(e.g. mountains). Default is
:data:`wrf.default_fill(numpy.float64)`.
meta (:obj:`bool`): Set to False to disable metadata and return meta (:obj:`bool`): Set to False to disable metadata and return
:class:`numpy.ndarray` instead of :class:`numpy.ndarray` instead of
:class:`xarray.DataArray`. Default is True. :class:`xarray.DataArray`. Default is True.
@ -1047,8 +1052,9 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh,
@set_alg_metadata(2, "pres_hpa", refvarndims=3, @set_alg_metadata(2, "pres_hpa", refvarndims=3,
description="cloud top temperature") description="cloud top temperature")
@convert_units("temp", "c") @convert_units("temp", "c")
def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None,
units="degC"): fill_nocloud=False, missing=default_fill(np.float64),
opt_thresh=1.0, meta=True, units="degC"):
"""Return the cloud top temperature. """Return the cloud top temperature.
This is the raw computational algorithm and does not extract any variables This is the raw computational algorithm and does not extract any variables
@ -1095,6 +1101,23 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True,
Ice mixing ratio in [kg/kg] with the same dimensionality as Ice mixing ratio in [kg/kg] with the same dimensionality as
*pres_hpa*. *pres_hpa*.
fill_nocloud (:obj:`bool`, optional): Set to True to use fill values in
regions where clouds are not detected (optical depth less than 1).
Otherwise, the output will contain the surface temperature for
areas without clouds. Default is False.
missing (:obj:`float`, optional): The fill value to use for areas
where no clouds are detected. Only used if *fill_nocloud* is
True. Default is
:data:`wrf.default_fill(numpy.float64)`.
opt_thresh (:obj:`float`, optional): The required amount of optical
depth (looking from top down) required to trigger a cloud top
temperature calculation. Values less than this threshold will be
treated as no cloud areas. For thin cirrus, a value of .1 might be
appropriate. For large CB clouds, 1000.0 might be appropriate.
Default is 1.0.
meta (:obj:`bool`): Set to False to disable metadata and return meta (:obj:`bool`): Set to False to disable metadata and return
:class:`numpy.ndarray` instead of :class:`numpy.ndarray` instead of
:class:`xarray.DataArray`. Default is True. :class:`xarray.DataArray`. Default is True.
@ -1129,7 +1152,12 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True,
else: else:
haveqci = 1 if qice.any() else 0 haveqci = 1 if qice.any() else 0
return _ctt(pres_hpa, tkel, qice, qcld, qv, height, terrain, haveqci) _fill_nocloud = 1 if fill_nocloud else 0
ctt = _ctt(pres_hpa, tkel, qice, qcld, qv, height, terrain, haveqci,
_fill_nocloud, missing, opt_thresh)
return ma.masked_values(ctt, missing)

10
src/wrf/extension.py

@ -767,10 +767,11 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1,
@check_args(0, 3, (3,3,3,3,3,3,2)) @check_args(0, 3, (3,3,3,3,3,3,2))
@left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,)) @left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,8,9,10))
@cast_type(arg_idxs=(0,1,2,3,4,5,6)) @cast_type(arg_idxs=(0,1,2,3,4,5,6))
@extract_and_transpose() @extract_and_transpose()
def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, fill_nocloud,
missing, opt_thresh, outview=None):
"""Wrapper for wrfcttcalc. """Wrapper for wrfcttcalc.
Located in wrf_fctt.f90. Located in wrf_fctt.f90.
@ -790,7 +791,10 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None):
ter, ter,
outview, outview,
pf, pf,
haveqci) haveqci,
fill_nocloud,
missing,
opt_thresh)
return result return result

38
src/wrf/g_ctt.py

@ -1,11 +1,12 @@
from __future__ import (absolute_import, division, print_function, from __future__ import (absolute_import, division, print_function,
unicode_literals) unicode_literals)
import numpy as n import numpy as np
import numpy.ma as ma
#from .extension import computectt, computetk #from .extension import computectt, computetk
from .extension import _ctt, _tk from .extension import _ctt, _tk
from .constants import Constants, ConversionFactors from .constants import Constants, ConversionFactors, default_fill
from .destag import destagger from .destag import destagger
from .decorators import convert_units from .decorators import convert_units
from .metadecorators import copy_and_set_metadata from .metadecorators import copy_and_set_metadata
@ -19,7 +20,8 @@ from .util import extract_vars
@convert_units("temp", "c") @convert_units("temp", "c")
def get_ctt(wrfin, timeidx=0, method="cat", def get_ctt(wrfin, timeidx=0, method="cat",
squeeze=True, cache=None, meta=True, _key=None, squeeze=True, cache=None, meta=True, _key=None,
units="degC"): fill_nocloud=False, missing=default_fill(np.float64),
opt_thresh=1.0, units="degC"):
"""Return the cloud top temperature. """Return the cloud top temperature.
This functions extracts the necessary variables from the NetCDF file This functions extracts the necessary variables from the NetCDF file
@ -64,6 +66,27 @@ def get_ctt(wrfin, timeidx=0, method="cat",
_key (:obj:`int`, optional): A caching key. This is used for internal _key (:obj:`int`, optional): A caching key. This is used for internal
purposes only. Default is None. purposes only. Default is None.
fill_nocloud (:obj:`bool`, optional): Set to True to use fill values in
regions where clouds are not detected (optical depth less than 1).
Otherwise, the output will contain the surface temperature for
areas without clouds. Default is False.
missing (:obj:`float`, optional): The fill value to use for areas
where no clouds are detected. Only used if *fill_nocloud* is
True. Default is
:data:`wrf.default_fill(numpy.float64)`.
opt_thresh (:obj:`float`, optional): The amount of optical
depth (integrated from top down) required to trigger a cloud top
temperature calculation. The cloud top temperature is calculated at
the vertical level where this threshold is met. Vertical columns
with less than this threshold will be treated as cloud free areas.
In general, the larger the value is for this
threshold, the lower the altitude will be for the cloud top
temperature calculation, and therefore higher cloud top
temperature values. Default is 1.0, which should be sufficient for
most users.
units (:obj:`str`): The desired units. Refer to the :meth:`getvar` units (:obj:`str`): The desired units. Refer to the :meth:`getvar`
product table for a list of available units for 'ctt'. Default product table for a list of available units for 'ctt'. Default
is 'degC'. is 'degC'.
@ -93,7 +116,7 @@ def get_ctt(wrfin, timeidx=0, method="cat",
method, squeeze, cache, meta=False, method, squeeze, cache, meta=False,
_key=_key) _key=_key)
except KeyError: except KeyError:
qice = n.zeros(qv.shape, qv.dtype) qice = np.zeros(qv.shape, qv.dtype)
haveqci = 0 haveqci = 0
else: else:
qice = icevars["QICE"] * 1000.0 #g/kg qice = icevars["QICE"] * 1000.0 #g/kg
@ -116,6 +139,9 @@ def get_ctt(wrfin, timeidx=0, method="cat",
geopt_unstag = destagger(geopt, -3) geopt_unstag = destagger(geopt, -3)
ght = geopt_unstag / Constants.G ght = geopt_unstag / Constants.G
ctt = _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci) _fill_nocloud = 1 if fill_nocloud else 0
ctt = _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, _fill_nocloud,
missing, opt_thresh)
return ctt return ma.masked_values(ctt, missing)

4
src/wrf/metadecorators.py

@ -2034,8 +2034,8 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"):
outcoords = {} outcoords = {}
# Left-most is always cape_cin or cape_cin_lcl_lfc # Left-most is always cape_cin or cape_cin_lcl_lfc
if is2d: if is2d:
outdims[0] = "cape_cin_lcl_lfc" outdims[0] = "mcape_mcin_lcl_lfc"
outcoords["cape_cin_lcl_lfc"] = ["cape", "cin", "lcl", "lfc"] outcoords["mcape_mcin_lcl_lfc"] = ["mcape", "mcin", "lcl", "lfc"]
else: else:
outdims[0] = "cape_cin" outdims[0] = "cape_cin"
outcoords["cape_cin"] = ["cape", "cin"] outcoords["cape_cin"] = ["cape", "cin"]

2
src/wrf/routines.py

@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"],
"wspd_wdir10" : ["units"], "wspd_wdir10" : ["units"],
"uvmet_wspd_wdir" : ["units"], "uvmet_wspd_wdir" : ["units"],
"uvmet10_wspd_wdir" : ["units"], "uvmet10_wspd_wdir" : ["units"],
"ctt" : ["fill_nocloud", "missing", "units"], "ctt" : ["fill_nocloud", "missing", "opt_thresh", "units"],
"cloudfrac" : ["vert_type", "low_thresh", "cloudfrac" : ["vert_type", "low_thresh",
"mid_thresh", "high_thresh"], "mid_thresh", "high_thresh"],
"geopt_stag" : [], "geopt_stag" : [],

2
src/wrf/version.py

@ -1,2 +1,2 @@
__version__ = "1.1.2" __version__ = "1.1.3"

2
test/ci_tests/make_test_file.py

@ -12,7 +12,7 @@ from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round,
VARS_TO_KEEP = ("Times", "XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", VARS_TO_KEEP = ("Times", "XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U",
"XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", "XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2",
"T2", "PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD", "T2", "PSFC", "U10", "V10", "XTIME", "QVAPOR", "QCLOUD",
"QGRAUP", "QRAIN", "QSNOW", "MAPFAC_M", "MAPFAC_U", "QGRAUP", "QRAIN", "QSNOW", "QICE", "MAPFAC_M", "MAPFAC_U",
"MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC") "MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC")
DIMS_TO_TRIM = ("west_east", "south_north", "bottom_top", "bottom_top_stag", DIMS_TO_TRIM = ("west_east", "south_north", "bottom_top", "bottom_top_stag",

4
test/comp_utest.py

@ -609,14 +609,14 @@ if __name__ == "__main__":
#varnames = ["helicity"] #varnames = ["helicity"]
varnames=["avo", "pvo", "eth", "dbz", "helicity", "updraft_helicity", varnames=["avo", "pvo", "eth", "dbz", "helicity", "updraft_helicity",
"omg", "pw", "rh", "slp", "td", "tk", "tv", "twb", "uvmet", "omg", "pw", "rh", "slp", "td", "tk", "tv", "twb", "uvmet",
"cloudfrac"] "cloudfrac", "ctt"]
omp_set_num_threads(omp_get_num_procs()-1) omp_set_num_threads(omp_get_num_procs()-1)
omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_schedule(OMP_SCHED_STATIC, 0)
omp_set_dynamic(False) omp_set_dynamic(False)
# Turn this one off when not needed, since it's slow # Turn this one off when not needed, since it's slow
#varnames += ["cape_2d", "cape_3d"] varnames += ["cape_2d", "cape_3d"]
for varname in varnames: for varname in varnames:
for i,wrfnc in enumerate((NCFILE, NCGROUP)): for i,wrfnc in enumerate((NCFILE, NCGROUP)):

214
test/ipynb/Test_CTT.ipynb

@ -0,0 +1,214 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from netCDF4 import Dataset\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import get_cmap\n",
"import cartopy.crs as crs\n",
"from cartopy.feature import NaturalEarthFeature\n",
"\n",
"from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n",
"\n",
"# Open the NetCDF file\n",
"ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n",
"\n",
"# Get the sea level pressure\n",
"ctt = getvar(ncfile, \"ctt\", fill_nocloud=False)\n",
"slp = getvar(ncfile, \"slp\")\n",
"\n",
"\n",
"# Get the latitude and longitude points\n",
"lats, lons = latlon_coords(ctt)\n",
"\n",
"# Get the cartopy mapping object\n",
"cart_proj = get_cartopy(ctt)\n",
"\n",
"# Create a figure\n",
"fig = plt.figure(figsize=(12,9))\n",
"# Set the GeoAxes to the projection used by WRF\n",
"ax = plt.axes(projection=cart_proj)\n",
"\n",
"# Download and add the states and coastlines\n",
"states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
" name='admin_1_states_provinces_shp')\n",
"ax.add_feature(states, linewidth=.5)\n",
"ax.coastlines('50m', linewidth=0.8)\n",
"\n",
"# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
"levels = np.arange(-80, 20, 5)\n",
"plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n",
" transform=crs.PlateCarree())\n",
"plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n",
" cmap=get_cmap(\"Greys\"))\n",
"\n",
"# Add a color bar\n",
"plt.colorbar(ax=ax, shrink=.88)\n",
"\n",
"# Set the map limits. Not really necessary, but used for demonstration.\n",
"ax.set_xlim(cartopy_xlim(ctt))\n",
"ax.set_ylim(cartopy_ylim(ctt))\n",
"\n",
"# Add the gridlines\n",
"ax.gridlines(color=\"black\", linestyle=\"dotted\")\n",
"\n",
"plt.title(\"Cloud Top Temperature (degC)\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from netCDF4 import Dataset\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import get_cmap\n",
"import cartopy.crs as crs\n",
"from cartopy.feature import NaturalEarthFeature\n",
"\n",
"from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n",
"\n",
"# Open the NetCDF file\n",
"ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n",
"\n",
"# Get the sea level pressure\n",
"ctt = getvar(ncfile, \"ctt\", fill_nocloud=True, opt_thresh=1.0)\n",
"slp = getvar(ncfile, \"slp\")\n",
"\n",
"\n",
"# Get the latitude and longitude points\n",
"lats, lons = latlon_coords(ctt)\n",
"\n",
"# Get the cartopy mapping object\n",
"cart_proj = get_cartopy(ctt)\n",
"\n",
"# Create a figure\n",
"fig = plt.figure(figsize=(12,9))\n",
"# Set the GeoAxes to the projection used by WRF\n",
"ax = plt.axes(projection=cart_proj)\n",
"\n",
"# Download and add the states and coastlines\n",
"states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
" name='admin_1_states_provinces_shp')\n",
"ax.add_feature(states, linewidth=.5)\n",
"ax.coastlines('50m', linewidth=0.8)\n",
"\n",
"# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
"levels = np.arange(-80, 20, 5)\n",
"plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n",
" transform=crs.PlateCarree())\n",
"plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n",
" cmap=get_cmap(\"Greys\"))\n",
"\n",
"# Add a color bar\n",
"plt.colorbar(ax=ax, shrink=.88)\n",
"\n",
"# Set the map limits. Not really necessary, but used for demonstration.\n",
"ax.set_xlim(cartopy_xlim(ctt))\n",
"ax.set_ylim(cartopy_ylim(ctt))\n",
"\n",
"# Add the gridlines\n",
"ax.gridlines(color=\"black\", linestyle=\"dotted\")\n",
"\n",
"plt.title(\"Cloud Top Temperature (degC)\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from netCDF4 import Dataset\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import get_cmap\n",
"import cartopy.crs as crs\n",
"from cartopy.feature import NaturalEarthFeature\n",
"\n",
"from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n",
"\n",
"# Open the NetCDF file\n",
"ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n",
"\n",
"# Get the sea level pressure\n",
"cfrac = getvar(ncfile, \"cfrac\")[2, :]\n",
"slp = getvar(ncfile, \"slp\")\n",
"\n",
"\n",
"# Get the latitude and longitude points\n",
"lats, lons = latlon_coords(ctt)\n",
"\n",
"# Get the cartopy mapping object\n",
"cart_proj = get_cartopy(ctt)\n",
"\n",
"# Create a figure\n",
"fig = plt.figure(figsize=(12,9))\n",
"# Set the GeoAxes to the projection used by WRF\n",
"ax = plt.axes(projection=cart_proj)\n",
"\n",
"# Download and add the states and coastlines\n",
"states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
" name='admin_1_states_provinces_shp')\n",
"ax.add_feature(states, linewidth=.5)\n",
"ax.coastlines('50m', linewidth=0.8)\n",
"\n",
"# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
"levels = np.arange(0, 1.1, .1)\n",
"plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n",
" transform=crs.PlateCarree())\n",
"plt.contourf(to_np(lons), to_np(lats), to_np(cfrac), levels=levels, transform=crs.PlateCarree(),\n",
" cmap=get_cmap(\"Greys_r\"))\n",
"\n",
"# Add a color bar\n",
"plt.colorbar(ax=ax, shrink=.88)\n",
"\n",
"# Set the map limits. Not really necessary, but used for demonstration.\n",
"ax.set_xlim(cartopy_xlim(ctt))\n",
"ax.set_ylim(cartopy_ylim(ctt))\n",
"\n",
"# Add the gridlines\n",
"ax.gridlines(color=\"black\", linestyle=\"dotted\")\n",
"\n",
"plt.title(\"Cloud Fraction\")\n",
"\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}

214
test/ipynb/Test_CTT_Prod.ipynb

@ -0,0 +1,214 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from netCDF4 import Dataset\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import get_cmap\n",
"import cartopy.crs as crs\n",
"from cartopy.feature import NaturalEarthFeature\n",
"\n",
"from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n",
"\n",
"# Open the NetCDF file\n",
"ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n",
"\n",
"# Get the sea level pressure\n",
"ctt = getvar(ncfile, \"ctt\")\n",
"slp = getvar(ncfile, \"slp\")\n",
"\n",
"\n",
"# Get the latitude and longitude points\n",
"lats, lons = latlon_coords(ctt)\n",
"\n",
"# Get the cartopy mapping object\n",
"cart_proj = get_cartopy(ctt)\n",
"\n",
"# Create a figure\n",
"fig = plt.figure(figsize=(12,9))\n",
"# Set the GeoAxes to the projection used by WRF\n",
"ax = plt.axes(projection=cart_proj)\n",
"\n",
"# Download and add the states and coastlines\n",
"states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
" name='admin_1_states_provinces_shp')\n",
"ax.add_feature(states, linewidth=.5)\n",
"ax.coastlines('50m', linewidth=0.8)\n",
"\n",
"# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
"levels = np.arange(-80, 20, 5)\n",
"plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n",
" transform=crs.PlateCarree())\n",
"plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n",
" cmap=get_cmap(\"Greys\"))\n",
"\n",
"# Add a color bar\n",
"plt.colorbar(ax=ax, shrink=.88)\n",
"\n",
"# Set the map limits. Not really necessary, but used for demonstration.\n",
"ax.set_xlim(cartopy_xlim(ctt))\n",
"ax.set_ylim(cartopy_ylim(ctt))\n",
"\n",
"# Add the gridlines\n",
"ax.gridlines(color=\"black\", linestyle=\"dotted\")\n",
"\n",
"plt.title(\"Cloud Top Temperature (degC)\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from netCDF4 import Dataset\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import get_cmap\n",
"import cartopy.crs as crs\n",
"from cartopy.feature import NaturalEarthFeature\n",
"\n",
"from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n",
"\n",
"# Open the NetCDF file\n",
"ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n",
"\n",
"# Get the sea level pressure\n",
"ctt = getvar(ncfile, \"ctt\")\n",
"slp = getvar(ncfile, \"slp\")\n",
"\n",
"\n",
"# Get the latitude and longitude points\n",
"lats, lons = latlon_coords(ctt)\n",
"\n",
"# Get the cartopy mapping object\n",
"cart_proj = get_cartopy(ctt)\n",
"\n",
"# Create a figure\n",
"fig = plt.figure(figsize=(12,9))\n",
"# Set the GeoAxes to the projection used by WRF\n",
"ax = plt.axes(projection=cart_proj)\n",
"\n",
"# Download and add the states and coastlines\n",
"states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
" name='admin_1_states_provinces_shp')\n",
"ax.add_feature(states, linewidth=.5)\n",
"ax.coastlines('50m', linewidth=0.8)\n",
"\n",
"# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
"levels = np.arange(-80, 20, 5)\n",
"plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n",
" transform=crs.PlateCarree())\n",
"plt.contourf(to_np(lons), to_np(lats), to_np(ctt), levels=levels, transform=crs.PlateCarree(),\n",
" cmap=get_cmap(\"Greys\"))\n",
"\n",
"# Add a color bar\n",
"plt.colorbar(ax=ax, shrink=.88)\n",
"\n",
"# Set the map limits. Not really necessary, but used for demonstration.\n",
"ax.set_xlim(cartopy_xlim(ctt))\n",
"ax.set_ylim(cartopy_ylim(ctt))\n",
"\n",
"# Add the gridlines\n",
"ax.gridlines(color=\"black\", linestyle=\"dotted\")\n",
"\n",
"plt.title(\"Cloud Top Temperature (degC)\")\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from netCDF4 import Dataset\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib.cm import get_cmap\n",
"import cartopy.crs as crs\n",
"from cartopy.feature import NaturalEarthFeature\n",
"\n",
"from wrf import to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords\n",
"\n",
"# Open the NetCDF file\n",
"ncfile = Dataset(\"/Users/ladwig/Documents/wrf_files/wrf_vortex_multi/wrfout_d01_2005-08-28_12:00:00\")\n",
"\n",
"# Get the sea level pressure\n",
"cfrac = getvar(ncfile, \"cfrac\")[2, :]\n",
"slp = getvar(ncfile, \"slp\")\n",
"\n",
"\n",
"# Get the latitude and longitude points\n",
"lats, lons = latlon_coords(ctt)\n",
"\n",
"# Get the cartopy mapping object\n",
"cart_proj = get_cartopy(ctt)\n",
"\n",
"# Create a figure\n",
"fig = plt.figure(figsize=(12,9))\n",
"# Set the GeoAxes to the projection used by WRF\n",
"ax = plt.axes(projection=cart_proj)\n",
"\n",
"# Download and add the states and coastlines\n",
"states = NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',\n",
" name='admin_1_states_provinces_shp')\n",
"ax.add_feature(states, linewidth=.5)\n",
"ax.coastlines('50m', linewidth=0.8)\n",
"\n",
"# Make the contour outlines and filled contours for the smoothed sea level pressure.\n",
"levels = np.arange(0, 1.1, .1)\n",
"plt.contour(to_np(lons), to_np(lats), to_np(slp), 10, colors=\"black\",\n",
" transform=crs.PlateCarree())\n",
"plt.contourf(to_np(lons), to_np(lats), to_np(cfrac), levels=levels, transform=crs.PlateCarree(),\n",
" cmap=get_cmap(\"Greys_r\"))\n",
"\n",
"# Add a color bar\n",
"plt.colorbar(ax=ax, shrink=.88)\n",
"\n",
"# Set the map limits. Not really necessary, but used for demonstration.\n",
"ax.set_xlim(cartopy_xlim(ctt))\n",
"ax.set_ylim(cartopy_ylim(ctt))\n",
"\n",
"# Add the gridlines\n",
"ax.gridlines(color=\"black\", linestyle=\"dotted\")\n",
"\n",
"plt.title(\"Cloud Fraction\")\n",
"\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Loading…
Cancel
Save