Browse Source

Merge branch 'release/1.0.4'

main 1.0.4
Bill Ladwig 8 years ago
parent
commit
7445c6358c
  1. 5
      conda_recipe/meta.yaml
  2. 9
      doc/source/new.rst
  3. 46
      fortran/rip_cape.f90
  4. 18
      fortran/wrf_cloud_fracf.f90
  5. 4
      fortran/wrf_rip_phys_routines.f90
  6. 26
      fortran/wrf_user.f90
  7. 4
      fortran/wrf_user_dbz.f90
  8. 18
      fortran/wrf_user_latlon_routines.f90
  9. 2
      fortran/wrf_vinterp.f90
  10. 45
      hpc_install/build-wrapt.ys
  11. 51
      hpc_install/build-wrf-python.ys
  12. 2
      src/wrf/computation.py
  13. 23
      src/wrf/decorators.py
  14. 70
      src/wrf/projection.py
  15. 20
      src/wrf/specialdec.py
  16. 17
      src/wrf/util.py
  17. 2
      src/wrf/version.py
  18. BIN
      test/ci_tests/ci_result_file.nc
  19. BIN
      test/ci_tests/ci_test_file.nc
  20. 8
      test/ci_tests/make_test_file.py

5
conda_recipe/meta.yaml

@ -1,4 +1,4 @@
{% set version = "1.0b3" %} {% set version = "1.0.2" %}
package: package:
name: wrf-python name: wrf-python
@ -7,7 +7,7 @@ package:
source: source:
fn: wrf-python-{{ version }}.tar.gz fn: wrf-python-{{ version }}.tar.gz
url: https://github.com/NCAR/wrf-python/archive/{{ version }}.tar.gz url: https://github.com/NCAR/wrf-python/archive/{{ version }}.tar.gz
sha256: 6fc6a268d91129ba3c16facb1b53194a551d6336582a1f2162ced584e646aa7e sha256: 42df1c46c146010e8d7a14cfa4daca3c30cb530e93d26e1fba34e2118540e0b2
build: build:
number: 0 number: 0
@ -29,6 +29,7 @@ requirements:
- mingwpy # [win] - mingwpy # [win]
- libgfortran # [unix] - libgfortran # [unix]
- libgcc # [unix] - libgcc # [unix]
- xarray
test: test:
requires: requires:

9
doc/source/new.rst

@ -4,6 +4,15 @@ What's New
Releases Releases
------------- -------------
v1.0.4
^^^^^^^^^^^^^^
- Release 1.0.4
- Fix warnings with CI tests which were caused by fill values being written
as NaN to the NetCDF result file.
- Added the __eq__ operator to the WrfProj projection base class.
- Fixed array order issue when using the raw CAPE routine with 1D arrays.
v1.0.3 v1.0.3
^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^

46
fortran/rip_cape.f90

