From ea7aba57b2c4703c524cbaea337258c8831971ad Mon Sep 17 00:00:00 2001 From: Bill Ladwig Date: Thu, 2 Nov 2017 16:36:41 -0600 Subject: [PATCH] Updated OpenMP directives. --- fortran/calc_uh.f90 | 13 +++++-- fortran/rip_cape.f90 | 58 +++++++++++++-------------- fortran/wrf_fctt.f90 | 45 +++++++++------------ fortran/wrf_pvo.f90 | 34 ++++++++-------- fortran/wrf_user.f90 | 86 +++++++++++++++++++---------------------- fortran/wrf_vinterp.f90 | 28 +++++++------- 6 files changed, 126 insertions(+), 138 deletions(-) diff --git a/fortran/calc_uh.f90 b/fortran/calc_uh.f90 index b99a48e..6eccb91 100644 --- a/fortran/calc_uh.f90 +++ b/fortran/calc_uh.f90 @@ -61,7 +61,10 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & twodx = 2.0*dx twody = 2.0*dy - !$OMP PARALLEL DO COLLAPSE(3) + + !$OMP PARALLEL + + !$OMP DO COLLAPSE(3) DO k=2,nz-2 DO j=2,ny-1 DO i=2,nx-1 @@ -73,13 +76,13 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & END DO END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO ! Integrate over depth uhminhgt to uhmxhgt AGL ! ! WRITE(6,'(a,f12.1,a,f12.1,a)') & ! 'Calculating UH from ',uhmnhgt,' to ',uhmxhgt,' m AGL' - !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, zbot, ztop, kbot, ktop, & + !$OMP DO COLLAPSE(2) PRIVATE(i, j, k, zbot, ztop, kbot, ktop, & !$OMP wgtlw, wbot, wtop, wsum, wmean, sum, helbot, heltop) DO j=2,ny-2 DO i=2,nx-2 @@ -147,7 +150,9 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, & END IF END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO + + !$OMP END PARALLEL uh = uh*1000. ! Scale according to Kain et al. (2008) diff --git a/fortran/rip_cape.f90 b/fortran/rip_cape.f90 index c4429b4..d474702 100644 --- a/fortran/rip_cape.f90 +++ b/fortran/rip_cape.f90 @@ -85,11 +85,11 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam rang1 = h1 - l1 mid1 = (h1 + l1) / 2 DO WHILE(rang1 .GT. 1) - if(thte .GE. psadithte(mid1)) then + IF (thte .GE. psadithte(mid1)) THEN l1 = mid1 - else + ELSE h1 = mid1 - end if + END IF rang1 = h1 - l1 mid1 = (h1 + l1) / 2 END DO @@ -103,17 +103,17 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam ! END IF ! END DO - ip = -1 + ip = -1 l2 = 1 h2 = 149 rang2 = h2 - l2 mid2 = (h2 + l2) / 2 DO WHILE(rang2 .GT. 1) - if(prs .LE. psadiprs(mid2)) then + IF (prs .LE. psadiprs(mid2)) THEN l2 = mid2 - else + ELSE h2 = mid2 - end if + END IF rang2 = h2 - l2 mid2 = (h2 + l2) / 2 END DO @@ -386,8 +386,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& !$OMP i,j,k,kpar) DO j = 1,mjy DO i = 1,mix - cape(i,j,1) = 0.d0 - cin(i,j,1) = 0.d0 + cape(i,j,1) = 0.D0 + cin(i,j,1) = 0.D0 !!$OMP SIMD DO kpar = 2, mkzh @@ -421,7 +421,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& klcl = 1 END IF -!!$OMP SIMD lastprivate(qvplift,tmklift,ghtlift,tvlift,tmkenv,qvpenv,tvenv,eslift,facden) + !!$OMP SIMD lastprivate(qvplift,tmklift,ghtlift,tvlift,tmkenv,qvpenv,tvenv,eslift,facden) DO k = kpar,1,-1 ! For arrays that go bottom to top kk = kk + 1 @@ -684,7 +684,7 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! kg/kg (should range from 0.000 to 0.025) ! -!$OMP PARALLEL DO COLLAPSE(3) + !$OMP PARALLEL DO COLLAPSE(3) DO j = 1,mjy DO i = 1,mix DO k = 1,mkzh @@ -695,7 +695,7 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& END DO END DO END DO -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO ! calculated the pressure at full sigma levels (a set of pressure @@ -715,29 +715,29 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& !nthreads = omp_get_num_threads() -!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, & -!$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, & -!$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, & -!$OMP benaccum, zrel, kmax, dz, elfound, & -!$OMP kel, klfc, pavg, p2, p1, totthe, totqvp, totprs, & -!$OMP i,j,k,kpar, kpar1, kpar2, qvppari, tmkpari,p, pup, pdn, q, th, & -!$OMP pp1, pp2, ethmax, eth_temp, klev) + !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, & + !$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, & + !$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, & + !$OMP benaccum, zrel, kmax, dz, elfound, & + !$OMP kel, klfc, pavg, p2, p1, totthe, totqvp, totprs, & + !$OMP i,j,k,kpar, kpar1, kpar2, qvppari, tmkpari,p, pup, pdn, q, th, & + !$OMP pp1, pp2, ethmax, eth_temp, klev) DO j = 1,mjy DO i = 1,mix - cape(i,j,1) = 0.d0 - cin(i,j,1) = 0.d0 + cape(i,j,1) = 0.D0 + cin(i,j,1) = 0.D0 ! find parcel with max theta-e in lowest 3 km agl. - ethmax = -1.d0 - eth_temp = -1.d0 + ethmax = -1.D0 + eth_temp = -1.D0 DO k = 1, mkzh - IF (ght_new(k,i,j)-ter(i,j) .LT. 3000.d0) THEN + IF (ght_new(k,i,j)-ter(i,j) .LT. 3000.D0) THEN tlcl = TLCLC1 / (LOG(tmk_new(k,i,j)**TLCLC2/& (MAX(qvp_new(k,i,j), 1.d-15)*prs_new(k,i,j)/(EPS+MAX(qvp_new(k,i,j), 1.d-15))))-TLCLC3)+& TLCLC4 - eth_temp(k) = tmk_new(k,i,j) * (1000.d0/prs_new(k,i,j))**& - (GAMMA*(1.d0 + GAMMAMD*(MAX(qvp_new(k,i,j), 1.d-15))))*& + eth_temp(k) = tmk_new(k,i,j) * (1000.D0/prs_new(k,i,j))**& + (GAMMA*(1.D0 + GAMMAMD*(MAX(qvp_new(k,i,j), 1.d-15))))*& EXP((THTECON1/tlcl - THTECON2)*(MAX(qvp_new(k,i,j), 1.d-15))*& - (1.d0 + THTECON3*(MAX(qvp_new(k,i,j), 1.d-15)))) + (1.D0 + THTECON3*(MAX(qvp_new(k,i,j), 1.d-15)))) END IF END DO klev = mkzh @@ -755,8 +755,8 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& ! Establish average properties of that parcel ! (over depth of approximately davg meters) - !davg = 500.d0 - pavg = 500.d0 * prs_new(kpar1,i,j)*& + !davg = 500.D0 + pavg = 500.D0 * prs_new(kpar1,i,j)*& G/(RD*tvirtual(tmk_new(kpar1,i,j), qvp_new(kpar1,i,j))) p2 = MIN(prs_new(kpar1,i,j)+.5d0*pavg, prsf(mkzh,i,j)) p1 = p2 - pavg diff --git a/fortran/wrf_fctt.f90 b/fortran/wrf_fctt.f90 index 97da1b9..a1c97da 100644 --- a/fortran/wrf_fctt.f90 +++ b/fortran/wrf_fctt.f90 @@ -18,23 +18,15 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew ! 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) :: 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 + REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix + REAL(KIND=8) :: arg2, agl_hgt, vt REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf REAL(KIND=8) :: p1, p2 - !mjx = ew - !miy = ns - !mkzh = nz + !$OMP PARALLEL - prsctt = 0 ! removes the warning - -! Calculate the surface pressure + ! Calculate the surface pressure + !$OMP DO COLLAPSE(2) DO j=1,ns DO i=1,ew ratmix = .001D0*qvp(i,j,1) @@ -46,7 +38,9 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew pf(i,j,nz) = prs(i,j,1)*(vt/(vt + USSALR*(agl_hgt)))**(arg1) END DO END DO + !$OMP END DO + !$OMP DO COLLAPSE(3) DO k=1,nz-1 DO j=1,ns DO i=1,ew @@ -55,26 +49,27 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew END DO END DO END DO + !$OMP END DO + !$OMP DO COLLAPSE(2) PRIVATE(i, j, k, ripk, opdepthd, opdepthu, & + !$OMP prsctt, dp, p1, p2, fac, arg1) DO j=1,ns DO i=1,ew opdepthd = 0.D0 k = 0 + prsctt = 0 -! Integrate downward from model top, calculating path at full -! model vertical levels. - -!20 opdepthu=opdepthd + ! Integrate downward from model top, calculating path at full + ! model vertical levels. DO k=1, nz opdepthu = opdepthd - !k=k+1 ripk = nz - k + 1 - IF (k .EQ. 1) THEN - dp = 200.D0*(pf(i,j,1) - prs(i,j,nz)) ! should be in Pa - ELSE + IF (k .NE. 1) THEN dp = 100.D0*(pf(i,j,k) - pf(i,j,k-1)) ! should be in Pa + ELSE + dp = 200.D0*(pf(i,j,1) - prs(i,j,nz)) ! should be in Pa END IF IF (haveqci .EQ. 0) then @@ -89,7 +84,6 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew END IF IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN - !GOTO 20 CYCLE ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN @@ -111,15 +105,14 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew 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 END IF END DO END DO END DO -! 30 CONTINUE -! 40 CONTINUE -! 190 CONTINUE + !$OMP END DO + + !$OMP END PARALLEL RETURN END SUBROUTINE wrfcttcalc diff --git a/fortran/wrf_pvo.f90 b/fortran/wrf_pvo.f90 index 209a257..cef3ed9 100644 --- a/fortran/wrf_pvo.f90 +++ b/fortran/wrf_pvo.f90 @@ -23,29 +23,29 @@ SUBROUTINE DCOMPUTEABSVORT(av, u, v, msfu, msfv, msft, cor, dx, dy, nx, ny, nz,& REAL(KIND=8) :: dsy, dsx, dudy, dvdx, avort REAL(KIND=8) :: mm - ! PRINT*,'nx,ny,nz,nxp1,nyp1' - ! PRINT*,nx,ny,nz,nxp1,nyp1 + !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i, j, k, jp1, jm1, ip1, im1, & + !$OMP dsx, dsy, mm, dudy, dvdx, avort) DO k = 1,nz DO j = 1,ny - jp1 = MIN(j+1, ny) - jm1 = MAX(j-1, 1) DO i = 1,nx + jp1 = MIN(j+1, ny) + jm1 = MAX(j-1, 1) ip1 = MIN(i+1, nx) im1 = MAX(i-1, 1) - ! PRINT *,jp1,jm1,ip1,im1 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) - & 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) - & 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 + END DO END DO END DO + !$OMP END PARALLEL DO RETURN @@ -80,22 +80,23 @@ SUBROUTINE DCOMPUTEPV(pv, u, v, theta, prs, msfu, msfv, msft, cor, dx, dy, nx, & REAL(KIND=8) :: dsy, dsx, dp, dudy, dvdx, dudp, dvdp, dthdp, avort REAL(KIND=8) :: dthdx, dthdy, mm - ! PRINT*,'nx,ny,nz,nxp1,nyp1' - ! PRINT*,nx,ny,nz,nxp1,nyp1 + !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i, j, k, kp1, km1, jp1, jm1, ip1, & + !$OMP im1, dsx, dsy, mm, dudy, dvdx, avort, & + !$OMP dp, dudp, dvdp, dthdp, dthdx, dthdy) DO k = 1,nz - kp1 = MIN(k+1, nz) - km1 = MAX(k-1, 1) DO J = 1,ny - jp1 = MIN(j+1, ny) - jm1 = MAX(j-1, 1) DO i = 1,nx + kp1 = MIN(k+1, nz) + km1 = MAX(k-1, 1) + jp1 = MIN(j+1, ny) + jm1 = MAX(j-1, 1) ip1 = MIN(i+1, nx) im1 = MAX(i-1, 1) - ! PRINT *,jp1,jm1,ip1,im1 + 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) - & 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) - & @@ -108,14 +109,11 @@ SUBROUTINE DCOMPUTEPV(pv, u, v, theta, prs, msfu, msfv, msft, cor, dx, dy, nx, & 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) - ! endif pv(i,j,k) = pv(i,j,k)*1.D2 END DO END DO END DO + !$OMP END PARALLEL DO RETURN diff --git a/fortran/wrf_user.f90 b/fortran/wrf_user.f90 index 78433e3..618b864 100644 --- a/fortran/wrf_user.f90 +++ b/fortran/wrf_user.f90 @@ -327,9 +327,9 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & INTEGER :: klo, khi INTEGER :: errcnt, bad_i, bad_j, bad_sfp - !REAL(KIND=8) :: plo, phi, tlo, thi, zlo, zhi - !REAL(KIND=8) :: p_at_pconst, t_at_pconst, z_at_pconst - !REAL(KIND=8) :: z_half_lowest + REAL(KIND=8) :: plo, phi, tlo, thi, zlo, zhi + REAL(KIND=8) :: p_at_pconst, t_at_pconst, z_at_pconst + REAL(KIND=8) :: z_half_lowest LOGICAL :: l1, l2, l3, found @@ -399,32 +399,18 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, & !$OMP END CRITICAL END IF - ! This is the readable version of the code below. Don't delete this! - !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)) - !zlo = z(i,j,klo) - !zhi = z(i,j,khi) - !p_at_pconst = p(i,j,1) - PCONST - !t_at_pconst = thi - (thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi) - !z_at_pconst = zhi - (zhi-zlo)*LOG(p_at_pconst/phi)*LOG(plo/phi) + 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)) + zlo = z(i,j,klo) + zhi = z(i,j,khi) + p_at_pconst = p(i,j,1) - PCONST + t_at_pconst = thi - (thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi) + z_at_pconst = zhi - (zhi-zlo)*LOG(p_at_pconst/phi)*LOG(plo/phi) ! - !t_surf(i,j) = t_at_pconst * (p(i,j,1)/p_at_pconst)**(USSALR*RD/G) - !t_sea_level(i,j) = t_at_pconst + USSALR*z_at_pconst - - ! The same code as above with temporaries removed to improve vectorization - t_surf(i,j) = ((t(i,j,khi)*(1.D0 + 0.608D0*q(i,j,khi))) - & - ((t(i,j,khi)*(1.D0 + 0.608D0*q(i,j,khi)))-(t(i,j,klo)*& - (1.D0 + 0.608D0*q(i,j,klo))))*LOG((p(i,j,1) - PCONST)/p(i,j,khi))& - *LOG(p(i,j,klo)/p(i,j,khi))) * (p(i,j,1)/(p(i,j,1) - PCONST))**(USSALR*RD/G) - - t_sea_level(i,j) = t(i,j,khi)*(1.D0 + 0.608D0*q(i,j,khi)) - & - ((t(i,j,khi)*(1.D0 + 0.608D0*q(i,j,khi)))-(t(i,j,klo)*& - (1.D0 + 0.608D0*q(i,j,klo))))*LOG((p(i,j,1) - PCONST)/& - p(i,j,khi))*LOG(p(i,j,klo)/p(i,j,khi)) + & - USSALR*(z(i,j,khi) - (z(i,j,khi)-z(i,j,klo))*& - LOG((p(i,j,1) - PCONST)/p(i,j,khi))*LOG(p(i,j,klo)/p(i,j,khi))) + t_surf(i,j) = t_at_pconst * (p(i,j,1)/p_at_pconst)**(USSALR*RD/G) + t_sea_level(i,j) = t_at_pconst + USSALR*z_at_pconst END DO END DO !$OMP END PARALLEL DO @@ -581,16 +567,18 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) INTEGER :: i, j, iter + !$OMP PARALLEL + DO iter=1,it - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j=1,ny DO i = 1,nx b(i,j) = a(i,j) END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j=2,ny-1 DO i=1,nx IF (b(i,j-1) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & @@ -601,9 +589,9 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) END IF END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j=1,ny DO i=2,nx-1 IF (b(i-1,j) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & @@ -614,7 +602,7 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) END IF END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO ! do j=1,ny ! do i=1,nx ! b(i,j) = a(i,j) @@ -632,6 +620,8 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing) ! enddo END DO + !$OMP END PARALLEL + RETURN END SUBROUTINE FILTER2D @@ -768,7 +758,9 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & ! msg stands for missing value in this code ! WRITE (6,FMT=*) ' in compute_uvmet ',NX,NY,NXP1,NYP1,ISTAG - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP PARALLEL + + !$OMP DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx @@ -790,23 +782,23 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO ! Note: Intentionally removed as many IF statements as possible from loops ! to improve vectorization. IF (istag .EQ. 0) THEN ! Not staggered IF (.NOT. is_msg_val) THEN ! No missing values used - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx uvmet(i,j,1) = v(i,j)*longcb(i,j) + u(i,j)*longca(i,j) uvmet(i,j,2) = v(i,j)*longca(i,j) - u(i,j)*longcb(i,j) END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO ELSE ! Missing values used - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx IF ((u(i,j) .NE. umsg .AND. v(i,j) .NE. vmsg)) THEN @@ -818,14 +810,14 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & END IF END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO END IF ELSE ! Staggered IF (.NOT. is_msg_val) THEN ! No missing values used - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx - ! This is the more readable version. Do not delete. + ! This is the more readable version. !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) @@ -838,13 +830,13 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO ELSE ! Missing values used - !$OMP PARALLEL DO COLLAPSE(2) + !$OMP DO COLLAPSE(2) DO j = 1,ny DO i = 1,nx IF (u(i,j) .NE. umsg .AND. v(i,j) .NE. vmsg .AND. u(i+1,j) .NE. umsg .AND. v(i,j+1) .NE. vmsg) THEN - ! This is the more readable version. Do not delete. + ! This is the more readable version. !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) @@ -860,10 +852,12 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, & END IF END DO END DO - !$OMP END PARALLEL DO + !$OMP END DO END IF END IF + !$OMP END PARALLEL + RETURN END SUBROUTINE DCOMPUTEUVMET diff --git a/fortran/wrf_vinterp.f90 b/fortran/wrf_vinterp.f90 index ee1835f..7f7c2f1 100644 --- a/fortran/wrf_vinterp.f90 +++ b/fortran/wrf_vinterp.f90 @@ -20,24 +20,24 @@ SUBROUTINE wrf_monotonic(out, in, lvprs, cor, idir, delta, ew, ns, nz, icorsw) !NCLEND - INTEGER :: i, j, k, ripk, k300 + INTEGER :: i, j, k, k300 - k300 = 1 ! removes the warning - - !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, ripk) FIRSTPRIVATE(k300) + !$OMP PARALLEL DO COLLAPSE(2) DO j=1,ns DO i=1,ew + IF (icorsw .EQ. 1 .AND. cor(i,j) .LT. 0.) THEN DO k=1,nz in(i,j,k) = -in(i,j,k) END DO END IF + k300 = nz + ! First find k index that is at or below (height-wise) ! the 300 hPa level. - DO k = 1,nz - ripk = nz-k+1 - IF (lvprs(i,j,k) .LE. 300.D0) THEN + DO k = nz,1,-1 + IF (lvprs(i,j,k) .GE. 300.D0) THEN k300 = k EXIT END IF @@ -184,18 +184,13 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& cvcord = 't' END IF - !miy = ns - !mjx = ew - !njx = ew - !niy = ns - -!$OMP PARALLEL DO + !$OMP PARALLEL DO COLLAPSE(2) DO j = 1,ns DO i = 1,ew tempout(i,j) = rmsg END DO END DO -!$OMP END PARALLEL DO + !$OMP END PARALLEL DO DO nreqlvs = 1,numlevels IF (cvcord .EQ. 'z') THEN @@ -257,6 +252,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& ifound = 1 END IF END IF + EXIT END IF END DO !end for the k loop @@ -415,6 +411,8 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& END DO !$OMP END PARALLEL DO + + IF (log_errcnt > 0) THEN errstat = ALGERR WRITE(errmsg, *) "Pres=0. Unable to take log of 0." @@ -427,7 +425,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& RETURN END IF - !$OMP PARALLEL DO + !$OMP PARALLEL DO COLLAPSE(2) DO j = 1,ns DO i = 1,ew dataout(i,j,nreqlvs) = tempout(i,j)