From 6564b3c0e3badb098bf7ebc45151c0f1ecb1a7ba Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 9 Mar 2018 11:40:33 -0700 Subject: [PATCH 1/5] Updated pyf file for new signutures --- fortran/wrffortran.pyf | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index 7211ff6..c4d9183 100644 --- a/fortran/wrffortran.pyf +++ b/fortran/wrffortran.pyf @@ -231,7 +231,7 @@ python module _wrffortran ! in integer, optional,intent(in),check(shape(prs,0)==mkzh),depend(prs) :: mkzh=shape(prs,0) integer intent(in) :: ter_follow end subroutine dpfcalc - subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 + subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,prsf,prs_new,tmk_new,qvp_new,ght_new,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,ezero,thtecon2 real(kind=8) dimension(mix,mjy,mkzh),intent(in) :: prs @@ -242,6 +242,11 @@ python module _wrffortran ! in real(kind=8) dimension(mix,mjy),intent(in),depend(mix,mjy) :: sfp real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cape real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cin + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prsf + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prs_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: tmk_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: qvp_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: ght_new real(kind=8) intent(in) :: cmsg integer, optional,intent(in),check(shape(prs,0)==mix),depend(prs) :: mix=shape(prs,0) integer, optional,intent(in),check(shape(prs,1)==mjy),depend(prs) :: mjy=shape(prs,1) @@ -251,7 +256,7 @@ python module _wrffortran ! in integer intent(inout) :: errstat character*(*) intent(inout) :: errmsg end subroutine dcapecalc3d - subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 + subroutine dcapecalc2d(prs,tmk,qvp,ght,ter,sfp,cape,cin,prsf,prs_new,tmk_new,qvp_new,ght_new,cmsg,mix,mjy,mkzh,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,ezero,thtecon2 real(kind=8) dimension(mix,mjy,mkzh),intent(in) :: prs @@ -262,6 +267,11 @@ python module _wrffortran ! in real(kind=8) dimension(mix,mjy),intent(in),depend(mix,mjy) :: sfp real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cape real(kind=8) dimension(mix,mjy,mkzh),intent(out,in),depend(mix,mjy,mkzh) :: cin + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prsf + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: prs_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: tmk_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: qvp_new + real(kind=8) dimension(mkzh,mix,mjy),intent(inout),depend(mkzh,mix,mjy) :: ght_new real(kind=8) intent(in) :: cmsg integer, optional,intent(in),check(shape(prs,0)==mix),depend(prs) :: mix=shape(prs,0) integer, optional,intent(in),check(shape(prs,1)==mjy),depend(prs) :: mjy=shape(prs,1) @@ -349,7 +359,7 @@ python module _wrffortran ! in real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285714285714 real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25864379523 end module wrf_constants - subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 + subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,pf,haveqci,nz,ns,ew) ! in :_wrffortran:wrf_fctt.f90 threadsafe use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel real(kind=8) dimension(ew,ns,nz),intent(in) :: prs @@ -360,6 +370,7 @@ python module _wrffortran ! in real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: ght real(kind=8) dimension(ew,ns),intent(in),depend(ew,ns) :: ter 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 integer intent(in) :: haveqci 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) @@ -730,7 +741,7 @@ python module _wrffortran ! in integer intent(inout) :: errstat real(kind=8) :: wrf_intrp_value end function wrf_intrp_value - subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90 + subroutine wrf_vintrp(datain,dataout,pres,tk,qvp,ght,terrain,sfp,smsfp,vcarray,interp_levels,numlevels,icase,ew,ns,nz,extrap,vcor,logp,tempout,rmsg,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90 threadsafe use wrf_constants, only: tlclc2,tlclc3,sclht,tlclc4,thtecon3,eps,tlclc1,celkel,exponi,gammamd,ussalr,expon,thtecon1,algerr,gamma,thtecon2 real(kind=8) dimension(ew,ns,nz),intent(in) :: datain @@ -752,6 +763,7 @@ python module _wrffortran ! in integer intent(in) :: extrap integer intent(in) :: vcor integer intent(in) :: logp + real(kind=8) dimension(ew,ns),intent(inout),depend(ew,ns) :: tempout real(kind=8) intent(in) :: rmsg integer intent(inout) :: errstat character*(*) intent(inout) :: errmsg From 71fe678e916ab632ad8fb2b99ddc19a34bf7aa9d Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 9 Mar 2018 11:42:19 -0700 Subject: [PATCH 2/5] Moved work arrays out of fortran. Fixes #47 --- fortran/rip_cape.f90 | 28 ++++++++++++++++--------- fortran/wrf_fctt.f90 | 6 ++++-- fortran/wrf_vinterp.f90 | 6 ++++-- src/wrf/extension.py | 45 +++++++++++++++++++++++++++++------------ 4 files changed, 58 insertions(+), 27 deletions(-) diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index 5e9cdac..1fc3b23 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -272,6 +272,7 @@ END SUBROUTINE DPFCALC ! NCLFORTSTART SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + prsf, prs_new, tmk_new, qvp_new, ght_new,& cmsg,mix,mjy,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, & @@ -293,6 +294,14 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin + + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new + + REAL(KIND=8), INTENT(IN) :: cmsg CHARACTER(LEN=*), INTENT(IN) :: psafile INTEGER, INTENT(INOUT) :: errstat @@ -308,15 +317,10 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat REAL(KIND=8) :: benamin, dz REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new ! To remove compiler warnings tmkpari = 0 @@ -597,6 +601,7 @@ END SUBROUTINE DCAPECALC3D ! NCLFORTSTART SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& + prsf, prs_new, tmk_new, qvp_new, ght_new,& cmsg,mix,mjy,mkzh,ter_follow,& psafile, errstat, errmsg) USE wrf_constants, ONLY : CELKEL, G, EZERO, ESLCON1, ESLCON2, & @@ -618,6 +623,13 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin + + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prsf + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: prs_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: tmk_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: qvp_new + REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(INOUT) :: ght_new + REAL(KIND=8), INTENT(IN) :: cmsg CHARACTER(LEN=*), INTENT(IN) :: psafile INTEGER, INTENT(INOUT) :: errstat @@ -635,15 +647,11 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat REAL(KIND=8) :: benamin, dz, pup, pdn REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs REAL(KIND=8), DIMENSION(150,150) :: psaditmk LOGICAL :: elfound REAL(KIND=8), DIMENSION(mkzh) :: eth_temp - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new - REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new + ! To remove compiler warnings errstat = 0 diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index 2ef7ecb..0be18ae 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -1,5 +1,5 @@ !NCLFORTSTART -SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew) +SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns, ew) USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL IMPLICIT NONE @@ -12,6 +12,8 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: ter REAL(KIND=8), DIMENSION(ew,ns), INTENT(OUT) :: ctt + REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: pf + !NCLEND ! REAL(KIND=8) :: znfac(nz) @@ -20,7 +22,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew INTEGER i,j,k,ripk REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix REAL(KIND=8) :: arg2, agl_hgt, vt - REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf + REAL(KIND=8) :: p1, p2 !$OMP PARALLEL diff --git a/fortran/wrf_vinterp.f90 b/fortran/wrf_vinterp.f90 index c34efd7..d81270b 100644 --- a/fortran/wrf_vinterp.f90 +++ b/fortran/wrf_vinterp.f90 @@ -127,7 +127,7 @@ END FUNCTION wrf_intrp_value !NCLFORTSTART SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& sfp, smsfp, vcarray, interp_levels, numlevels,& - icase, ew, ns, nz, extrap, vcor, logp, rmsg,& + icase, ew, ns, nz, extrap, vcor, logp, tempout, rmsg,& errstat, errmsg) USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, & @@ -146,6 +146,9 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels + + REAL(KIND=8), DIMENSION(ew,ns), INTENT(INOUT) :: tempout + REAL(KIND=8), INTENT(IN) :: rmsg INTEGER, INTENT(INOUT) :: errstat CHARACTER(LEN=*), INTENT(INOUT) :: errmsg @@ -156,7 +159,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& INTEGER :: i, j, k, kupper !itriv INTEGER :: ifound, isign !miy,mjx INTEGER :: log_errcnt, interp_errcnt, interp_errstat - REAL(KIND=8), DIMENSION(ew,ns) :: tempout REAL(KIND=8) :: rlevel, vlev, diff REAL(KIND=8) :: tmpvlev REAL(KIND=8) :: vcp1, vcp0, valp0, valp1 diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 167d27c..5bcccfa 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -618,20 +618,33 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, else: cape_routine = dcapecalc2d + # Work arrays + k_left_shape = (p_hpa.shape[2], p_hpa.shape[0], p_hpa.shape[1]) + prsf = np.empty(k_left_shape, np.float64, order="F") + prs_new = np.empty(k_left_shape, np.float64, order="F") + tmk_new = np.empty(k_left_shape, np.float64, order="F") + qvp_new = np.empty(k_left_shape, np.float64, order="F") + ght_new = np.empty(k_left_shape, np.float64, order="F") + # note that p_hpa, tk, qv, and ht have the vertical flipped result = cape_routine(p_hpa, - tk, - qv, - ht, - ter, - sfp, - capeview, - cinview, - missing, - ter_follow, - psafile, - errstat, - errmsg) + tk, + qv, + ht, + ter, + sfp, + capeview, + cinview, + prsf, + prs_new, + tmk_new, + qvp_new, + ght_new, + missing, + ter_follow, + psafile, + errstat, + errmsg) if int(errstat) != 0: raise DiagnosticError("".join(npbytes_to_str(errmsg)).strip()) @@ -766,6 +779,8 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): if outview is None: outview = np.empty_like(ter) + pf = np.empty(p_hpa.shape[0:3], np.float64, order="F") + result = wrfcttcalc(p_hpa, tk, qice, @@ -774,6 +789,7 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): ght, ter, outview, + pf, haveqci) return result @@ -858,7 +874,9 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, if outview is None: outdims = field.shape[0:2] + interp_levels.shape outview = np.empty(outdims, field.dtype, order="F") - + + tempout = np.zeros(field.shape[0:2], np.float64, order="F") + errstat = np.array(0) errmsg = np.zeros(Constants.ERRLEN, "c") @@ -877,6 +895,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, extrap, vcor, logp, + tempout, missing, errstat, errmsg) From e6b4f7382b4cb485aab9e8182444e78af543418e Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Wed, 14 Mar 2018 15:54:20 -0600 Subject: [PATCH 3/5] Added checks for _key=None in extract_vars. Fix typo with updraft helicity cache key. Fixes #51. --- src/wrf/g_helicity.py | 2 +- src/wrf/util.py | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/wrf/g_helicity.py b/src/wrf/g_helicity.py index 02ae780..cb3beb5 100755 --- a/src/wrf/g_helicity.py +++ b/src/wrf/g_helicity.py @@ -176,7 +176,7 @@ def get_uh(wrfin, timeidx=0, method="cat", squeeze=True, """ ncvars = extract_vars(wrfin, timeidx, ("W", "PH", "PHB", "MAPFAC_M"), - method, squeeze, cache, meta=False, _key=None) + method, squeeze, cache, meta=False, _key=_key) wstag = ncvars["W"] ph = ncvars["PH"] diff --git a/src/wrf/util.py b/src/wrf/util.py index 3aa9583..7bf9927 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -973,10 +973,12 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): is_moving = is_moving_domain(wrfdict, varname, _key=_key) + _cache_key = _key[first_key] if _key is not None else None + first_array = _extract_var(wrfdict[first_key], varname, timeidx, is_moving=is_moving, method=method, squeeze=False, cache=None, meta=meta, - _key=_key[first_key]) + _key=_cache_key) # Create the output data numpy array based on the first array outdims = [numkeys] @@ -992,10 +994,11 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): break else: keynames.append(key) + _cache_key = _key[key] if _key is not None else None vardata = _extract_var(wrfdict[key], varname, timeidx, is_moving=is_moving, method=method, squeeze=False, cache=None, meta=meta, - _key=_key[key]) + _key=_cache_key) if outdata.shape[1:] != vardata.shape: raise ValueError("data sequences must have the " From 196e28760045e5763f034eeb7b5196f9006083c9 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 16 Mar 2018 10:05:18 -0600 Subject: [PATCH 4/5] Added th alias for theta. Fixes #52. --- src/wrf/routines.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/wrf/routines.py b/src/wrf/routines.py index 5442956..f3665fe 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], "wspd_wdir10" : ["units"], "uvmet_wspd_wdir" : ["units"], "uvmet10_wspd_wdir" : ["units"], - "ctt" : [], + "ctt" : ["fill_nocloud", "missing", "units"], "cloudfrac" : ["vert_type", "low_thresh", "mid_thresh", "high_thresh"], "geopt_stag" : [], @@ -138,7 +138,8 @@ _ALIASES = {"cape_2d" : "cape2d", "td2" : "dp2m", "cfrac" : "cloudfrac", "wspd_wdir_uvmet" : "uvmet_wspd_wdir", - "wspd_wdir_uvmet10" : "uvmet10_wspd_wdir" + "wspd_wdir_uvmet10" : "uvmet10_wspd_wdir", + "th" : "theta" } class ArgumentError(Exception): From cc86d68a35efb5de440426e595c4286c875d1108 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 16 Mar 2018 16:34:44 -0600 Subject: [PATCH 5/5] 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. --- doc/source/_templates/product_table.txt | 6 +- doc/source/new.rst | 20 +++ fortran/wrf_fctt.f90 | 54 +++--- fortran/wrffortran.pyf | 5 +- src/wrf/computation.py | 34 +++- src/wrf/extension.py | 10 +- src/wrf/g_ctt.py | 38 ++++- src/wrf/metadecorators.py | 4 +- src/wrf/routines.py | 2 +- src/wrf/version.py | 2 +- test/ci_tests/make_test_file.py | 2 +- test/comp_utest.py | 4 +- test/ipynb/Test_CTT.ipynb | 214 ++++++++++++++++++++++++ test/ipynb/Test_CTT_Prod.ipynb | 214 ++++++++++++++++++++++++ 14 files changed, 565 insertions(+), 44 deletions(-) create mode 100644 test/ipynb/Test_CTT.ipynb create mode 100644 test/ipynb/Test_CTT_Prod.ipynb diff --git a/doc/source/_templates/product_table.txt b/doc/source/_templates/product_table.txt index c49c3ac..b2ee337 100644 --- a/doc/source/_templates/product_table.txt +++ b/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 | +--------------------+---------------------------------------------------------------+-----------------------------+---------------------------------------------------------------------------------------------------------------------------------------------------------+ -| 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 | | +| | | | **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'. | | | | | | diff --git a/doc/source/new.rst b/doc/source/new.rst index 2962534..c07e5a5 100644 --- a/doc/source/new.rst +++ b/doc/source/new.rst @@ -4,6 +4,26 @@ What's New 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 ^^^^^^^^^^^^^^ diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index 0be18ae..2e61632 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -1,5 +1,6 @@ !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 IMPLICIT NONE @@ -7,13 +8,15 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns !f2py threadsafe !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), INTENT(IN) :: ter 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 + !NCLEND ! 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 opdepthd = 0.D0 k = 0 - prsctt = 0 + prsctt = -1 ! Integrate downward from model top, calculating path at full ! model vertical levels. - DO k=1, nz + DO k=2,nz opdepthu = opdepthd ripk = nz - k + 1 @@ -76,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, pf, haveqci, nz, ns END IF 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 - opdepthd = opdepthu + ABSCOEFI*qcw(i,j,k) * dp/G + opdepthd = opdepthu + ABSCOEFI*qcw(i,j,ripk) * dp/G ELSE - opdepthd = opdepthu + ABSCOEF*qcw(i,j,k) * dp/G + opdepthd = opdepthu + ABSCOEF*qcw(i,j,ripk) * dp/G END IF ELSE opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G END IF - IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN + IF (opdepthd .LT. opt_thresh .AND. k .LT. nz) THEN CYCLE - ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN - prsctt = prs(i,j,1) + ELSE IF (opdepthd .LT. opt_thresh .AND. k .EQ. nz) THEN + IF (fill_nocloud .EQ. 0) THEN + prsctt = prs(i,j,1) + ENDIF EXIT ELSE 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 DO - DO k=2,nz - ripk = nz - k + 1 - p1 = prs(i,j,ripk+1) - p2 = prs(i,j,ripk) - IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN - fac = (prsctt - p1)/(p2 - p1) - arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL - ctt(i,j) = tk(i,j,ripk+1) + arg1 - EXIT - END IF - END DO + ! prsctt should only be 0 if fill values are used + IF (prsctt .GT. -1) THEN + DO k=2,nz + ripk = nz - k + 1 + p1 = prs(i,j,ripk+1) + p2 = prs(i,j,ripk) + IF (prsctt .GE. p1 .AND. prsctt .LE. p2) THEN + fac = (prsctt - p1)/(p2 - p1) + arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL + ctt(i,j) = tk(i,j,ripk+1) + arg1 + EXIT + END IF + END DO + ELSE + ctt(i,j) = missing + END IF END DO END DO !$OMP END DO diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index c4d9183..be27452 100644 --- a/fortran/wrffortran.pyf +++ b/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(expon,rd,ussalr,g) :: exponi=5.25864379523 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 use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel 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,nz),intent(inout),depend(ew,ns,nz) :: pf 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,1)==ns),depend(prs) :: ns=shape(prs,1) integer, optional,intent(in),check(shape(prs,0)==ew),depend(prs) :: ew=shape(prs,0) diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 3973709..164d905 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -1014,6 +1014,11 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, high_thresh (:obj:`float`): The bottom vertical threshold for what is 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 :class:`numpy.ndarray` instead of @@ -1047,8 +1052,9 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, @set_alg_metadata(2, "pres_hpa", refvarndims=3, description="cloud top temperature") @convert_units("temp", "c") -def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, - units="degC"): +def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, + fill_nocloud=False, missing=default_fill(np.float64), + opt_thresh=1.0, meta=True, units="degC"): """Return the cloud top temperature. This is the raw computational algorithm and does not extract any variables @@ -1094,6 +1100,23 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, qice (:class:`xarray.DataArray` or :class:`numpy.ndarray`, optional): Ice mixing ratio in [kg/kg] with the same dimensionality as *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 :class:`numpy.ndarray` instead of @@ -1128,8 +1151,13 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, haveqci = 0 else: haveqci = 1 if qice.any() else 0 + + _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 _ctt(pres_hpa, tkel, qice, qcld, qv, height, terrain, haveqci) + return ma.masked_values(ctt, missing) diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 5bcccfa..42c4734 100755 --- a/src/wrf/extension.py +++ b/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)) -@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)) @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. Located in wrf_fctt.f90. @@ -790,7 +791,10 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): ter, outview, pf, - haveqci) + haveqci, + fill_nocloud, + missing, + opt_thresh) return result diff --git a/src/wrf/g_ctt.py b/src/wrf/g_ctt.py index c3dff09..1c59c01 100644 --- a/src/wrf/g_ctt.py +++ b/src/wrf/g_ctt.py @@ -1,11 +1,12 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) -import numpy as n +import numpy as np +import numpy.ma as ma #from .extension import computectt, computetk from .extension import _ctt, _tk -from .constants import Constants, ConversionFactors +from .constants import Constants, ConversionFactors, default_fill from .destag import destagger from .decorators import convert_units from .metadecorators import copy_and_set_metadata @@ -19,7 +20,8 @@ from .util import extract_vars @convert_units("temp", "c") def get_ctt(wrfin, timeidx=0, method="cat", 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. 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 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` product table for a list of available units for 'ctt'. Default is 'degC'. @@ -93,7 +116,7 @@ def get_ctt(wrfin, timeidx=0, method="cat", method, squeeze, cache, meta=False, _key=_key) except KeyError: - qice = n.zeros(qv.shape, qv.dtype) + qice = np.zeros(qv.shape, qv.dtype) haveqci = 0 else: qice = icevars["QICE"] * 1000.0 #g/kg @@ -116,6 +139,9 @@ def get_ctt(wrfin, timeidx=0, method="cat", geopt_unstag = destagger(geopt, -3) 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) diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 27425c1..b35d57d 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -2034,8 +2034,8 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): outcoords = {} # Left-most is always cape_cin or cape_cin_lcl_lfc if is2d: - outdims[0] = "cape_cin_lcl_lfc" - outcoords["cape_cin_lcl_lfc"] = ["cape", "cin", "lcl", "lfc"] + outdims[0] = "mcape_mcin_lcl_lfc" + outcoords["mcape_mcin_lcl_lfc"] = ["mcape", "mcin", "lcl", "lfc"] else: outdims[0] = "cape_cin" outcoords["cape_cin"] = ["cape", "cin"] diff --git a/src/wrf/routines.py b/src/wrf/routines.py index f3665fe..df6d533 100644 --- a/src/wrf/routines.py +++ b/src/wrf/routines.py @@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], "wspd_wdir10" : ["units"], "uvmet_wspd_wdir" : ["units"], "uvmet10_wspd_wdir" : ["units"], - "ctt" : ["fill_nocloud", "missing", "units"], + "ctt" : ["fill_nocloud", "missing", "opt_thresh", "units"], "cloudfrac" : ["vert_type", "low_thresh", "mid_thresh", "high_thresh"], "geopt_stag" : [], diff --git a/src/wrf/version.py b/src/wrf/version.py index e39ded1..18b2dba 100644 --- a/src/wrf/version.py +++ b/src/wrf/version.py @@ -1,2 +1,2 @@ -__version__ = "1.1.2" +__version__ = "1.1.3" diff --git a/test/ci_tests/make_test_file.py b/test/ci_tests/make_test_file.py index e6ec29b..fcdfac3 100644 --- a/test/ci_tests/make_test_file.py +++ b/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", "XLONG_V", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", "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") DIMS_TO_TRIM = ("west_east", "south_north", "bottom_top", "bottom_top_stag", diff --git a/test/comp_utest.py b/test/comp_utest.py index 55448a8..958cd75 100644 --- a/test/comp_utest.py +++ b/test/comp_utest.py @@ -609,14 +609,14 @@ if __name__ == "__main__": #varnames = ["helicity"] varnames=["avo", "pvo", "eth", "dbz", "helicity", "updraft_helicity", "omg", "pw", "rh", "slp", "td", "tk", "tv", "twb", "uvmet", - "cloudfrac"] + "cloudfrac", "ctt"] omp_set_num_threads(omp_get_num_procs()-1) omp_set_schedule(OMP_SCHED_STATIC, 0) omp_set_dynamic(False) # 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 i,wrfnc in enumerate((NCFILE, NCGROUP)): diff --git a/test/ipynb/Test_CTT.ipynb b/test/ipynb/Test_CTT.ipynb new file mode 100644 index 0000000..533ea13 --- /dev/null +++ b/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 +} diff --git a/test/ipynb/Test_CTT_Prod.ipynb b/test/ipynb/Test_CTT_Prod.ipynb new file mode 100644 index 0000000..1bf8707 --- /dev/null +++ b/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 +}