@ -374,17 +374,17 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
G/(RD*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1))) G/(RD*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1)))
p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh)) p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh))
p1 = p2 - pavg p1 = p2 - pavg
totthe = 0.d0 totthe = 0.D0
totqvp = 0.d0 totqvp = 0.D0
totprs = 0.d0 totprs = 0.D0
DO k = mkzh,2,-1 DO k = mkzh,2,-1
IF (prsf(i,j,k) .LE. p1) EXIT !GOTO 35 IF (prsf(i,j,k) .LE. p1) EXIT !GOTO 35
IF (prsf(i,j,k-1) .GE. p2) CYCLE !GOTO 34 IF (prsf(i,j,k-1) .GE. p2) CYCLE !GOTO 34
p = prs(i,j,k) p = prs(i,j,k)
pup = prsf(i,j,k) pup = prsf(i,j,k)
pdn = prsf(i,j,k-1) pdn = prsf(i,j,k-1)
q = MAX(qvp(i,j,k),1.d-15) q = MAX(qvp(i,j,k),1.D-15)
th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q)) th = tmk(i,j,k)*(1000.D0/p)**(GAMMA*(1.D0 + GAMMAMD*q))
pp1 = MAX(p1,pdn) pp1 = MAX(p1,pdn)
pp2 = MIN(p2,pup) pp2 = MIN(p2,pup)
IF (pp2 .GT. pp1) THEN IF (pp2 .GT. pp1) THEN
@ -398,7 +398,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! 35 CONTINUE ! 35 CONTINUE
qvppari = totqvp/totprs qvppari = totqvp/totprs
tmkpari = (totthe/totprs)*& tmkpari = (totthe/totprs)*&
(prs(i,j,kpar1)/1000.d0)**(GAMMA*(1.d0+GAMMAMD*qvp(i,j,kpar1))) (prs(i,j,kpar1)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp(i,j,kpar1)))
END IF END IF
DO kpar = kpar1, kpar2 DO kpar = kpar1, kpar2
@ -412,13 +412,13 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
END IF END IF
prspari = prs(i,j,kpar) prspari = prs(i,j,kpar)
ghtpari = ght(i,j,kpar) ghtpari = ght(i,j,kpar)
gammam = GAMMA * (1.d0 + GAMMAMD*qvppari) gammam = GAMMA * (1.D0 + GAMMAMD*qvppari)
cpm = CP * (1.d0 + CPMD*qvppari) cpm = CP * (1.D0 + CPMD*qvppari)
e = MAX(1.d-20,qvppari*prspari/(EPS + qvppari)) e = MAX(1.D-20,qvppari*prspari/(EPS + qvppari))
tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4 tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4
ethpari = tmkpari*(1000.d0/prspari)**(GAMMA*(1.d0 + GAMMAMD*qvppari))*& ethpari = tmkpari*(1000.D0/prspari)**(GAMMA*(1.D0 + GAMMAMD*qvppari))*&
EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.d0 + THTECON3*qvppari)) EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.D0 + THTECON3*qvppari))
zlcl = ghtpari + (tmkpari - tlcl)/(G/cpm) zlcl = ghtpari + (tmkpari - tlcl)/(G/cpm)
! Calculate buoyancy and relative height of lifted parcel at ! Calculate buoyancy and relative height of lifted parcel at
@ -476,13 +476,13 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
buoy(kk) = G*(tvlift - tvenv)/tvenv buoy(kk) = G*(tvlift - tvenv)/tvenv
zrel(kk) = ghtlift - ghtpari zrel(kk) = ghtlift - ghtpari
IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0d0)) THEN IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN
! Parcel ascent curve crosses sounding curve, so create a new level ! Parcel ascent curve crosses sounding curve, so create a new level
! in the bottom-up array at the crossing. ! in the bottom-up array at the crossing.
kk = kk + 1 kk = kk + 1
buoy(kk) = buoy(kk-1) buoy(kk) = buoy(kk-1)
zrel(kk) = zrel(kk-1) zrel(kk) = zrel(kk-1)
buoy(kk-1) = 0.d0 buoy(kk-1) = 0.D0
zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/&
(buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) (buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2))
END IF END IF
@ -511,11 +511,11 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! Get the accumulated buoyant energy from the parcel's starting ! Get the accumulated buoyant energy from the parcel's starting
! point, at all levels up to the top level. ! point, at all levels up to the top level.
benaccum(1) = 0.0d0 benaccum(1) = 0.0D0
benamin = 9d9 benamin = 9d9
DO k = 2,kmax DO k = 2,kmax
dz = zrel(k) - zrel(k-1) dz = zrel(k) - zrel(k-1)
benaccum(k) = benaccum(k-1) + .5d0*dz*(buoy(k-1) + buoy(k)) benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k))
IF (benaccum(k) .LT. benamin) THEN IF (benaccum(k) .LT. benamin) THEN
benamin = benaccum(k) benamin = benaccum(k)
END IF END IF
@ -527,7 +527,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
elfound = .FALSE. elfound = .FALSE.
DO k = kmax,klcl,-1 DO k = kmax,klcl,-1
IF (buoy(k) .GE. 0.d0) THEN IF (buoy(k) .GE. 0.D0) THEN
! k of equilibrium level ! k of equilibrium level
kel = k kel = k
elfound = .TRUE. elfound = .TRUE.
@ -549,8 +549,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! to a more appropriate missing value, which is passed into this ! to a more appropriate missing value, which is passed into this
! routine as cmsg. ! routine as cmsg.
! cape(i,j,kpar) = -0.1d0 ! cape(i,j,kpar) = -0.1D0
! cin(i,j,kpar) = -0.1d0 ! cin(i,j,kpar) = -0.1D0
IF (.NOT. elfound) THEN IF (.NOT. elfound) THEN
cape(i,j,kpar) = cmsg cape(i,j,kpar) = cmsg
cin(i,j,kpar) = cmsg cin(i,j,kpar) = cmsg
@ -572,7 +572,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! inhibition (cin). ! inhibition (cin).
! First get the lfc according to the above definition. ! First get the lfc according to the above definition.
benamin = 9d9 benamin = 9D9
klfc = kmax klfc = kmax
DO k = klcl,kel DO k = klcl,kel
IF (benaccum(k) .LT. benamin) THEN IF (benaccum(k) .LT. benamin) THEN
@ -583,8 +583,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! Now we can assign values to cape and cin ! Now we can assign values to cape and cin
cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1d0) cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0)
cin(i,j,kpar) = MAX(-benamin, 0.1d0) cin(i,j,kpar) = MAX(-benamin, 0.1D0)
! cin is uninteresting when cape is small (< 100 j/kg), so set ! cin is uninteresting when cape is small (< 100 j/kg), so set
! cin to -0.1 (see note about missing values in v6.1.0) in ! cin to -0.1 (see note about missing values in v6.1.0) in
@ -595,8 +595,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! to a more appropriate missing value, which is passed into this ! to a more appropriate missing value, which is passed into this
! routine as cmsg. ! routine as cmsg.
! IF (cape(i,j,kpar).lt.100.d0) cin(i,j,kpar) = -0.1d0 ! IF (cape(i,j,kpar).lt.100.D0) cin(i,j,kpar) = -0.1D0
IF (cape(i,j,kpar) .LT. 100.d0) cin(i,j,kpar) = cmsg IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg
! 102 CONTINUE ! 102 CONTINUE
END DO END DO

18
fortran/wrf_cloud_fracf.f90

@ -29,11 +29,11 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew)
DO k = 1,nz-1 DO k = 1,nz-1
IF (k .GE. kclo .AND. k .LT. kcmi) THEN IF (k .GE. kclo .AND. k .LT. kcmi) THEN
lowc(i,j) = AMAX1(rh(i,j,k), lowc(i,j)) lowc(i,j) = MAX(rh(i,j,k), lowc(i,j))
ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud ELSE IF (k .GE. kcmi .AND. k .LT. kchi) THEN ! mid cloud
midc(i,j) = AMAX1(rh(i,j,k), midc(i,j)) midc(i,j) = MAX(rh(i,j,k), midc(i,j))
ELSE if (k .GE. kchi) THEN ! high cloud ELSE if (k .GE. kchi) THEN ! high cloud
highc(i,j) = AMAX1(rh(i,j,k), highc(i,j)) highc(i,j) = MAX(rh(i,j,k), highc(i,j))
END IF END IF
END DO END DO
@ -41,12 +41,12 @@ SUBROUTINE DCLOUDFRAC(pres, rh, lowc, midc, highc, nz, ns, ew)
midc(i,j) = 4.0*midc(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 highc(i,j) = 2.5*highc(i,j)/100. - 1.5
lowc(i,j) = amin1(lowc(i,j), 1.0) lowc(i,j) = MIN(lowc(i,j), 1.0)
lowc(i,j) = amax1(lowc(i,j), 0.0) lowc(i,j) = MAX(lowc(i,j), 0.0)
midc(i,j) = amin1(midc(i,j), 1.0) midc(i,j) = MIN(midc(i,j), 1.0)
midc(i,j) = amax1(midc(i,j), 0.0) midc(i,j) = MAX(midc(i,j), 0.0)
highc(i,j) = amin1(highc(i,j), 1.0) highc(i,j) = MIN(highc(i,j), 1.0)
highc(i,j) = amax1(highc(i,j), 0.0) highc(i,j) = MAX(highc(i,j), 0.0)
END DO END DO
END DO END DO

4
fortran/wrf_rip_phys_routines.f90

@ -70,11 +70,11 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg)
DO k=1,nz DO k=1,nz
DO j=1,ny DO j=1,ny
DO i=1,nx DO i=1,nx
q = DMAX1(qvp(i,j,k), 1.D-15) q = MAX(qvp(i,j,k), 1.D-15)
t = tmk(i,j,k) t = tmk(i,j,k)
p = prs(i,j,k)/100. p = prs(i,j,k)/100.
e = q*p/(EPS + q) e = q*p/(EPS + q)
tlcl = TLCLC1/(DLOG(t**TLCLC2/e) - TLCLC3) + TLCLC4 tlcl = TLCLC1/(LOG(t**TLCLC2/e) - TLCLC3) + TLCLC4
eth = t*(1000./p)**(GAMMA*(1. + GAMMAMD*q))*& eth = t*(1000./p)**(GAMMA*(1. + GAMMAMD*q))*&
EXP((THTECON1/tlcl - THTECON2)*q*(1. + THTECON3*q)) EXP((THTECON1/tlcl - THTECON2)*q*(1. + THTECON3*q))

26
fortran/wrf_user.f90

