Browse Source

Merge branch 'feature/dyn_temp' into develop

lon0
Bill Ladwig 7 years ago
parent
commit
2495b51944
  1. 6
      doc/source/_templates/product_table.txt
  2. 20
      doc/source/new.rst
  3. 28
      fortran/rip_cape.f90
  4. 56
      fortran/wrf_fctt.f90
  5. 6
      fortran/wrf_vinterp.f90
  6. 23
      fortran/wrffortran.pyf
  7. 34
      src/wrf/computation.py
  8. 53
      src/wrf/extension.py
  9. 38
      src/wrf/g_ctt.py
  10. 2
      src/wrf/g_helicity.py
  11. 4
      src/wrf/metadecorators.py
  12. 5
      src/wrf/routines.py
  13. 7
      src/wrf/util.py
  14. 2
      src/wrf/version.py
  15. 2
      test/ci_tests/make_test_file.py
  16. 4
      test/comp_utest.py
  17. 214
      test/ipynb/Test_CTT.ipynb
  18. 214
      test/ipynb/Test_CTT_Prod.ipynb

6
doc/source/_templates/product_table.txt

@ -13,11 +13,13 @@ @@ -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'. |
| | | | |

20
doc/source/new.rst

@ -4,6 +4,26 @@ What's New @@ -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
^^^^^^^^^^^^^^

28
fortran/rip_cape.f90

@ -272,6 +272,7 @@ END SUBROUTINE DPFCALC @@ -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,& @@ -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,& @@ -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 @@ -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,& @@ -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,& @@ -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

56
fortran/wrf_fctt.f90

@ -1,5 +1,6 @@ @@ -1,5 +1,6 @@
!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,&
fill_nocloud, missing, opt_thresh, nz, ns, ew)
USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL
IMPLICIT NONE
@ -7,10 +8,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -7,10 +8,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
!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
@ -20,7 +25,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -20,7 +25,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
@ -58,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -58,12 +63,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
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
@ -74,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -74,21 +79,23 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
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)
@ -98,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew @@ -98,17 +105,22 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
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

6
fortran/wrf_vinterp.f90

@ -127,7 +127,7 @@ END FUNCTION wrf_intrp_value @@ -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,& @@ -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,& @@ -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

23
fortran/wrffortran.pyf

@ -231,7 +231,7 @@ python module _wrffortran ! in @@ -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 @@ -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 @@ -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 @@ -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 @@ -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,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
@ -360,7 +370,11 @@ python module _wrffortran ! in @@ -360,7 +370,11 @@ 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 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)
@ -730,7 +744,7 @@ python module _wrffortran ! in @@ -730,7 +744,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 +766,7 @@ python module _wrffortran ! in @@ -752,6 +766,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

34
src/wrf/computation.py

@ -1015,6 +1015,11 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, @@ -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
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
:class:`xarray.DataArray`. Default is True.
@ -1047,8 +1052,9 @@ def cloudfrac(vert, relh, vert_inc_w_height, low_thresh, mid_thresh, @@ -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
@ -1095,6 +1101,23 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, @@ -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
*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
:class:`xarray.DataArray`. Default is True.
@ -1129,7 +1152,12 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True, @@ -1129,7 +1152,12 @@ def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True,
else:
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)

53
src/wrf/extension.py

@ -618,20 +618,33 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, @@ -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())
@ -754,10 +767,11 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1, @@ -754,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.
@ -766,6 +780,8 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): @@ -766,6 +780,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,7 +790,11 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): @@ -774,7 +790,11 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None):
ght,
ter,
outview,
haveqci)
pf,
haveqci,
fill_nocloud,
missing,
opt_thresh)
return result
@ -859,6 +879,8 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, @@ -859,6 +879,8 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp,
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 +899,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, @@ -877,6 +899,7 @@ def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp,
extrap,
vcor,
logp,
tempout,
missing,
errstat,
errmsg)

38
src/wrf/g_ctt.py

@ -1,11 +1,12 @@ @@ -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 @@ -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", @@ -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", @@ -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", @@ -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)

2
src/wrf/g_helicity.py

@ -176,7 +176,7 @@ def get_uh(wrfin, timeidx=0, method="cat", squeeze=True, @@ -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"]

4
src/wrf/metadecorators.py

@ -2034,8 +2034,8 @@ def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): @@ -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"]

5
src/wrf/routines.py

@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"], @@ -112,7 +112,7 @@ _VALID_KARGS = {"cape2d" : ["missing"],
"wspd_wdir10" : ["units"],
"uvmet_wspd_wdir" : ["units"],
"uvmet10_wspd_wdir" : ["units"],
"ctt" : [],
"ctt" : ["fill_nocloud", "missing", "opt_thresh", "units"],
"cloudfrac" : ["vert_type", "low_thresh",
"mid_thresh", "high_thresh"],
"geopt_stag" : [],
@ -138,7 +138,8 @@ _ALIASES = {"cape_2d" : "cape2d", @@ -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):

7
src/wrf/util.py

@ -973,10 +973,12 @@ def _combine_dict(wrfdict, varname, timeidx, method, meta, _key): @@ -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): @@ -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 "

2
src/wrf/version.py

@ -1,2 +1,2 @@ @@ -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, @@ -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",

4
test/comp_utest.py

@ -609,14 +609,14 @@ if __name__ == "__main__": @@ -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)):

214
test/ipynb/Test_CTT.ipynb

@ -0,0 +1,214 @@ @@ -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 @@ @@ -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