diff --git a/fortran/wrf_vinterp.f90 b/fortran/wrf_vinterp.f90 index c9bebb0..b7a5f64 100644 --- a/fortran/wrf_vinterp.f90 +++ b/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 + IF (vlev .NE. 0.0D0) THEN + tmpvlev = LOG(vlev) + ELSE + log_errcnt = log_errcnt + 1 + tmpvlev = rmsg END IF - tmpvlev = LOG(vlev) ELSE tmpvlev = vlev END IF - tempout(i,j) = wrf_intrp_value(valp0, valp1, tmpvlev, vcp0, & - vcp1, icase, errstat, errmsg) - IF (errstat .NE. 0) THEN - RETURN - END IF - ! print *,"one ",i,j,tempout(i,j) - ifound = 1 + IF (tmpvlev .NE. rmsg) THEN + tempout(i,j) = wrf_intrp_value(valp0, valp1, tmpvlev, vcp0, & + vcp1, icase, interp_errstat) + + IF (interp_errstat .NE. 0) THEN + tempout(i,j) = rmsg + interp_errcnt = interp_errcnt + 1 + END IF + + ifound = 1 + END IF END IF - !GOTO 115 ! EXIT 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