@ -137,8 +137,8 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain)
DO k = 1,nz DO k = 1,nz
DO j = 1,ny DO j = 1,ny
DO i = 1,nx DO i = 1,nx
ii = MIN0(i,nxz) ii = MIN(i,nxz)
im1 = MAX0(i-1,1) im1 = MAX(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 END DO
@ -147,8 +147,8 @@ SUBROUTINE DZSTAG(znew, nx, ny, nz, z, nxz, nyz ,nzz, terrain)
ELSE IF (ny .GT. nyz) THEN ELSE IF (ny .GT. nyz) THEN
DO k = 1,nz DO k = 1,nz
DO j = 1,NY DO j = 1,NY
jj = MIN0(j,nyz) jj = MIN(j,nyz)
jm1 = MAX0(j-1,1) jm1 = MAX(j-1,1)
DO i = 1,nx 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
@ -196,8 +196,8 @@ SUBROUTINE DINTERP2DXY(v3d, v2d, xy, nx, ny, nz, nxy)
REAL(KIND=8) :: w11, w12, w21, w22, wx, wy REAL(KIND=8) :: w11, w12, w21, w22, wx, wy
DO ij = 1,nxy DO ij = 1,nxy
i = MAX0(1,MIN0(nx-1,INT(xy(1,ij)+1))) i = MAX(1,MIN(nx-1,INT(xy(1,ij)+1)))
j = MAX0(1,MIN0(ny-1,INT(xy(2,ij)+1))) j = MAX(1,MIN(ny-1,INT(xy(2,ij)+1)))
wx = DBLE(i+1) - (xy(1,ij)+1) wx = DBLE(i+1) - (xy(1,ij)+1)
wy = DBLE(j+1) - (xy(2,ij)+1) wy = DBLE(j+1) - (xy(2,ij)+1)
w11 = wx*wy w11 = wx*wy
@ -255,7 +255,7 @@ SUBROUTINE DINTERP1D(v_in, v_out, z_in, z_out, vmsg, nz_in, nz_out)
DO WHILE ((.NOT. interp) .AND. (kp .GE. 2)) DO WHILE ((.NOT. interp) .AND. (kp .GE. 2))
IF (((z_in(kp-im) .LE. height) .AND. (z_in(kp-ip) .GT. height))) THEN IF (((z_in(kp-im) .LE. height) .AND. (z_in(kp-ip) .GT. height))) THEN
w2 = (height - z_in(kp-im))/(z_in(kp-ip) - z_in(kp-im)) w2 = (height - z_in(kp-im))/(z_in(kp-ip) - z_in(kp-im))
w1 = 1.d0 - w2 w1 = 1.D0 - w2
v_out(k) = w1*v_in(kp-im) + w2*v_in(kp-ip) v_out(k) = w1*v_in(kp-im) + w2*v_in(kp-ip)
interp = .TRUE. interp = .TRUE.
END IF END IF
@ -608,9 +608,9 @@ SUBROUTINE DCOMPUTERH(qv, p, t, rh, nx)
es = EZERO*EXP(ESLCON1*(temperature - CELKEL)/(temperature - ESLCON2)) es = EZERO*EXP(ESLCON1*(temperature - CELKEL)/(temperature - ESLCON2))
! qvs = ep_2*es/(pressure-es) ! 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 = 100*MAX(1., qv(i)/qvs)
! rh(i) = 100.*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*MAX(MIN(qv(i)/qvs, 1.0D0), 0.0D0)
END DO END DO
RETURN RETURN
@ -645,7 +645,7 @@ SUBROUTINE DGETIJLATLONG(lat_array, long_array, lat, longitude, ii, jj, nx, ny,
DO i = 1,nx DO i = 1,nx
latd = (lat_array(i,j) - lat)**2 latd = (lat_array(i,j) - lat)**2
longd = (long_array(i,j) - longitude)**2 longd = (long_array(i,j) - longitude)**2
! longd = dmin1((long_array(i,j)-longitude)**2, & ! longd = MIN((long_array(i,j)-longitude)**2, &
! (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 IF (dist_min .GT. dist) THEN
@ -792,12 +792,12 @@ SUBROUTINE DCOMPUTETD(td, pressure, qv_in, nx)
INTEGER :: i INTEGER :: i
DO i = 1,nx DO i = 1,nx
qv = DMAX1(qv_in(i), 0.D0) qv = MAX(qv_in(i), 0.D0)
! vapor pressure ! vapor pressure
tdc = qv*pressure(i)/(.622D0 + qv) tdc = qv*pressure(i)/(.622D0 + qv)
! avoid problems near zero ! avoid problems near zero
tdc = DMAX1(tdc, 0.001D0) tdc = MAX(tdc, 0.001D0)
td(i) = (243.5D0*LOG(tdc) - 440.8D0)/(19.48D0 - LOG(tdc)) td(i) = (243.5D0*LOG(tdc) - 440.8D0)/(19.48D0 - LOG(tdc))
END DO END DO
@ -834,7 +834,7 @@ SUBROUTINE DCOMPUTEICLW(iclw, pressure, qc_in, nx, ny, nz)
DO j = 3,ny - 2 DO j = 3,ny - 2
DO i = 3,nx - 2 DO i = 3,nx - 2
DO k = 1,nz DO k = 1,nz
qclw = DMAX1(qc_in(i,j,k), 0.D0) qclw = MAX(qc_in(i,j,k), 0.D0)
IF (k.EQ.1) THEN IF (k.EQ.1) THEN
dp = pressure(i,j,k-1) - pressure(i,j,k) dp = pressure(i,j,k-1) - pressure(i,j,k)
ELSE IF (k.EQ.nz) then ELSE IF (k.EQ.nz) then

4
fortran/wrf_user_dbz.f90

@ -136,8 +136,8 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx
! Calculate variable intercept parameters ! Calculate variable intercept parameters
IF (ivarint .EQ. 1) THEN IF (ivarint .EQ. 1) THEN
temp_c = DMIN1(-0.001D0, tmk(i,j,k)-CELKEL) temp_c = MIN(-0.001D0, tmk(i,j,k)-CELKEL)
sonv = DMIN1(2.0D8, 2.0D6*EXP(-0.12D0*temp_c)) sonv = MIN(2.0D8, 2.0D6*EXP(-0.12D0*temp_c))
gonv = gon gonv = gon
IF (qgr(i,j,k) .GT. R1) THEN IF (qgr(i,j,k) .GT. R1) THEN

18
fortran/wrf_user_latlon_routines.f90

@ -142,14 +142,14 @@ SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
! the rsw tag. ! the rsw tag.
rsw = 0.D0 rsw = 0.D0
IF (lat1 .NE. 0.D0) THEN IF (lat1 .NE. 0.D0) THEN
rsw = (DLOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon rsw = (LOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon
END IF END IF
deltalon = lon - lon1 deltalon = lon - lon1
IF (deltalon .LT. -180.D0) deltalon = deltalon + 360.D0 IF (deltalon .LT. -180.D0) deltalon = deltalon + 360.D0
IF (deltalon .GT. 180.D0) deltalon = deltalon - 360.D0 IF (deltalon .GT. 180.D0) deltalon = deltalon - 360.D0
i = knowni + (deltalon/(dlon*DEG_PER_RAD)) i = knowni + (deltalon/(dlon*DEG_PER_RAD))
j = knownj + (DLOG(TAN(0.5D0*((lat + 90.D0)*RAD_PER_DEG))))/dlon - rsw j = knownj + (LOG(TAN(0.5D0*((lat + 90.D0)*RAD_PER_DEG))))/dlon - rsw
! ps ! ps
ELSE IF (map_proj .EQ. 2) THEN ELSE IF (map_proj .EQ. 2) THEN
@ -183,9 +183,9 @@ SUBROUTINE DLLTOIJ(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
END IF END IF
IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN
cone = (DLOG(COS(truelat1*RAD_PER_DEG))-DLOG(COS(truelat2*RAD_PER_DEG)))/& cone = (LOG(COS(truelat1*RAD_PER_DEG))-LOG(COS(truelat2*RAD_PER_DEG)))/&
(DLOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0))-& (LOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0))-&
DLOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0))) LOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0)))
ELSE ELSE
cone = SIN(ABS(truelat1)*RAD_PER_DEG) cone = SIN(ABS(truelat1)*RAD_PER_DEG)
END IF END IF
@ -358,7 +358,7 @@ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
! the rsw tag. ! the rsw tag.
rsw = 0.D0 rsw = 0.D0
IF (lat1 .NE. 0.D0) THEN IF (lat1 .NE. 0.D0) THEN
rsw = (DLOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon rsw = (LOG(TAN(0.5D0*((lat1 + 90.D0)*RAD_PER_DEG))))/dlon
END IF END IF
lat = 2.0D0*ATAN(EXP(dlon*(rsw + aj - knownj)))*DEG_PER_RAD - 90.D0 lat = 2.0D0*ATAN(EXP(dlon*(rsw + aj - knownj)))*DEG_PER_RAD - 90.D0
@ -418,9 +418,9 @@ SUBROUTINE DIJTOLL(map_proj, truelat1, truelat2, stdlon, lat1, lon1,&
END IF END IF
IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN IF (ABS(truelat1 - truelat2) .GT. 0.1D0) THEN
cone = (DLOG(COS(truelat1*RAD_PER_DEG)) - DLOG(COS(truelat2*RAD_PER_DEG)))/& cone = (LOG(COS(truelat1*RAD_PER_DEG)) - LOG(COS(truelat2*RAD_PER_DEG)))/&
(DLOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0)) - & (LOG(TAN((90.D0 - ABS(truelat1))*RAD_PER_DEG*0.5D0)) - &
DLOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0))) LOG(TAN((90.D0 - ABS(truelat2))*RAD_PER_DEG*0.5D0)))
ELSE ELSE
cone = SIN(ABS(truelat1)*RAD_PER_DEG) cone = SIN(ABS(truelat1)*RAD_PER_DEG)
END IF END IF

