From a07297636a3aed242ecabe92144488d506a01b64 Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Fri, 19 Aug 2016 15:58:23 -0600 Subject: [PATCH] Added computational unit tests. Most, if not all, computational unit tests passed. Removed obsolute routines. More code cleanup for fortran. Changed constants module to wrf_constants.f90. Fixed bugs. --- fortran/calc_uh.f90 | 50 +- fortran/cloud_fracf.f90 | 6 +- fortran/eqthecalc.f90 | 2 +- fortran/rip_cape.f90 | 12 +- fortran/{constants.f90 => wrf_constants.f90} | 6 +- fortran/wrf_fctt.f90 | 38 +- fortran/wrf_pvo.f90 | 24 +- fortran/wrf_pw.f90 | 4 +- fortran/wrf_relhl.f90 | 60 +- fortran/wrf_rip_phys_routines.f90 | 6 +- fortran/wrf_testfunc.f90 | 32 + fortran/wrf_user.f90 | 126 ++-- fortran/wrf_user_dbz.f90 | 49 +- fortran/wrf_user_latlon_routines.f90 | 6 +- fortran/wrf_vinterp.f90 | 8 +- fortran/wrffortran.pyf | 161 ++--- setup.py | 3 +- src/wrf/api.py | 12 +- src/wrf/computation.py | 123 ++-- src/wrf/constants.py | 11 +- src/wrf/dbz.py | 10 +- src/wrf/decorators.py | 188 +----- src/wrf/destag.py | 4 +- src/wrf/extension.py | 675 +++---------------- src/wrf/helicity.py | 14 +- src/wrf/interp.py | 9 +- src/wrf/interputils.py | 10 +- src/wrf/latlon.py | 29 +- src/wrf/latlonutils.py | 16 +- src/wrf/metadecorators.py | 89 +-- src/wrf/projection.py | 11 +- src/wrf/pw.py | 5 +- src/wrf/util.py | 9 +- src/wrf/uvmet.py | 13 +- test/comp_utest.py | 580 ++++++++++++++++ test/ipynb/WRF_python_demo.ipynb | 22 +- test/utests.py | 64 +- 37 files changed, 1287 insertions(+), 1200 deletions(-) rename fortran/{constants.f90 => wrf_constants.f90} (96%) create mode 100644 fortran/wrf_testfunc.f90 create mode 100644 test/comp_utest.py diff --git a/fortran/calc_uh.f90 b/fortran/calc_uh.f90 index 262d954..a893189 100644 --- a/fortran/calc_uh.f90 +++ b/fortran/calc_uh.f90 @@ -1,3 +1,27 @@ +!################################################################## +!################################################################## +!###### ###### +!###### SUBROUTINE CALC_UH ###### +!###### ###### +!###### Developed by ###### +!###### Center for Analysis and Prediction of Storms ###### +!###### University of Oklahoma ###### +!###### ###### +!################################################################## +!################################################################## +! +! Calculates updraft helicity (UH) to detect rotating updrafts. +! Formula follows Kain et al, 2008, Wea. and Forecasting, 931-952, +! but this version has controls for the limits of integration +! uhminhgt to uhmxhgt, in m AGL. Kain et al used 2000 to 5000 m. +! Units of UH are m^2/s^2. +! +! Note here that us and vs are at ARPS scalar points. +! w is at w-point (scalar pt in horiz, staggered vertical) +! +! Keith Brewster, CAPS/Univ. of Oklahoma +! March, 2010 + !NCLFORTSTART SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & vs, w, uh, tem1, tem2) @@ -41,8 +65,8 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & DO j=2,ny-1 DO i=2,nx-1 wavg = 0.5*(w(i,j,k)+w(i,j,k+1)) - tem1(i,j,k) = wavg * ((vs(i+1,j,k) - vs(i-1,j,k))/(twodx * mapfct(i,j)) - & - (us(i,j+1,k) - us(i,j-1,k))/(twody * mapfct(i,j))) + tem1(i,j,k) = wavg*((vs(i+1,j,k) - vs(i-1,j,k))/(twodx*mapfct(i,j)) - & + (us(i,j+1,k) - us(i,j-1,k))/(twody*mapfct(i,j))) tem2(i,j,k) = 0.5*(zp(i,j,k) + zp(i,j,k+1)) END DO END DO @@ -65,7 +89,7 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & END DO kbot = k wgtlw = (zp(i,j,kbot) - zbot)/(zp(i,j,kbot) - zp(i,j,kbot-1)) - wbot = (wgtlw*w(i,j,kbot-1)) + ((1.-wgtlw)*w(i,j,kbot)) + wbot = (wgtlw*w(i,j,kbot-1)) + ((1. - wgtlw)*w(i,j,kbot)) ! Find w at uhmxhgt AGL (top) DO k=2,nz-3 @@ -73,18 +97,18 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & END DO ktop = k wgtlw = (zp(i,j,ktop) - ztop)/(zp(i,j,ktop) - zp(i,j,ktop-1)) - wtop = (wgtlw*w(i,j,ktop-1)) + ((1.-wgtlw)*w(i,j,ktop)) + wtop = (wgtlw*w(i,j,ktop-1)) + ((1. - wgtlw)*w(i,j,ktop)) ! First part, uhmnhgt to kbot - wsum = 0.5*(w(i,j,kbot) + wbot) * (zp(i,j,kbot) - zbot) + wsum = 0.5*(w(i,j,kbot) + wbot)*(zp(i,j,kbot) - zbot) ! Integrate up through column DO k=(kbot+1),(ktop-1) - wsum = wsum + 0.5*(w(i,j,k) + w(i,j,k-1)) * (zp(i,j,k) - zp(i,j,k-1)) + wsum = wsum + 0.5*(w(i,j,k) + w(i,j,k-1))*(zp(i,j,k) - zp(i,j,k-1)) END DO ! Last part, ktop-1 to uhmxhgt - wsum = wsum + 0.5*(wtop + w(i,j,ktop-1)) * (ztop - zp(i,j,ktop-1)) + wsum = wsum + 0.5*(wtop + w(i,j,ktop-1))*(ztop - zp(i,j,ktop-1)) wmean = wsum/(uhmxhgt - uhmnhgt) IF (wmean > 0.) THEN ! column updraft, not downdraft @@ -95,7 +119,7 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & END DO kbot = k wgtlw = (tem2(i,j,kbot) - zbot)/(tem2(i,j,kbot) - tem2(i,j,kbot-1)) - helbot = (wgtlw*tem1(i,j,kbot-1)) + ((1.-wgtlw)*tem1(i,j,kbot)) + helbot = (wgtlw*tem1(i,j,kbot-1)) + ((1. - wgtlw)*tem1(i,j,kbot)) ! Find helicity at uhmxhgt AGL (top) DO k=2,nz-3 @@ -103,23 +127,23 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & END DO ktop = k wgtlw = (tem2(i,j,ktop) - ztop)/(tem2(i,j,ktop) - tem2(i,j,ktop-1)) - heltop = (wgtlw*tem1(i,j,ktop-1)) + ((1.-wgtlw)*tem1(i,j,ktop)) + heltop = (wgtlw*tem1(i,j,ktop-1)) + ((1. - wgtlw)*tem1(i,j,ktop)) ! First part, uhmnhgt to kbot - sum = 0.5*(tem1(i,j,kbot) + helbot) * (tem2(i,j,kbot) - zbot) + sum = 0.5*(tem1(i,j,kbot) + helbot)*(tem2(i,j,kbot) - zbot) ! Integrate up through column DO k=(kbot+1),(ktop-1) - sum = sum + 0.5*(tem1(i,j,k) + tem1(i,j,k-1)) * (tem2(i,j,k) - tem2(i,j,k-1)) + sum = sum + 0.5*(tem1(i,j,k) + tem1(i,j,k-1))*(tem2(i,j,k) - tem2(i,j,k-1)) END DO ! Last part, ktop-1 to uhmxhgt - uh(i,j) = sum + 0.5*(heltop + tem1(i,j,ktop-1)) * (ztop - tem2(i,j,ktop-1)) + uh(i,j) = sum + 0.5*(heltop + tem1(i,j,ktop-1))*(ztop - tem2(i,j,ktop-1)) END IF END DO END DO - uh = uh * 1000. ! Scale according to Kain et al. (2008) + uh = uh*1000. ! Scale according to Kain et al. (2008) RETURN diff --git a/fortran/cloud_fracf.f90 b/fortran/cloud_fracf.f90 index 99a8c63..f3b23ea 100644 --- a/fortran/cloud_fracf.f90 +++ b/fortran/cloud_fracf.f90 @@ -37,9 +37,9 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew) END IF END DO - lowc(i,j) = 4.0 * lowc(i,j)/100. - 3.0 - midc(i,j) = 4.0 * midc(i,j)/100. - 3.0 - highc(i,j) = 2.5 * highc(i,j)/100. - 1.5 + lowc(i,j) = 4.0*lowc(i,j)/100. - 3.0 + midc(i,j) = 4.0*midc(i,j)/100. - 3.0 + highc(i,j) = 2.5*highc(i,j)/100. - 1.5 lowc(i,j) = amin1(lowc(i,j), 1.0) lowc(i,j) = amax1(lowc(i,j), 0.0) diff --git a/fortran/eqthecalc.f90 b/fortran/eqthecalc.f90 index 9bbd03c..a4cb6ba 100644 --- a/fortran/eqthecalc.f90 +++ b/fortran/eqthecalc.f90 @@ -1,7 +1,7 @@ ! Theta-e !NCLFORTSTART SUBROUTINE DEQTHECALC(qvp, tmk, prs, eth, miy, mjx, mkzh) - USE constants, ONLY : EPS, GAMMA, GAMMAMD, TLCLC1, TLCLC2, TLCLC3, TLCLC4, & + USE wrf_constants, ONLY : EPS, GAMMA, GAMMAMD, TLCLC1, TLCLC2, TLCLC3, TLCLC4, & THTECON1, THTECON2, THTECON3 IMPLICIT NONE diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index b97adb9..0778fe4 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -18,7 +18,7 @@ ! NCLFORTSTART REAL(KIND=8) FUNCTION TVIRTUAL(temp, ratmix) - USE constants, ONLY : EPS + USE wrf_constants, ONLY : EPS !f2py threadsafe @@ -34,9 +34,9 @@ REAL(KIND=8) FUNCTION TVIRTUAL(temp, ratmix) END FUNCTION TVIRTUAL ! NCLFORTSTART -REAL(KIND=8) FUNCTION TONPSADIABAT(thte,prs,psadithte,psadiprs,psaditmk,gamma,& +REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gamma,& errstat, errmsg) - USE constants, ONLY : ALGERR + USE wrf_constants, ONLY : ALGERR !f2py threadsafe !f2py intent(in,out) :: cape, cin @@ -127,7 +127,7 @@ END FUNCTION TONPSADIABAT !NCLFORTSTART SUBROUTINE DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, fname, errstat, errmsg) - USE constants, ONLY : ALGERR + USE wrf_constants, ONLY : ALGERR !f2py threadsafe @@ -250,7 +250,7 @@ END SUBROUTINE DPFCALC SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& cmsg,miy,mjx,mkzh,i3dflag,ter_follow,& psafile, errstat, errmsg) - USE constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & + USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 @@ -259,7 +259,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& !f2py threadsafe !f2py intent(in,out) :: cape, cin - INTEGER, INTENT(IN) :: miy,mjx,mkzh,i3dflag,ter_follow + INTEGER, INTENT(IN) :: miy, mjx, mkzh, i3dflag, ter_follow REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: prs REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: tmk REAL(KIND=8), DIMENSION(miy,mjx,mkzh), INTENT(IN) :: qvp diff --git a/fortran/constants.f90 b/fortran/wrf_constants.f90 similarity index 96% rename from fortran/constants.f90 rename to fortran/wrf_constants.f90 index 81594f1..617e63f 100644 --- a/fortran/constants.f90 +++ b/fortran/wrf_constants.f90 @@ -1,11 +1,11 @@ ! These are chosen to match the wrf module_model_constants.F where ! applicable -MODULE constants +MODULE wrf_constants INTEGER :: ERRLEN=512 INTEGER :: ALGERR=64 REAL(KIND=8), PARAMETER :: WRF_EARTH_RADIUS = 6370000.D0 - REAL(KIND=8), PARAMETER :: T_BASE = 300.0 + REAL(KIND=8), PARAMETER :: T_BASE = 300.0D0 REAL(KIND=8), PARAMETER :: PI = 3.1415926535897932384626433D0 REAL(KIND=8), PARAMETER :: RAD_PER_DEG = PI/180.D0 REAL(KIND=8), PARAMETER :: DEG_PER_RAD = 180.D0/PI @@ -62,5 +62,5 @@ MODULE constants REAL(KIND=8), PARAMETER :: EXPONI = 1./EXPON -END MODULE constants +END MODULE wrf_constants diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index ebcb07f..97da1b9 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -1,6 +1,6 @@ !NCLFORTSTART -SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz) - USE constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL +SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew) + USE wrf_constants, ONLY : EPS, USSALR, RD, G, ABSCOEFI, ABSCOEF, CELKEL IMPLICIT NONE @@ -19,9 +19,9 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz ! LOCAL VARIABLES INTEGER i,j,k,ripk !INTEGER :: mjx,miy,mkzh - REAL(KIND=8) :: vt,opdepthu,opdepthd,dp - REAL(KIND=8) :: ratmix,arg1,arg2,agl_hgt - REAL(KIND=8) :: fac,prsctt + REAL(KIND=8) :: vt, opdepthu, opdepthd, dp + REAL(KIND=8) :: ratmix, arg1, arg2, agl_hgt + REAL(KIND=8) :: fac, prsctt !REAL(KIND=8) :: eps,ussalr,rgas,grav,abscoefi,abscoef,celkel,wrfout !REAL(KIND=8) :: ght(ew,ns,nz),stuff(ew,ns) !REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf(ns,ew,nz),p1,p2 @@ -37,13 +37,13 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz ! Calculate the surface pressure DO j=1,ns DO i=1,ew - ratmix = .001d0*qvp(i,j,1) + ratmix = .001D0*qvp(i,j,1) arg1 = EPS + ratmix - arg2 = EPS * (1. + ratmix) - vt = tk(i,j,1) * arg1/arg2 !Virtual temperature + arg2 = EPS*(1. + ratmix) + vt = tk(i,j,1)*arg1/arg2 !Virtual temperature agl_hgt = ght(i,j,nz) - ter(i,j) - arg1 = -G / (RD * USSALR) - pf(i,j,nz) = prs(i,j,1) * (vt / (vt + USSALR*(agl_hgt)))**(arg1) + arg1 = -G/(RD*USSALR) + pf(i,j,nz) = prs(i,j,1)*(vt/(vt + USSALR*(agl_hgt)))**(arg1) END DO END DO @@ -51,14 +51,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz DO j=1,ns DO i=1,ew ripk = nz-k+1 - pf(i,j,k) = .5d0 * (prs(i,j,ripk) + prs(i,j,ripk-1)) + pf(i,j,k) = .5D0*(prs(i,j,ripk) + prs(i,j,ripk-1)) END DO END DO END DO DO j=1,ns DO i=1,ew - opdepthd = 0.d0 + opdepthd = 0.D0 k = 0 ! Integrate downward from model top, calculating path at full @@ -72,9 +72,9 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz ripk = nz - k + 1 IF (k .EQ. 1) THEN - dp = 200.d0 * (pf(i,j,1) - prs(i,j,nz)) ! should be in Pa + dp = 200.D0*(pf(i,j,1) - prs(i,j,nz)) ! should be in Pa ELSE - dp = 100.d0 * (pf(i,j,k) - pf(i,j,k-1)) ! should be in Pa + dp = 100.D0*(pf(i,j,k) - pf(i,j,k-1)) ! should be in Pa END IF IF (haveqci .EQ. 0) then @@ -85,7 +85,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz opdepthd = opdepthu + ABSCOEF*qcw(i,j,k) * dp/G END IF ELSE - opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk)) * dp/G + opdepthd = opdepthd + (ABSCOEF*qcw(i,j,ripk) + ABSCOEFI*qci(i,j,ripk))*dp/G END IF IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN @@ -96,7 +96,7 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz prsctt = prs(i,j,1) EXIT ELSE - fac = (1. - opdepthu) / (opdepthd - opdepthu) + fac = (1. - opdepthu)/(opdepthd - opdepthu) prsctt = pf(i,j,k-1) + fac*(pf(i,j,k) - pf(i,j,k-1)) prsctt = MIN(prs(i,j,1), MAX(prs(i,j,nz), prsctt)) EXIT @@ -104,12 +104,12 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, ew, ns, nz END DO DO k=2,nz - ripk = nz-k+1 + 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 + 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 !GOTO 40 EXIT diff --git a/fortran/wrf_pvo.f90 b/fortran/wrf_pvo.f90 index a810edc..209a257 100644 --- a/fortran/wrf_pvo.f90 +++ b/fortran/wrf_pvo.f90 @@ -37,9 +37,9 @@ SUBROUTINE DCOMPUTEABSVORT(av, u, v, msfu, msfv, msft, cor, dx, dy, nx, ny, nz,& dsy = (jp1 - jm1) * dy mm = msft(i,j)*msft(i,j) ! PRINT *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1) - dudy = 0.5D0 * (u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - & + dudy = 0.5D0*(u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - & u(i,jm1,k)/msfu(i,jm1) - u(i+1,jm1,k)/msfu(i+1,jm1))/dsy*mm - dvdx = 0.5D0 * (v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - & + dvdx = 0.5D0*(v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - & v(im1,j,k)/msfv(im1,j) - v(im1,j+1,k)/msfv(im1,j+1))/dsx*mm avort = dvdx - dudy + cor(i,j) av(i,j,k) = avort*1.D5 @@ -55,7 +55,7 @@ END SUBROUTINE DCOMPUTEABSVORT !NCLFORTSTART SUBROUTINE DCOMPUTEPV(pv, u, v, theta, prs, msfu, msfv, msft, cor, dx, dy, nx, & ny, nz, nxp1, nyp1) - USE constants, ONLY : G + USE wrf_constants, ONLY : G IMPLICIT NONE @@ -92,22 +92,22 @@ SUBROUTINE DCOMPUTEPV(pv, u, v, theta, prs, msfu, msfv, msft, cor, dx, dy, nx, & ip1 = MIN(i+1, nx) im1 = MAX(i-1, 1) ! PRINT *,jp1,jm1,ip1,im1 - dsx = (ip1 - im1) * dx - dsy = (jp1 - jm1) * dy + dsx = (ip1 - im1)*dx + dsy = (jp1 - jm1)*dy mm = msft(i,j)*msft(i,j) ! PRINT *,j,i,u(i,jp1,k),msfu(i,jp1),u(i,jp1,k)/msfu(i,jp1) - dudy = 0.5D0 * (u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - & + dudy = 0.5D0*(u(i,jp1,k)/msfu(i,jp1) + u(i+1,jp1,k)/msfu(i+1,jp1) - & u(i,jm1,k)/msfu(i,jm1) - u(i+1,jm1,k)/msfu(i+1,jm1))/dsy*mm - dvdx = 0.5D0 * (v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - & + dvdx = 0.5D0*(v(ip1,j,k)/msfv(ip1,j) + v(ip1,j+1,k)/msfv(ip1,j+1) - & v(im1,j,k)/msfv(im1,j) - v(im1,j+1,k)/msfv(im1,j+1))/dsx*mm avort = dvdx - dudy + cor(i,j) dp = prs(i,j,kp1) - prs(i,j,km1) - dudp = 0.5D0 * (u(i,j,kp1) + u(i+1,j,kp1) - u(i,j,km1) - u(i+1,j,km1))/dp - dvdp = 0.5D0 * (v(i,j,kp1) + v(i,j+1,kp1) - v(i,j,km1) - v(i,J+1,km1))/dp + dudp = 0.5D0*(u(i,j,kp1) + u(i+1,j,kp1) - u(i,j,km1) - u(i+1,j,km1))/dp + dvdp = 0.5D0*(v(i,j,kp1) + v(i,j+1,kp1) - v(i,j,km1) - v(i,J+1,km1))/dp dthdp = (theta(i,j,kp1) - theta(i,j,km1))/dp - dthdx = (theta(ip1,j,k) - theta(im1,j,k))/dsx * msft(i,j) - dthdy = (theta(i,jp1,k) - theta(i,jm1,k))/dsy * msft(i,j) - pv(i,j,k) = -G * (dthdp*avort - dvdp*dthdx + dudp*dthdy)*10000.D0 + dthdx = (theta(ip1,j,k) - theta(im1,j,k))/dsx*msft(i,j) + dthdy = (theta(i,jp1,k) - theta(i,jm1,k))/dsy*msft(i,j) + pv(i,j,k) = -G*(dthdp*avort - dvdp*dthdx + dudp*dthdy)*10000.D0 ! if(i.eq.300 .and. j.eq.300) then ! PRINT*,'avort,dudp,dvdp,dthdp,dthdx,dthdy,pv' ! PRINT*,avort,dudp,dvdp,dthdp,dthdx,dthdy,pv(i,j,k) diff --git a/fortran/wrf_pw.f90 b/fortran/wrf_pw.f90 index 6a7e8d3..ba19ffb 100644 --- a/fortran/wrf_pw.f90 +++ b/fortran/wrf_pw.f90 @@ -1,6 +1,6 @@ !NCLFORTSTART SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh) - USE constants, ONLY : RD + USE wrf_constants, ONLY : RD IMPLICIT NONE @@ -21,7 +21,7 @@ SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh) DO k=1,nz DO j=1,ny DO i=1,nx - pw(i,j) = pw(i,j) + ((p(i,j,k)/(RD*tv(i,j,k))) * qv(i,j,k) * (ht(i,j,k+1) - ht(i,j,k))) + pw(i,j) = pw(i,j) + ((p(i,j,k)/(RD*tv(i,j,k)))*qv(i,j,k)*(ht(i,j,k+1) - ht(i,j,k))) END DO END DO END DO diff --git a/fortran/wrf_relhl.f90 b/fortran/wrf_relhl.f90 index b72acab..9efeec8 100644 --- a/fortran/wrf_relhl.f90 +++ b/fortran/wrf_relhl.f90 @@ -1,6 +1,38 @@ +! *************************************************************** +! * Storm Relative Helicity (SRH) is a measure of the * +! * streamwise vorticity within the inflow environment of a * +! * convective storm. It is calculated by multiplying the * +! * storm-relative inflow velocity vector (Vh-C) by the * +! * streamwise vorticity (Zh) and integrating this quantity * +! * over the inflow depth (lowest 1-3 km layers above ground * +! * level). It describes the extent to which corkscrew-like * +! * motion occurs (similar to the spiraling motion of an * +! * American football). SRH corresponds to the transfer of * +! * vorticity from the environment to an air parcel in * +! * convective motion and is used to predict the potential * +! * for tornadic development (cyclonic updraft rotation) in * +! * right-moving supercells. * +! * * +! * There is no clear threshold value for SRH when forecasting * +! * supercells, since the formation of supercells appears to be * +! * related more strongly to the deeper layer vertical shear. * +! * Larger values of 0-3-km SRH (greater than 250 m**2/s**2) * +! * and 0-1-km SRH (greater than 100 m**2/s**2), suggest an * +! * increased threat of tornadoes with supercells. For SRH, * +! * larger values are generally better, but there are no clear * +! * "boundaries" between non-tornadic and significant tornadic * +! * supercells. * +! * * +! * SRH < 100 (lowest 1 km): cutoff value * +! * SRH = 150-299: supercells possible with weak tornadoes * +! * SRH = 300-499: very favorable to supercell development and * +! * strong tornadoes * +! * SRH > 450 : violent tornadoes * +! *************************************************************** + ! NCLFORTSTART SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) - USE constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD + USE wrf_constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD IMPLICIT NONE @@ -24,11 +56,13 @@ SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) INTEGER :: i, j, k, k10, k3, ktop !REAL(KIND=8), PARAMETER :: DTR=PI/180.d0, DPR=180.d0/PI - DO j = 1, mjx-1 - DO i = 1, miy-1 - sdh = 0.d0 - su = 0.d0 - sv = 0.d0 + !DO j = 1, mjx-1 + DO j=1, mjx + DO i=1, miy + !DO i=1, miy-1 + sdh = 0.D0 + su = 0.D0 + sv = 0.D0 k3 = 0 k10 = 0 ktop = 0 @@ -54,25 +88,25 @@ SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) su = su + 0.5D0*dh*(u(i,j,k-1) + u(i,j,k)) sv = sv + 0.5D0*dh*(v(i,j,k-1) + v(i,j,k)) END DO - ua = su / sdh - va = sv / sdh + ua = su/sdh + va = sv/sdh asp = SQRT(ua*ua + va*va) IF (ua .EQ. 0.D0 .AND. va .EQ. 0.D0) THEN adr = 0.D0 ELSE adr = DEG_PER_RAD * (PI + ATAN2(ua,va)) ENDIF - bsp = 0.75D0 * asp + bsp = 0.75D0*asp bdr = adr + 30.D0 IF (bdr .GT. 360.D0) THEN bdr = bdr - 360.D0 ENDIF - cu = -bsp * SIN(bdr * RAD_PER_DEG) - cv = -bsp * COS(bdr * RAD_PER_DEG) + cu = -bsp*SIN(bdr * RAD_PER_DEG) + cv = -bsp*COS(bdr * RAD_PER_DEG) sum = 0.D0 DO k = mkzh-1, ktop, -1 - x = ((u(i,j,k) - cu) * (v(i,j,k) - v(i,j,k+1))) - & - ((v(i,j,k) - cv) * (u(i,j,k) - u(i,j,k+1))) + x = ((u(i,j,k) - cu)*(v(i,j,k) - v(i,j,k+1))) - & + ((v(i,j,k) - cv)*(u(i,j,k) - u(i,j,k+1))) sum = sum + x END DO sreh(i,j) = -sum diff --git a/fortran/wrf_rip_phys_routines.f90 b/fortran/wrf_rip_phys_routines.f90 index 3476fbd..ed40d5e 100644 --- a/fortran/wrf_rip_phys_routines.f90 +++ b/fortran/wrf_rip_phys_routines.f90 @@ -30,7 +30,7 @@ ! NCLFORTSTART SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) - USE constants, ONLY : ALGERR, GAMMA, GAMMAMD, TLCLC1, TLCLC2, TLCLC3, & + USE wrf_constants, ONLY : ALGERR, GAMMA, GAMMAMD, TLCLC1, TLCLC2, TLCLC3, & EPS, TLCLC4, THTECON1, THTECON2, THTECON3 IMPLICIT NONE @@ -175,7 +175,7 @@ END SUBROUTINE WETBULBCALC !NCLFORTSTART SUBROUTINE OMGCALC(qvp, tmk, www, prs, omg, mx, my, mz) - USE constants, ONLY : G, RD, EPS + USE wrf_constants, ONLY : G, RD, EPS IMPLICIT NONE @@ -238,7 +238,7 @@ END SUBROUTINE OMGCALC ! ------------------------------------------------------------------ !NCLFORTSTART SUBROUTINE VIRTUAL_TEMP(temp, ratmix, tv, nx, ny, nz) - USE constants, ONLY : EPS + USE wrf_constants, ONLY : EPS IMPLICIT NONE diff --git a/fortran/wrf_testfunc.f90 b/fortran/wrf_testfunc.f90 new file mode 100644 index 0000000..dd3c915 --- /dev/null +++ b/fortran/wrf_testfunc.f90 @@ -0,0 +1,32 @@ +SUBROUTINE testfunc(a, b, c, nx, ny, nz, errstat, errmsg) + USE wrf_constants, ONLY : ALGERR + IMPLICIT NONE + + !threadsafe + !f2py intent(in,out) :: b + + INTEGER, INTENT(IN) :: nx, ny, nz + + REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: a + REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: b + CHARACTER(LEN=*), INTENT(IN) :: c + INTEGER, INTENT(INOUT) :: errstat + REAL(KIND=8), PARAMETER :: blah=123.45 + CHARACTER(LEN=*), INTENT(INOUT) :: errmsg + + INTEGER :: i,j,k + + DO k=1,nz + DO j=1,ny + DO i=1,nx + b(i,j,k) = a(i,j,k) + END DO + END DO + END DO + + errstat = ALGERR + WRITE(errmsg, *) c(1:20), blah + + RETURN + +END SUBROUTINE testfunc diff --git a/fortran/wrf_user.f90 b/fortran/wrf_user.f90 index b3081d7..cee1632 100644 --- a/fortran/wrf_user.f90 +++ b/fortran/wrf_user.f90 @@ -1,6 +1,6 @@ ! NCLFORTSTART SUBROUTINE DCOMPUTEPI(pi, pressure, nx, ny, nz) - USE constants, ONLY : P1000MB, RD, CP + USE wrf_constants, ONLY : P1000MB, RD, CP IMPLICIT NONE @@ -29,7 +29,7 @@ END SUBROUTINE DCOMPUTEPI ! Temperature from potential temperature in kelvin. !NCLFORTSTART SUBROUTINE DCOMPUTETK(tk, pressure, theta, nx) - USE constants, ONLY : P1000MB, RD, CP + USE wrf_constants, ONLY : P1000MB, RD, CP IMPLICIT NONE @@ -65,7 +65,7 @@ SUBROUTINE DINTERP3DZ(data3d, out2d, zdata, desiredloc, nx, ny, nz, missingval) !f2py threadsafe !f2py intent(in,out) :: out2d - INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: nx, ny, nz REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: data3d REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: out2d REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: zdata @@ -99,7 +99,7 @@ SUBROUTINE DINTERP3DZ(data3d, out2d, zdata, desiredloc, nx, ny, nz, missingval) DO WHILE ((.NOT. dointerp) .AND. (kp >= 2)) IF (((zdata(i,j,kp-im) < height) .AND. (zdata(i,j,kp-ip) > height))) THEN - w2 = (height-zdata(i,j,kp-im))/(zdata(i,j,kp-ip)-zdata(i,j,kp-im)) + w2 = (height - zdata(i,j,kp-im))/(zdata(i,j,kp-ip) - zdata(i,j,kp-im)) w1 = 1.D0 - w2 out2d(i,j) = w1*data3d(i,j,kp-im) + w2*data3d(i,j,kp-ip) dointerp = .TRUE. @@ -139,7 +139,7 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain) DO i = 1,nx ii = MIN0(i,nxz) im1 = MAX0(i-1,1) - znew(i,j,k) = 0.5D0 * (z(ii,j,k) + z(im1,j,k)) + znew(i,j,k) = 0.5D0*(z(ii,j,k) + z(im1,j,k)) END DO END DO END DO @@ -150,7 +150,7 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain) jj = MIN0(j,nyz) jm1 = MAX0(j-1,1) DO i = 1,nx - znew(i,j,k) = 0.5D0 * (z(i,jj,k) + z(i,jm1,k)) + znew(i,j,k) = 0.5D0*(z(i,jj,k) + z(i,jm1,k)) END DO END DO END DO @@ -166,7 +166,7 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain) DO k = 2,nz DO j = 1,ny DO i = 1,nx - znew(i,j,k) = znew(i,j,k-1) + 2.D0 * (z(i,j,k-1) - znew(i,j,k-1)) + znew(i,j,k) = znew(i,j,k-1) + 2.D0*(z(i,j,k-1) - znew(i,j,k-1)) END DO END DO END DO @@ -282,7 +282,7 @@ END SUBROUTINE DINTERP1D ! NCLFORTSTART SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & t_sea_level, t_surf, level, errstat, errmsg) - USE constants, ONLY : ALGERR, RD, G, USSALR + USE wrf_constants, ONLY : ALGERR, RD, G, USSALR IMPLICIT NONE @@ -364,8 +364,8 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & DO j = 1,ny DO i = 1,nx - klo = MAX(level(i,j)-1,1) - khi = MIN(klo+1,nz-1) + klo = MAX(level(i,j)-1, 1) + khi = MIN(klo+1, nz-1) IF (klo == khi) THEN errstat = ALGERR @@ -379,8 +379,8 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & plo = p(i,j,klo) phi = p(i,j,khi) - tlo = t(i,j,klo) * (1.D0+0.608D0*q(i,j,klo)) - thi = t(i,j,khi) * (1.D0+0.608D0*q(i,j,khi)) + tlo = t(i,j,klo)*(1.D0 + 0.608D0*q(i,j,klo)) + thi = t(i,j,khi)*(1.D0 + 0.608D0*q(i,j,khi)) ! zlo = zetahalf(klo)/ztop*(ztop-terrain(i,j))+terrain(i,j) ! zhi = zetahalf(khi)/ztop*(ztop-terrain(i,j))+terrain(i,j) zlo = z(i,j,klo) @@ -424,8 +424,8 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & ! Convert to hPa in this step, by multiplying by 0.01. The original ! Fortran routine didn't do this, but the NCL script that called it ! did, so we moved it here. - sea_level_pressure(i,j) = 0.01 * (p(i,j,1)*EXP((2.D0*G*z_half_lowest)/& - (RD*(t_sea_level(i,j)+t_surf(i,j))))) + sea_level_pressure(i,j) = 0.01*(p(i,j,1)*EXP((2.D0*G*z_half_lowest)/& + (RD*(t_sea_level(i,j) + t_surf(i,j))))) END DO END DO @@ -546,7 +546,7 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) b(i,j+1) .EQ. missing) THEN a(i,j) = a(i,j) ELSE - a(i,j) = a(i,j) + COEF*(b(i,j-1)-2*b(i,j) + b(i,j+1)) + a(i,j) = a(i,j) + COEF*(b(i,j-1) - 2*b(i,j) + b(i,j+1)) END IF END DO END DO @@ -557,7 +557,7 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) b(i+1,j) .EQ. missing) THEN a(i,j) = a(i,j) ELSE - a(i,j) = a(i,j) + COEF*(b(i-1,j)-2*b(i,j)+b(i+1,j)) + a(i,j) = a(i,j) + COEF*(b(i-1,j) - 2*b(i,j)+b(i+1,j)) END IF END DO END DO @@ -585,7 +585,7 @@ END SUBROUTINE FILTER2D ! NCLFORTSTART SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx) - USE constants, ONLY : EZERO, ESLCON1, ESLCON2, CELKEL, RD, RV, EPS + USE wrf_constants, ONLY : EZERO, ESLCON1, ESLCON2, CELKEL, RD, RV, EPS IMPLICIT NONE @@ -593,7 +593,7 @@ SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx) !f2py intent(in,out) :: rh INTEGER, INTENT(IN) :: nx - REAL(KIND=8), DIMENSION(nx), INTENT(IN) :: qv,p,t + REAL(KIND=8), DIMENSION(nx), INTENT(IN) :: qv, p, t REAL(KIND=8), DIMENSION(nx), INTENT(OUT) :: rh ! NCLEND @@ -605,12 +605,12 @@ SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx) pressure = p(i) temperature = t(i) ! es = 1000.*svp1* - es = EZERO * EXP(ESLCON1 * (temperature-CELKEL) / (temperature-ESLCON2)) + es = EZERO*EXP(ESLCON1*(temperature - CELKEL)/(temperature - ESLCON2)) ! qvs = ep_2*es/(pressure-es) - qvs = EPS * es / (0.01D0 * pressure - (1.D0 - EPS) * es) + qvs = EPS*es/(0.01D0*pressure - (1.D0 - EPS)*es) ! rh = 100*amax1(1., qv(i)/qvs) ! rh(i) = 100.*qv(i)/qvs - rh(i) = 100.D0 * DMAX1(DMIN1(qv(i) / qvs, 1.0D0), 0.0D0) + rh(i) = 100.D0*DMAX1(DMIN1(qv(i)/qvs, 1.0D0), 0.0D0) END DO RETURN @@ -625,10 +625,10 @@ SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny, !f2py threadsafe !f2py intent(in,out) :: ii,jj - INTEGER, INTENT(IN) :: nx,ny,imsg - INTEGER, INTENT(OUT) :: ii,jj - REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: lat_array,long_array - REAL(KIND=8) :: lat,longitude + INTEGER, INTENT(IN) :: nx, ny, imsg + INTEGER, INTENT(OUT) :: ii, jj + REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: lat_array, long_array + REAL(KIND=8) :: lat, longitude ! NCLEND REAL(KIND=8) :: longd,latd @@ -640,14 +640,14 @@ SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny, ir = imsg jr = imsg - dist_min = 1.d+20 + dist_min = 1.D+20 DO j = 1,ny DO i = 1,nx - latd = (lat_array(i,j)-lat)**2 - longd = (long_array(i,j)-longitude)**2 + latd = (lat_array(i,j) - lat)**2 + longd = (long_array(i,j) - longitude)**2 ! longd = dmin1((long_array(i,j)-longitude)**2, & ! (long_array(i,j)+longitude)**2) - dist = SQRT(latd+longd) + dist = SQRT(latd + longd) IF (dist_min .GT. dist) THEN dist_min = dist ir = DBLE(i) @@ -664,8 +664,8 @@ SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny, ! script which has 0-based indexing. IF (ir .NE. imsg .AND. jr .NE. imsg) THEN - ii = NINT(ir)-1 - jj = NINT(jr)-1 + ii = NINT(ir) - 1 + jj = NINT(jr) - 1 ELSE ii = imsg jj = imsg @@ -691,7 +691,7 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & !f2py threadsafe !f2py intent(in,out) :: uvmet - INTEGER,INTENT(IN) :: nx,ny,nxp1,nyp1,istag + INTEGER,INTENT(IN) :: nx, ny, nxp1, nyp1, istag LOGICAL,INTENT(IN) :: is_msg_val REAL(KIND=8), DIMENSION(nxp1,ny), INTENT(IN):: u REAL(KIND=8), DIMENSION(nx,nyp1), INTENT(IN) :: v @@ -699,15 +699,15 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: flat REAL(KIND=8), DIMENSION(nx,ny), INTENT(INOUT) :: longca REAL(KIND=8), DIMENSION(nx,ny), INTENT(INOUT) :: longcb - REAL(KIND=8), INTENT(IN) :: cen_long,cone,rpd - REAL(KIND=8), INTENT(IN) :: umsg,vmsg,uvmetmsg + REAL(KIND=8), INTENT(IN) :: cen_long, cone, rpd + REAL(KIND=8), INTENT(IN) :: umsg, vmsg, uvmetmsg REAL(KIND=8), DIMENSION(nx,ny,2), INTENT(OUT) :: uvmet ! NCLEND INTEGER :: i,j - REAL(KIND=8) :: uk,vk + REAL(KIND=8) :: uk, vk ! msg stands for missing value in this code ! WRITE (6,FMT=*) ' in compute_uvmet ',NX,NY,NXP1,NYP1,ISTAG @@ -743,8 +743,8 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & uvmet(i,j,1) = uvmetmsg uvmet(i,j,2) = uvmetmsg ELSE - uk = 0.5D0* (u(i,j)+u(i+1,j)) - vk = 0.5D0* (v(i,j)+v(i,j+1)) + uk = 0.5D0*(u(i,j) + u(i+1,j)) + vk = 0.5D0*(v(i,j) + v(i,j+1)) uvmet(i,j,1) = vk*longcb(i,j) + uk*longca(i,j) uvmet(i,j,2) = vk*longca(i,j) - uk*longcb(i,j) END IF @@ -792,13 +792,13 @@ SUBROUTINE DCOMPUTETD(td, pressure, qv_in, nx) INTEGER :: i DO i = 1,nx - qv = DMAX1(qv_in(i),0.D0) + qv = DMAX1(qv_in(i), 0.D0) ! vapor pressure - tdc = qv*pressure(i)/ (.622D0 + qv) + tdc = qv*pressure(i)/(.622D0 + qv) ! avoid problems near zero - tdc = DMAX1(tdc,0.001D0) - td(i) = (243.5D0*LOG(tdc)-440.8D0) / (19.48D0-LOG(tdc)) + tdc = DMAX1(tdc, 0.001D0) + td(i) = (243.5D0*LOG(tdc) - 440.8D0)/(19.48D0 - LOG(tdc)) END DO RETURN @@ -807,6 +807,8 @@ END SUBROUTINE DCOMPUTETD ! NCLFORTSTART SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz) + USE wrf_constants, ONLY : G + IMPLICIT NONE !f2py threadsafe @@ -819,26 +821,26 @@ SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz) ! NCLEND - REAL(KIND=8) :: qclw,dp - REAL(KIND=8), PARAMETER :: GG = 1000.d0/9.8d0 + REAL(KIND=8) :: qclw, dp + REAL(KIND=8), PARAMETER :: GG = 1000.D0/G INTEGER i,j,k DO j = 1,ny DO i = 1,nx - iclw(i,j) = 0.d0 + iclw(i,j) = 0.D0 END DO END DO DO j = 3,ny - 2 DO i = 3,nx - 2 DO k = 1,nz - qclw = DMAX1(qc_in(i,j,k),0.d0) + qclw = DMAX1(qc_in(i,j,k), 0.D0) IF (k.EQ.1) THEN dp = pressure(i,j,k-1) - pressure(i,j,k) ELSE IF (k.EQ.nz) then dp = pressure(i,j,k) - pressure(i,j,k+1) ELSE - dp = (pressure(i,j,k-1) - pressure(i,j,k+1)) / 2.d0 + dp = (pressure(i,j,k-1) - pressure(i,j,k+1))/2.D0 END IF iclw(i,j) = iclw(i,j) + qclw*dp*GG END DO @@ -849,38 +851,6 @@ SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz) END SUBROUTINE DCOMPUTEICLW -SUBROUTINE testfunc(a, b, c, nx, ny, nz, errstat, errmsg) - USE constants, ONLY : ALGERR - IMPLICIT NONE - - !threadsafe - !f2py intent(in,out) :: b - - INTEGER, INTENT(IN) :: nx, ny, nz - - REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(IN) :: a - REAL(KIND=8), DIMENSION(nx,ny,nz), INTENT(OUT) :: b - CHARACTER(LEN=*), INTENT(IN) :: c - INTEGER, INTENT(INOUT) :: errstat - REAL(KIND=8), PARAMETER :: blah=123.45 - CHARACTER(LEN=*), INTENT(INOUT) :: errmsg - - INTEGER :: i,j,k - - DO k=1,nz - DO j=1,ny - DO i=1,nx - b(i,j,k) = a(i,j,k) - END DO - END DO - END DO - - errstat = ALGERR - WRITE(errmsg, *) c(1:20), blah - - RETURN - -END SUBROUTINE testfunc diff --git a/fortran/wrf_user_dbz.f90 b/fortran/wrf_user_dbz.f90 index 93ad9af..0ff3470 100644 --- a/fortran/wrf_user_dbz.f90 +++ b/fortran/wrf_user_dbz.f90 @@ -1,6 +1,37 @@ +! This routine computes equivalent reflectivity factor (in dBZ) at +! each model grid point. In calculating Ze, the RIP algorithm makes +! assumptions consistent with those made in an early version +! (ca. 1996) of the bulk mixed-phase microphysical scheme in the MM5 +! model (i.e., the scheme known as "Resiner-2"). For each species: +! +! 1. Particles are assumed to be spheres of constant density. The +! densities of rain drops, snow particles, and graupel particles are +! taken to be rho_r = rho_l = 1000 kg m^-3, rho_s = 100 kg m^-3, and +! rho_g = 400 kg m^-3, respectively. (l refers to the density of +! liquid water.) +! +! 2. The size distribution (in terms of the actual diameter of the +! particles, rather than the melted diameter or the equivalent solid +! ice sphere diameter) is assumed to follow an exponential +! distribution of the form N(D) = N_0 * exp( lambda*D ). +! +! 3. If ivarint=0, the intercept parameters are assumed constant +! (as in early Reisner-2), with values of 8x10^6, 2x10^7, +! and 4x10^6 m^-4, for rain, snow, and graupel, respectively. +! If ivarint=1, variable intercept parameters are used, as +! calculated in Thompson, Rasmussen, and Manning (2004, Monthly +! Weather Review, Vol. 132, No. 2, pp. 519-542.) +! +! 4. If iliqskin=1, frozen particles that are at a temperature above +! freezing are assumed to scatter as a liquid particle. +! +! More information on the derivation of simulated reflectivity in +! RIP can be found in Stoelinga (2005, unpublished write-up). +! Contact Mark Stoelinga (stoeling@atmos.washington.edu) for a copy. + !NCLFORTSTART SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx, ny, nz) - USE constants, ONLY : GAMMA_SEVEN, RHOWAT, RHO_R, RHO_S, RHO_G, ALPHA, & + USE wrf_constants, ONLY : GAMMA_SEVEN, RHOWAT, RHO_R, RHO_S, RHO_G, ALPHA, & CELKEL, PI, RD IMPLICIT NONE @@ -50,16 +81,16 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx DO k = 1,nz DO j = 1,ny DO i = 1,nx - IF (qvp(i,j,k).LT.0.0) THEN + IF (qvp(i,j,k) .LT. 0.0) THEN qvp(i,j,k) = 0.0 END IF - IF (qra(i,j,k).LT.0.0) THEN + IF (qra(i,j,k) .LT. 0.0) THEN qra(i,j,k) = 0.0 END IF - IF (qsn(i,j,k).LT.0.0) THEN + IF (qsn(i,j,k) .LT. 0.0) THEN qsn(i,j,k) = 0.0 END IF - IF (qgr(i,j,k).LT.0.0) THEN + IF (qgr(i,j,k) .LT. 0.0) THEN qgr(i,j,k) = 0.0 END IF END DO @@ -89,7 +120,7 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx DO j = 1,ny DO i = 1,nx virtual_t = tmk(i,j,k)*(0.622D0 + qvp(i,j,k))/(0.622D0*(1.D0 + qvp(i,j,k))) - rhoair = prs(i,j,k) / (RD*virtual_t) + rhoair = prs(i,j,k)/(RD*virtual_t) ! Adjust factor for brightband, where snow or graupel particle ! scatters like liquid water (alpha=1.0) because it is assumed to @@ -116,7 +147,7 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx ronv = RON2 IF (qra(i,j,k) .GT. R1) THEN - ronv = RON_CONST1R*TANH((RON_QR0-qra(i,j,k))/RON_DELQR0) + RON_CONST2R + ronv = RON_CONST1R*TANH((RON_QR0 - qra(i,j,k))/RON_DELQR0) + RON_CONST2R END IF ELSE @@ -129,8 +160,8 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx ! the sum of z_e for each hydrometeor species: z_e = factor_r*(rhoair*qra(i,j,k))**1.75D0/ronv**.75D0 + & - factorb_s*(rhoair*qsn(i,j,k))**1.75D0/sonv**.75D0 + & - factorb_g* (rhoair*qgr(i,j,k))**1.75D0/gonv**.75D0 + factorb_s*(rhoair*qsn(i,j,k))**1.75D0/sonv**.75D0 + & + factorb_g*(rhoair*qgr(i,j,k))**1.75D0/gonv**.75D0 ! Adjust small values of Z_e so that dBZ is no lower than -30 z_e = MAX(z_e, .001D0) diff --git a/fortran/wrf_user_latlon_routines.f90 b/fortran/wrf_user_latlon_routines.f90 index 86e638f..42eb7fd 100644 --- a/fortran/wrf_user_latlon_routines.f90 +++ b/fortran/wrf_user_latlon_routines.f90 @@ -1,6 +1,6 @@ !NCLFORTSTART SUBROUTINE ROTATECOORDS(ilat, ilon, olat, olon, lat_np, lon_np, lon_0, direction) - USE constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD + USE wrf_constants, ONLY : PI, RAD_PER_DEG, DEG_PER_RAD IMPLICIT NONE @@ -61,7 +61,7 @@ END SUBROUTINE ROTATECOORDS SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& pole_lat, pole_lon, knowni, knownj, dx, dy, latinc,& loninc, lat, lon, loc, errstat, errmsg) - USE constants, ONLY : ALGERR, PI, RAD_PER_DEG, DEG_PER_RAD, WRF_EARTH_RADIUS + USE wrf_constants, ONLY : ALGERR, PI, RAD_PER_DEG, DEG_PER_RAD, WRF_EARTH_RADIUS ! Converts input lat/lon values to the cartesian (i,j) value ! for the given projection. @@ -278,7 +278,7 @@ END SUBROUTINE DLLTOIJ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,& pole_lat, pole_lon, knowni, knownj, dx, dy, latinc,& loninc, ai, aj, loc, errstat, errmsg) - USE constants, ONLY : ALGERR, PI, RAD_PER_DEG, DEG_PER_RAD, WRF_EARTH_RADIUS + USE wrf_constants, ONLY : ALGERR, PI, RAD_PER_DEG, DEG_PER_RAD, WRF_EARTH_RADIUS ! converts input lat/lon values to the cartesian (i,j) value ! for the given projection. diff --git a/fortran/wrf_vinterp.f90 b/fortran/wrf_vinterp.f90 index 6fb2328..995ebe6 100644 --- a/fortran/wrf_vinterp.f90 +++ b/fortran/wrf_vinterp.f90 @@ -1,3 +1,7 @@ +! The subroutines in this file were taken directly from RIP code written +! by Dr. Mark Stoelinga. They were modified by Sherrie +! Fredrick(NCAR/MMM) to work with NCL February 2015. + !NCLFORTSTART SUBROUTINE wrf_monotonic(out, in, lvprs, cor, idir, delta, ew, ns, nz, icorsw) @@ -62,7 +66,7 @@ END SUBROUTINE wrf_monotonic !NCLFORTSTART FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat, errmsg) - USE constants, ONLY : ALGERR, SCLHT + USE wrf_constants, ONLY : ALGERR, SCLHT IMPLICIT NONE @@ -122,7 +126,7 @@ 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) - USE constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & + USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, & CELKEL, EPS, USSALR diff --git a/fortran/wrffortran.pyf b/fortran/wrffortran.pyf index fe3f9ad..831eea3 100644 --- a/fortran/wrffortran.pyf +++ b/fortran/wrffortran.pyf @@ -33,53 +33,9 @@ python module _wrffortran ! in integer, optional,check(shape(pres,1)==ns),depend(pres) :: ns=shape(pres,1) integer, optional,check(shape(pres,0)==ew),depend(pres) :: ew=shape(pres,0) end subroutine dcloudfrac - module constants ! in :_wrffortran:constants.f90 - real(kind=8), parameter,optional :: wrf_earth_radius=6370000.d0 - real(kind=8), parameter,optional :: rhowat=1000.d0 - real(kind=8), parameter,optional :: t_base=300.0 - real(kind=8), parameter,optional :: cp=1004.d0 - real(kind=8), parameter,optional :: thtecon1=3376.d0 - real(kind=8), parameter,optional :: thtecon3=.81d0 - real(kind=8), parameter,optional :: thtecon2=2.54d0 - real(kind=8), parameter,optional :: p1000mb=100000.d0 - real(kind=8), parameter,optional :: rv=461.5d0 - real(kind=8), parameter,optional,depend(pi) :: rad_per_deg=pi/180.d0 - real(kind=8), parameter,optional :: rd=287.04d0 - real(kind=8), parameter,optional :: abscoef=.145d0 - real(kind=8), parameter,optional :: celkel=273.15d0 - real(kind=8), parameter,optional :: abscoefi=.272d0 - real(kind=8), parameter,optional :: eslcon2=29.65d0 - real(kind=8), parameter,optional :: eslcon1=17.67d0 - real(kind=8), parameter,optional :: pi=3.141592653589793d0 - real(kind=8), parameter,optional :: tlclc2=3.5d0 - real(kind=8), parameter,optional :: tlclc3=4.805d0 - real(kind=8), parameter,optional :: rho_g=400.d0 - real(kind=8), parameter,optional :: tlclc1=2840.d0 - real(kind=8), parameter,optional :: tlclc4=55.d0 - real(kind=8), parameter,optional,depend(pi) :: deg_per_rad=180.d0/pi - real(kind=8), parameter,optional :: cpmd=.887d0 - real(kind=8), parameter,optional,depend(rd,g) :: sclht=rd*256.d0/g - real(kind=8), parameter,optional :: ussalr=0.0065d0 - real(kind=8), parameter,optional :: default_fill=9.9692099683868690d36 - real(kind=8), parameter,optional :: rho_s=100.d0 - real(kind=8), parameter,optional,depend(rhowat) :: rho_r=1000.0 - real(kind=8), parameter,optional :: alpha=0.224d0 - real(kind=8), parameter,optional :: celkel_triple=273.16d0 - real(kind=8), parameter,optional :: ezero=6.112d0 - real(kind=8), parameter,optional,depend(rd,ussalr,g) :: expon=0.190189602446 - real(kind=8), parameter,optional :: rgasmd=.608d0 - real(kind=8), parameter,optional :: g=9.81d0 - integer, optional :: errlen=512 - real(kind=8), parameter,optional :: eps=0.622d0 - real(kind=8), parameter,optional :: gamma_seven=720.d0 - real(kind=8), parameter,optional,depend(cpmd,rgasmd) :: gammamd=-0.279 - integer, optional :: algerr=64 - real(kind=8), parameter,optional,depend(cp,rd) :: gamma=0.285896414343 - real(kind=8), parameter,optional,depend(expon,rd,ussalr,g) :: exponi=5.25791098534 - end module constants subroutine deqthecalc(qvp,tmk,prs,eth,miy,mjx,mkzh) ! in :_wrffortran:eqthecalc.f90 threadsafe - use constants, only: tlclc2,tlclc3,tlclc1,tlclc4,eps,thtecon3,gammamd,thtecon1,gamma,thtecon2 + use wrf_constants, only: tlclc2,tlclc3,tlclc1,tlclc4,eps,thtecon3,gammamd,thtecon1,gamma,thtecon2 real(kind=8) dimension(miy,mjx,mkzh),intent(in) :: qvp real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: tmk real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: prs @@ -90,14 +46,14 @@ python module _wrffortran ! in end subroutine deqthecalc function tvirtual(temp,ratmix) ! in :_wrffortran:rip_cape.f90 threadsafe - use constants, only: eps + use wrf_constants, only: eps real(kind=8) intent(in) :: temp real(kind=8) intent(in) :: ratmix real(kind=8) :: tvirtual end function tvirtual function tonpsadiabat(thte,prs,psadithte,psadiprs,psaditmk,gamma,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe - use constants, only: algerr + use wrf_constants, only: algerr real(kind=8) intent(in) :: thte real(kind=8) intent(in) :: prs real(kind=8) dimension(150),intent(in) :: psadithte @@ -110,7 +66,7 @@ python module _wrffortran ! in end function tonpsadiabat subroutine dlookup_table(psadithte,psadiprs,psaditmk,fname,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe - use constants, only: algerr + use wrf_constants, only: algerr real(kind=8) dimension(150),intent(inout) :: psadithte real(kind=8) dimension(150),intent(inout) :: psadiprs real(kind=8) dimension(150,150),intent(inout) :: psaditmk @@ -129,7 +85,7 @@ python module _wrffortran ! in end subroutine dpfcalc subroutine dcapecalc3d(prs,tmk,qvp,ght,ter,sfp,cape,cin,cmsg,miy,mjx,mkzh,i3dflag,ter_follow,psafile,errstat,errmsg) ! in :_wrffortran:rip_cape.f90 threadsafe - use constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,algerr,ezero,thtecon2 + use wrf_constants, only: tlclc2,gamma,tlclc1,rgasmd,tlclc4,g,tlclc3,thtecon3,eps,rd,cpmd,celkel,gammamd,eslcon2,eslcon1,cp,thtecon1,algerr,ezero,thtecon2 real(kind=8) dimension(miy,mjx,mkzh),intent(in) :: prs real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: tmk real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: qvp @@ -148,9 +104,53 @@ python module _wrffortran ! in integer intent(inout) :: errstat character*(*) intent(inout) :: errmsg end subroutine dcapecalc3d - subroutine wrfcttcalc(prs,tk,qci,qcw,qvp,ght,ter,ctt,haveqci,ew,ns,nz) ! in :_wrffortran:wrf_fctt.f90 + module wrf_constants ! in :_wrffortran:wrf_constants.f90 + real(kind=8), parameter,optional :: wrf_earth_radius=6370000.d0 + real(kind=8), parameter,optional :: rhowat=1000.d0 + real(kind=8), parameter,optional :: t_base=300.0d0 + real(kind=8), parameter,optional :: cp=1004.5d0 + real(kind=8), parameter,optional :: thtecon1=3376.d0 + real(kind=8), parameter,optional :: thtecon3=.81d0 + real(kind=8), parameter,optional :: thtecon2=2.54d0 + real(kind=8), parameter,optional :: p1000mb=100000.d0 + real(kind=8), parameter,optional :: rv=461.6d0 + real(kind=8), parameter,optional,depend(pi) :: rad_per_deg=pi/180.d0 + real(kind=8), parameter,optional :: rd=287.d0 + real(kind=8), parameter,optional :: abscoef=.145d0 + real(kind=8), parameter,optional :: celkel=273.15d0 + real(kind=8), parameter,optional :: abscoefi=.272d0 + real(kind=8), parameter,optional :: eslcon2=29.65d0 + real(kind=8), parameter,optional :: eslcon1=17.67d0 + real(kind=8), parameter,optional :: pi=3.1415926535897932384626433d0 + real(kind=8), parameter,optional :: tlclc2=3.5d0 + real(kind=8), parameter,optional :: tlclc3=4.805d0 + real(kind=8), parameter,optional :: rho_g=400.d0 + real(kind=8), parameter,optional :: tlclc1=2840.d0 + real(kind=8), parameter,optional :: tlclc4=55.d0 + real(kind=8), parameter,optional,depend(pi) :: deg_per_rad=180.d0/pi + real(kind=8), parameter,optional :: cpmd=.887d0 + real(kind=8), parameter,optional,depend(rd,g) :: sclht=rd*256.d0/g + real(kind=8), parameter,optional :: ussalr=0.0065d0 + real(kind=8), parameter,optional :: default_fill=9.9692099683868690d36 + real(kind=8), parameter,optional :: rho_s=100.d0 + real(kind=8), parameter,optional,depend(rhowat) :: rho_r=1000.0 + real(kind=8), parameter,optional :: alpha=0.224d0 + real(kind=8), parameter,optional :: celkel_triple=273.16d0 + real(kind=8), parameter,optional :: ezero=6.112d0 + real(kind=8), parameter,optional,depend(rd,ussalr,g) :: expon=0.190163098879 + real(kind=8), parameter,optional :: rgasmd=.608d0 + real(kind=8), parameter,optional :: g=9.81d0 + integer, optional :: errlen=512 + real(kind=8), parameter,optional :: eps=0.622d0 + real(kind=8), parameter,optional :: gamma_seven=720.d0 + real(kind=8), parameter,optional,depend(cpmd,rgasmd) :: gammamd=-0.279 + integer, optional :: algerr=64 + 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 threadsafe - use constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel + use wrf_constants, only: g,eps,rd,ussalr,abscoefi,abscoef,celkel real(kind=8) dimension(ew,ns,nz),intent(in) :: prs real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: tk real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: qci @@ -160,9 +160,9 @@ python module _wrffortran ! in 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 integer intent(in) :: haveqci - integer, optional,intent(in),check(shape(prs,0)==ew),depend(prs) :: ew=shape(prs,0) - integer, optional,intent(in),check(shape(prs,1)==ns),depend(prs) :: ns=shape(prs,1) 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) end subroutine wrfcttcalc subroutine dcomputeabsvort(av,u,v,msfu,msfv,msft,cor,dx,dy,nx,ny,nz,nxp1,nyp1) ! in :_wrffortran:wrf_pvo.f90 threadsafe @@ -183,7 +183,7 @@ python module _wrffortran ! in end subroutine dcomputeabsvort subroutine dcomputepv(pv,u,v,theta,prs,msfu,msfv,msft,cor,dx,dy,nx,ny,nz,nxp1,nyp1) ! in :_wrffortran:wrf_pvo.f90 threadsafe - use constants, only: g + use wrf_constants, only: g real(kind=8) dimension(nx,ny,nz),intent(out,in) :: pv real(kind=8) dimension(nxp1,ny,nz),intent(in),depend(ny,nz) :: u real(kind=8) dimension(nx,nyp1,nz),intent(in),depend(nx,nz) :: v @@ -203,7 +203,7 @@ python module _wrffortran ! in end subroutine dcomputepv subroutine dcomputepw(p,tv,qv,ht,pw,nx,ny,nz,nzh) ! in :_wrffortran:wrf_pw.f90 threadsafe - use constants, only: rd + use wrf_constants, only: rd real(kind=8) dimension(nx,ny,nz),intent(in) :: p real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: tv real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: qv @@ -216,7 +216,7 @@ python module _wrffortran ! in end subroutine dcomputepw subroutine dcalrelhl(u,v,ght,ter,top,sreh,miy,mjx,mkzh) ! in :_wrffortran:wrf_relhl.f90 threadsafe - use constants, only: rad_per_deg,deg_per_rad,pi + use wrf_constants, only: rad_per_deg,deg_per_rad,pi real(kind=8) dimension(miy,mjx,mkzh),intent(in) :: u real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: v real(kind=8) dimension(miy,mjx,mkzh),intent(in),depend(miy,mjx,mkzh) :: ght @@ -229,7 +229,7 @@ python module _wrffortran ! in end subroutine dcalrelhl subroutine wetbulbcalc(prs,tmk,qvp,twb,nx,ny,nz,psafile,errstat,errmsg) ! in :_wrffortran:wrf_rip_phys_routines.f90 threadsafe - use constants, only: tlclc2,tlclc3,tlclc1,tlclc4,eps,thtecon3,gammamd,thtecon1,algerr,gamma,thtecon2 + use wrf_constants, only: tlclc2,tlclc3,tlclc1,tlclc4,eps,thtecon3,gammamd,thtecon1,algerr,gamma,thtecon2 real(kind=8) dimension(nx,ny,nz),intent(in) :: prs real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: tmk real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: qvp @@ -243,7 +243,7 @@ python module _wrffortran ! in end subroutine wetbulbcalc subroutine omgcalc(qvp,tmk,www,prs,omg,mx,my,mz) ! in :_wrffortran:wrf_rip_phys_routines.f90 threadsafe - use constants, only: rd,eps,g + use wrf_constants, only: rd,eps,g real(kind=8) dimension(mx,my,mz),intent(in) :: qvp real(kind=8) dimension(mx,my,mz),intent(in),depend(mx,my,mz) :: tmk real(kind=8) dimension(mx,my,mz),intent(in),depend(mx,my,mz) :: www @@ -255,7 +255,7 @@ python module _wrffortran ! in end subroutine omgcalc subroutine virtual_temp(temp,ratmix,tv,nx,ny,nz) ! in :_wrffortran:wrf_rip_phys_routines.f90 threadsafe - use constants, only: eps + use wrf_constants, only: eps real(kind=8) dimension(nx,ny,nz),intent(in) :: temp real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: ratmix real(kind=8) dimension(nx,ny,nz),intent(out,in),depend(nx,ny,nz) :: tv @@ -263,9 +263,20 @@ python module _wrffortran ! in integer, optional,intent(in),check(shape(temp,1)==ny),depend(temp) :: ny=shape(temp,1) integer, optional,intent(in),check(shape(temp,2)==nz),depend(temp) :: nz=shape(temp,2) end subroutine virtual_temp + subroutine testfunc(a,b,c,nx,ny,nz,errstat,errmsg) ! in :_wrffortran:wrf_testfunc.f90 + use wrf_constants, only: algerr + real(kind=8) dimension(nx,ny,nz),intent(in) :: a + real(kind=8) dimension(nx,ny,nz),intent(out,in),depend(nx,ny,nz) :: b + character*(*) intent(in) :: c + integer, optional,intent(in),check(shape(a,0)==nx),depend(a) :: nx=shape(a,0) + integer, optional,intent(in),check(shape(a,1)==ny),depend(a) :: ny=shape(a,1) + integer, optional,intent(in),check(shape(a,2)==nz),depend(a) :: nz=shape(a,2) + integer intent(inout) :: errstat + character*(*) intent(inout) :: errmsg + end subroutine testfunc subroutine dcomputepi(pi,pressure,nx,ny,nz) ! in :_wrffortran:wrf_user.f90 threadsafe - use constants, only: rd,p1000mb,cp + use wrf_constants, only: rd,p1000mb,cp real(kind=8) dimension(nx,ny,nz),intent(out,in) :: pi real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: pressure integer, optional,intent(in),check(shape(pi,0)==nx),depend(pi) :: nx=shape(pi,0) @@ -274,7 +285,7 @@ python module _wrffortran ! in end subroutine dcomputepi subroutine dcomputetk(tk,pressure,theta,nx) ! in :_wrffortran:wrf_user.f90 threadsafe - use constants, only: rd,p1000mb,cp + use wrf_constants, only: rd,p1000mb,cp real(kind=8) dimension(nx),intent(out,in) :: tk real(kind=8) dimension(nx),intent(in),depend(nx) :: pressure real(kind=8) dimension(nx),intent(in),depend(nx) :: theta @@ -325,7 +336,7 @@ python module _wrffortran ! in end subroutine dinterp1d subroutine dcomputeseaprs(nx,ny,nz,z,t,p,q,sea_level_pressure,t_sea_level,t_surf,level,errstat,errmsg) ! in :_wrffortran:wrf_user.f90 threadsafe - use constants, only: rd,ussalr,algerr,g + use wrf_constants, only: rd,ussalr,algerr,g integer, optional,intent(in),check(shape(z,0)==nx),depend(z) :: nx=shape(z,0) integer, optional,intent(in),check(shape(z,1)==ny),depend(z) :: ny=shape(z,1) integer, optional,intent(in),check(shape(z,2)==nz),depend(z) :: nz=shape(z,2) @@ -360,7 +371,7 @@ python module _wrffortran ! in end subroutine filter2d subroutine dcomputerh(qv,p,t,rh,nx) ! in :_wrffortran:wrf_user.f90 threadsafe - use constants, only: rv,eps,rd,ezero,eslcon2,eslcon1,celkel + use wrf_constants, only: rv,eps,rd,ezero,eslcon2,eslcon1,celkel real(kind=8) dimension(nx),intent(in) :: qv real(kind=8) dimension(nx),intent(in),depend(nx) :: p real(kind=8) dimension(nx),intent(in),depend(nx) :: t @@ -410,6 +421,7 @@ python module _wrffortran ! in end subroutine dcomputetd subroutine dcomputeiclw(iclw,pressure,qc_in,nx,ny,nz) ! in :_wrffortran:wrf_user.f90 threadsafe + use wrf_constants, only: g real(kind=8) dimension(nx,ny),intent(out,in) :: iclw real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny) :: pressure real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: qc_in @@ -417,20 +429,9 @@ python module _wrffortran ! in integer, optional,intent(in),check(shape(iclw,1)==ny),depend(iclw) :: ny=shape(iclw,1) integer, optional,intent(in),check(shape(pressure,2)==nz),depend(pressure) :: nz=shape(pressure,2) end subroutine dcomputeiclw - subroutine testfunc(a,b,c,nx,ny,nz,errstat,errmsg) ! in :_wrffortran:wrf_user.f90 - use constants, only: algerr - real(kind=8) dimension(nx,ny,nz),intent(in) :: a - real(kind=8) dimension(nx,ny,nz),intent(out,in),depend(nx,ny,nz) :: b - character*(*) intent(in) :: c - integer, optional,intent(in),check(shape(a,0)==nx),depend(a) :: nx=shape(a,0) - integer, optional,intent(in),check(shape(a,1)==ny),depend(a) :: ny=shape(a,1) - integer, optional,intent(in),check(shape(a,2)==nz),depend(a) :: nz=shape(a,2) - integer intent(inout) :: errstat - character*(*) intent(inout) :: errmsg - end subroutine testfunc subroutine calcdbz(prs,tmk,qvp,qra,qsn,qgr,sn0,ivarint,iliqskin,dbz,nx,ny,nz) ! in :_wrffortran:wrf_user_dbz.f90 threadsafe - use constants, only: rho_g,gamma_seven,rho_r,rd,rho_s,alpha,rhowat,pi,celkel + use wrf_constants, only: rho_g,gamma_seven,rho_r,rd,rho_s,alpha,rhowat,pi,celkel real(kind=8) dimension(nx,ny,nz),intent(in) :: prs real(kind=8) dimension(nx,ny,nz),intent(in),depend(nx,ny,nz) :: tmk real(kind=8) dimension(nx,ny,nz),intent(inout),depend(nx,ny,nz) :: qvp @@ -447,7 +448,7 @@ python module _wrffortran ! in end subroutine calcdbz subroutine rotatecoords(ilat,ilon,olat,olon,lat_np,lon_np,lon_0,direction) ! in :_wrffortran:wrf_user_latlon_routines.f90 threadsafe - use constants, only: rad_per_deg,deg_per_rad,pi + use wrf_constants, only: rad_per_deg,deg_per_rad,pi real(kind=8) intent(in) :: ilat real(kind=8) intent(in) :: ilon real(kind=8) intent(out,in) :: olat @@ -459,7 +460,7 @@ python module _wrffortran ! in end subroutine rotatecoords subroutine dlltoij(map_proj,truelat1,truelat2,stdlon,lat1,lon1,pole_lat,pole_lon,knowni,knownj,dx,dy,latinc,loninc,lat,lon,loc,errstat,errmsg) ! in :_wrffortran:wrf_user_latlon_routines.f90 threadsafe - use constants, only: rad_per_deg,deg_per_rad,pi,algerr,wrf_earth_radius + use wrf_constants, only: rad_per_deg,deg_per_rad,pi,algerr,wrf_earth_radius integer intent(in) :: map_proj real(kind=8) intent(inout) :: truelat1 real(kind=8) intent(inout) :: truelat2 @@ -482,7 +483,7 @@ python module _wrffortran ! in end subroutine dlltoij subroutine dijtoll(map_proj,truelat1,truelat2,stdlon,lat1,lon1,pole_lat,pole_lon,knowni,knownj,dx,dy,latinc,loninc,ai,aj,loc,errstat,errmsg) ! in :_wrffortran:wrf_user_latlon_routines.f90 threadsafe - use constants, only: rad_per_deg,deg_per_rad,pi,algerr,wrf_earth_radius + use wrf_constants, only: rad_per_deg,deg_per_rad,pi,algerr,wrf_earth_radius integer intent(in) :: map_proj real(kind=8) intent(inout) :: truelat1 real(kind=8) intent(inout) :: truelat2 @@ -518,7 +519,7 @@ python module _wrffortran ! in end subroutine wrf_monotonic function wrf_intrp_value(wvalp0,wvalp1,vlev,vcp0,vcp1,icase,errstat,errmsg) ! in :_wrffortran:wrf_vinterp.f90 threadsafe - use constants, only: sclht,algerr + use wrf_constants, only: sclht,algerr real(kind=8) intent(in) :: wvalp0 real(kind=8) intent(in) :: wvalp1 real(kind=8) intent(in) :: vlev @@ -531,7 +532,7 @@ python module _wrffortran ! in 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 threadsafe - use constants, only: tlclc2,tlclc3,sclht,tlclc4,thtecon3,eps,tlclc1,celkel,exponi,gammamd,ussalr,expon,thtecon1,algerr,gamma,thtecon2 + 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 real(kind=8) dimension(ew,ns,numlevels),intent(out,in),depend(ew,ns) :: dataout real(kind=8) dimension(ew,ns,nz),intent(in),depend(ew,ns,nz) :: pres diff --git a/setup.py b/setup.py index dc12688..25634cc 100755 --- a/setup.py +++ b/setup.py @@ -15,7 +15,8 @@ import numpy.distutils.core ext1 = numpy.distutils.core.Extension( name = "wrf._wrffortran", - sources = ["fortran/constants.f90", + sources = ["fortran/wrf_constants.f90", + "fortran/wrf_testfunc.f90", "fortran/wrf_user.f90", "fortran/rip_cape.f90", "fortran/cloud_fracf.f90", diff --git a/src/wrf/api.py b/src/wrf/api.py index 88f4315..c247da4 100644 --- a/src/wrf/api.py +++ b/src/wrf/api.py @@ -2,11 +2,13 @@ from .config import (xarray_enabled, disable_xarray, enable_xarray, cartopy_enabled, enable_cartopy, enable_cartopy, basemap_enabled, disable_basemap, enable_basemap, pyngl_enabled, enable_pyngl, disable_pyngl) -from .constants import ALL_TIMES, Constants, ConversionFactors +from .constants import ALL_TIMES, Constants, ConversionFactors, ProjectionTypes from .destag import destagger from .routines import getvar from .computation import (xy, interp1d, interp2dxy, interpz3d, slp, tk, td, rh, - uvmet, smooth2d, cape_2d, cape_3d, cloudfrac) + uvmet, smooth2d, cape_2d, cape_3d, cloudfrac, ctt, + dbz, srhel, udhel, avo, pvo, eth, wetbulb, tvirtual, + omega, pw) from .interp import (interplevel, vertcross, interpline, vinterp) from .latlon import (xy_to_ll, ll_to_xy) from .py3compat import (viewitems, viewkeys, viewvalues, isstr, py2round, @@ -21,11 +23,13 @@ __all__ += ["xarray_enabled", "disable_xarray", "enable_xarray", "cartopy_enabled", "enable_cartopy", "enable_cartopy", "basemap_enabled", "disable_basemap", "enable_basemap", "pyngl_enabled", "enable_pyngl", "disable_pyngl"] -__all__ += ["ALL_TIMES", "Constants", "ConversionFactors"] +__all__ += ["ALL_TIMES", "Constants", "ConversionFactors", "ProjectionTypes"] __all__ += ["destagger"] __all__ += ["getvar"] __all__ += ["xy", "interp1d", "interp2dxy", "interpz3d", "slp", "tk", "td", - "rh", "uvmet", "smooth2d", "cape_2d", "cape_3d", "cloudfrac"] + "rh", "uvmet", "smooth2d", "cape_2d", "cape_3d", "cloudfrac", + "ctt", "dbz", "srhel", "udhel", "avo", "pvo", "eth", "wetbulb", + "tvirtual", "omega", "pw"] __all__ += ["interplevel", "vertcross", "interpline", "vinterp"] __all__ += ["xy_to_ll", "ll_to_xy"] __all__ += ["viewitems", "viewkeys", "viewvalues", "isstr", "py2round", diff --git a/src/wrf/computation.py b/src/wrf/computation.py index 5dc015c..a76a014 100644 --- a/src/wrf/computation.py +++ b/src/wrf/computation.py @@ -7,7 +7,8 @@ import numpy.ma as ma from .constants import Constants from .extension import (_interpz3d, _interp2dxy, _interp1d, _slp, _tk, _td, _rh, _uvmet, _smooth2d, _cape, _cloudfrac, _ctt, _dbz, - _srhel, _udhel, _eth, _wetbulb, _tv, _omega) + _srhel, _udhel, _avo, _pvo, _eth, _wetbulb, _tv, + _omega, _pw) from .util import from_var from .metadecorators import (set_alg_metadata, set_uvmet_alg_metadata, set_interp_metadata, set_cape_alg_metadata, @@ -21,9 +22,9 @@ def xy(field, pivot_point=None, angle=None, start_point=None, end_point=None, @set_interp_metadata("1d") -def interp1d(v_in, z_in, z_out, missingval=Constants.DEFAULT_FILL, +def interp1d(field, z_in, z_out, missingval=Constants.DEFAULT_FILL, meta=True): - return _interp1d(v_in, z_in, z_out, missingval) + return _interp1d(field, z_in, z_out, missingval) @set_interp_metadata("2dxy") @@ -37,31 +38,31 @@ def interpz3d(field3d, z, desiredloc, missingval=Constants.DEFAULT_FILL, return _interpz3d(field3d, z, desiredloc, missingval) -@set_alg_metadata(2, "z", refvarndims=3, units="hpa", +@set_alg_metadata(2, "pres", refvarndims=3, units="hpa", description="sea level pressure") -def slp(z, t, p, q, meta=True): - return _slp(z, t, p, q) +def slp(height, tkel, pres, qv, meta=True): + return _slp(height, tkel, pres, qv) -@set_alg_metadata(3, "pressure", units="K", +@set_alg_metadata(3, "pres", units="K", description="temperature") -def tk(pressure, theta, meta=True): - return _tk(pressure, theta) +def tk(pres, theta, meta=True): + return _tk(pres, theta) -@set_alg_metadata(3, "pressure", units="degC", +@set_alg_metadata(3, "pres", units="degC", description="dew point temperature") -def td(pressure, qv_in, meta=True): - return _td(pressure, qv_in) +def td(pres, qv, meta=True): + return _td(pres, qv) -@set_alg_metadata(3, "pressure", +@set_alg_metadata(3, "pres", description="relative humidity", units=None) -def rh(qv, q, t, meta=True): - return _rh(qv, q, t, meta) +def rh(qv, pres, tkel, meta=True): + return _rh(qv, pres, tkel) -@set_uvmet_alg_metadata() +@set_uvmet_alg_metadata(latarg="lat", windarg="u") def uvmet(u, v, lat, lon, cen_long, cone, meta=True): return _uvmet(u, v, lat, lon, cen_long, cone) @@ -73,15 +74,15 @@ def smooth2d(field, passes, meta=True): return _smooth2d(field, passes) -@set_cape_alg_metadata(is2d=True) -def cape_2d(p_hpa, tk, qvapor, height, terrain, psfc_hpa, ter_follow, +@set_cape_alg_metadata(is2d=True, copyarg="pres_hpa") +def cape_2d(pres_hpa, tkel, qv, height, terrain, psfc_hpa, ter_follow, missing=Constants.DEFAULT_FILL, meta=True): if isinstance(ter_follow, bool): ter_follow = 1 if ter_follow else 0 i3dflag = 0 - cape_cin = _cape(p_hpa, tk, qvapor, height, terrain, psfc_hpa, + cape_cin = _cape(pres_hpa, tkel, qvapor, height, terrain, psfc_hpa, missing, i3dflag, ter_follow) left_dims = cape_cin.shape[1:-3] @@ -102,28 +103,28 @@ def cape_2d(p_hpa, tk, qvapor, height, terrain, psfc_hpa, ter_follow, return ma.masked_values(result, missing) -@set_cape_alg_metadata(is2d=False) -def cape_3d(p_hpa, tk, qvapor, height, terrain, psfc_hpa, ter_follow, +@set_cape_alg_metadata(is2d=False, copyarg="pres_hpa") +def cape_3d(pres_hpa, tkel, qvapor, height, terrain, psfc_hpa, ter_follow, missing=Constants.DEFAULT_FILL, meta=True): if isinstance(ter_follow, bool): ter_follow = 1 if ter_follow else 0 i3dflag = 1 - cape_cin = _cape(p_hpa, tk, qvapor, height, terrain, psfc_hpa, + cape_cin = _cape(pres_hpa, tkel, qvapor, height, terrain, psfc_hpa, missing, i3dflag, ter_follow) return ma.masked_values(cape_cin, missing) -@set_cloudfrac_alg_metadata() -def cloudfrac(p, rh, meta=True): - return _cloudfrac(p, rh) +@set_cloudfrac_alg_metadata(copyarg="pres") +def cloudfrac(pres, relh, meta=True): + return _cloudfrac(pres, relh) -@set_alg_metadata(2, "p_hpa", refvarndims=3, units="degC", +@set_alg_metadata(2, "pres_hpa", refvarndims=3, units="degC", description="cloud top temperature") -def ctt(p_hpa, tk, qv, qcld, ght, ter, qice=None, meta=True): +def ctt(pres_hpa, tkel, qv, qcld, height, terrain, qice=None, meta=True): # Qice and QCLD need to be in g/kg if qice is None: @@ -132,13 +133,14 @@ def ctt(p_hpa, tk, qv, qcld, ght, ter, qice=None, meta=True): else: haveqci = 1 if qice.any() else 0 - ctt = _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci) - return _ctt(p, rh) + return _ctt(pres_hpa, tkel, qice, qcld, qv, height, terrain, haveqci) + -@set_alg_metadata(3, "p", units="dBZ", +@set_alg_metadata(3, "pres", units="dBZ", description="radar reflectivity") -def dbz(p, tk, qv, qr, ivarint, iliqskin, qs=None, qg=None, meta=True): +def dbz(pres, tkel, qv, qr, qs=None, qg=None, use_varint=False, use_liqskin=False, + meta=True): if qs is None: qs = np.zeros(qv.shape, qv.dtype) @@ -147,67 +149,74 @@ def dbz(p, tk, qv, qr, ivarint, iliqskin, qs=None, qg=None, meta=True): qg = np.zeros(qv.shape, qv.dtype) sn0 = 1 if qs.any() else 0 - ivarint = 1 if do_varint else 0 - iliqskin = 1 if do_liqskin else 0 + ivarint = 1 if use_varint else 0 + iliqskin = 1 if use_liqskin else 0 - return _dbz(p, tk, qv, qr, qs, qg, sn0, ivarint, iliqskin) + return _dbz(pres, tkel, qv, qr, qs, qg, sn0, ivarint, iliqskin) -@set_alg_metadata(2, "ter", units="m-2/s-2", +@set_alg_metadata(2, "terrain", units="m-2/s-2", description="storm relative helicity") -def srhel(u, v, z, ter, top): - return _srhel(u, v, z, ter, top=3000.0) +def srhel(u, v, z, terrain, top=3000.0, meta=True): + # u, v get swapped in vertical + _u = np.ascontiguousarray(u[...,::-1,:,:]) + _v = np.ascontiguousarray(v[...,::-1,:,:]) + _z = np.ascontiguousarray(z[...,::-1,:,:]) + + return _srhel(_u, _v, _z, terrain, top) @set_alg_metadata(2, "u", refvarndims=3, units="m-2/s-2", description="updraft helicity") -def udhel(zstag, mapfct, u, v, wstag, dx, dy, bottom=2000.0, top=5000.0): +def udhel(zstag, mapfct, u, v, wstag, dx, dy, bottom=2000.0, top=5000.0, + meta=True): return _udhel(zstag, mapfct, u, v, wstag, dx, dy, bottom, top) # Requires both u an v for dimnames -@set_alg_metadata(3, "u", units="10-5 s-1", +@set_alg_metadata(3, "ustag", units="10-5 s-1", stagdim=-1, stagsubvar="vstag", description="absolute vorticity") -def avo(ustag, vstag, msfu, msfv, msfm, cor, dx, dy): +def avo(ustag, vstag, msfu, msfv, msfm, cor, dx, dy, meta=True): return _avo(ustag, vstag, msfu, msfv, msfm, cor, dx, dy) -@set_alg_metadata(3, "t", units="PVU", +@set_alg_metadata(3, "tkel", units="PVU", description="potential vorticity") -def pvo(ustag, vstag, t, p, msfu, msfv, msfm, cor, dx, dy): - return _pvo(ustag, vstag, t, p, msfu, msfv, msfm, cor, dx, dy) +def pvo(ustag, vstag, tkel, pres, msfu, msfv, msfm, cor, dx, dy, meta=True): + return _pvo(ustag, vstag, tkel, pres, msfu, msfv, msfm, cor, dx, dy) @set_alg_metadata(3, "qv", units="K", description="equivalent potential temperature") -def eth(qv, tk, p): - return _eth(qv, tk, p) +def eth(qv, tkel, pres, meta=True): + return _eth(qv, tkel, pres) -@set_alg_metadata(3, "p", units="K", +@set_alg_metadata(3, "pres", units="K", description="wetbulb temperature") -def wetbulb(p, tk, qv): - return _wetbulb(p, tk, qv) +def wetbulb(pres, tkel, qv, meta=True): + return _wetbulb(pres, tkel, qv) -@set_alg_metadata(3, "tk", units="K", +@set_alg_metadata(3, "tkel", units="K", description="virtual temperature") -def tvirtual(tk, qv): - return _tv(tk, qv) +def tvirtual(tkel, qv, meta=True): + return _tv(tkel, qv) @set_alg_metadata(3, "qv", units="Pa/s", description="omega") -def omega(qv, tk, w, p): - return _omega(qv, tk, w, p) +def omega(qv, tkel, w, pres, meta=True): + return _omega(qv, tkel, w, pres) -@set_alg_metadata(2, "p", refvarndims=3, units="kg m-2", +@set_alg_metadata(2, "pres", refvarndims=3, units="kg m-2", description="precipitable water") -def pw(p, tk, qv, ht): - tv = _tv(tk, qv) - return _pw(full_p, tv, qv, ht) +def pw(pres, tkel, qv, height, meta=True): + tv = _tv(tkel, qv) + + return _pw(pres, tv, qv, height) diff --git a/src/wrf/constants.py b/src/wrf/constants.py index b364acc..50d0d9f 100755 --- a/src/wrf/constants.py +++ b/src/wrf/constants.py @@ -4,14 +4,14 @@ from __future__ import (absolute_import, division, print_function, import numpy as np from .py3compat import viewitems -from wrf._wrffortran import constants +from wrf._wrffortran import wrf_constants ALL_TIMES = None class Constants(object): pass -for key,val in viewitems(constants.__dict__): +for key,val in viewitems(wrf_constants.__dict__): setattr(Constants, key.upper(), np.asscalar(val)) class ConversionFactors(object): @@ -28,4 +28,11 @@ class ConversionFactors(object): M_TO_FT = 3.28084 M_TO_MILES = .000621371 +class ProjectionTypes(object): + ZERO = 0 + LAMBERT_CONFORMAL = 1 + POLAR_STEREOGRAPHIC = 2 + MERCATOR = 3 + LAT_LON = 6 + \ No newline at end of file diff --git a/src/wrf/dbz.py b/src/wrf/dbz.py index dd30da4..a374cb6 100755 --- a/src/wrf/dbz.py +++ b/src/wrf/dbz.py @@ -15,14 +15,14 @@ from .metadecorators import copy_and_set_metadata units="dBz") def get_dbz(wrfnc, timeidx=0, method="cat", squeeze=True, cache=None, meta=True, - do_varint=False, do_liqskin=False): + use_varint=False, use_liqskin=False): """ Return the dbz - do_varint - do variable intercept (if False, constants are used. Otherwise, + use_varint - do variable intercept (if False, constants are used. Otherwise, intercepts are calculated using a technique from Thompson, Rasmussen, and Manning (2004, Monthly Weather Review, Vol. 132, No. 2, pp. 519-542.) - do_liqskin - do liquid skin for snow (frozen particles above freezing scatter + use_liqskin - do liquid skin for snow (frozen particles above freezing scatter as liquid) """ @@ -57,8 +57,8 @@ def get_dbz(wrfnc, timeidx=0, method="cat", # If qsnow is not all 0, set sn0 to 1 sn0 = 1 if qs.any() else 0 - ivarint = 1 if do_varint else 0 - iliqskin = 1 if do_liqskin else 0 + ivarint = 1 if use_varint else 0 + iliqskin = 1 if use_liqskin else 0 return _dbz(full_p, tk, qv, qr, qs, qg, sn0, ivarint, iliqskin) diff --git a/src/wrf/decorators.py b/src/wrf/decorators.py index bde6102..d9dc645 100644 --- a/src/wrf/decorators.py +++ b/src/wrf/decorators.py @@ -42,168 +42,18 @@ def _calc_out_dims(outvar, left_dims): left_dims = [x for x in left_dims] right_dims = [x for x in outvar.shape] return left_dims + right_dims - -def handle_left_iter(ref_var_expected_dims, ref_var_idx=-1, - ref_var_name=None, - ignore_args=None, ignore_kargs=None): - """Decorator to handle iterating over leftmost dimensions when using - multiple files and/or multiple times. - - Arguments: - - ref_var_expected_dims - the number of dimensions that the Fortran - algorithm is expecting for the reference variable - - ref_var_idx - the index in args used as the reference variable for - calculating leftmost dimensions - - ref_var_name - the keyword argument name for kwargs used as the - reference varible for calculating leftmost dimensions - - alg_out_fixed_dims - additional fixed dimension sizes for the - numerical algorithm (e.g. uvmet has a fixed left dimsize of 2) - - ignore_args - indexes of any arguments which should be passed - directly without any slicing - - ignore_kargs - keys of any keyword arguments which should be passed - directly without slicing - - ref_stag_dim - in some cases the reference variable dimensions - have a staggered dimension that needs to be corrected when calculating - the output dimensions (e.g. avo) - - """ - @wrapt.decorator - def func_wrapper(wrapped, instance, args, kwargs): - _ignore_args = ignore_args if ignore_args is not None else () - _ignore_kargs = ignore_kargs if ignore_kargs is not None else () - - if ref_var_idx >= 0: - ref_var = args[ref_var_idx] - else: - ref_var = kwargs[ref_var_name] - - ref_var_shape = ref_var.shape - extra_dim_num = ref_var.ndim - ref_var_expected_dims - - # No special left side iteration, return the function result - if (extra_dim_num == 0): - return wrapped(*args, **kwargs) - - # Start by getting the left-most 'extra' dims - extra_dims = ref_var_shape[0:extra_dim_num] - - mask_output = False - out_inited = False - for left_idxs in iter_left_indexes(extra_dims): - # Make the left indexes plus a single slice object - # The single slice will handle all the dimensions to - # the right (e.g. [1,1,:]) - left_and_slice_idxs = left_idxs + (slice(None), ) - - # Slice the args if applicable - new_args = [arg[left_and_slice_idxs] - if i not in _ignore_args else arg - for i,arg in enumerate(args)] - - # Slice the kwargs if applicable - new_kargs = {key:(val[left_and_slice_idxs] - if key not in _ignore_kargs else val) - for key,val in viewitems(kwargs)} - - # Skip the possible empty/missing arrays for the join method - skip_missing = False - if out_inited: - for i,arg in enumerate(new_args): - if isinstance(arg, DataArray): - arr = npvalues(arg) - elif isinstance(arg, np.ndarray): - arr = arg - else: - continue - - if isinstance(arr, np.ma.MaskedArray): - if arr.mask.all(): - if isinstance(output, np.ndarray): - output[left_and_slice_idxs] = ( - Constants.DEFAULT_FILL) - else: - for arr in output: - arr[left_and_slice_idxs] = ( - Constants.DEFAULT_FILL) - skip_missing = True - mask_output = True - - if skip_missing: - continue - - # Call the numerical routine - result = wrapped(*new_args, **new_kargs) - - if isinstance(result, np.ndarray): - # Output array - if not out_inited: - outdims = _calc_out_dims(result, extra_dims) - if not isinstance(result, ma.MaskedArray): - output = np.empty(outdims, ref_var.dtype) - masked = False - else: - output = ma.MaskedArray( - np.zeros(outdims, ref_var.dtype), - mask=np.zeros(outdims, np.bool_), - fill_value=result.fill_value) - masked = True - - out_inited = True - - if not masked: - output[left_and_slice_idxs] = result[:] - else: - output.data[left_and_slice_idxs] = result.data[:] - output.mask[left_and_slice_idxs] = result.mask[:] - - else: # This should be a list or a tuple (cape) - if not out_inited: - outdims = _calc_out_dims(result[0], extra_dims) - if not isinstance(result[0], ma.MaskedArray): - output = [np.empty(outdims, ref_var.dtype) - for i in py3range(len(result))] - masked = False - else: - output = [ma.MaskedArray( - np.zeros(outdims, ref_var.dtype), - mask=np.zeros(outdims, np.bool_), - fill_value=result[0].fill_value) - for i in py3range(len(result))] - masked = True - - out_inited = True - - for i,outarr in enumerate(result): - if not masked: - output[i][left_and_slice_idxs] = outarr[:] - else: - output[i].data[left_and_slice_idxs] = outarr.data[:] - output[i].mask[left_and_slice_idxs] = outarr.mask[:] - - - if mask_output: - if isinstance(output, np.ndarray): - output = ma.masked_values(output, Constants.DEFAULT_FILL) - else: - output = tuple(ma.masked_values(arr, Constants.DEFAULT_FILL) - for arr in output) - - return output - - return func_wrapper - -def left_iter_nocopy(ref_var_expected_dims, - ref_var_right_ndims, - insert_dims=None, - ref_var_idx=None, - ref_var_name=None, - ignore_args=None, - ignore_kargs=None, - outviews="outview", - alg_dtype=np.float64, - cast_output=True): +def left_iteration(ref_var_expected_dims, + ref_var_right_ndims, + insert_dims=None, + ref_var_idx=None, + ref_var_name=None, + ignore_args=None, + ignore_kargs=None, + outviews="outview", + alg_dtype=np.float64, + cast_output=True): """Decorator to handle iterating over leftmost dimensions when using multiple files and/or multiple times. @@ -352,7 +202,7 @@ def left_iter_nocopy(ref_var_expected_dims, return func_wrapper -def handle_casting(ref_idx=0, arg_idxs=None, karg_names=None, +def cast_type(ref_idx=0, arg_idxs=None, karg_names=None, alg_dtype=np.float64, outviews="outview"): """Decorator to handle casting to/from required dtype used in @@ -417,7 +267,7 @@ def _extract_and_transpose(arg, do_transpose): return arg -def handle_extract_transpose(do_transpose=True, outviews="outview"): +def extract_and_transpose(do_transpose=True, outviews="outview"): """Decorator to extract the data array from a DataArray and also transposes the view of the data if the data is not fortran contiguous. @@ -475,7 +325,13 @@ def check_args(refvaridx, refvarndim, rightdims, stagger=None, def func_wrapper(wrapped, instance, args, kwargs): refvar = args[refvaridx] - extra_dims = refvar.ndim - refvarndim + try: + _ndim = refvar.ndim + except AttributeError: + raise ValueError("argument {} is not an arraylike " + "object".format(refvaridx)) + else: + extra_dims = refvar.ndim - refvarndim # Always use unstaggered as the basis of comparison if refstagdim is not None: @@ -496,6 +352,12 @@ def check_args(refvaridx, refvarndim, rightdims, stagger=None, var = args[i] + try: + _ = var.ndim + except AttributeError: + raise ValueError("argument {} is not an arraylike " + "object".format(i)) + right_var_ndims = rightdims[i] # Check that the number of dims is correct diff --git a/src/wrf/destag.py b/src/wrf/destag.py index a53811f..e8e4c45 100755 --- a/src/wrf/destag.py +++ b/src/wrf/destag.py @@ -1,12 +1,12 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) -from .decorators import handle_extract_transpose +from .decorators import extract_and_transpose from .metadecorators import set_destag_metadata @set_destag_metadata() -@handle_extract_transpose(do_transpose=False) +@extract_and_transpose(do_transpose=False) def destagger(var, stagger_dim, meta=False): """ De-stagger the variable. diff --git a/src/wrf/extension.py b/src/wrf/extension.py index 03e6858..d4ea1cb 100755 --- a/src/wrf/extension.py +++ b/src/wrf/extension.py @@ -4,19 +4,7 @@ from __future__ import (absolute_import, division, print_function, import numpy as np from .constants import Constants -#from .psadlookup import get_lookup_tables -# Old way -#from ._wrfext import (f_interpz3d, f_interp2dxy,f_interp1d, -# f_computeslp, f_computetk, f_computetd, f_computerh, -# f_computeabsvort,f_computepvo, f_computeeth, -# f_computeuvmet, -# f_computeomega, f_computetv, f_computewetbulb, -# f_computesrh, f_computeuh, f_computepw, f_computedbz, -# f_lltoij, f_ijtoll, f_converteta, f_computectt, -# f_monotonic, f_filter2d, f_vintrp) -#from ._wrfcape import f_computecape - -# New way + from ._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d, dcomputeseaprs, dfilter2d, dcomputerh, dcomputeuvmet, dcomputetd, dcapecalc3d, dcloudfrac, wrfcttcalc, @@ -25,24 +13,13 @@ from ._wrffortran import (dcomputetk, dinterp3dz, dinterp2dxy, dinterp1d, omgcalc, virtual_temp, wetbulbcalc, dcomputepw, wrf_monotonic, wrf_vintrp) -from .decorators import (handle_left_iter, handle_casting, - handle_extract_transpose) -from .decorators import (left_iter_nocopy, check_args) +from .decorators import (left_iteration, cast_type, + extract_and_transpose, check_args) from .util import combine_dims, npbytes_to_str, psafilepath from .py3compat import py3range from .specialdec import (uvmet_left_iter_nocopy, cape_left_iter, cloudfrac_left_iter) -#__all__ = [] -# __all__ += ["FortranException", "computeslp", "computetk", "computetd", -# "computerh", "computeavo", "computepvo", "computeeth", -# "computeuvmet","computeomega", "computetv", "computesrh", -# "computeuh", "computepw","computedbz","computecape", -# "computeij", "computell", "computeeta", "computectt", -# "interp2dxy", "interpz3d", "interp1d", "computeinterpline", -# "computevertcross"] -# __all__ += ["", ""] - class FortranError(Exception): def __init__(self, message=None): self._msg = message @@ -55,24 +32,13 @@ class FortranError(Exception): # IMPORTANT! Unless otherwise noted, all variables used in the routines -# below assume that fortran-ordering views are being used. - - - -# @handle_left_iter(3,0, ignore_args=(2,3)) -# @handle_casting(arg_idxs=(0,1)) -# @handle_extract_transpose() -# def interpz3d(field3d, z, desiredloc, missingval): -# result = f_interpz3d(field3d, -# z, -# desiredloc, -# missingval) -# return result +# below assume that fortran-ordered views are being used. This allows +# f2py to pass the array pointers directly to the fortran routine. @check_args(0, 3, (3, 3, None, None)) -@left_iter_nocopy(3, 2, ref_var_idx=0, ignore_args=(2,3)) -@handle_casting(arg_idxs=(0,1)) -@handle_extract_transpose() +@left_iteration(3, 2, ref_var_idx=0, ignore_args=(2,3)) +@cast_type(arg_idxs=(0,1)) +@extract_and_transpose() def _interpz3d(field3d, z, desiredloc, missingval, outview=None): if outview is None: outview = np.empty(field3d.shape[0:2], np.float64, order="F") @@ -84,19 +50,12 @@ def _interpz3d(field3d, z, desiredloc, missingval, outview=None): missingval) return result -# @handle_left_iter(3,0, ignore_args=(1,)) -# @handle_casting(arg_idxs=(0,1)) -# @handle_extract_transpose() -# def interp2dxy(field3d, xy): -# result = f_interp2dxy(field3d, -# xy) -# return result @check_args(0, 3, (3,)) -@left_iter_nocopy(3, combine_dims([(0,-3),(1,-2)]), ref_var_idx=0, +@left_iteration(3, combine_dims([(0,-3),(1,-2)]), ref_var_idx=0, ignore_args=(1,)) -@handle_casting(arg_idxs=(0,1)) -@handle_extract_transpose() +@cast_type(arg_idxs=(0,1)) +@extract_and_transpose() def _interp2dxy(field3d, xy, outview=None): if outview is None: outview = np.empty((xy.shape[-1], field3d.shape[-1]), np.float64, @@ -107,21 +66,11 @@ def _interp2dxy(field3d, xy, outview=None): xy) return result -# @handle_left_iter(1, 0, ignore_args=(2,3)) -# @handle_casting(arg_idxs=(0,1,2)) -# @handle_extract_transpose() -# def interp1d(v_in, z_in, z_out, missingval): -# result = f_interp1d(v_in, -# z_in, -# z_out, -# missingval) -# -# return result @check_args(0, 1, (1,1,None,None)) -@left_iter_nocopy(1, combine_dims([(2,0)]), ref_var_idx=0, ignore_args=(2,3)) -@handle_casting(arg_idxs=(0,1,2)) -@handle_extract_transpose() +@left_iteration(1, combine_dims([(2,0)]), ref_var_idx=0, ignore_args=(2,3)) +@cast_type(arg_idxs=(0,1,2)) +@extract_and_transpose() def _interp1d(v_in, z_in, z_out, missingval, outview=None): if outview is None: @@ -135,22 +84,11 @@ def _interp1d(v_in, z_in, z_out, missingval, outview=None): return result -# @handle_left_iter(3, 0, ignore_args=(1,4,3)) -# @handle_casting(arg_idxs=(0,)) -# @handle_extract_transpose(do_transpose=False) -# def computevertcross(field3d, xy, var2dz, z_var2d, missingval): -# var2d = np.empty((z_var2d.size, xy.shape[0]), dtype=var2dz.dtype) -# var2dtmp = interp2dxy(field3d, xy) -# -# for i in py3range(xy.shape[0]): -# var2d[:,i] = interp1d(var2dtmp[:,i], var2dz[:,i], z_var2d, missingval) -# -# return var2d - -@left_iter_nocopy(3, combine_dims([(3,0), (1,0)]), + +@left_iteration(3, combine_dims([(3,0), (1,0)]), ref_var_idx=0, ignore_args=(1,3,4)) -@handle_casting(arg_idxs=(0,)) -@handle_extract_transpose(do_transpose=False) +@cast_type(arg_idxs=(0,)) +@extract_and_transpose(do_transpose=False) def _vertcross(field3d, xy, var2dz, z_var2d, missingval, outview=None): # Note: This is using C-ordering if outview is None: @@ -164,23 +102,10 @@ def _vertcross(field3d, xy, var2dz, z_var2d, missingval, outview=None): return outview -# @handle_left_iter(2, 0, ignore_args=(1,)) -# @handle_casting(arg_idxs=(0,)) -# @handle_extract_transpose(do_transpose=False) -# def interpline(field2d, xy): -# -# tmp_shape = [1] + [x for x in field2d.shape] -# var2dtmp = np.empty(tmp_shape, field2d.dtype) -# var2dtmp[0,:,:] = field2d[:,:] -# -# var1dtmp = interp2dxy(var2dtmp, xy) -# -# return var1dtmp[0,:] - - -@left_iter_nocopy(2, combine_dims([(1,0)]), ref_var_idx=0, ignore_args=(1,)) -@handle_casting(arg_idxs=(0,)) -@handle_extract_transpose(do_transpose=False) + +@left_iteration(2, combine_dims([(1,0)]), ref_var_idx=0, ignore_args=(1,)) +@cast_type(arg_idxs=(0,)) +@extract_and_transpose(do_transpose=False) def _interpline(field2d, xy, outview=None): # Note: This is using C-ordering if outview is None: @@ -196,31 +121,11 @@ def _interpline(field2d, xy, outview=None): return outview -# @handle_left_iter(3,0) -# @handle_casting(arg_idxs=(0,1,2,3)) -# @handle_extract_transpose() -# def computeslp(z, t, p, q): -# -# t_surf = np.zeros(z.shape[0:2], np.float64, order="F") -# t_sea_level = np.zeros(z.shape[0:2], np.float64, order="F") -# level = np.zeros(z.shape[0:2], np.int32, order="F") -# -# result = f_computeslp(z, -# t, -# p, -# q, -# t_sea_level, -# t_surf, -# level, -# FortranException()) -# -# return result - @check_args(0, 3, (3,3,3,3)) -@left_iter_nocopy(3, 2, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1,2,3)) -@handle_extract_transpose() +@left_iteration(3, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0,1,2,3)) +@extract_and_transpose() def _slp(z, t, p, q, outview=None): t_surf = np.zeros(z.shape[0:2], np.float64, order="F") @@ -249,24 +154,11 @@ def _slp(z, t, p, q, outview=None): return result -# @handle_left_iter(3,0) -# @handle_casting(arg_idxs=(0,1)) -# @handle_extract_transpose() -# def computetk(pres, theta): -# # No need to transpose here since operations on 1D array -# shape = pres.shape -# result = f_computetk(pres.ravel(order="A"), -# theta.ravel(order="A")) -# result = np.reshape(result, shape, order="F") -# -# return result - - @check_args(0, 3, (3,3)) -@left_iter_nocopy(3, 3, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0) +@cast_type(arg_idxs=(0,1)) +@extract_and_transpose() def _tk(pressure, theta, outview=None): # No need to transpose here since operations on 1D array shape = pressure.shape @@ -279,20 +171,11 @@ def _tk(pressure, theta, outview=None): return result -# Note: No left iteration decorator needed with 1D arrays -# @handle_casting(arg_idxs=(0,1)) -# @handle_extract_transpose() -# def computetd(pressure, qv_in): -# shape = pressure.shape -# result = f_computetd(pressure.ravel(order="A"), -# qv_in.ravel(order="A")) -# result = np.reshape(result, shape, order="F") -# return result @check_args(0, 2, (2,2)) -@left_iter_nocopy(2, 2, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1)) -@handle_extract_transpose() +@left_iteration(2, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0,1)) +@extract_and_transpose() def _td(pressure, qv_in, outview=None): shape = pressure.shape if outview is None: @@ -305,22 +188,11 @@ def _td(pressure, qv_in, outview=None): return result -# # Note: No decorator needed with 1D arrays -# @handle_casting(arg_idxs=(0,1,2)) -# @handle_extract_transpose() -# def computerh(qv, q, t): -# shape = qv.shape -# result = f_computerh(qv.ravel(order="A"), -# q.ravel(order="A"), -# t.ravel(order="A")) -# result = np.reshape(result, shape, order="F") -# return result - -# Note: No left iteration decorator needed with 1D arrays + @check_args(0, 2, (2,2,2)) -@left_iter_nocopy(2, 2, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1,2)) -@handle_extract_transpose() +@left_iteration(2, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0,1,2)) +@extract_and_transpose() def _rh(qv, q, t, outview=None): shape = qv.shape if outview is None: @@ -333,30 +205,16 @@ def _rh(qv, q, t, outview=None): return result -#@handle_left_iter(3,0, ignore_args=(6,7)) -#@handle_casting(arg_idxs=(0,1,2,3,4,5)) -#@handle_extract_transpose() -#def computeavo(u, v, msfu, msfv, msfm, cor, dx, dy): -# result = f_computeabsvort(u, -# v, -# msfu, -# msfv, -# msfm, -# cor, -# dx, -# dy) -# -# return result # Note: combining the -3 and -2 dimensions from u, then the -1 dimension # from v @check_args(0, 3, (3,3,2,2,2,2), stagger=(-1,-2,-1,-2,None,None), refstagdim=-1) -@left_iter_nocopy(3, combine_dims([(0, (-3,-2)), +@left_iteration(3, combine_dims([(0, (-3,-2)), (1, (-1,))]), ref_var_idx=0, ignore_args=(6,7)) -@handle_casting(arg_idxs=(0,1,2,3,4,5)) -@handle_extract_transpose() +@cast_type(arg_idxs=(0,1,2,3,4,5)) +@extract_and_transpose() def _avo(u, v, msfu, msfv, msfm, cor, dx, dy, outview=None): if outview is None: outshape = (v.shape[0],) + u.shape[1:] @@ -374,30 +232,13 @@ def _avo(u, v, msfu, msfv, msfm, cor, dx, dy, outview=None): return result -#@handle_left_iter(3,2, ignore_args=(8,9)) -#@handle_casting(arg_idxs=(0,1,2,3,4,5,6,7)) -#@handle_extract_transpose() -#def computepvo(u, v, theta, prs, msfu, msfv, msfm, cor, dx, dy): -# -# result = f_computepvo(u, -# v, -# theta, -# prs, -# msfu, -# msfv, -# msfm, -# cor, -# dx, -# dy) -# -# return result @check_args(0, 3, (3,3,3,3,2,2,2,2), stagger=(-1,-2,None,None,-1,-2,None, None), refstagdim=-1) -@left_iter_nocopy(3, 3, ref_var_idx=2, ignore_args=(8,9)) -@handle_casting(arg_idxs=(0,1,2,3,4,5,6,7)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=2, ignore_args=(8,9)) +@cast_type(arg_idxs=(0,1,2,3,4,5,6,7)) +@extract_and_transpose() def _pvo(u, v, theta, prs, msfu, msfv, msfm, cor, dx, dy, outview=None): if outview is None: outview = np.empty_like(prs) @@ -416,22 +257,11 @@ def _pvo(u, v, theta, prs, msfu, msfv, msfm, cor, dx, dy, outview=None): return result -#@handle_left_iter(3,0) -#@handle_casting(arg_idxs=(0,1,2)) -#@handle_extract_transpose() -#def computeeth(qv, tk, p): -# -# result = f_computeeth(qv, -# tk, -# p) -# -# return result - @check_args(0, 3, (3,3,3)) -@left_iter_nocopy(3, 3, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1,2)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0) +@cast_type(arg_idxs=(0,1,2)) +@extract_and_transpose() def _eth(qv, tk, p, outview=None): if outview is None: outview = np.empty_like(qv) @@ -443,43 +273,10 @@ def _eth(qv, tk, p, outview=None): return result -# @uvmet_left_iter() -# @handle_casting(arg_idxs=(0,1,2,3)) -# @handle_extract_transpose() -# def computeuvmet(u, v, lat, lon, cen_long, cone): -# longca = np.zeros((lat.shape[0], lat.shape[1]), np.float64, order="F") -# longcb = np.zeros((lon.shape[0], lon.shape[1]), np.float64, order="F") -# rpd = Constants.PI/180. -# -# -# # Make the 2D array a 3D array with 1 dimension -# if u.ndim < 3: -# u = u[:, :, np.newaxis] -# v = v[:, :, np.newaxis] -# -# # istag will always be false since winds are destaggered already -# # Missing values don't appear to be used, so setting to 0 -# result = f_computeuvmet(u, -# v, -# longca, -# longcb, -# lon, -# lat, -# cen_long, -# cone, -# rpd, -# 0, -# False, -# 0, -# 0, -# 0) -# -# -# return np.squeeze(result) @uvmet_left_iter_nocopy() -@handle_casting(arg_idxs=(0,1,2,3)) -@handle_extract_transpose() +@cast_type(arg_idxs=(0,1,2,3)) +@extract_and_transpose() def _uvmet(u, v, lat, lon, cen_long, cone, isstag=0, has_missing=False, umissing=Constants.DEFAULT_FILL, vmissing=Constants.DEFAULT_FILL, uvmetmissing=Constants.DEFAULT_FILL, outview=None): @@ -509,22 +306,11 @@ def _uvmet(u, v, lat, lon, cen_long, cone, isstag=0, has_missing=False, return result -#@handle_left_iter(3,0) -#@handle_casting(arg_idxs=(0,1,2,3)) -#@handle_extract_transpose() -#def computeomega(qv, tk, w, p): -# -# result = f_computeomega(qv, -# tk, -# w, -# p) -# -# return result @check_args(0, 3, (3,3,3,3)) -@left_iter_nocopy(3, 3, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1,2,3)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0) +@cast_type(arg_idxs=(0,1,2,3)) +@extract_and_transpose() def _omega(qv, tk, w, p, outview=None): if outview is None: outview = np.empty_like(qv) @@ -537,20 +323,11 @@ def _omega(qv, tk, w, p, outview=None): return result -#@handle_left_iter(3,0) -#@handle_casting(arg_idxs=(0,1)) -#@handle_extract_transpose() -#def computetv(tk, qv): -# result = f_computetv(tk, -# qv) -# -# return result - @check_args(0, 3, (3,3)) -@left_iter_nocopy(3, 3, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0) +@cast_type(arg_idxs=(0,1)) +@extract_and_transpose() def _tv(tk, qv, outview=None): if outview is None: outview = np.empty_like(tk) @@ -561,28 +338,11 @@ def _tv(tk, qv, outview=None): return result -#@handle_left_iter(3,0) -#@handle_casting(arg_idxs=(0,1,2)) -#@handle_extract_transpose() -#def computewetbulb(p, tk, qv): -# PSADITHTE, PSADIPRS, PSADITMK = get_lookup_tables() -# PSADITMK = PSADITMK.T -# -# result = f_computewetbulb(p, -# tk, -# qv, -# PSADITHTE, -# PSADIPRS, -# PSADITMK, -# FortranError()) -# -# return result - @check_args(0, 3, (3,3,3)) -@left_iter_nocopy(3, 3, ref_var_idx=0, ignore_args=(3,)) -@handle_casting(arg_idxs=(0,1,2)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0, ignore_args=(3,)) +@cast_type(arg_idxs=(0,1,2)) +@extract_and_transpose() def _wetbulb(p, tk, qv, psafile=psafilepath(), outview=None): if outview is None: outview = np.empty_like(p) @@ -603,23 +363,11 @@ def _wetbulb(p, tk, qv, psafile=psafilepath(), outview=None): return result -#@handle_left_iter(3,0, ignore_args=(4,)) -#@handle_casting(arg_idxs=(0,1,2,3)) -#@handle_extract_transpose() -#def computesrh(u, v, z, ter, top): -# -# result = f_computesrh(u, -# v, -# z, -# ter, -# top) -# -# return result @check_args(0, 3, (3,3,3,2)) -@left_iter_nocopy(3, 2, ref_var_idx=0, ignore_args=(4,)) -@handle_casting(arg_idxs=(0,1,2,3)) -@handle_extract_transpose() +@left_iteration(3, 2, ref_var_idx=0, ignore_args=(4,)) +@cast_type(arg_idxs=(0,1,2,3)) +@extract_and_transpose() def _srhel(u, v, z, ter, top, outview=None): if outview is None: outview = np.empty_like(ter) @@ -633,34 +381,11 @@ def _srhel(u, v, z, ter, top, outview=None): return result -#@handle_left_iter(3,2, ignore_args=(5,6,7,8)) -#@handle_casting(arg_idxs=(0,1,2,3,4)) -#@handle_extract_transpose() -#def computeuh(zp, mapfct, u, v, wstag, dx, dy, bottom, top): -# -# tem1 = np.zeros((u.shape[0], u.shape[1], u.shape[2]), np.float64, -# order="F") -# tem2 = np.zeros((u.shape[0], u.shape[1], u.shape[2]), np.float64, -# order="F") -# -# result = f_computeuh(zp, -# mapfct, -# dx, -# dy, -# bottom, -# top, -# u, -# v, -# wstag, -# tem1, -# tem2) -# -# return result @check_args(2, 3, (3,2,3,3,3), stagger=(-3,None,None,None,-3)) -@left_iter_nocopy(3, 2, ref_var_idx=2, ignore_args=(5,6,7,8)) -@handle_casting(arg_idxs=(0,1,2,3,4)) -@handle_extract_transpose() +@left_iteration(3, 2, ref_var_idx=2, ignore_args=(5,6,7,8)) +@cast_type(arg_idxs=(0,1,2,3,4)) +@extract_and_transpose() def _udhel(zstag, mapfct, u, v, wstag, dx, dy, bottom, top, outview=None): if outview is None: outview = np.empty_like(mapfct) @@ -685,25 +410,11 @@ def _udhel(zstag, mapfct, u, v, wstag, dx, dy, bottom, top, outview=None): return result -#@handle_left_iter(3,0) -#@handle_casting(arg_idxs=(0,1,2,3)) -#@handle_extract_transpose() -#def computepw(p, tv, qv, ht): -# -# zdiff = np.zeros((p.shape[0], p.shape[1]), np.float64, order="F") -# result = f_computepw(p, -# tv, -# qv, -# ht, -# zdiff) -# -# return result - @check_args(0, 3, (3,3,3,3), stagger=(None, None, None, -3)) -@left_iter_nocopy(3, 2, ref_var_idx=0) -@handle_casting(arg_idxs=(0,1,2,3)) -@handle_extract_transpose() +@left_iteration(3, 2, ref_var_idx=0) +@cast_type(arg_idxs=(0,1,2,3)) +@extract_and_transpose() def _pw(p, tv, qv, ht, outview=None): if outview is None: @@ -717,27 +428,11 @@ def _pw(p, tv, qv, ht, outview=None): return result -#@handle_left_iter(3,0, ignore_args=(6,7,8)) -#@handle_casting(arg_idxs=(0,1,2,3,4,5)) -#@handle_extract_transpose() -#def computedbz(p, tk, qv, qr, qs, qg, sn0, ivarint, iliqskin): -# -# result = f_computedbz(p, -# tk, -# qv, -# qr, -# qs, -# qg, -# sn0, -# ivarint, -# iliqskin) -# -# return result @check_args(0, 3, (3,3,3,3,3,3)) -@left_iter_nocopy(3, 3, ref_var_idx=0, ignore_args=(6,7,8)) -@handle_casting(arg_idxs=(0,1,2,3,4,5)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0, ignore_args=(6,7,8)) +@cast_type(arg_idxs=(0,1,2,3,4,5)) +@extract_and_transpose() def _dbz(p, tk, qv, qr, qs, qg, sn0, ivarint, iliqskin, outview=None): if outview is None: outview = np.empty_like(p) @@ -756,56 +451,10 @@ def _dbz(p, tk, qv, qr, qs, qg, sn0, ivarint, iliqskin, outview=None): return result -# TODO: Make a new decorator to handle the flipping of the vertical -# The output arrays can be ignored and passed directly without flipping -# Then flip the final result. The copy is unavoidable, but at least it is -# only happening once, not every time. - -# @handle_left_iter(3,0,ignore_args=(6,7,8)) -# @handle_casting(arg_idxs=(0,1,2,3,4,5)) -# @handle_extract_transpose() -# def computecape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow): -# flip_cape = np.zeros(p_hpa.shape[0:3], np.float64, order="F") -# flip_cin = np.zeros(p_hpa.shape[0:3], np.float64, order="F") -# PSADITHTE, PSADIPRS, PSADITMK = get_lookup_tables() -# PSADITMK = PSADITMK.T -# -# # The fortran routine needs pressure to be ascending in z-direction, -# # along with tk,qv,and ht. -# # The extra mumbo-jumbo is so that the view created by numpy is fortran -# # contiguous. 'ascontiguousarray' only works in C ordering, hence the -# # extra transposes. Note, this is probably making a copy -# flip_p = np.ascontiguousarray(p_hpa[:,:,::-1].T).T -# flip_tk = np.ascontiguousarray(tk[:,:,::-1].T).T -# flip_qv = np.ascontiguousarray(qv[:,:,::-1].T).T -# flip_ht = np.ascontiguousarray(ht[:,:,::-1].T).T -# -# f_computecape(flip_p, -# flip_tk, -# flip_qv, -# flip_ht, -# ter, -# sfp, -# flip_cape, -# flip_cin, -# PSADITHTE, -# PSADIPRS, -# PSADITMK, -# missing, -# i3dflag, -# ter_follow, -# FortranException()) -# -# # Need to flip the vertical back to decending pressure with height. -# cape = np.ascontiguousarray(flip_cape[:,:,::-1].T).T -# cin = np.ascontiguousarray(flip_cin[:,:,::-1].T).T -# -# return (cape, cin) - @check_args(0, 3, (3,3,3,3,2,2)) @cape_left_iter() -@handle_casting(arg_idxs=(0,1,2,3,4,5), outviews=("capeview", "cinview")) -@handle_extract_transpose(outviews=("capeview", "cinview")) +@cast_type(arg_idxs=(0,1,2,3,4,5), outviews=("capeview", "cinview")) +@extract_and_transpose(outviews=("capeview", "cinview")) def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, psafile=psafilepath(), capeview=None, cinview=None): @@ -841,8 +490,8 @@ def _cape(p_hpa, tk, qv, ht, ter, sfp, missing, i3dflag, ter_follow, @check_args(0, 3, (3,3)) @cloudfrac_left_iter() -@handle_casting(arg_idxs=(0, 1), outviews=("lowview", "medview", "highview")) -@handle_extract_transpose(outviews=("lowview", "medview", "hightview")) +@cast_type(arg_idxs=(0, 1), outviews=("lowview", "medview", "highview")) +@extract_and_transpose(outviews=("lowview", "medview", "hightview")) def _cloudfrac(p, rh, lowview=None, medview=None, highview=None): if lowview is None: @@ -863,55 +512,13 @@ def _cloudfrac(p, rh, lowview=None, medview=None, highview=None): return result -#def computeij(map_proj, truelat1, truelat2, stdlon, -# lat1, lon1, pole_lat, pole_lon, -# knowni, knownj, dx, latinc, loninc, lat, lon): -# -# result = f_lltoij(map_proj, -# truelat1, -# truelat2, -# stdlon, -# lat1, -# lon1, -# pole_lat, -# pole_lon, -# knowni, -# knownj, -# dx, -# latinc, -# loninc, -# lat, -# lon, -# FortranError()) -# -# return result - -#def computell(map_proj, truelat1, truelat2, stdlon, lat1, lon1, -# pole_lat, pole_lon, knowni, knownj, dx, latinc, -# loninc, i, j): -# -# result = f_ijtoll(map_proj, -# truelat1, -# truelat2, -# stdlon, -# lat1, -# lon1, -# pole_lat, -# pole_lon, -# knowni, -# knownj, -# dx, -# latinc, -# loninc, -# i, -# j, -# FortranError()) -# -# return result - def _lltoxy(map_proj, truelat1, truelat2, stdlon, lat1, lon1, pole_lat, pole_lon, - known_x, known_y, dx, dy, latinc, loninc, lat, lon): + known_x, known_y, dx, dy, latinc, loninc, lat, lon, + outview=None): + + if outview is None: + outview = np.zeros((2), dtype=np.float64, order="F") errstat = np.array(0) errmsg = np.zeros(Constants.ERRLEN, "c") @@ -932,6 +539,7 @@ def _lltoxy(map_proj, truelat1, truelat2, stdlon, loninc, lat, lon, + outview, errstat, errmsg) @@ -943,7 +551,10 @@ def _lltoxy(map_proj, truelat1, truelat2, stdlon, def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1, pole_lat, pole_lon, known_x, known_y, dx, dy, latinc, - loninc, x, y): + loninc, x, y, outview=None): + + if outview is None: + outview = np.zeros((2), dtype=np.float64, order="F") errstat = np.array(0) errmsg = np.zeros(Constants.ERRLEN, "c") @@ -964,6 +575,7 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1, loninc, x, y, + outview, errstat, errmsg) @@ -973,43 +585,10 @@ def _xytoll(map_proj, truelat1, truelat2, stdlon, lat1, lon1, return result -#@handle_left_iter(3,0, ignore_args=(3,)) -#@handle_casting(arg_idxs=(0,1,2)) -#@handle_extract_transpose() -#def computeeta(full_t, znu, psfc, ptop): -# pcalc = np.zeros(full_t.shape, np.float64, order="F") -# mean_t = np.zeros(full_t.shape, np.float64, order="F") -# temp_t = np.zeros(full_t.shape, np.float64, order="F") -# -# result = f_converteta(full_t, -# znu, -# psfc, -# ptop, -# pcalc, -# mean_t, -# temp_t) -# -# return result - -#@handle_left_iter(3,0,ignore_args=(7,)) -#@handle_casting(arg_idxs=(0,1,2,3,4,5,6)) -#@handle_extract_transpose() -#def computectt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci): -# result = f_computectt(p_hpa, -# tk, -# qice, -# qcld, -# qv, -# ght, -# ter, -# haveqci) -# -# return result - @check_args(0, 3, (3,3,3,3,3,3,2)) -@left_iter_nocopy(3, 2, ref_var_idx=0, ignore_args=(7,)) -@handle_casting(arg_idxs=(0,1,2,3,4,5,6)) -@handle_extract_transpose() +@left_iteration(3, 2, ref_var_idx=0, ignore_args=(7,)) +@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): if outview is None: outview = np.empty_like(ter) @@ -1027,34 +606,10 @@ def _ctt(p_hpa, tk, qice, qcld, qv, ght, ter, haveqci, outview=None): return result -# @handle_left_iter(2,0,ignore_args=(1,)) -# @handle_casting(arg_idxs=(0,)) -# @handle_extract_transpose() -# def smooth2d(field, passes): -# # Unlike NCL, this routine will not modify the values in place, but -# # copies the original data before modifying it. This allows the decorator -# # to work properly and also behaves like the other methods. -# -# if isinstance(field, np.ma.MaskedArray): -# missing = field.fill_value -# else: -# missing = Constants.DEFAULT_FILL -# -# field_copy = field.copy(order="A") -# field_tmp = np.zeros(field_copy.shape, field_copy.dtype, order="F") -# -# f_filter2d(field_copy, -# field_tmp, -# missing, -# passes) -# -# # Don't transpose here since the fortran routine is not returning an -# # array. It's only modifying the existing array. -# return field_copy - -@left_iter_nocopy(2, 2, ref_var_idx=0, ignore_args=(1,)) -@handle_casting(arg_idxs=(0,)) -@handle_extract_transpose() +@check_args(0, 2, (2,)) +@left_iteration(2, 2, ref_var_idx=0, ignore_args=(1,)) +@cast_type(arg_idxs=(0,)) +@extract_and_transpose() def _smooth2d(field, passes, outview=None): # Unlike NCL, this routine will not modify the values in place, but # copies the original data before modifying it. @@ -1076,28 +631,13 @@ def _smooth2d(field, passes, outview=None): passes, missing) - # Don't transpose here since the fortran routine is not returning an - # array. It's only modifying the existing array. return outview - -#@handle_left_iter(3,0,ignore_args=(3,4,5)) -#@handle_casting(arg_idxs=(0,1,2)) -#@handle_extract_transpose() -#def monotonic(var, lvprs, coriolis, idir, delta, icorsw): -# result = f_monotonic(var, -# lvprs, -# coriolis, -# idir, -# delta, -# icorsw) -# -# return result @check_args(0, 3, (3,3,2)) -@left_iter_nocopy(3, 3, ref_var_idx=0, ignore_args=(3,4,5)) -@handle_casting(arg_idxs=(0,1,2)) -@handle_extract_transpose() +@left_iteration(3, 3, ref_var_idx=0, ignore_args=(3,4,5)) +@cast_type(arg_idxs=(0,1,2)) +@extract_and_transpose() def _monotonic(var, lvprs, coriolis, idir, delta, icorsw, outview=None): # If icorsw is not 0, then the input variable might get modified by the # fortran routine. We don't want this, so make a copy and pass that on. @@ -1116,39 +656,14 @@ def _monotonic(var, lvprs, coriolis, idir, delta, icorsw, outview=None): return result -#@handle_left_iter(3,0,ignore_args=(9,10,11,12,13,14)) -#@handle_casting(arg_idxs=(0,1,2,3,4,5,6,7,8,9)) -#@handle_extract_transpose() -#def vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, -# vcarray, interp_levels, icase, extrap, vcor, logp, -# missing): -# -# result = f_vintrp(field, -# pres, -# tk, -# qvp, -# ght, -# terrain, -# sfp, -# smsfp, -# vcarray, -# interp_levels, -# icase, -# extrap, -# vcor, -# logp, -# missing) -# -# return result - # Output shape is interp_levels.shape + field.shape[-2:] @check_args(0, 3, (3,3,3,3,3,2,2,2,3)) -@left_iter_nocopy(3, combine_dims([(9, (-1,)), +@left_iteration(3, combine_dims([(9, (-1,)), (0, (-2,-1))]), ref_var_idx=0, ignore_args=(9,10,11,12,13,14)) -@handle_casting(arg_idxs=(0,1,2,3,4,5,6,7,8,9)) -@handle_extract_transpose() +@cast_type(arg_idxs=(0,1,2,3,4,5,6,7,8,9)) +@extract_and_transpose() def _vintrp(field, pres, tk, qvp, ght, terrain, sfp, smsfp, vcarray, interp_levels, icase, extrap, vcor, logp, missing, outview=None): diff --git a/src/wrf/helicity.py b/src/wrf/helicity.py index dcb3dff..70ad456 100755 --- a/src/wrf/helicity.py +++ b/src/wrf/helicity.py @@ -1,8 +1,9 @@ from __future__ import (absolute_import, division, print_function, unicode_literals) -from .constants import Constants +import numpy as np +from .constants import Constants from .extension import _srhel, _udhel from .destag import destagger from .util import extract_vars, extract_global_attrs, either @@ -12,8 +13,7 @@ from .metadecorators import copy_and_set_metadata description="storm relative helicity", units="m-2/s-2") def get_srh(wrfnc, timeidx=0, method="cat", squeeze=True, - cache=None, meta=True, - top=3000.0): + cache=None, meta=True, top=3000.0): # Top can either be 3000 or 1000 (for 0-1 srh or 0-3 srh) ncvars = extract_vars(wrfnc, timeidx, ("HGT", "PH", "PHB"), @@ -40,9 +40,9 @@ def get_srh(wrfnc, timeidx=0, method="cat", squeeze=True, z = geopt_unstag / Constants.G # Re-ordering from high to low - u1 = u[...,::-1,:,:] - v1 = v[...,::-1,:,:] - z1 = z[...,::-1,:,:] + u1 = np.ascontiguousarray(u[...,::-1,:,:]) + v1 = np.ascontiguousarray(v[...,::-1,:,:]) + z1 = np.ascontiguousarray(z[...,::-1,:,:]) srh = _srhel(u1, v1, z1, ter, top) @@ -56,7 +56,7 @@ def get_uh(wrfnc, timeidx=0, method="cat", squeeze=True, bottom=2000.0, top=5000.0): ncvars = extract_vars(wrfnc, timeidx, ("W", "PH", "PHB", "MAPFAC_M"), - method, squeeze, cache) + method, squeeze, cache, meta=False) wstag = ncvars["W"] ph = ncvars["PH"] diff --git a/src/wrf/interp.py b/src/wrf/interp.py index 13abb1e..0391c70 100755 --- a/src/wrf/interp.py +++ b/src/wrf/interp.py @@ -67,8 +67,7 @@ def interplevel(field3d, z, desiredlev, missing=Constants.DEFAULT_FILL, def vertcross(field3d, z, missing=Constants.DEFAULT_FILL, pivot_point=None, angle=None, start_point=None, end_point=None, - include_latlon=False, - cache=None, meta=True): + latlon=False, cache=None, meta=True): """Return the vertical cross section for a 3D field, interpolated to a vertical plane defined by a horizontal line. @@ -105,7 +104,7 @@ def vertcross(field3d, z, missing=Constants.DEFAULT_FILL, (or [west_east, south_north]), which indicates the end x,y location through which the plane will pass. - include_latlon : {True, False}, optional + latlon : {True, False}, optional Set to True to also interpolate the two-dimensional latitude and longitude coordinates along the same horizontal line and include this information in the metadata (if enabled). This can be @@ -147,7 +146,7 @@ def vertcross(field3d, z, missing=Constants.DEFAULT_FILL, @set_interp_metadata("line") def interpline(field2d, pivot_point=None, angle=None, start_point=None, - end_point=None, include_latlon=False, + end_point=None, latlon=False, cache=None, meta=True): """Return the two-dimensional field interpolated along a line. @@ -176,7 +175,7 @@ def interpline(field2d, pivot_point=None, (or [west_east, south_north]), which indicates the end x,y location through which the plane will pass. - include_latlon : {True, False}, optional + latlon : {True, False}, optional Set to True to also interpolate the two-dimensional latitude and longitude coordinates along the same horizontal line and include this information in the metadata (if enabled). This can be diff --git a/src/wrf/interputils.py b/src/wrf/interputils.py index a9b9e56..1ed89cd 100644 --- a/src/wrf/interputils.py +++ b/src/wrf/interputils.py @@ -138,18 +138,18 @@ def get_xy_z_params(z, pivot_point=None, angle=None, # interp to constant z grid if(var2dz[idx1] > var2dz[idx2]): # monotonically decreasing coordinate - z_max = floor(np.amax(z) / 10) * 10 # bottom value - z_min = ceil(np.amin(z) / 10) * 10 # top value + z_max = floor(np.amax(z)/10) * 10 # bottom value + z_min = ceil(np.amin(z)/10) * 10 # top value dz = 10 - nlevels = int((z_max-z_min) / dz) + nlevels = int((z_max - z_min)/dz) z_var2d = np.zeros((nlevels), dtype=z.dtype) z_var2d[0] = z_max dz = -dz else: z_max = np.amax(z) z_min = 0. - dz = 0.01 * z_max - nlevels = int(z_max / dz) + dz = 0.01*z_max + nlevels = int(z_max/dz) z_var2d = np.zeros((nlevels), dtype=z.dtype) z_var2d[0] = z_min diff --git a/src/wrf/latlon.py b/src/wrf/latlon.py index 641936e..da231d1 100755 --- a/src/wrf/latlon.py +++ b/src/wrf/latlon.py @@ -31,19 +31,28 @@ def get_lon(wrfnc, timeidx=0, method="cat", squeeze=True, # Can either use wrfnc as a single file or sequence, or provide # projection parameters (which don't allow for moving domains) @set_latlon_metadata(xy=True) -def ll_to_xy(latitude, longitude, wrfnc=None, timeidx=0, - stagger=None, method="cat", squeeze=True, cache=None, meta=True, - **projparms): - return _ll_to_xy(latitude, longitude, wrfnc, timeidx, stagger, - method, squeeze, cache, **projparams) +def ll_to_xy(wrfnc, latitude, longitude, timeidx=0, stagger=None, method="cat", + squeeze=True, cache=None, meta=True): + return _ll_to_xy(latitude, longitude, wrfnc, timeidx, stagger, method, + squeeze, cache, **{}) + + +@set_latlon_metadata(xy=True) +def ll_to_xy_proj(latitude, longitude, meta=True, squeeze=True, **projparams): + return _ll_to_xy(latitude, longitude, None, 0, squeeze, "cat", True, None, + **projparams) @set_latlon_metadata(xy=False) -def xy_to_ll(x, y, wrfnc=None, timeidx=0, - stagger=None, method="cat", squeeze=True, cache=None, meta=True, - **projparms): - return _xy_to_ll(x, y, wrfnc, timeidx, stagger, - method, squeeze, cache, **projparams) +def xy_to_ll(wrfnc, x, y, timeidx=0, stagger=None, method="cat", squeeze=True, + cache=None, meta=True): + return _xy_to_ll(x, y, wrfnc, timeidx, stagger, method, squeeze, cache, + **{}) + + +@set_latlon_metadata(xy=False) +def xy_to_ll_proj(x, y, meta=True, squeeze=True, **projparams): + return _xy_to_ll(x, y, None, 0, None, "cat", squeeze, None, **projparams) \ No newline at end of file diff --git a/src/wrf/latlonutils.py b/src/wrf/latlonutils.py index febf60e..583c965 100644 --- a/src/wrf/latlonutils.py +++ b/src/wrf/latlonutils.py @@ -5,7 +5,7 @@ from collections import Iterable import numpy as np -from .constants import Constants +from .constants import Constants, ProjectionTypes from .extension import _lltoxy, _xytoll from .util import (extract_vars, extract_global_attrs, either, is_moving_domain, is_multi_time_req, @@ -46,7 +46,7 @@ def _get_proj_params(wrfnc, timeidx, stagger, method, squeeze, cache): dx = attrs["DX"] dy = attrs["DY"] - if map_proj == 6: + if map_proj == ProjectionTypes.LAT_LON: pole_attrs = extract_global_attrs(wrfnc, attrs=("POLE_LAT", "POLE_LON")) pole_lat = pole_attrs["POLE_LAT"] @@ -133,7 +133,9 @@ def _kwarg_proj_params(projparams): known_x = known_x + 1 known_y = known_y + 1 - if map_proj in (1, 2, 3): + if map_proj in (ProjectionTypes.LAMBERT_CONFORMAL, + ProjectionTypes.POLAR_STEREOGRAPHIC, + ProjectionTypes.MERCATOR): if truelat1 is None: raise ValueError("'TRUELAT1' argument required") else: @@ -141,7 +143,7 @@ def _kwarg_proj_params(projparams): truelat1 = 0.0 # Map projection 6 (lat lon) required latinc, loninc, and dy - if map_proj == 6: + if map_proj == ProjectionTypes.LAT_LON: if latinc is None: raise ValueError("'LATINC' argument required") @@ -227,9 +229,6 @@ def _ll_to_xy(latitude, longitude, wrfnc=None, timeidx=0, # Make indexes 0-based result = result - 1 - - if squeeze: - result = result.squeeze() return result @@ -303,8 +302,5 @@ def _xy_to_ll(x, y, wrfnc=None, timeidx=0, stagger=None, result = _xytoll(map_proj, truelat1, truelat2, stdlon, ref_lat, ref_lon, pole_lat, pole_lon, known_x, known_y, dx, dy, latinc, loninc, x_val, y_val) - - if squeeze: - result = result.squeeze() return result diff --git a/src/wrf/metadecorators.py b/src/wrf/metadecorators.py index 29f4a2e..4fe0acc 100644 --- a/src/wrf/metadecorators.py +++ b/src/wrf/metadecorators.py @@ -423,11 +423,11 @@ def set_latlon_metadata(xy=False): outname = "latlon" if not xy else "xy" if result.ndim == 2: - dimnames = (["xy", "lat_lon"] if not xy - else ["latlon", "x_y"]) + dimnames = (["idx", "lat_lon"] if not xy + else ["idx", "x_y"]) else: - dimnames = (["xy", "domain", "lat_lon"] if not xy - else ["latlon", "domain", "x_y"]) + dimnames = (["idx", "domain", "lat_lon"] if not xy + else ["idx", "domain", "x_y"]) argvars = from_args(wrapped, argnames, *args, **kwargs) @@ -439,12 +439,12 @@ def set_latlon_metadata(xy=False): arr2 = np.asarray(var2).ravel() coords = {} - if not ij: - coords["coord_pair"] = (dimnames[0], [CoordPair(x=x[0], y=x[1]) + if not xy: + coords["xy_coord"] = (dimnames[0], [CoordPair(x=x[0], y=x[1]) for x in zip(arr1, arr2)]) coords[dimnames[-1]] = ["lat", "lon"] else: - coords["coord_pair"] = (dimnames[0], [CoordPair(lat=x[0], lon=x[1]) + coords["latlon_coord"] = (dimnames[0], [CoordPair(lat=x[0], lon=x[1]) for x in zip(arr1, arr2)]) coords[dimnames[-1]] = ["x", "y"] @@ -581,7 +581,7 @@ def _set_horiz_meta(wrapped, instance, args, kwargs): coords=outcoords, attrs=outattrs) def _set_cross_meta(wrapped, instance, args, kwargs): - argvars = from_args(wrapped, ("field3d", "z", "include_latlon", "missing", + argvars = from_args(wrapped, ("field3d", "z", "latlon", "missing", "pivot_point", "angle", "start_point", "end_point", "cache"), @@ -589,7 +589,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): field3d = argvars["field3d"] z = argvars["z"] - inc_latlon = argvars["include_latlon"] + inc_latlon = argvars["latlon"] missingval = argvars["missing"] pivot_point = argvars["pivot_point"] angle = argvars["angle"] @@ -650,7 +650,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): del outcoords[key] outdimnames.append("vertical") - outdimnames.append("xy") + outdimnames.append("idx") outattrs.update(field3d.attrs) outname = "{0}_cross".format(field3d.name) @@ -673,13 +673,13 @@ def _set_cross_meta(wrapped, instance, args, kwargs): lats = _interpline(latcoord, xy) lons = _interpline(loncoord, xy) - outcoords["xy_loc"] = ("xy", + outcoords["xy_loc"] = ("idx", np.asarray(tuple( CoordPair(x=xy[i,0], y=xy[i,1], lat=lats[i], lon=lons[i]) for i in py3range(xy.shape[-2]))) ) - + # Moving domain else: extra_dims = latcoord.shape[0:-2] outdims = extra_dims + xy.shape[-2:-1] @@ -698,16 +698,16 @@ def _set_cross_meta(wrapped, instance, args, kwargs): extra_dimnames = latcoord.dims[0:-2] - loc_dimnames = extra_dimnames + ("xy",) + loc_dimnames = extra_dimnames + ("idx",) outcoords["xy_loc"] = (loc_dimnames, latlon_loc) else: - outcoords["xy_loc"] = ("xy", np.asarray(tuple( + outcoords["xy_loc"] = ("idx", np.asarray(tuple( CoordPair(xy[i,0], xy[i,1]) for i in py3range(xy.shape[-2])))) else: - outcoords["xy_loc"] = ("xy", np.asarray(tuple( + outcoords["xy_loc"] = ("idx", np.asarray(tuple( CoordPair(xy[i,0], xy[i,1]) for i in py3range(xy.shape[-2])))) @@ -728,7 +728,7 @@ def _set_cross_meta(wrapped, instance, args, kwargs): def _set_line_meta(wrapped, instance, args, kwargs): argvars = from_args(wrapped, ("field2d", "pivot_point", "angle", - "start_point", "end_point", "include_latlon", + "start_point", "end_point", "latlon", "cache"), *args, **kwargs) @@ -737,7 +737,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): angle = argvars["angle"] start_point = argvars["start_point"] end_point = argvars["end_point"] - inc_latlon = argvars["include_latlon"] + inc_latlon = argvars["latlon"] cache = argvars["cache"] if cache is None: @@ -788,7 +788,7 @@ def _set_line_meta(wrapped, instance, args, kwargs): for key in delkeys: del outcoords[key] - outdimnames.append("xy") + outdimnames.append("line_idx") outattrs.update(field2d.attrs) outname = "{0}_line".format(field2d.name) @@ -811,12 +811,14 @@ def _set_line_meta(wrapped, instance, args, kwargs): lats = _interpline(latcoord, xy) lons = _interpline(loncoord, xy) - outcoords["xy_loc"] = np.asarray(tuple( + outcoords["xy_loc"] = ("line_idx", + np.asarray(tuple( CoordPair(x=xy[i,0], y=xy[i,1], lat=lats[i], lon=lons[i]) - for i in py3range(xy.shape[-2]))) + for i in py3range(xy.shape[-2]))) + ) - + # Moving domain else: extra_dims = latcoord.shape[0:-2] outdims = extra_dims + xy.shape[-2:-1] @@ -835,16 +837,16 @@ def _set_line_meta(wrapped, instance, args, kwargs): extra_dimnames = latcoord.dims[0:-2] - loc_dimnames = extra_dimnames + ("xy",) + loc_dimnames = extra_dimnames + ("line_idx",) outcoords["xy_loc"] = (loc_dimnames, latlon_loc) else: - outcoords["xy_loc"] = ("xy", np.asarray(tuple( + outcoords["xy_loc"] = ("line_idx", np.asarray(tuple( CoordPair(xy[i,0], xy[i,1]) for i in py3range(xy.shape[-2])))) else: - outcoords["xy_loc"] = ("xy", np.asarray(tuple( + outcoords["xy_loc"] = ("line_idx", np.asarray(tuple( CoordPair(xy[i,0], xy[i,1]) for i in py3range(xy.shape[-2])))) @@ -905,8 +907,7 @@ def _set_vinterp_meta(wrapped, instance, args, kwargs): def _set_2dxy_meta(wrapped, instance, args, kwargs): - argvars = from_args(wrapped, ("field3d", "xy"), - *args, **kwargs) + argvars = from_args(wrapped, ("field3d", "xy"), *args, **kwargs) field3d = argvars["field3d"] xy = argvars["xy"] @@ -941,7 +942,7 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs): for key in delkeys: del outcoords[key] - outdimnames.append("xy") + outdimnames.append("line_idx") #outattrs.update(field3d.attrs) desc = field3d.attrs.get("description", None) @@ -954,7 +955,7 @@ def _set_2dxy_meta(wrapped, instance, args, kwargs): outname = "{0}_2dxy".format(field3d.name) - outcoords["xy_loc"] = ("xy", [CoordPair(xy[i,0], xy[i,1]) + outcoords["xy_loc"] = ("line_idx", [CoordPair(xy[i,0], xy[i,1]) for i in py3range(xy.shape[-2])]) for key in ("MemoryOrder",): @@ -1040,7 +1041,7 @@ def _set_xy_meta(wrapped, instance, args, kwargs): else: outname = "xy" - outdimnames = ["idx", "x_y"] + outdimnames = ["line_idx", "x_y"] outcoords = OrderedDict() outattrs = OrderedDict() @@ -1098,7 +1099,8 @@ def set_alg_metadata(alg_ndims, refvarname, refvarname: argument name for the reference variable missingarg: argument name for the missing value stagdim: staggered dimension in reference - stagsubvar: the variable to use to supply the staggered dimension name + stagsubvar: the variable name to use to supply the unstaggered + dimension size """ @@ -1170,7 +1172,8 @@ def set_alg_metadata(alg_ndims, refvarname, outdims[-alg_ndims:] = refvar.dims[-alg_ndims:] # Use the stagsubvar if applicable - outdims[stagdim] = stagvar.dims[stagdim] + if stagvar is not None and stagdim is not None: + outdims[stagdim] = stagvar.dims[stagdim] # Left dims if refvarndims is None: @@ -1191,11 +1194,10 @@ def set_alg_metadata(alg_ndims, refvarname, ref_left_dimnames = refvar.dims[0:ref_extra] for i,dimname in enumerate(ref_left_dimnames[::-1], 1): - idx = -i - if -idx <= result.shape: - outdims[idx] = dimname + if i <= result.ndim: + outdims[-alg_ndims - i] = dimname else: - continute + continue out = DataArray(result, name=outname, dims=outdims, attrs=outattrs) @@ -1206,7 +1208,8 @@ def set_alg_metadata(alg_ndims, refvarname, return func_wrapper -def set_uvmet_alg_metadata(units="mps", description="earth rotated u,v"): +def set_uvmet_alg_metadata(units="mps", description="earth rotated u,v", + latarg="lat", windarg="u"): @wrapt.decorator def func_wrapper(wrapped, instance, args, kwargs): @@ -1232,8 +1235,8 @@ def set_uvmet_alg_metadata(units="mps", description="earth rotated u,v"): if description is not None: outattrs["description"] = description - latvar = from_args(wrapped, ("lat",), *args, **kwargs)["lat"] - uvar = from_args(wrapped, ("u",), *args, **kwargs)["u"] + latvar = from_args(wrapped, latarg, *args, **kwargs)[latarg] + uvar = from_args(wrapped, windarg, *args, **kwargs)[windarg] if isinstance(uvar, DataArray): # Right dims come from latvar @@ -1254,7 +1257,7 @@ def set_uvmet_alg_metadata(units="mps", description="earth rotated u,v"): return func_wrapper -def set_cape_alg_metadata(is2d): +def set_cape_alg_metadata(is2d, copyarg="pres_hpa"): @wrapt.decorator def func_wrapper(wrapped, instance, args, kwargs): @@ -1285,8 +1288,8 @@ def set_cape_alg_metadata(is2d): outattrs["MemoryOrder"] = "XYZ" - argvals = from_args(wrapped, ("p_hpa","missing"), *args, **kwargs) - p = argvals["p_hpa"] + argvals = from_args(wrapped, (copyarg,"missing"), *args, **kwargs) + p = argvals[copyarg] missing = argvals["missing"] if isinstance(p, DataArray): @@ -1322,7 +1325,7 @@ def set_cape_alg_metadata(is2d): return func_wrapper -def set_cloudfrac_alg_metadata(): +def set_cloudfrac_alg_metadata(copyarg="pres"): @wrapt.decorator def func_wrapper(wrapped, instance, args, kwargs): @@ -1347,7 +1350,7 @@ def set_cloudfrac_alg_metadata(): outattrs["MemoryOrder"] = "XY" - argvals = from_args(wrapped, "p", *args, **kwargs)["p"] + p = from_args(wrapped, copyarg, *args, **kwargs)[copyarg] if isinstance(p, DataArray): # Right dims diff --git a/src/wrf/projection.py b/src/wrf/projection.py index 85c76f8..726d362 100644 --- a/src/wrf/projection.py +++ b/src/wrf/projection.py @@ -4,7 +4,7 @@ import numpy as np import math from .config import basemap_enabled, cartopy_enabled, pyngl_enabled -from .constants import Constants +from .constants import Constants, ProjectionTypes if cartopy_enabled(): from cartopy import crs @@ -719,16 +719,17 @@ def getproj(bottom_left=None, top_right=None, lats=None, lons=None, **proj_params): proj_type = proj_params.get("MAP_PROJ", 0) - if proj_type == 1: + if proj_type == ProjectionTypes.LAMBERT_CONFORMAL: return LambertConformal(bottom_left, top_right, lats, lons, **proj_params) - elif proj_type == 2: + elif proj_type == ProjectionTypes.POLAR_STEREOGRAPHIC: return PolarStereographic(bottom_left, top_right, lats, lons, **proj_params) - elif proj_type == 3: + elif proj_type == ProjectionTypes.MERCATOR: return Mercator(bottom_left, top_right, lats, lons, **proj_params) - elif proj_type == 0 or proj_type == 6: + elif (proj_type == ProjectionTypes.ZERO or + proj_type == ProjectionTypes.LAT_LON): if (proj_params.get("POLE_LAT", None) == 90. and proj_params.get("POLE_LON", None) == 0.): return LatLon(bottom_left, top_right, diff --git a/src/wrf/pw.py b/src/wrf/pw.py index b4e3e80..e5b4c34 100755 --- a/src/wrf/pw.py +++ b/src/wrf/pw.py @@ -26,10 +26,9 @@ def get_pw(wrfnc, timeidx=0, method="cat", squeeze=True, cache=None, phb = ncvars["PHB"] qv = ncvars["QVAPOR"] - # Change this to use real virtual temperature! - full_p = p + pb + full_p = p + pb ht = (ph + phb)/Constants.G - full_t = t + Constants.T_BASE + full_t = t + Constants.T_BASE tk = _tk(full_p, full_t) tv = _tv(tk, qv) diff --git a/src/wrf/util.py b/src/wrf/util.py index 4a70020..4e3eca5 100644 --- a/src/wrf/util.py +++ b/src/wrf/util.py @@ -656,7 +656,9 @@ def _find_forward(wrfseq, varname, timeidx, is_moving, meta): return _build_data_array(wrfnc, varname, filetimeidx, is_moving) else: - return wrfnc.variables[varname][filetimeidx, :] + result = wrfnc.variables[varname][filetimeidx, :] + return result[np.newaxis, :] # So that nosqueee works + else: comboidx += numtimes @@ -692,7 +694,8 @@ def _find_reverse(wrfseq, varname, timeidx, is_moving, meta): return _build_data_array(wrfnc, varname, filetimeidx, is_moving) else: - return wrfnc.variables[varname][filetimeidx, :] + result = wrfnc.variables[varname][filetimeidx, :] + return result[np.newaxis, :] # So that nosqueeze works else: comboidx += numtimes @@ -1090,7 +1093,7 @@ def _extract_var(wrfnc, varname, timeidx, is_moving, # Squeeze handled in this routine, so just return it return combine_files(wrfnc, varname, timeidx, is_moving, method, squeeze, meta) - + return result.squeeze() if squeeze else result diff --git a/src/wrf/uvmet.py b/src/wrf/uvmet.py index 301853c..e14b9ee 100755 --- a/src/wrf/uvmet.py +++ b/src/wrf/uvmet.py @@ -63,17 +63,17 @@ def _get_uvmet(wrfnc, timeidx=0, method="cat", squeeze=True, resdim = (2,) + u.shape[0:end_idx] + u.shape[end_idx:] # Make a new output array for the result - res = np.empty(resdim, u.dtype) + result = np.empty(resdim, u.dtype) # For 2D array, this makes (0,...,:,:) and (1,...,:,:) # For 3D array, this makes (0,...,:,:,:) and (1,...,:,:,:) idx0 = (0,) + (Ellipsis,) + (slice(None),)*(-end_idx) idx1 = (1,) + (Ellipsis,) + (slice(None),)*(-end_idx) - res[idx0] = u[:] - res[idx1] = v[:] + result[idx0] = u[:] + result[idx1] = v[:] - return res + return result elif map_proj in (1,2): lat_attrs = extract_global_attrs(wrfnc, attrs=("TRUELAT1", "TRUELAT2")) @@ -120,7 +120,10 @@ def _get_uvmet(wrfnc, timeidx=0, method="cat", squeeze=True, result = _uvmet(u, v, lat, lon, cen_lon, cone) - return result.squeeze() + if squeeze: + result = result.squeeze() + + return result @set_wind_metadata(copy_varname=either("P", "PRES"), diff --git a/test/comp_utest.py b/test/comp_utest.py new file mode 100644 index 0000000..2ba98bf --- /dev/null +++ b/test/comp_utest.py @@ -0,0 +1,580 @@ +from math import fabs, log, tan, sin, cos + +import unittest as ut +import numpy.testing as nt +import numpy as np +import numpy.ma as ma +import os, sys +import subprocess + +from netCDF4 import Dataset as nc + +from wrf import * + +NCL_EXE = "/Users/ladwig/nclbuild/6.3.0/bin/ncl" +TEST_FILE = "/Users/ladwig/Documents/wrf_files/wrfout_d01_2010-06-13_21:00:00" +OUT_NC_FILE = "/tmp/wrftest.nc" +NCFILE = nc(TEST_FILE) +NCGROUP = [NCFILE, NCFILE, NCFILE] + +# Python 3 +if sys.version_info > (3,): + xrange = range + + +ROUTINE_MAP = {"avo" : avo, + "eth" : eth, + "cape_2d" : cape_2d, + "cape_3d" : cape_3d, + "ctt" : ctt, + "dbz" : dbz, + "helicity" : srhel, + "omg" : omega, + "pvo" : pvo, + "pw" : pw, + "rh" : rh, + "slp" : slp, + "td" : td, + "tk" : tk, + "tv" : tvirtual, + "twb" : wetbulb, + "updraft_helicity" : udhel, + "uvmet" : uvmet, + "cloudfrac" : cloudfrac} + +class ProjectionError(RuntimeError): + pass + +def get_args(varname, wrfnc, timeidx, method, squeeze): + if varname == "avo": + ncvars = extract_vars(wrfnc, timeidx, ("U", "V", "MAPFAC_U", + "MAPFAC_V", "MAPFAC_M", + "F"), + method, squeeze, cache=None, meta=True) + + attrs = extract_global_attrs(wrfnc, attrs=("DX", "DY")) + u = ncvars["U"] + v = ncvars["V"] + msfu = ncvars["MAPFAC_U"] + msfv = ncvars["MAPFAC_V"] + msfm = ncvars["MAPFAC_M"] + cor = ncvars["F"] + + dx = attrs["DX"] + dy = attrs["DY"] + + return (u, v, msfu, msfv, msfm, cor, dx, dy) + + if varname == "pvo": + ncvars = extract_vars(wrfnc, timeidx, ("U", "V", "T", "P", + "PB", "MAPFAC_U", + "MAPFAC_V", "MAPFAC_M", + "F"), + method, squeeze, cache=None, meta=True) + attrs = extract_global_attrs(wrfnc, attrs=("DX", "DY")) + + u = ncvars["U"] + v = ncvars["V"] + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + msfu = ncvars["MAPFAC_U"] + msfv = ncvars["MAPFAC_V"] + msfm = ncvars["MAPFAC_M"] + cor = ncvars["F"] + + dx = attrs["DX"] + dy = attrs["DY"] + + full_t = t + 300 + full_p = p + pb + + return (u, v, full_t, full_p, msfu, msfv, msfm, cor, dx, dy) + + if varname == "eth": + varnames=("T", "P", "PB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + return (qv, tkel, full_p) + + if varname == "cape_2d": + varnames = ("T", "P", "PB", "QVAPOR", "PH","PHB", "HGT", "PSFC") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + ter = ncvars["HGT"] + psfc = ncvars["PSFC"] + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + geopt = ph + phb + geopt_unstag = destagger(geopt, -3) + z = geopt_unstag/Constants.G + + # Convert pressure to hPa + p_hpa = ConversionFactors.PA_TO_HPA * full_p + psfc_hpa = ConversionFactors.PA_TO_HPA * psfc + + i3dflag = 0 + ter_follow = 1 + + return (p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow) + + if varname == "cape_3d": + varnames = ("T", "P", "PB", "QVAPOR", "PH", "PHB", "HGT", "PSFC") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + ter = ncvars["HGT"] + psfc = ncvars["PSFC"] + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + geopt = ph + phb + geopt_unstag = destagger(geopt, -3) + z = geopt_unstag/Constants.G + + # Convert pressure to hPa + p_hpa = ConversionFactors.PA_TO_HPA * full_p + psfc_hpa = ConversionFactors.PA_TO_HPA * psfc + + i3dflag = 1 + ter_follow = 1 + + return (p_hpa, tkel, qv, z, ter, psfc_hpa, ter_follow) + + if varname == "ctt": + varnames = ("T", "P", "PB", "PH", "PHB", "HGT", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + ter = ncvars["HGT"] + qv = ncvars["QVAPOR"] * 1000.0 # g/kg + + haveqci = 1 + try: + icevars = extract_vars(wrfnc, timeidx, "QICE", + method, squeeze, cache=None, meta=False) + except KeyError: + qice = np.zeros(qv.shape, qv.dtype) + haveqci = 0 + else: + qice = icevars["QICE"] * 1000.0 #g/kg + + try: + cldvars = extract_vars(wrfnc, timeidx, "QCLOUD", + method, squeeze, cache=None, meta=False) + except KeyError: + raise RuntimeError("'QCLOUD' not found in NetCDF file") + else: + qcld = cldvars["QCLOUD"] * 1000.0 #g/kg + + full_p = p + pb + p_hpa = full_p * ConversionFactors.PA_TO_HPA + full_t = t + Constants.T_BASE + tkel = tk(full_p, full_t, meta=False) + + geopt = ph + phb + geopt_unstag = destagger(geopt, -3) + ght = geopt_unstag / Constants.G + + return (p_hpa, tkel, qv, qcld, ght, ter, qice) + + if varname == "dbz": + varnames = ("T", "P", "PB", "QVAPOR", "QRAIN") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + qr = ncvars["QRAIN"] + + try: + snowvars = extract_vars(wrfnc, timeidx, "QSNOW", + method, squeeze, cache=None, meta=False) + except KeyError: + qs = np.zeros(qv.shape, qv.dtype) + else: + qs = snowvars["QSNOW"] + + try: + graupvars = extract_vars(wrfnc, timeidx, "QGRAUP", + method, squeeze, cache=None, meta=False) + except KeyError: + qg = np.zeros(qv.shape, qv.dtype) + else: + qg = graupvars["QGRAUP"] + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + return (full_p, tkel, qv, qr, qs, qg) + + if varname == "helicity": + # Top can either be 3000 or 1000 (for 0-1 srh or 0-3 srh) + + ncvars = extract_vars(wrfnc, timeidx, ("HGT", "PH", "PHB"), + method, squeeze, cache=None, meta=True) + + ter = ncvars["HGT"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + + # As coded in NCL, but not sure this is possible + varname = "U" + u_vars = extract_vars(wrfnc, timeidx, varname, method, squeeze, + cache=None, meta=False) + u = destagger(u_vars[varname], -1) + + varname = "V" + v_vars = extract_vars(wrfnc, timeidx, varname, method, squeeze, + cache=None, meta=False) + v = destagger(v_vars[varname], -2) + + geopt = ph + phb + geopt_unstag = destagger(geopt, -3) + + z = geopt_unstag / Constants.G + + return (u, v, z, ter) + + if varname == "updraft_helicity": + ncvars = extract_vars(wrfnc, timeidx, ("W", "PH", "PHB", "MAPFAC_M"), + method, squeeze, cache=None, meta=True) + + wstag = ncvars["W"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + mapfct = ncvars["MAPFAC_M"] + + attrs = extract_global_attrs(wrfnc, attrs=("DX", "DY")) + dx = attrs["DX"] + dy = attrs["DY"] + + # As coded in NCL, but not sure this is possible + varname = "U" + u_vars = extract_vars(wrfnc, timeidx, varname, method, squeeze, + cache=None, meta=True) + u = destagger(u_vars[varname], -1, meta=True) + + varname = "V" + v_vars = extract_vars(wrfnc, timeidx, varname, method, squeeze, + cache=None, meta=True) + v = destagger(v_vars[varname], -2, meta=True) + + zstag = ph + phb + + return (zstag, mapfct, u, v, wstag, dx, dy) + + if varname == "omg": + varnames=("T", "P", "W", "PB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + w = ncvars["W"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + + wa = destagger(w, -3) + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t, meta=False) + + return (qv, tkel, wa, full_p) + + if varname == "pw": + varnames=("T", "P", "PB", "PH", "PHB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + ph = ncvars["PH"] + phb = ncvars["PHB"] + qv = ncvars["QVAPOR"] + + # Change this to use real virtual temperature! + full_p = p + pb + ht = (ph + phb)/Constants.G + full_t = t + Constants.T_BASE + + tkel = tk(full_p, full_t, meta=False) + + return (full_p, tkel, qv, ht) + + if varname == "rh": + varnames=("T", "P", "PB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qvapor = npvalues(ncvars["QVAPOR"]) + + full_t = t + Constants.T_BASE + full_p = p + pb + qvapor[qvapor < 0] = 0 + tkel = tk(full_p, full_t, meta=False) + + return (qvapor, full_p, tkel) + + if varname == "slp": + varnames=("T", "P", "PB", "QVAPOR", "PH", "PHB") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qvapor = npvalues(ncvars["QVAPOR"]) + ph = ncvars["PH"] + phb = ncvars["PHB"] + + full_t = t + Constants.T_BASE + full_p = p + pb + qvapor[qvapor < 0] = 0. + + full_ph = (ph + phb) / Constants.G + + destag_ph = destagger(full_ph, -3) + + tkel = tk(full_p, full_t, meta=False) + + return (destag_ph, tkel, full_p, qvapor) + + if varname == "td": + varnames=("P", "PB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + + p = ncvars["P"] + pb = ncvars["PB"] + qvapor = npvalues(ncvars["QVAPOR"]) + + # Algorithm requires hPa + full_p = .01*(p + pb) + qvapor[qvapor < 0] = 0 + + return (full_p, qvapor) + + if varname == "tk": + varnames=("T", "P", "PB") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + + full_t = t + Constants.T_BASE + full_p = p + pb + + return (full_p, full_t) + + if varname == "tv": + varnames=("T", "P", "PB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + + full_t = t + Constants.T_BASE + full_p = p + pb + tkel = tk(full_p, full_t) + + return (tkel, qv) + + if varname == "twb": + varnames=("T", "P", "PB", "QVAPOR") + ncvars = extract_vars(wrfnc, timeidx, varnames, method, squeeze, + cache=None, meta=True) + t = ncvars["T"] + p = ncvars["P"] + pb = ncvars["PB"] + qv = ncvars["QVAPOR"] + + full_t = t + Constants.T_BASE + full_p = p + pb + + tkel = tk(full_p, full_t) + + return (full_p, tkel, qv) + + if varname == "uvmet": + varname = "U" + u_vars = extract_vars(wrfnc, timeidx, varname, method, squeeze, + cache=None, meta=True) + + u = destagger(u_vars[varname], -1, meta=True) + + varname = "V" + v_vars = extract_vars(wrfnc, timeidx, varname, method, squeeze, + cache=None, meta=True) + v = destagger(v_vars[varname], -2, meta=True) + + map_proj_attrs = extract_global_attrs(wrfnc, attrs="MAP_PROJ") + map_proj = map_proj_attrs["MAP_PROJ"] + + if map_proj in (0,3,6): + raise ProjectionError("Map projection does not need rotation") + elif map_proj in (1,2): + lat_attrs = extract_global_attrs(wrfnc, attrs=("TRUELAT1", + "TRUELAT2")) + radians_per_degree = Constants.PI/180.0 + # Rotation needed for Lambert and Polar Stereographic + true_lat1 = lat_attrs["TRUELAT1"] + true_lat2 = lat_attrs["TRUELAT2"] + + try: + lon_attrs = extract_global_attrs(wrfnc, attrs="STAND_LON") + except AttributeError: + try: + cen_lon_attrs = extract_global_attrs(wrfnc, attrs="CEN_LON") + except AttributeError: + raise RuntimeError("longitude attributes not found in NetCDF") + else: + cen_lon = cen_lon_attrs["CEN_LON"] + else: + cen_lon = lon_attrs["STAND_LON"] + + + varname = "XLAT" + xlat_var = extract_vars(wrfnc, timeidx, varname, + method, squeeze, cache=None, meta=True) + lat = xlat_var[varname] + + varname = "XLONG" + xlon_var = extract_vars(wrfnc, timeidx, varname, + method, squeeze, cache=None, meta=True) + lon = xlon_var[varname] + + if map_proj == 1: + if((fabs(true_lat1 - true_lat2) > 0.1) and + (fabs(true_lat2 - 90.) > 0.1)): + cone = (log(cos(true_lat1*radians_per_degree)) + - log(cos(true_lat2*radians_per_degree))) + cone = (cone / + (log(tan((45.-fabs(true_lat1/2.))*radians_per_degree)) + - log(tan((45.-fabs(true_lat2/2.))*radians_per_degree)))) + else: + cone = sin(fabs(true_lat1)*radians_per_degree) + else: + cone = 1 + + return (u, v, lat, lon, cen_lon, cone) + + if varname == "cloudfrac": + vars = extract_vars(wrfnc, timeidx, ("P", "PB", "QVAPOR", "T"), + method, squeeze, cache=None, meta=True) + + p = vars["P"] + pb = vars["PB"] + qv = vars["QVAPOR"] + t = vars["T"] + + full_p = p + pb + full_t = t + Constants.T_BASE + + tkel = tk(full_p, full_t) + relh = rh(qv, full_p, tkel) + + return (full_p, relh) + + + +class WRFVarsTest(ut.TestCase): + longMessage = True + +def make_func(varname, wrfnc, timeidx, method, squeeze, meta): + def func(self): + + try: + args = get_args(varname, wrfnc, timeidx, method, squeeze) + except ProjectionError: # Don't fail for this + return + + routine = ROUTINE_MAP[varname] + + kwargs = {"meta" : meta} + result = routine(*args, **kwargs) + + ref = getvar(wrfnc, varname, timeidx, method, squeeze, cache=None, + meta=meta) + + nt.assert_allclose(npvalues(result), npvalues(ref)) + + if meta: + self.assertEqual(result.dims, ref.dims) + + return func + +if __name__ == "__main__": + varnames = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", + "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", + "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", + "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", + "wa", "uvmet10", "uvmet", "z", "cloudfrac"] + + #varnames = ["helicity"] + varnames=["avo", "pvo", "eth", "dbz", "helicity", "updraft_helicity", + "omg", "pw", "rh", "slp", "td", "tk", "tv", "twb", "uvmet", + "cloudfrac"] + + # Turn this one off when not needed, since it's slow + #varnames += ["cape_2d", "cape_3d"] + + for varname in varnames: + for i,wrfnc in enumerate((NCFILE, NCGROUP)): + for j,timeidx in enumerate((0, ALL_TIMES)): + for method in ("cat", "join"): + for squeeze in (True, False): + for meta in (True, False): + func = make_func(varname, wrfnc, timeidx, method, + squeeze, meta) + ncname = "single" if i == 0 else "multi" + timename = "t0" if j == 0 else "all" + squeeze_name = "squeeze" if squeeze else "nosqueeze" + meta_name = "meta" if meta else "nometa" + test_name = "test_{}_{}_{}_{}_{}_{}".format(varname, + ncname, timename, method, + squeeze_name, meta_name) + + setattr(WRFVarsTest, test_name, func) + + ut.main() + + \ No newline at end of file diff --git a/test/ipynb/WRF_python_demo.ipynb b/test/ipynb/WRF_python_demo.ipynb index 0fbdd4c..f9de36e 100644 --- a/test/ipynb/WRF_python_demo.ipynb +++ b/test/ipynb/WRF_python_demo.ipynb @@ -405,7 +405,7 @@ "pivot_point = (z.shape[-1] / 2, z.shape[-2] / 2) \n", "angle = 90.0\n", "\n", - "p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle)\n", + "p_vert = vertcross(p, z, pivot_point=pivot_point, angle=angle, latlon=True)\n", "print(p_vert)\n", "print (\"\\n\")\n", "del p_vert\n", @@ -414,7 +414,7 @@ "start_point = (0, z.shape[-2]/2)\n", "end_point = (-1, z.shape[-2]/2)\n", "\n", - "p_vert = vertcross(p, z, start_point=start_point, end_point=end_point)\n", + "p_vert = vertcross(p, z, start_point=start_point, end_point=end_point, latlon=True)\n", "print(p_vert)\n", "del p_vert, p, z" ] @@ -441,7 +441,7 @@ "pivot_point = (t2.shape[-2] / 2, t2.shape[-1] / 2) \n", "angle = 90.0\n", "\n", - "t2_line = interpline(t2, pivot_point=pivot_point, angle=angle)\n", + "t2_line = interpline(t2, pivot_point=pivot_point, angle=angle, latlon=True)\n", "print(t2_line, \"\\n\")\n", "\n", "del t2_line\n", @@ -450,7 +450,7 @@ "start_point = (t2.shape[-2]/2, 0)\n", "end_point = (t2.shape[-2]/2, -1)\n", "\n", - "t2_line = interpline(t2, start_point=start_point, end_point=end_point)\n", + "t2_line = interpline(t2, start_point=start_point, end_point=end_point, latlon=True)\n", "print(t2_line, \"\\n\")\n", "\n", "del t2_line, t2" @@ -563,16 +563,16 @@ }, "outputs": [], "source": [ - "from wrf.latlon import get_ll, get_ij # These names are going to change\n", + "from wrf.latlon import xy_to_ll, ll_to_xy \n", "\n", - "a = get_ll(ncfile, [400,105], [200,205])\n", - "b = get_ij(ncfile, 45.5, -110.8)\n", + "a = xy_to_ll(ncfile, [400,105], [200,205])\n", + "b = ll_to_xy(ncfile, 45.5, -110.8)\n", "\n", "# Note: Lists/Dictionaries of files will add a new dimension ('domain') only if the domain is moving\n", - "c = get_ll([ncfile, ncfile, ncfile], [400,105], [200,205])\n", - "d = get_ll({\"label1\" : [ncfile, ncfile],\n", - " \"label2\" : [ncfile, ncfile]},\n", - " [400,105], [200,205])\n", + "c = xy_to_ll([ncfile, ncfile, ncfile], [400,105], [200,205])\n", + "d = xy_to_ll({\"label1\" : [ncfile, ncfile],\n", + " \"label2\" : [ncfile, ncfile]}, \n", + " [400,105], [200,205])\n", "\n", "print(a)\n", "print(\"\\n\")\n", diff --git a/test/utests.py b/test/utests.py index fbf8d50..03c3efe 100644 --- a/test/utests.py +++ b/test/utests.py @@ -1,6 +1,6 @@ import unittest as ut import numpy.testing as nt -import numpy as n +import numpy as np import numpy.ma as ma import os, sys import subprocess @@ -111,9 +111,9 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): masked=True if not masked: - ref_vals = n.zeros(new_dims, data.dtype) + ref_vals = np.zeros(new_dims, data.dtype) else: - ref_vals = ma.asarray(n.zeros(new_dims, data.dtype)) + ref_vals = ma.asarray(np.zeros(new_dims, data.dtype)) for i in xrange(repeat): if (varname != "uvmet" and varname != "uvmet10" @@ -173,13 +173,13 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): tol = 0/100. atol = 200.0 - #print n.amax(n.abs(npvalues(cape_3d[0,:]) - ref_vals[0,:])) + #print np.amax(np.abs(npvalues(cape_3d[0,:]) - ref_vals[0,:])) nt.assert_allclose(npvalues(cape_3d), ref_vals, tol, atol) else: my_vals = getvar(in_wrfnc, varname, timeidx=timeidx) tol = 2/100. atol = 0.1 - #print (n.amax(n.abs(npvalues(my_vals) - ref_vals))) + #print (np.amax(np.abs(npvalues(my_vals) - ref_vals))) nt.assert_allclose(npvalues(my_vals), ref_vals, tol, atol) @@ -203,9 +203,9 @@ def _get_refvals(referent, varname, repeat, multi): masked=True if not masked: - ref_vals = n.zeros(new_dims, data.dtype) + ref_vals = np.zeros(new_dims, data.dtype) else: - ref_vals = ma.asarray(n.zeros(new_dims, data.dtype)) + ref_vals = ma.asarray(np.zeros(new_dims, data.dtype)) for i in xrange(repeat): ref_vals[i,:] = data[:] @@ -311,7 +311,7 @@ def make_interp_test(varname, wrf_in, referent, multi=False, elif (varname == "vinterp"): # Tk to theta fld_tk_theta = _get_refvals(referent, "fld_tk_theta", repeat, multi) - fld_tk_theta = n.squeeze(fld_tk_theta) + fld_tk_theta = np.squeeze(fld_tk_theta) tk = getvar(in_wrfnc, "temp", timeidx=timeidx, units="k") @@ -329,13 +329,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, tol = 5/100. atol = 0.0001 - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_tk_theta))) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_tk_theta))) nt.assert_allclose(npvalues(field), fld_tk_theta, tol, atol) # Tk to theta-e fld_tk_theta_e = _get_refvals(referent, "fld_tk_theta_e", repeat, multi) - fld_tk_theta_e = n.squeeze(fld_tk_theta_e) + fld_tk_theta_e = np.squeeze(fld_tk_theta_e) interp_levels = [200,300,500,1000] @@ -351,13 +351,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, tol = 3/100. atol = 50.0001 - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_tk_theta_e)/fld_tk_theta_e)*100) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_tk_theta_e)/fld_tk_theta_e)*100) nt.assert_allclose(npvalues(field), fld_tk_theta_e, tol, atol) # Tk to pressure fld_tk_pres = _get_refvals(referent, "fld_tk_pres", repeat, multi) - fld_tk_pres = n.squeeze(fld_tk_pres) + fld_tk_pres = np.squeeze(fld_tk_pres) interp_levels = [850,500] @@ -370,14 +370,14 @@ def make_interp_test(varname, wrf_in, referent, multi=False, timeidx=timeidx, log_p=True) - field = n.squeeze(field) + field = np.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_tk_pres))) + #print (np.amax(np.abs(npvalues(field) - fld_tk_pres))) nt.assert_allclose(npvalues(field), fld_tk_pres, tol, atol) # Tk to geoht_msl fld_tk_ght_msl = _get_refvals(referent, "fld_tk_ght_msl", repeat, multi) - fld_tk_ght_msl = n.squeeze(fld_tk_ght_msl) + fld_tk_ght_msl = np.squeeze(fld_tk_ght_msl) interp_levels = [1,2] field = vinterp(in_wrfnc, @@ -389,13 +389,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, timeidx=timeidx, log_p=True) - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_tk_ght_msl))) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_tk_ght_msl))) nt.assert_allclose(npvalues(field), fld_tk_ght_msl, tol, atol) # Tk to geoht_agl fld_tk_ght_agl = _get_refvals(referent, "fld_tk_ght_agl", repeat, multi) - fld_tk_ght_agl = n.squeeze(fld_tk_ght_agl) + fld_tk_ght_agl = np.squeeze(fld_tk_ght_agl) interp_levels = [1,2] field = vinterp(in_wrfnc, @@ -407,13 +407,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, timeidx=timeidx, log_p=True) - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_tk_ght_agl))) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_tk_ght_agl))) nt.assert_allclose(npvalues(field), fld_tk_ght_agl, tol, atol) # Hgt to pressure fld_ht_pres = _get_refvals(referent, "fld_ht_pres", repeat, multi) - fld_ht_pres = n.squeeze(fld_ht_pres) + fld_ht_pres = np.squeeze(fld_ht_pres) z = getvar(in_wrfnc, "height", timeidx=timeidx, units="m") interp_levels = [500,50] @@ -426,13 +426,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, timeidx=timeidx, log_p=True) - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_ht_pres))) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_ht_pres))) nt.assert_allclose(npvalues(field), fld_ht_pres, tol, atol) # Pressure to theta fld_pres_theta = _get_refvals(referent, "fld_pres_theta", repeat, multi) - fld_pres_theta = n.squeeze(fld_pres_theta) + fld_pres_theta = np.squeeze(fld_pres_theta) p = getvar(in_wrfnc, "pressure", timeidx=timeidx) interp_levels = [200,300,500,1000] @@ -445,13 +445,13 @@ def make_interp_test(varname, wrf_in, referent, multi=False, timeidx=timeidx, log_p=True) - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_pres_theta))) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_pres_theta))) nt.assert_allclose(npvalues(field), fld_pres_theta, tol, atol) # Theta-e to pres fld_thetae_pres = _get_refvals(referent, "fld_thetae_pres", repeat, multi) - fld_thetae_pres = n.squeeze(fld_thetae_pres) + fld_thetae_pres = np.squeeze(fld_thetae_pres) eth = getvar(in_wrfnc, "eth", timeidx=timeidx) interp_levels = [850,500,5] @@ -464,8 +464,8 @@ def make_interp_test(varname, wrf_in, referent, multi=False, timeidx=timeidx, log_p=True) - field = n.squeeze(field) - #print (n.amax(n.abs(npvalues(field) - fld_thetae_pres))) + field = np.squeeze(field) + #print (np.amax(np.abs(npvalues(field) - fld_thetae_pres))) nt.assert_allclose(npvalues(field), fld_thetae_pres, tol, atol) return test @@ -483,7 +483,7 @@ if __name__ == "__main__": "geopt", "helicity", "lat", "lon", "omg", "p", "pressure", "pvo", "pw", "rh2", "rh", "slp", "ter", "td2", "td", "tc", "theta", "tk", "tv", "twb", "updraft_helicity", "ua", "va", - "wa", "uvmet10", "uvmet", "z", "ctt", "cloudfrac"] + "wa", "uvmet10", "uvmet", "z", "cloudfrac"] interp_methods = ["interplevel", "vertcross", "interpline", "vinterp"] try: