nilyin
/
wrf-python
forked from 3rdparty/wrf-python
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Updated OpenMP directives.

lon0
Bill Ladwig 8 years ago
parent
4135d06edb
commit
ea7aba57b2
6 changed files with 126 additions and 138 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 13
      fortran/calc_uh.f90
  2. 58
      fortran/rip_cape.f90
  3. 45
      fortran/wrf_fctt.f90
  4. 34
      fortran/wrf_pvo.f90
  5. 86
      fortran/wrf_user.f90
  6. 28
      fortran/wrf_vinterp.f90

13
fortran/calc_uh.f90
Unescape View File

@ -61,7 +61,10 @@ SUBROUTINE DCALCUH(nx, ny, nz, nzp1, zp, mapfct, dx, dy, uhmnhgt, uhmxhgt, us, &
twodx = 2.0*dx twodx = 2.0*dx
twody = 2.0*dy twody = 2.0*dy
!$OMP PARALLEL DO COLLAPSE(3)
!$OMP PARALLEL
!$OMP DO COLLAPSE(3)
DO k=2,nz-2 DO k=2,nz-2
DO j=2,ny-1 DO j=2,ny-1
DO i=2,nx-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 END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
! Integrate over depth uhminhgt to uhmxhgt AGL ! Integrate over depth uhminhgt to uhmxhgt AGL
! !
! WRITE(6,'(a,f12.1,a,f12.1,a)') & ! WRITE(6,'(a,f12.1,a,f12.1,a)') &
! 'Calculating UH from ',uhmnhgt,' to ',uhmxhgt,' m AGL' ! '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) !$OMP wgtlw, wbot, wtop, wsum, wmean, sum, helbot, heltop)
DO j=2,ny-2 DO j=2,ny-2
DO i=2,nx-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 IF
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
!$OMP END PARALLEL
uh = uh*1000. ! Scale according to Kain et al. (2008) uh = uh*1000. ! Scale according to Kain et al. (2008)

58
fortran/rip_cape.f90
Unescape View File

@ -85,11 +85,11 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam
rang1 = h1 - l1 rang1 = h1 - l1
mid1 = (h1 + l1) / 2 mid1 = (h1 + l1) / 2
DO WHILE(rang1 .GT. 1) DO WHILE(rang1 .GT. 1)
if(thte .GE. psadithte(mid1)) then IF (thte .GE. psadithte(mid1)) THEN
l1 = mid1 l1 = mid1
else ELSE
h1 = mid1 h1 = mid1
end if END IF
rang1 = h1 - l1 rang1 = h1 - l1
mid1 = (h1 + l1) / 2 mid1 = (h1 + l1) / 2
END DO END DO
@ -103,17 +103,17 @@ REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gam
! END IF ! END IF
! END DO ! END DO
ip = -1 ip = -1
l2 = 1 l2 = 1
h2 = 149 h2 = 149
rang2 = h2 - l2 rang2 = h2 - l2
mid2 = (h2 + l2) / 2 mid2 = (h2 + l2) / 2
DO WHILE(rang2 .GT. 1) DO WHILE(rang2 .GT. 1)
if(prs .LE. psadiprs(mid2)) then IF (prs .LE. psadiprs(mid2)) THEN
l2 = mid2 l2 = mid2
else ELSE
h2 = mid2 h2 = mid2
end if END IF
rang2 = h2 - l2 rang2 = h2 - l2
mid2 = (h2 + l2) / 2 mid2 = (h2 + l2) / 2
END DO END DO
@ -386,8 +386,8 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
!$OMP i,j,k,kpar) !$OMP i,j,k,kpar)
DO j = 1,mjy DO j = 1,mjy
DO i = 1,mix DO i = 1,mix
cape(i,j,1) = 0.d0 cape(i,j,1) = 0.D0
cin(i,j,1) = 0.d0 cin(i,j,1) = 0.D0
!!$OMP SIMD !!$OMP SIMD
DO kpar = 2, mkzh DO kpar = 2, mkzh
@ -421,7 +421,7 @@ SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
klcl = 1 klcl = 1
END IF 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 DO k = kpar,1,-1
! For arrays that go bottom to top ! For arrays that go bottom to top
kk = kk + 1 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) ! 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 j = 1,mjy
DO i = 1,mix DO i = 1,mix
DO k = 1,mkzh DO k = 1,mkzh
@ -695,7 +695,7 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
END DO END DO
END DO 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 ! 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() !nthreads = omp_get_num_threads()
!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, & !$OMP PARALLEL DO COLLAPSE(2) PRIVATE(tlcl, ethpari, &
!$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, & !$OMP zlcl, kk, ilcl, klcl, tmklift, tvenv, tvlift, ghtlift, &
!$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, & !$OMP facden, tmkenv, qvpenv, eslift, qvplift, buoy, benamin, &
!$OMP benaccum, zrel, kmax, dz, elfound, & !$OMP benaccum, zrel, kmax, dz, elfound, &
!$OMP kel, klfc, pavg, p2, p1, totthe, totqvp, totprs, & !$OMP kel, klfc, pavg, p2, p1, totthe, totqvp, totprs, &
!$OMP i,j,k,kpar, kpar1, kpar2, qvppari, tmkpari,p, pup, pdn, q, th, & !$OMP i,j,k,kpar, kpar1, kpar2, qvppari, tmkpari,p, pup, pdn, q, th, &
!$OMP pp1, pp2, ethmax, eth_temp, klev) !$OMP pp1, pp2, ethmax, eth_temp, klev)
DO j = 1,mjy DO j = 1,mjy
DO i = 1,mix DO i = 1,mix
cape(i,j,1) = 0.d0 cape(i,j,1) = 0.D0
cin(i,j,1) = 0.d0 cin(i,j,1) = 0.D0
! find parcel with max theta-e in lowest 3 km agl. ! find parcel with max theta-e in lowest 3 km agl.
ethmax = -1.d0 ethmax = -1.D0
eth_temp = -1.d0 eth_temp = -1.D0
DO k = 1, mkzh 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/& 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)+& (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 TLCLC4
eth_temp(k) = tmk_new(k,i,j) * (1000.d0/prs_new(k,i,j))**& 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))))*& (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))*& 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 IF
END DO END DO
klev = mkzh klev = mkzh
@ -755,8 +755,8 @@ SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,&
! Establish average properties of that parcel ! Establish average properties of that parcel
! (over depth of approximately davg meters) ! (over depth of approximately davg meters)
!davg = 500.d0 !davg = 500.D0
pavg = 500.d0 * prs_new(kpar1,i,j)*& pavg = 500.D0 * prs_new(kpar1,i,j)*&
G/(RD*tvirtual(tmk_new(kpar1,i,j), qvp_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)) p2 = MIN(prs_new(kpar1,i,j)+.5d0*pavg, prsf(mkzh,i,j))
p1 = p2 - pavg p1 = p2 - pavg