2
fortran/wrf_vinterp.f90

@ -369,7 +369,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
ELSE !else for checking above ground ELSE !else for checking above ground
ptarget = psurfsm - 150.D0 ptarget = psurfsm - 150.D0
dpmin = 1.e4 dpmin = 1.E4
DO k=1,nz DO k=1,nz
ripk = nz-k+1 ripk = nz-k+1
dp = ABS((pres(i,j,ripk) * 0.01D0) - ptarget) dp = ABS((pres(i,j,ripk) * 0.01D0) - ptarget)

45
hpc_install/build-wrapt.ys

@ -0,0 +1,45 @@
#!/bin/bash
#
#HPCI -n wrapt
#HPCI -v 1.10.10
#HPCI -a wrapt-1.10.10.tar.gz
#HPCI -p gnu/4.8.2
#HPCI -p python/2.7.7
set -e
PKGSRC=$HPCI_SW_NAME"-"${HPCI_SW_VERSION}
tar -xf $PKGSRC".tar.gz"
cd $PKGSRC
python setup.py build
# Create the module directory
mkdir -p $HPCI_MOD_DIR/pythonpkgs/$HPCI_SW_NAME/
# Create the module with the following template
cat << EOF > $HPCI_MOD_DIR/pythonpkgs/$HPCI_SW_NAME/${HPCI_SW_VERSION}.lua
require("posix")
whatis("wrapt v$HPCI_SW_VERSION")
help([[
This module loads wrapt, a Python package that provides fully functional
decorators. See https://wrapt.readthedocs.io/en/latest/ for details.
]])
local verpath = "$HPCI_SW_DIR" -- specific version path
local pytestpath = pathJoin(verpath, "lib/python2.7/site-packages/") -- internal python libs
prepend_path("PYTHONPATH", pytestpath)
conflict("all-python-libs")
EOF
mkdir -p $HPCI_SW_DIR/lib/python2.7/site-packages/ # create the install directory (python does not install in not existing dirs)
# Need to set pythonpath in order to dump the .pth files
ml use $HPCI_MOD_DIR/pythonpkgs
ml $HPCI_SW_NAME/${HPCI_SW_VERSION}
python setup.py install --prefix=$HPCI_SW_DIR

