Added OpenMP directives for wrf_vintrp

fortran/wrf_vinterp.f90

@@ -12,10 +12,11 @@ SUBROUTINE wrf_monotonic(out, in, lvprs, cor, idir, delta, ew, ns, nz, icorsw)
 
     INTEGER, INTENT(IN) :: idir, ew, ns, nz, icorsw
     REAL(KIND=8), INTENT(IN) :: delta
-    REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(INOUT) :: in
+    REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: in
     REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(OUT) :: out
-    REAL(KIND=8), DIMENSION(ew,ns,nz) :: lvprs
-    REAL(KIND=8), DIMENSION(ew,ns) :: cor
+    REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: lvprs
+    REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: cor
+
 
 !NCLEND
 
@@ -65,7 +66,7 @@ END SUBROUTINE wrf_monotonic
 
 
 !NCLFORTSTART
-FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat, errmsg)
+FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat)
     USE wrf_constants, ONLY : ALGERR, SCLHT
 
     IMPLICIT NONE
@@ -75,7 +76,6 @@ FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat, errms
     INTEGER, INTENT(IN) :: icase
     REAL(KIND=8), INTENT(IN) :: wvalp0, wvalp1, vlev, vcp0, vcp1
     INTEGER, INTENT(INOUT) :: errstat
-    CHARACTER(LEN=*), INTENT(INOUT) :: errmsg
     REAL(KIND=8) :: wrf_intrp_value
 
 !NCLEND
@@ -83,10 +83,6 @@ FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat, errms
     REAL(KIND=8) :: valp0, valp1, rvalue
     REAL(KIND=8) :: chkdiff
 
-    !REAL(KIND=8), PARAMETER :: RGAS=287.04d0
-    !REAL(KIND=8), PARAMETER :: USSALR=0.0065d0
-    !REAL(KIND=8), PARAMETER :: SCLHT=RGAS*256.d0/9.81d0
-
     errstat = 0
 
     valp0 = wvalp0
@@ -99,7 +95,7 @@ FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat, errms
     chkdiff = vcp1 - vcp0
     IF(chkdiff .EQ. 0) THEN
         errstat = ALGERR
-        errmsg = "bad difference in vcp's"
+        !errmsg = "bad difference in vcp's"
         wrf_intrp_value = 0
         RETURN
         !PRINT *,"bad difference in vcp's"
@@ -129,6 +125,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
     USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, &
                           TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, &
                           CELKEL, EPS, USSALR
+    USE omp_lib
 
     IMPLICIT NONE
 
@@ -152,6 +149,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
     INTEGER :: nreqlvs, ripk !njx,niy
     INTEGER :: i, j, k, kupper !itriv
     INTEGER :: ifound, isign !miy,mjx
+    INTEGER :: log_errcnt, interp_errcnt, interp_errstat
     REAL(KIND=8), DIMENSION(ew,ns) :: tempout
     REAL(KIND=8) :: rlevel, vlev, diff
     REAL(KIND=8) :: tmpvlev
@@ -165,27 +163,7 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
     REAL(KIND=8) :: pbot, zbot, tbotextrap, e
     REAL(KIND=8) :: tlcl, gammam
     CHARACTER(LEN=1) :: cvcord
-
-    !REAL(KIND=8), PARAMETER :: RGAS    = 287.04d0     !J/K/kg
-    !REAL(KIND=8), PARAMETER :: RGASMD  = .608d0
-    !REAL(KIND=8), PARAMETER :: USSALR  = .0065d0      ! deg C per m
-    !REAL(KIND=8), PARAMETER :: SCLHT   = RGAS*256.d0/9.81d0
-    !REAL(KIND=8), PARAMETER :: EPS     = 0.622d0
-    !REAL(KIND=8), PARAMETER :: RCONST  = -9.81d0/(RGAS * USSALR)
-    !REAL(KIND=8), PARAMETER :: EXPON   =  RGAS*USSALR/9.81d0
-    !REAL(KIND=8), PARAMETER :: EXPONI  =  1./EXPON
-    !REAL(KIND=8), PARAMETER :: TLCLC1   = 2840.d0
-    !REAL(KIND=8), PARAMETER :: TLCLC2   = 3.5d0
-    !REAL(KIND=8), PARAMETER :: TLCLC3   = 4.805d0
-    !REAL(KIND=8), PARAMETER :: TLCLC4   = 55.d0
-    !REAL(KIND=8), PARAMETER :: THTECON1 = 3376.d0 ! K
-    !REAL(KIND=8), PARAMETER :: THTECON2 = 2.54d0
-    !REAL(KIND=8), PARAMETER :: THTECON3 = 0.81d0
-    !REAL(KIND=8), PARAMETER :: CP       = 1004.d0
-    !REAL(KIND=8), PARAMETER :: CPMD     = 0.887d0
-    !REAL(KIND=8), PARAMETER :: GAMMA    = RGAS/CP
-    !REAL(KIND=8), PARAMETER :: GAMMAMD  = RGASMD-CPMD
-    !REAL(KIND=8), PARAMETER :: CELKEL   = 273.16d0
+    INTEGER :: thd
 
     ! Removes the warnings for uninitialized variables
     cvcord = ''
@@ -193,6 +171,9 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
     zlev = 0
     vlev = 0
     errstat = 0
+    interp_errcnt = 0
+    interp_errstat = 0
+    log_errcnt = 0
 
     IF (vcor .EQ. 1) THEN
         cvcord = 'p'
@@ -207,11 +188,13 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
     !njx = ew
     !niy = ns
 
+!$OMP PARALLEL DO
     DO j = 1,ns
         DO i = 1,ew
             tempout(i,j) = rmsg
         END DO
     END DO
+!$OMP END PARALLEL DO
 
     DO nreqlvs = 1,numlevels
         IF (cvcord .EQ. 'z') THEN
@@ -224,9 +207,16 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
             vlev = interp_levels(nreqlvs)
         END IF
 
+!$OMP PARALLEL DO COLLAPSE(2) PRIVATE(i, j, k, ifound, &
+!$OMP ripk, vcp1, vcp0, valp0, valp1, tmpvlev, interp_errstat, &
+!$OMP vclhsl, vctophsl, diff, isign, plhsl, zlhsl, ezlhsl, tlhsl, &
+!$OMP zsurf, qvapor, psurf, psurfsm, ezsurf, plev, ezlev, zlev, &
+!$OMP ptarget, dpmin, kupper, pbot, zbot, pratio, tbotextrap, &
+!$OMP vt, tlev, gammam, e, tlcl) REDUCTION (+:log_errcnt, interp_errcnt)
         DO j=1,ns
             DO i=1,ew
                 ! Get the interpolated value that is within the model domain
+                thd = omp_get_thread_num()
                 ifound = 0
                 DO k = 1,nz-1
                     ripk = nz-k+1
@@ -245,35 +235,33 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                             IF (logp .EQ. 1) THEN
                                 vcp1 = LOG(vcp1)
                                 vcp0 = LOG(vcp0)
-                                IF (vlev .EQ. 0.0D0) THEN
-                                    errstat = ALGERR
-                                    WRITE(errmsg, *) "Pres=0.  Unable to take log of 0."
-                                    RETURN
-                                    !PRINT *,"Pressure value = 0"
-                                    !PRINT *,"Unable to take log of 0"
-                                    !STOP
-                                END IF
+                                IF (vlev .NE. 0.0D0) THEN
                                     tmpvlev = LOG(vlev)
+                                ELSE
+                                    log_errcnt = log_errcnt + 1
+                                    tmpvlev = rmsg
+                                END IF
                             ELSE
                                 tmpvlev = vlev
                             END IF
+
+                            IF (tmpvlev .NE. rmsg) THEN
                                 tempout(i,j) = wrf_intrp_value(valp0, valp1, tmpvlev, vcp0, &
-                                                           vcp1, icase, errstat, errmsg)
-                            IF (errstat .NE. 0) THEN
-                                RETURN
+                                                               vcp1, icase, interp_errstat)
+
+                                IF (interp_errstat .NE. 0) THEN
+                                    tempout(i,j) = rmsg
+                                    interp_errcnt = interp_errcnt + 1
                                 END IF
 
-                            ! print *,"one ",i,j,tempout(i,j)
                                 ifound = 1
                             END IF
-                        !GOTO 115 ! EXIT
+                        END IF
                         EXIT
                     END IF
                 END DO !end for the k loop
- !115  CONTINUE
 
                 IF (ifound .EQ. 1) THEN !Grid point is in the model domain
-                    !GOTO 333 ! CYCLE
                     CYCLE
                 END IF
 
@@ -281,7 +269,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                 !all values above or below the model level to rmsg.
                 IF (extrap .EQ. 0) THEN
                     tempout(i,j) = rmsg
-                    !GOTO 333 ! CYCLE
                     CYCLE
                 END IF
 
@@ -296,8 +283,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                 IF (isign*vlev .GE. isign*vctophsl) THEN
                     ! Assign the highest model level to the out array
                     tempout(i,j) = datain(i,j,nz)
-                    ! print *,"at warn",i,j,tempout(i,j)
-                    !GOTO 333 ! CYCLE
                     CYCLE
                 END IF
 
@@ -307,7 +292,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
 
                 IF (datain(i,j,1) .EQ. rmsg) THEN
                     tempout(i,j) = rmsg
-                    !GOTO 333 ! CYCLE
                     CYCLE
                 END IF
 
@@ -351,7 +335,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                             zlev = -SCLHT*LOG(ezlev)
                             IF (icase .EQ. 2) THEN
                                 tempout(i,j) = zlev
-                                !GOTO 333 ! CYCLE
                                 CYCLE
                             END IF
 
@@ -362,7 +345,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                                    psurf + (ezsurf - ezlev)*plhsl)/(ezsurf - ezlhsl)
                             IF (icase .EQ. 1) THEN
                                 tempout(i,j) = plev
-                                !GOTO 333 ! CYCLE
                                 CYCLE
                             END IF
                         END IF
@@ -374,12 +356,11 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                             ripk = nz-k+1
                             dp = ABS((pres(i,j,ripk) * 0.01D0) - ptarget)
                             IF (dp .GT. dpmin) THEN
-                                !GOTO 334 ! EXIT
                                 EXIT
                             END IF
                             dpmin = MIN(dpmin, dp)
                         END DO
-         !334
+
                         kupper = k-1
 
                         ripk = nz - kupper + 1
@@ -394,7 +375,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                             zlev = zbot + vt/USSALR*(1. - (vlev/pbot)**EXPON)
                             IF (icase .EQ. 2) THEN
                                 tempout(i,j) = zlev
-                                !GOTO 333 ! CYCLE
                                 CYCLE
                             END IF
                         ELSE IF (cvcord .EQ. 'z') THEN
@@ -402,7 +382,6 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
                             plev = pbot*(1. + USSALR/vt*(zbot - zlev))**EXPONI
                             IF (icase .EQ. 1) THEN
                                 tempout(i,j) = plev
-                                !GOTO 333 ! CYCLE
                                 CYCLE
                             END IF
                         END IF
@@ -434,13 +413,27 @@ SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,&
 
             END DO
         END DO
+!$OMP END PARALLEL DO
+
+        IF (log_errcnt > 0) THEN
+            errstat = ALGERR
+            WRITE(errmsg, *) "Pres=0.  Unable to take log of 0."
+            RETURN
+        END IF
+
+        IF (interp_errcnt > 0) THEN
+            errstat = ALGERR
+            WRITE(errmsg, *) "bad difference in vcp's"
+            RETURN
+        END IF
 
-        ! print *,"----done----",interp_levels(nreqlvs)
+!$OMP PARALLEL DO
         DO j = 1,ns
             DO i = 1,ew
                 dataout(i,j,nreqlvs) = tempout(i,j)
             END DO
         END DO
+!$OMP END PARALLEL DO
 
     END DO !end for the nreqlvs