45
fortran/wrf_fctt.f90
Unescape View File

@ -18,23 +18,15 @@ SUBROUTINE wrfcttcalc(prs, tk, qci, qcw, qvp, ght, ter, ctt, haveqci, nz, ns, ew
! LOCAL VARIABLES ! LOCAL VARIABLES
INTEGER i,j,k,ripk INTEGER i,j,k,ripk
!INTEGER :: mjx,miy,mkzh REAL(KIND=8) :: opdepthu, opdepthd, dp, arg1, fac, prsctt, ratmix
REAL(KIND=8) :: vt, opdepthu, opdepthd, dp REAL(KIND=8) :: arg2, agl_hgt, vt
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), DIMENSION(ew,ns,nz) :: pf REAL(KIND=8), DIMENSION(ew,ns,nz) :: pf
REAL(KIND=8) :: p1, p2 REAL(KIND=8) :: p1, p2
!mjx = ew !$OMP PARALLEL
!miy = ns
!mkzh = nz
prsctt = 0 ! removes the warning ! Calculate the surface pressure
!$OMP DO COLLAPSE(2)
! Calculate the surface pressure
DO j=1,ns DO j=1,ns
DO i=1,ew DO i=1,ew
ratmix = .001D0*qvp(i,j,1) 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) pf(i,j,nz) = prs(i,j,1)*(vt/(vt + USSALR*(agl_hgt)))**(arg1)
END DO END DO
END DO END DO
!$OMP END DO
!$OMP DO COLLAPSE(3)
DO k=1,nz-1 DO k=1,nz-1
DO j=1,ns DO j=1,ns
DO i=1,ew 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 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 j=1,ns
DO i=1,ew DO i=1,ew
opdepthd = 0.D0 opdepthd = 0.D0
k = 0 k = 0
prsctt = 0
! Integrate downward from model top, calculating path at full ! Integrate downward from model top, calculating path at full
! model vertical levels. ! model vertical levels.
!20 opdepthu=opdepthd
DO k=1, nz DO k=1, nz
opdepthu = opdepthd opdepthu = opdepthd
!k=k+1
ripk = nz - k + 1 ripk = nz - k + 1
IF (k .EQ. 1) THEN IF (k .NE. 1) THEN
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
ELSE
dp = 200.D0*(pf(i,j,1) - prs(i,j,nz)) ! should be in Pa
END IF END IF
IF (haveqci .EQ. 0) then 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 END IF
IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN IF (opdepthd .LT. 1. .AND. k .LT. nz) THEN
!GOTO 20
CYCLE CYCLE
ELSE IF (opdepthd .LT. 1. .AND. k .EQ. nz) THEN 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) fac = (prsctt - p1)/(p2 - p1)
arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL arg1 = fac*(tk(i,j,ripk) - tk(i,j,ripk+1)) - CELKEL
ctt(i,j) = tk(i,j,ripk+1) + arg1 ctt(i,j) = tk(i,j,ripk+1) + arg1
!GOTO 40
EXIT EXIT
END IF END IF
END DO END DO
END DO END DO
END DO END DO
! 30 CONTINUE !$OMP END DO
! 40 CONTINUE
! 190 CONTINUE !$OMP END PARALLEL
RETURN RETURN
END SUBROUTINE wrfcttcalc END SUBROUTINE wrfcttcalc