51
hpc_install/build-wrf-python.ys

@ -0,0 +1,51 @@
#!/bin/bash
#
#HPCI -n wrf-python
#HPCI -v 1.0.1
#HPCI -a wrf-python-1.0.1.tar.gz
#HPCI -p gnu/4.8.2
#HPCI -p python/2.7.7
#HPCI -p numpy/1.11.0
#HPCI -p wrapt/1.10.10
#HPCI -p scipy/0.17.1 bottleneck/1.1.0 numexpr/2.6.0 pyside/1.1.2 matplotlib/1.5.1 pandas/0.18.1 netcdf4python/1.2.4 xarray/0.8.2
set -e
PKGSRC=$HPCI_SW_NAME"-"${HPCI_SW_VERSION}
tar -xf $PKGSRC".tar.gz"
cd $PKGSRC
python setup.py build
# Create the module directory
mkdir -p $HPCI_MOD_DIR/pythonpkgs/$HPCI_SW_NAME/
# Create the module with the following template
cat << EOF > $HPCI_MOD_DIR/pythonpkgs/$HPCI_SW_NAME/${HPCI_SW_VERSION}.lua
require("posix")
whatis("wrf-python v$HPCI_SW_VERSION")
help([[
This module loads wrf-python, a Python package that provides diagnostic and
interpolation routines for users of the WRF-ARW model.
See http://wrf-python.rtfd.org for details.
]])
local verpath = "$HPCI_SW_DIR" -- specific version path
local pypath = pathJoin(verpath, "lib/python2.7/site-packages/") -- internal python libs
prepend_path("PYTHONPATH", pypath)
conflict("all-python-libs")
prereq($HPCI_MOD_PREREQ)
EOF
mkdir -p $HPCI_SW_DIR/lib/python2.7/site-packages/ # create the install directory (python does not install in not existing dirs)
# Need to set pythonpath in order to dump the .pth files
ml use $HPCI_MOD_DIR/pythonpkgs
ml $HPCI_SW_NAME/${HPCI_SW_VERSION}
python setup.py install --prefix=$HPCI_SW_DIR

2
src/wrf/computation.py

