diff --git a/fortran/wrf_pw.f90 b/fortran/wrf_pw.f90 index ba19ffb..275deb4 100644 --- a/fortran/wrf_pw.f90 +++ b/fortran/wrf_pw.f90 @@ -18,14 +18,21 @@ SUBROUTINE DCOMPUTEPW(p, tv, qv, ht, pw, nx, ny, nz, nzh) !REAL(KIND=8),PARAMETER :: R=287.06 pw = 0 + + !$OMP PARALLEL + DO k=1,nz + !$OMP DO COLLAPSE(2) 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))) END DO END DO + !$OMP END DO END DO + !$OMP END PARALLEL + RETURN END SUBROUTINE DCOMPUTEPW diff --git a/fortran/wrf_relhl.f90 b/fortran/wrf_relhl.f90 index 9efeec8..43d4950 100644 --- a/fortran/wrf_relhl.f90 +++ b/fortran/wrf_relhl.f90 @@ -56,10 +56,10 @@ 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 + !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i,j,k,k10,k3,ktop, cu, cv, x, & + !$OMP sum, dh, sdh, su, sv, ua, va, asp, adr, bsp, bdr) DO j=1, mjx DO i=1, miy - !DO i=1, miy-1 sdh = 0.D0 su = 0.D0 sv = 0.D0 @@ -82,12 +82,14 @@ SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) IF (k10 .EQ. 0) THEN k10 = 2 ENDIF + DO k = k3, k10, -1 dh = ght(i,j,k-1) - ght(i,j,k) sdh = sdh + dh 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 asp = SQRT(ua*ua + va*va) @@ -96,11 +98,13 @@ SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) ELSE adr = DEG_PER_RAD * (PI + ATAN2(ua,va)) ENDIF + 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) sum = 0.D0 @@ -112,6 +116,7 @@ SUBROUTINE DCALRELHL(u, v, ght, ter, top, sreh, miy, mjx, mkzh) sreh(i,j) = -sum END DO END DO + !$OMP END PARALLEL DO RETURN diff --git a/fortran/wrf_rip_phys_routines.f90 b/fortran/wrf_rip_phys_routines.f90 index 0b307f2..ad1306f 100644 --- a/fortran/wrf_rip_phys_routines.f90 +++ b/fortran/wrf_rip_phys_routines.f90 @@ -58,8 +58,16 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) REAL(KIND=8) :: tonpsadiabat INTEGER :: l1, h1, mid1, rang1, l2, h2, mid2, rang2 + INTEGER :: errcnt1, errcnt2, bad_i, bad_j, bad_k + REAL(KIND=8) :: bad_p, bad_eth !INTEGER :: ip, ipch, jt, jtch + errcnt1 = 0 + errcnt2 = 0 + bad_i = -1 + bad_j = -1 + bad_k = -1 + ! Before looping, set lookup table for getting temperature on ! a pseudoadiabat. @@ -69,6 +77,9 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) RETURN END IF + !$OMP PARALLEL DO COLLAPSE(3) PRIVATE (i, j, k, jt, ip, q, t, p, e, tlcl, & + !$OMP eth, fracip, fracip2, fracjt, fracjt2, l1, h1, mid1, rang1, l2, h2, & + !$OMP mid2, rang2, tonpsadiabat) REDUCTION(+:errcnt1, errcnt2) DO k=1,nz DO j=1,ny DO i=1,nx @@ -100,52 +111,45 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) rang1 = h1 - l1 mid1 = (h1 + l1) / 2 DO WHILE(rang1 .GT. 1) - if(eth .GE. psadithte(mid1)) then + IF (eth .GE. psadithte(mid1)) THEN l1 = mid1 - else + ELSE h1 = mid1 - end if + END IF rang1 = h1 - l1 mid1 = (h1 + l1) / 2 END DO jt = l1 -! jt=-1 -! DO jtch=1,150-1 -! IF (eth .GE. PSADITHTE(jtch) .AND. eth .LT. PSADITHTE(jtch+1)) THEN -! jt = jtch -! EXIT -! ENDIF -! END DO - - ip = -1 - l2 = 1 - h2 = 149 - rang2 = h2 - l2 - mid2 = (h2 + l2) / 2 - DO WHILE(rang2 .GT. 1) - if(p .LE. psadiprs(mid2)) then - l2 = mid2 - else - h2 = mid2 - end if - rang2 = h2 - l2 - mid2 = (h2 + l2) / 2 - END DO - ip = l2 - -! ip=-1 -! DO ipch=1,150-1 -! IF (p .LE. PSADIPRS(ipch) .AND. p .GT. PSADIPRS(ipch+1)) THEN -! ip = ipch -! EXIT -! ENDIF -! END DO + ip = -1 + l2 = 1 + h2 = 149 + rang2 = h2 - l2 + mid2 = (h2 + l2) / 2 + DO WHILE(rang2 .GT. 1) + IF (p .LE. psadiprs(mid2)) THEN + l2 = mid2 + ELSE + h2 = mid2 + END IF + rang2 = h2 - l2 + mid2 = (h2 + l2) / 2 + END DO + ip = l2 IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN - errstat = ALGERR - WRITE(errmsg, *) "Outside of lookup table bounds. prs,thte=", p, eth - RETURN + errcnt1 = errcnt1 + 1 + !$OMP CRITICAL + ! Only do this the first time + IF (bad_i .EQ. -1) THEN + bad_i = i + bad_j = j + bad_k = k + bad_p = p + bad_eth = eth + END IF + !$OMP END CRITICAL + CYCLE ENDIF fracjt = (eth - PSADITHTE(jt))/(PSADITHTE(jt+1) - PSADITHTE(jt)) @@ -153,12 +157,21 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) fracip = (PSADIPRS(ip) - p)/(PSADIPRS(ip) - PSADIPRS(ip+1)) fracip2 = 1. - fracip + IF (PSADITMK(ip,jt) .GT. 1e9 .OR. PSADITMK(ip+1,jt) .GT. 1e9 .OR. & PSADITMK(ip,jt+1) .GT. 1e9 .OR. PSADITMK(ip+1,jt+1) .GT. 1e9) THEN - !PRINT*,'Tried to access missing tmperature in lookup table.' - errstat = ALGERR - WRITE(errmsg, *) "Prs and Thte probably unreasonable. prs, thte=", p, eth - RETURN + errcnt2 = errcnt2 + 1 + !$OMP CRITICAL + ! Only do this the first time + IF (bad_i .EQ. -1) THEN + bad_i = i + bad_j = j + bad_k = k + bad_p = p + bad_eth = eth + END IF + !$OMP END CRITICAL + CYCLE ENDIF tonpsadiabat = fracip2*fracjt2*PSADITMK(ip,jt) + & @@ -173,6 +186,20 @@ SUBROUTINE WETBULBCALC(prs, tmk, qvp, twb, nx, ny, nz, psafile, errstat, errmsg) END DO END DO + !$OMP END PARALLEL DO + + IF (errcnt1 > 0) THEN + errstat = ALGERR + WRITE(errmsg, *) "Outside of lookup table bounds. i=", bad_i, ",j=", bad_j, ",k=", bad_k, ",p=", bad_p, ",eth=", bad_eth + RETURN + END IF + + IF (errcnt2 > 0) THEN + errstat = ALGERR + WRITE(errmsg, *) "PRS and THTE unreasonable. i=", bad_i, ",j=", bad_j, ",k=", bad_k, ",p=", bad_p, ",eth=", bad_eth + RETURN + END IF + RETURN END SUBROUTINE WETBULBCALC @@ -229,6 +256,7 @@ SUBROUTINE OMGCALC(qvp, tmk, www, prs, omg, mx, my, mz) INTEGER :: i, j, k !REAL(KIND=8), PARAMETER :: GRAV=9.81, RGAS=287.04, EPS=0.622 + !$OMP PARALLEL DO COLLAPSE(3) DO k=1,mz DO j=1,my DO i=1,mx @@ -238,6 +266,7 @@ SUBROUTINE OMGCALC(qvp, tmk, www, prs, omg, mx, my, mz) END DO END DO END DO + !$OMP END PARALLEL DO RETURN @@ -289,6 +318,7 @@ SUBROUTINE VIRTUAL_TEMP(temp, ratmix, tv, nx, ny, nz) INTEGER :: i,j,k !REAL(KIND=8),PARAMETER :: EPS = 0.622D0 + !$OMP PARALLEL DO COLLAPSE(3) DO k=1,nz DO j=1,ny DO i=1,nx @@ -296,6 +326,7 @@ SUBROUTINE VIRTUAL_TEMP(temp, ratmix, tv, nx, ny, nz) END DO END DO END DO + !$OMP END PARALLEL DO RETURN diff --git a/fortran/wrf_user_dbz.f90 b/fortran/wrf_user_dbz.f90 index d1151f0..d2c6098 100644 --- a/fortran/wrf_user_dbz.f90 +++ b/fortran/wrf_user_dbz.f90 @@ -77,7 +77,10 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx REAL(KIND=8), PARAMETER :: RN0_S = 2.D7 REAL(KIND=8), PARAMETER :: RN0_G = 4.D6 + !$OMP PARALLEL + ! Force all Q arrays to be 0.0 or greater. + !$OMP DO COLLAPSE(3) DO k = 1,nz DO j = 1,ny DO i = 1,nx @@ -96,10 +99,12 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx END DO END DO END DO + !$OMP END DO ! Input pressure is Pa, but we need hPa in calculations IF (sn0 .EQ. 0) THEN + !$OMP DO COLLAPSE(3) DO k = 1,nz DO j = 1,ny DO i = 1,nx @@ -110,12 +115,17 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx END DO END DO END DO + !$OMP END DO END IF factor_r = GAMMA_SEVEN*1.D18*(1.D0/(PI*RHO_R))**1.75D0 factor_s = GAMMA_SEVEN*1.D18*(1.D0/(PI*RHO_S))**1.75D0*(RHO_S/RHOWAT)**2*ALPHA factor_g = GAMMA_SEVEN*1.D18*(1.D0/(PI*RHO_G))**1.75D0*(RHO_G/RHOWAT)**2*ALPHA + + !$OMP DO COLLAPSE(3) PRIVATE(i, j, k, temp_c, virtual_t, gonv, ronv, sonv, & + !$OMP factorb_g, factorb_s, rhoair, z_e) & + !$OMP FIRSTPRIVATE(factor_r, factor_s, factor_g) DO k = 1,nz DO j = 1,ny DO i = 1,nx @@ -171,6 +181,9 @@ SUBROUTINE CALCDBZ(prs, tmk, qvp, qra, qsn, qgr, sn0, ivarint, iliqskin, dbz, nx END DO END DO END DO + !$OMP END DO + + !$OMP END PARALLEL RETURN diff --git a/fortran/wrf_wind.f90 b/fortran/wrf_wind.f90 index b08c019..88d022e 100644 --- a/fortran/wrf_wind.f90 +++ b/fortran/wrf_wind.f90 @@ -13,11 +13,13 @@ SUBROUTINE DCOMPUTEWSPD(wspd, u, v, nx, ny) INTEGER i, j + !$OMP PARALLEL DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx wspd(i,j) = SQRT(u(i,j)*u(i,j) + v(i,j)*v(i,j)) END DO END DO + !$OMP END PARALLEL DO END SUBROUTINE DCOMPUTEWSPD @@ -38,11 +40,13 @@ SUBROUTINE DCOMPUTEWDIR(wdir, u, v, nx, ny) INTEGER i, j + !$OMP PARALLEL DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx wdir(i,j) = MOD(270.0 - ATAN2(v(i,j), u(i,j)) * DEG_PER_RAD, 360.) END DO END DO + !$OMP END PARALLEL DO END SUBROUTINE DCOMPUTEWDIR diff --git a/test/utests.py b/test/utests.py index b7b26f1..51f0a97 100644 --- a/test/utests.py +++ b/test/utests.py @@ -143,7 +143,11 @@ def make_test(varname, wrf_in, referent, multi=False, repeat=3, pynio=False): tol = .5/100.0 atol = 0 # NCL uses different constants and doesn't use same # handrolled virtual temp in method - nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol) + try: + nt.assert_allclose(to_np(my_vals), ref_vals, tol, atol) + except AssertionError: + print (np.amax(np.abs(to_np(my_vals) - ref_vals))) + raise elif (varname == "cape_2d"): cape_2d = getvar(in_wrfnc, varname, timeidx=timeidx) tol = 0/100. @@ -264,22 +268,25 @@ def make_interp_test(varname, wrf_in, referent, multi=False, pivot_point = CoordPair(hts.shape[-1] / 2, hts.shape[-2] / 2) ht_cross = vertcross(hts, p, pivot_point=pivot_point, angle=90.) - - nt.assert_allclose(to_np(ht_cross), ref_ht_cross, rtol=.01) + + # Note: Until the bug is fixed in NCL, the wrf-python cross + # sections will contain one extra point + nt.assert_allclose(to_np(ht_cross)[...,0:-1], ref_ht_cross, rtol=.01) # Test opposite p_cross1 = vertcross(p,hts,pivot_point=pivot_point, angle=90.0) - - nt.assert_allclose(to_np(p_cross1), + + nt.assert_allclose(to_np(p_cross1)[...,0:-1], ref_p_cross, rtol=.01) # Test point to point start_point = CoordPair(0, hts.shape[-2]/2) end_point = CoordPair(-1,hts.shape[-2]/2) - + + p_cross2 = vertcross(p,hts,start_point=start_point, end_point=end_point) - + nt.assert_allclose(to_np(p_cross1), to_np(p_cross2)) @@ -292,7 +299,8 @@ def make_interp_test(varname, wrf_in, referent, multi=False, t2_line1 = interpline(t2, pivot_point=pivot_point, angle=90.0) - nt.assert_allclose(to_np(t2_line1), ref_t2_line) + # Note: After NCL is fixed, remove the slice. + nt.assert_allclose(to_np(t2_line1)[...,0:-1], ref_t2_line) # Test point to point start_point = CoordPair(0, t2.shape[-2]/2)