34
fortran/wrf_pvo.f90
Unescape View File

@ -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) :: dsy, dsx, dudy, dvdx, avort
REAL(KIND=8) :: mm REAL(KIND=8) :: mm
! PRINT*,'nx,ny,nz,nxp1,nyp1' !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i, j, k, jp1, jm1, ip1, im1, &
! PRINT*,nx,ny,nz,nxp1,nyp1 !$OMP dsx, dsy, mm, dudy, dvdx, avort)
DO k = 1,nz DO k = 1,nz
DO j = 1,ny DO j = 1,ny
jp1 = MIN(j+1, ny)
jm1 = MAX(j-1, 1)
DO i = 1,nx DO i = 1,nx
jp1 = MIN(j+1, ny)
jm1 = MAX(j-1, 1)
ip1 = MIN(i+1, nx) ip1 = MIN(i+1, nx)
im1 = MAX(i-1, 1) im1 = MAX(i-1, 1)
! PRINT *,jp1,jm1,ip1,im1
dsx = (ip1 - im1) * dx dsx = (ip1 - im1) * dx
dsy = (jp1 - jm1) * dy dsy = (jp1 - jm1) * dy
mm = msft(i,j)*msft(i,j) 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 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 v(im1,j,k)/msfv(im1,j) - v(im1,j+1,k)/msfv(im1,j+1))/dsx*mm
avort = dvdx - dudy + cor(i,j) avort = dvdx - dudy + cor(i,j)
av(i,j,k) = avort*1.D5 av(i,j,k) = avort*1.D5
END DO END DO
END DO END DO
END DO END DO
!$OMP END PARALLEL DO
RETURN 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) :: dsy, dsx, dp, dudy, dvdx, dudp, dvdp, dthdp, avort
REAL(KIND=8) :: dthdx, dthdy, mm REAL(KIND=8) :: dthdx, dthdy, mm
! PRINT*,'nx,ny,nz,nxp1,nyp1' !$OMP PARALLEL DO COLLAPSE(3) PRIVATE(i, j, k, kp1, km1, jp1, jm1, ip1, &
! PRINT*,nx,ny,nz,nxp1,nyp1 !$OMP im1, dsx, dsy, mm, dudy, dvdx, avort, &
!$OMP dp, dudp, dvdp, dthdp, dthdx, dthdy)
DO k = 1,nz DO k = 1,nz
kp1 = MIN(k+1, nz)
km1 = MAX(k-1, 1)
DO J = 1,ny DO J = 1,ny
jp1 = MIN(j+1, ny)
jm1 = MAX(j-1, 1)
DO i = 1,nx 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) ip1 = MIN(i+1, nx)
im1 = MAX(i-1, 1) im1 = MAX(i-1, 1)
! PRINT *,jp1,jm1,ip1,im1
dsx = (ip1 - im1)*dx dsx = (ip1 - im1)*dx
dsy = (jp1 - jm1)*dy dsy = (jp1 - jm1)*dy
mm = msft(i,j)*msft(i,j) 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 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) - &
@ -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) 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) 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 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 pv(i,j,k) = pv(i,j,k)*1.D2
END DO END DO
END DO END DO
END DO END DO
!$OMP END PARALLEL DO
RETURN RETURN