@ -26,7 +26,7 @@ def xy(field, pivot_point=None, angle=None, start_point=None, end_point=None,
Args: Args:
field3d (:class:`xarray.DataArray` or :class:`numpy.ndarray`): A field (:class:`xarray.DataArray` or :class:`numpy.ndarray`): A
field with at least two dimensions. field with at least two dimensions.
pivot_point (:obj:`tuple` or :obj:`list`, optional): A pivot_point (:obj:`tuple` or :obj:`list`, optional): A

23
src/wrf/decorators.py

@ -186,26 +186,29 @@ def left_iteration(ref_var_expected_dims,
# Skip the possible empty/missing arrays for the join method # Skip the possible empty/missing arrays for the join method
skip_missing = False skip_missing = False
for arg in new_args: for arg in new_args:
if isinstance(arg, DataArray): try:
arr = to_np(arg) _ = arg.ndim
elif isinstance(arg, np.ndarray): except AttributeError:
arr = arg continue # Not an array object
else: else:
continue arr = to_np(arg)
if isinstance(arr, np.ma.MaskedArray): try:
if arr.mask.all(): all_masked = arr.mask.all()
except AttributeError:
pass # Not a masked array
else:
if all_masked:
for output in viewvalues(outd): for output in viewvalues(outd):
output[left_and_slice_idxs] = ( output[left_and_slice_idxs] = (
Constants.DEFAULT_FILL) Constants.DEFAULT_FILL)
skip_missing = True skip_missing = True
mask_output = True mask_output = True
break
if skip_missing: if skip_missing:
continue continue
# Insert the output views if one hasn't been provided # Insert the output views if one hasn't been provided
if "outview" not in new_kargs: if "outview" not in new_kargs:
for outkey,output in viewitems(outd): for outkey,output in viewitems(outd):

70
src/wrf/projection.py

@ -2,6 +2,7 @@ from __future__ import (absolute_import, division, print_function,
unicode_literals) unicode_literals)
import numpy as np import numpy as np
import math import math
from decimal import Decimal, Context, ROUND_HALF_UP
from .config import basemap_enabled, cartopy_enabled, pyngl_enabled from .config import basemap_enabled, cartopy_enabled, pyngl_enabled
from .constants import Constants, ProjectionTypes from .constants import Constants, ProjectionTypes
@ -176,6 +177,73 @@ class WrfProj(object):
if self.stand_lon is None: if self.stand_lon is None:
self.stand_lon = self._cen_lon self.stand_lon = self._cen_lon
@staticmethod
def _context_equal(x, y, ctx):
"""Return True if both objects are equal based on the provided context.
Args:
x (numeric): A numeric value.
y (numeric): A numeric value.
ctx (:class:`decimal.Context`): A decimal Context object.
Returns:
:obj:`bool`: True if the values are equal based on the provided
context, False otherwise.
"""
if x is not None:
if y is None:
return False
# Note: The float conversion is because these may come in as
# numpy.float32 or numpy.float64, which Decimal does not know
# how to handle.
if (Decimal(float(x)).normalize(ctx) !=
Decimal(float(y)).normalize(ctx)):
return False
else:
if y is not None:
return False
return True
def __eq__(self, other):
"""Return True if this projection object is the same as *other*.
Note: WRF can use either floats or doubles, so only going to
guarantee floating point precision equivalence, in case the different
types are ever compared (or a projection is hand constructed). For WRF
domains, 7-digit equivalence should be good enough.
"""
if self.map_proj is not None:
if self.map_proj != other.map_proj:
return False
else:
if other.map_proj is not None:
return False
# Only checking for floating point equal.
ctx = Context(prec=7, rounding=ROUND_HALF_UP)
return (WrfProj._context_equal(self.truelat1, other.truelat1, ctx) and
WrfProj._context_equal(self.truelat2, other.truelat2, ctx) and
WrfProj._context_equal(self.moad_cen_lat, other.moad_cen_lat,
ctx) and
WrfProj._context_equal(self.stand_lon, other.stand_lon,
ctx) and
WrfProj._context_equal(self.pole_lat, other.pole_lat, ctx) and
WrfProj._context_equal(self.pole_lon, other.pole_lon, ctx) and
WrfProj._context_equal(self.dx, other.dx, ctx) and
WrfProj._context_equal(self.dy, other.dy, ctx))
def _basemap(self, geobounds, **kwargs): def _basemap(self, geobounds, **kwargs):
@ -291,7 +359,7 @@ class WrfProj(object):
"""Return a :class:`matplotlib.mpl_toolkits.basemap.Basemap` object """Return a :class:`matplotlib.mpl_toolkits.basemap.Basemap` object
for the map projection. for the map projection.
Arguments: Args:
geobounds (:class:`wrf.GeoBounds`, optional): The geobounds to geobounds (:class:`wrf.GeoBounds`, optional): The geobounds to
get the extents. If set to None and using the *var* parameter, get the extents. If set to None and using the *var* parameter,

20
src/wrf/specialdec.py

@ -70,13 +70,9 @@ def uvmet_left_iter(alg_dtype=np.float64):
mode = 2 # probably 3D with 2D lat/lon plus Time mode = 2 # probably 3D with 2D lat/lon plus Time
has_missing = False has_missing = False
u_arr = u u_arr = to_np(u)
if isinstance(u, DataArray):
u_arr = to_np(u)
v_arr = v v_arr = to_np(v)
if isinstance(v, DataArray):
v_arr = to_np(v)
umissing = Constants.DEFAULT_FILL umissing = Constants.DEFAULT_FILL
if isinstance(u_arr, np.ma.MaskedArray): if isinstance(u_arr, np.ma.MaskedArray):
@ -262,7 +258,6 @@ def cape_left_iter(alg_dtype=np.float64):
flip = False flip = False
if p_hpa[0] > p_hpa[-1]: if p_hpa[0] > p_hpa[-1]:
flip = True flip = True
p_hpa = np.ascontiguousarray(p_hpa[::-1]) p_hpa = np.ascontiguousarray(p_hpa[::-1])
tk = np.ascontiguousarray(tk[::-1]) tk = np.ascontiguousarray(tk[::-1])
@ -270,10 +265,13 @@ def cape_left_iter(alg_dtype=np.float64):
ht = np.ascontiguousarray(ht[::-1]) ht = np.ascontiguousarray(ht[::-1])
# Need to make 3D views for the fortran code. # Need to make 3D views for the fortran code.
new_args[0] = p_hpa[:, np.newaxis, np.newaxis] # Going to make these fortran ordered, since the f_contiguous and
new_args[1] = tk[:, np.newaxis, np.newaxis] # c_contiguous flags are broken in numpy 1.11 (always false). This
new_args[2] = qv[:, np.newaxis, np.newaxis] # should work across all numpy versions.
new_args[3] = ht[:, np.newaxis, np.newaxis] new_args[0] = p_hpa.reshape((1, 1, p_hpa.shape[0]), order='F')
new_args[1] = tk.reshape((1, 1, tk.shape[0]), order='F')
new_args[2] = qv.reshape((1, 1, qv.shape[0]), order='F')
new_args[3] = ht.reshape((1, 1, ht.shape[0]), order='F')
new_args[4] = np.full((1,1), ter, orig_dtype) new_args[4] = np.full((1,1), ter, orig_dtype)
new_args[5] = np.full((1,1), sfp, orig_dtype) new_args[5] = np.full((1,1), sfp, orig_dtype)

17
src/wrf/util.py

@ -3021,6 +3021,9 @@ def geo_bounds(var=None, wrfin=None, varname=None, timeidx=0, method="cat",
# Getting lat/lon from xarray coordinates # Getting lat/lon from xarray coordinates
if var is not None: if var is not None:
if not xarray_enabled():
raise ValueError("xarray is not installed or is disabled")
is_moving = None is_moving = None
try: try:
var_coords = var.coords var_coords = var.coords
@ -3142,8 +3145,15 @@ def _get_wrf_proj_geobnds(var, wrfin, varname, timeidx, method, squeeze,
""" """
# Using a variable # Using a variable
if var is not None: if var is not None:
if not xarray_enabled():
raise ValueError("xarray is not installed or is disabled")
geobnds = geo_bounds(var) geobnds = geo_bounds(var)
wrf_proj = var.attrs["projection"] try:
wrf_proj = var.attrs["projection"]
except AttributeError:
raise ValueError("variable does not contain projection "
"information")
else: else:
geobnds = geo_bounds(wrfin=wrfin, varname=varname, timeidx=timeidx, geobnds = geo_bounds(wrfin=wrfin, varname=varname, timeidx=timeidx,
method=method, cache=cache) method=method, cache=cache)
@ -3260,6 +3270,9 @@ def latlon_coords(var, as_np=False):
""" """
if not xarray_enabled():
raise ValueError("xarray is not installed or is disabled")
try: try:
var_coords = var.coords var_coords = var.coords
except AttributeError: except AttributeError:
@ -3284,8 +3297,6 @@ def latlon_coords(var, as_np=False):
return lats, lons return lats, lons
def get_cartopy(var=None, wrfin=None, varname=None, timeidx=0, method="cat", def get_cartopy(var=None, wrfin=None, varname=None, timeidx=0, method="cat",
squeeze=True, cache=None): squeeze=True, cache=None):
"""Return a :class:`cartopy.crs.Projection` subclass for the """Return a :class:`cartopy.crs.Projection` subclass for the

2
src/wrf/version.py

@ -1,2 +1,2 @@
__version__ = "1.0.3" __version__ = "1.0.4"

BIN
test/ci_tests/ci_result_file.nc

Binary file not shown.

BIN
test/ci_tests/ci_test_file.nc

Binary file not shown.

8
test/ci_tests/make_test_file.py

@ -7,7 +7,7 @@ import numpy as np
from netCDF4 import Dataset from netCDF4 import Dataset
from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round, from wrf import (getvar, interplevel, interpline, vertcross, vinterp, py2round,
CoordPair, ll_to_xy, xy_to_ll) CoordPair, ll_to_xy, xy_to_ll, to_np)
VARS_TO_KEEP = ("XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", "XLONG_V", VARS_TO_KEEP = ("XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", "XLONG_V",
"U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", "T2", "U", "V", "W", "PH", "PHB", "T", "P", "PB", "Q2", "T2",
@ -15,7 +15,7 @@ VARS_TO_KEEP = ("XLAT", "XLONG", "XLAT_U", "XLAT_V", "XLONG_U", "XLONG_V",
"QGRAUP", "QRAIN", "QSNOW", "MAPFAC_M", "MAPFAC_U", "QGRAUP", "QRAIN", "QSNOW", "MAPFAC_M", "MAPFAC_U",
"MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC") "MAPFAC_V", "F", "HGT", "RAINC", "RAINSH", "RAINNC")
DIMS_TO_TRIM = ("west_east", "south_north", "bottom_top", "bottom_top_stag", DIMS_TO_TRIM = ("west_east", "south_north", #"bottom_top", "bottom_top_stag",
"west_east_stag", "south_north_stag") "west_east_stag", "south_north_stag")
WRF_DIAGS = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz", WRF_DIAGS = ["avo", "eth", "cape_2d", "cape_3d", "ctt", "dbz", "mdbz",
@ -81,7 +81,9 @@ def add_to_ncfile(outfile, var, varname):
ncvar = outfile.createVariable(varname, var.dtype, var.dims, ncvar = outfile.createVariable(varname, var.dtype, var.dims,
zlib=True, fill_value=fill_value) zlib=True, fill_value=fill_value)
ncvar[:] = var[:]
var_vals = to_np(var)
ncvar[:] = var_vals[:]
def make_result_file(opts): def make_result_file(opts):

Loading…
Cancel
Save