86
fortran/wrf_user.f90
Unescape View File

@ -327,9 +327,9 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, &
INTEGER :: klo, khi INTEGER :: klo, khi
INTEGER :: errcnt, bad_i, bad_j, bad_sfp INTEGER :: errcnt, bad_i, bad_j, bad_sfp
!REAL(KIND=8) :: plo, phi, tlo, thi, zlo, zhi REAL(KIND=8) :: plo, phi, tlo, thi, zlo, zhi
!REAL(KIND=8) :: p_at_pconst, t_at_pconst, z_at_pconst REAL(KIND=8) :: p_at_pconst, t_at_pconst, z_at_pconst
!REAL(KIND=8) :: z_half_lowest REAL(KIND=8) :: z_half_lowest
LOGICAL :: l1, l2, l3, found LOGICAL :: l1, l2, l3, found
@ -399,32 +399,18 @@ SUBROUTINE DCOMPUTESEAPRS(nx, ny, nz, z, t, p, q, sea_level_pressure, &
!$OMP END CRITICAL !$OMP END CRITICAL
END IF END IF
! This is the readable version of the code below. Don't delete this! plo = p(i,j,klo)
!plo = p(i,j,klo) phi = p(i,j,khi)
!phi = p(i,j,khi) tlo = t(i,j,klo)*(1.D0 + 0.608D0*q(i,j,klo))
!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))
!thi = t(i,j,khi)*(1.D0 + 0.608D0*q(i,j,khi)) zlo = z(i,j,klo)
!zlo = z(i,j,klo) zhi = z(i,j,khi)
!zhi = z(i,j,khi) p_at_pconst = p(i,j,1) - PCONST
!p_at_pconst = p(i,j,1) - PCONST t_at_pconst = thi - (thi-tlo)*LOG(p_at_pconst/phi)*LOG(plo/phi)
!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)
!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_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 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)))
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -581,16 +567,18 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing)
INTEGER :: i, j, iter INTEGER :: i, j, iter
!$OMP PARALLEL
DO iter=1,it DO iter=1,it
!$OMP PARALLEL DO COLLAPSE(2) !$OMP DO COLLAPSE(2)
DO j=1,ny DO j=1,ny
DO i = 1,nx DO i = 1,nx
b(i,j) = a(i,j) b(i,j) = a(i,j)
END DO END DO
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 j=2,ny-1
DO i=1,nx DO i=1,nx
IF (b(i,j-1) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & 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 IF
END DO END DO
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 j=1,ny
DO i=2,nx-1 DO i=2,nx-1
IF (b(i-1,j) .EQ. missing .OR. b(i,j) .EQ. missing .OR. & 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 IF
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
! do j=1,ny ! do j=1,ny
! do i=1,nx ! do i=1,nx
! b(i,j) = a(i,j) ! b(i,j) = a(i,j)
@ -632,6 +620,8 @@ SUBROUTINE FILTER2D(a, b, nx, ny, it, missing)
! enddo ! enddo
END DO END DO
!$OMP END PARALLEL
RETURN RETURN
END SUBROUTINE FILTER2D END SUBROUTINE FILTER2D
@ -768,7 +758,9 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, &
! msg stands for missing value in this code ! msg stands for missing value in this code
! WRITE (6,FMT=*) ' in compute_uvmet ',NX,NY,NXP1,NYP1,ISTAG ! 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 j = 1,ny
DO i = 1,nx DO i = 1,nx
@ -790,23 +782,23 @@ SUBROUTINE DCOMPUTEUVMET(u, v, uvmet, longca,longcb,flong,flat, &
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
! Note: Intentionally removed as many IF statements as possible from loops ! Note: Intentionally removed as many IF statements as possible from loops
! to improve vectorization. ! to improve vectorization.
IF (istag .EQ. 0) THEN ! Not staggered IF (istag .EQ. 0) THEN ! Not staggered
IF (.NOT. is_msg_val) THEN ! No missing values used IF (.NOT. is_msg_val) THEN ! No missing values used
!$OMP PARALLEL DO COLLAPSE(2) !$OMP DO COLLAPSE(2)
DO j = 1,ny DO j = 1,ny
DO i = 1,nx DO i = 1,nx
uvmet(i,j,1) = v(i,j)*longcb(i,j) + u(i,j)*longca(i,j) 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) uvmet(i,j,2) = v(i,j)*longca(i,j) - u(i,j)*longcb(i,j)
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
ELSE ! Missing values used ELSE ! Missing values used
!$OMP PARALLEL DO COLLAPSE(2) !$OMP DO COLLAPSE(2)
DO j = 1,ny DO j = 1,ny
DO i = 1,nx DO i = 1,nx
IF ((u(i,j) .NE. umsg .AND. v(i,j) .NE. vmsg)) THEN 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 IF
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
END IF END IF
ELSE ! Staggered ELSE ! Staggered
IF (.NOT. is_msg_val) THEN ! No missing values used IF (.NOT. is_msg_val) THEN ! No missing values used
!$OMP PARALLEL DO COLLAPSE(2) !$OMP DO COLLAPSE(2)
DO j = 1,ny DO j = 1,ny
DO i = 1,nx 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)) !uk = 0.5D0*(u(i,j) + u(i+1,j))
!vk = 0.5D0*(v(i,j) + v(i,j+1)) !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,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
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
ELSE ! Missing values used ELSE ! Missing values used
!$OMP PARALLEL DO COLLAPSE(2) !$OMP DO COLLAPSE(2)
DO j = 1,ny DO j = 1,ny
DO i = 1,nx 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 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)) !uk = 0.5D0*(u(i,j) + u(i+1,j))
!vk = 0.5D0*(v(i,j) + v(i,j+1)) !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,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 IF
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END DO
END IF END IF
END IF END IF
!$OMP END PARALLEL
RETURN RETURN
END SUBROUTINE DCOMPUTEUVMET END SUBROUTINE DCOMPUTEUVMET

28
fortran/wrf_vinterp.f90
Unescape View File

@ -20,24 +20,24 @@ SUBROUTINE wrf_monotonic(out, in, lvprs, cor, idir, delta, ew, ns, nz, icorsw)
!NCLEND !NCLEND
INTEGER :: i, j, k, ripk, k300 INTEGER :: i, j, k, k300
k300 = 1 ! removes the warning !$OMP PARALLEL DO COLLAPSE(2)
!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, ripk) FIRSTPRIVATE(k300)
DO j=1,ns DO j=1,ns
DO i=1,ew DO i=1,ew
IF (icorsw .EQ. 1 .AND. cor(i,j) .LT. 0.) THEN IF (icorsw .EQ. 1 .AND. cor(i,j) .LT. 0.) THEN
DO k=1,nz DO k=1,nz
in(i,j,k) = -in(i,j,k) in(i,j,k) = -in(i,j,k)
END DO END DO
END IF END IF
k300 = nz
! First find k index that is at or below (height-wise) ! First find k index that is at or below (height-wise)
! the 300 hPa level. ! the 300 hPa level.
DO k = 1,nz DO k = nz,1,-1
ripk = nz-k+1 IF (lvprs(i,j,k) .GE. 300.D0) THEN
IF (lvprs(i,j,k) .LE. 300.D0) THEN
k300 = k k300 = k
EXIT EXIT
END IF END IF
@ -184,18 +184,13 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
cvcord = 't' cvcord = 't'
END IF END IF
!miy = ns !$OMP PARALLEL DO COLLAPSE(2)
!mjx = ew
!njx = ew
!niy = ns
!$OMP PARALLEL DO
DO j = 1,ns DO j = 1,ns
DO i = 1,ew DO i = 1,ew
tempout(i,j) = rmsg tempout(i,j) = rmsg
END DO END DO
END DO END DO
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
DO nreqlvs = 1,numlevels DO nreqlvs = 1,numlevels
IF (cvcord .EQ. 'z') THEN IF (cvcord .EQ. 'z') THEN
@ -257,6 +252,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
ifound = 1 ifound = 1
END IF END IF
END IF END IF
EXIT EXIT
END IF END IF
END DO !end for the k loop END DO !end for the k loop
@ -415,6 +411,8 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
END DO END DO
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
IF (log_errcnt > 0) THEN IF (log_errcnt > 0) THEN
errstat = ALGERR errstat = ALGERR
WRITE(errmsg, *) "Pres=0. Unable to take log of 0." 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 RETURN
END IF END IF
!$OMP PARALLEL DO !$OMP PARALLEL DO COLLAPSE(2)
DO j = 1,ns DO j = 1,ns
DO i = 1,ew DO i = 1,ew
dataout(i,j,nreqlvs) = tempout(i,j) dataout(i,j,nreqlvs) = tempout(i,j)

Write Preview
Loading…
Cancel
Save