forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
440 lines
16 KiB
440 lines
16 KiB
! The subroutines in this file were taken directly from RIP code written |
|
! by Dr. Mark Stoelinga. They were modified by Sherrie |
|
! Fredrick(NCAR/MMM) to work with NCL February 2015. |
|
|
|
!NCLFORTSTART |
|
SUBROUTINE wrf_monotonic(out, in, lvprs, cor, idir, delta, ew, ns, nz, icorsw) |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: out |
|
|
|
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(OUT) :: out |
|
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: lvprs |
|
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: cor |
|
|
|
|
|
!NCLEND |
|
|
|
INTEGER :: i, j, k, 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 = nz,1,-1 |
|
IF (lvprs(i,j,k) .GE. 300.D0) THEN |
|
k300 = k |
|
EXIT |
|
END IF |
|
END DO |
|
|
|
DO k = k300,1,-1 |
|
IF (idir .EQ. 1) THEN |
|
out(i,j,k) = MIN(in(i,j,k), in(i,j,k+1) + delta) |
|
ELSE IF (idir .EQ. -1) THEN |
|
out(i,j,k) = MAX(in(i,j,k), in(i,j,k+1) - delta) |
|
END IF |
|
END DO |
|
|
|
DO k = k300+1, nz |
|
IF (idir .EQ. 1) THEN |
|
out(i,j,k) = MAX(in(i,j,k), in(i,j,k-1) - delta) |
|
ELSE IF (idir .EQ. -1) THEN |
|
out(i,j,k) = MIN(in(i,j,k), in(i,j,k-1) + delta) |
|
END IF |
|
END DO |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE wrf_monotonic |
|
|
|
|
|
!NCLFORTSTART |
|
FUNCTION wrf_intrp_value(wvalp0, wvalp1, vlev, vcp0, vcp1, icase, errstat) |
|
USE wrf_constants, ONLY : ALGERR, SCLHT |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
|
|
INTEGER, INTENT(IN) :: icase |
|
REAL(KIND=8), INTENT(IN) :: wvalp0, wvalp1, vlev, vcp0, vcp1 |
|
INTEGER, INTENT(INOUT) :: errstat |
|
REAL(KIND=8) :: wrf_intrp_value |
|
|
|
!NCLEND |
|
|
|
REAL(KIND=8) :: valp0, valp1, rvalue |
|
REAL(KIND=8) :: chkdiff |
|
|
|
errstat = 0 |
|
|
|
valp0 = wvalp0 |
|
valp1 = wvalp1 |
|
IF ( icase .EQ. 2) THEN !GHT |
|
valp0=EXP(-wvalp0/SCLHT) |
|
valp1=EXP(-wvalp1/SCLHT) |
|
END IF |
|
|
|
chkdiff = vcp1 - vcp0 |
|
IF(chkdiff .EQ. 0) THEN |
|
errstat = ALGERR |
|
!errmsg = "bad difference in vcp's" |
|
wrf_intrp_value = 0 |
|
RETURN |
|
!PRINT *,"bad difference in vcp's" |
|
!STOP |
|
END IF |
|
|
|
rvalue = (vlev - vcp0)*(valp1 - valp0)/(vcp1 - vcp0) + valp0 |
|
IF (icase .EQ. 2) THEN !GHT |
|
wrf_intrp_value = -SCLHT*LOG(rvalue) |
|
ELSE |
|
wrf_intrp_value = rvalue |
|
END IF |
|
|
|
RETURN |
|
|
|
END FUNCTION wrf_intrp_value |
|
|
|
! NOTES: |
|
! vcarray is the array holding the values for the vertical coordinate. |
|
! It will always come in with the dimensions of the staggered U and V grid. |
|
|
|
!NCLFORTSTART |
|
SUBROUTINE wrf_vintrp(datain, dataout, pres, tk, qvp, ght, terrain,& |
|
sfp, smsfp, vcarray, interp_levels, numlevels,& |
|
icase, ew, ns, nz, extrap, vcor, logp, rmsg,& |
|
errstat, errmsg) |
|
USE wrf_constants, ONLY : ALGERR, SCLHT, EXPON, EXPONI, GAMMA, GAMMAMD, TLCLC1, & |
|
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3, & |
|
CELKEL, EPS, USSALR |
|
|
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: dataout |
|
|
|
INTEGER, INTENT(IN) :: ew, ns, nz, icase, extrap |
|
INTEGER, INTENT(IN) :: vcor, numlevels, logp |
|
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: datain, pres, tk, qvp |
|
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: ght |
|
REAL(KIND=8), DIMENSION(ew,ns), INTENT(IN) :: terrain, sfp, smsfp |
|
REAL(KIND=8), DIMENSION(ew,ns,numlevels), INTENT(OUT) :: dataout |
|
REAL(KIND=8), DIMENSION(ew,ns,nz), INTENT(IN) :: vcarray |
|
REAL(KIND=8), DIMENSION(numlevels), INTENT(IN) :: interp_levels |
|
REAL(KIND=8), INTENT(IN) :: rmsg |
|
INTEGER, INTENT(INOUT) :: errstat |
|
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
|
|
|
!NCLEND |
|
|
|
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 |
|
REAL(KIND=8) :: vcp1, vcp0, valp0, valp1 |
|
! REAL(KIND=8) :: cvc |
|
REAL(KIND=8) :: vclhsl, vctophsl !qvlhsl,ttlhsl |
|
REAL(KIND=8) :: wrf_intrp_value |
|
REAL(KIND=8) :: plhsl, zlhsl, ezlhsl, tlhsl, psurf, pratio, tlev |
|
REAL(KIND=8) :: ezsurf, psurfsm, zsurf, qvapor, vt |
|
REAL(KIND=8) :: ezlev, plev, zlev, ptarget, dpmin, dp |
|
REAL(KIND=8) :: pbot, zbot, tbotextrap, e |
|
REAL(KIND=8) :: tlcl, gammam |
|
CHARACTER(LEN=1) :: cvcord |
|
INTEGER :: thd |
|
|
|
! Removes the warnings for uninitialized variables |
|
cvcord = '' |
|
plev = 0 |
|
zlev = 0 |
|
vlev = 0 |
|
errstat = 0 |
|
interp_errcnt = 0 |
|
interp_errstat = 0 |
|
log_errcnt = 0 |
|
|
|
IF (vcor .EQ. 1) THEN |
|
cvcord = 'p' |
|
ELSE IF ((vcor .EQ. 2) .OR. (vcor .EQ. 3)) THEN |
|
cvcord = 'z' |
|
ELSE IF ((vcor .EQ. 4) .OR. (vcor .EQ. 5)) THEN |
|
cvcord = 't' |
|
END IF |
|
|
|
!$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 |
|
|
|
DO nreqlvs = 1,numlevels |
|
IF (cvcord .EQ. 'z') THEN |
|
! Convert rlevel to meters from km |
|
rlevel = interp_levels(nreqlvs) * 1000.D0 |
|
vlev = EXP(-rlevel/SCLHT) |
|
ELSE IF (cvcord .EQ. 'p') THEN |
|
vlev = interp_levels(nreqlvs) |
|
ELSE IF (cvcord .EQ. 't') THEN |
|
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 |
|
ifound = 0 |
|
DO k = 1,nz-1 |
|
ripk = nz-k+1 |
|
vcp1 = vcarray(i,j,ripk-1) |
|
vcp0 = vcarray(i,j,ripk) |
|
valp0 = datain(i,j,ripk) |
|
valp1 = datain(i,j,ripk-1) |
|
|
|
IF ((vlev .GE. vcp0 .AND. vlev .LE. vcp1) .OR. & |
|
(vlev .LE. vcp0 .AND. vlev .GE. vcp1)) THEN |
|
! print *,i,j,valp0,valp1 |
|
IF ((valp0 .EQ. rmsg) .OR. (valp1 .EQ. rmsg)) THEN |
|
tempout(i,j) = rmsg |
|
ifound = 1 |
|
ELSE |
|
IF (logp .EQ. 1) THEN |
|
vcp1 = LOG(vcp1) |
|
vcp0 = LOG(vcp0) |
|
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, 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 |
|
|
|
EXIT |
|
END IF |
|
END DO !end for the k loop |
|
|
|
IF (ifound .EQ. 1) THEN !Grid point is in the model domain |
|
CYCLE |
|
END IF |
|
|
|
!If the user has requested no extrapolatin then just assign |
|
!all values above or below the model level to rmsg. |
|
IF (extrap .EQ. 0) THEN |
|
tempout(i,j) = rmsg |
|
CYCLE |
|
END IF |
|
|
|
! The grid point is either above or below the model domain |
|
! First we will check to see if the grid point is above the |
|
! model domain. |
|
vclhsl = vcarray(i,j,1) !lowest model level |
|
vctophsl = vcarray(i,j,nz) !highest model level |
|
diff = vctophsl - vclhsl |
|
isign = NINT(diff/ABS(diff)) |
|
|
|
IF (isign*vlev .GE. isign*vctophsl) THEN |
|
! Assign the highest model level to the out array |
|
tempout(i,j) = datain(i,j,nz) |
|
CYCLE |
|
END IF |
|
|
|
! Only remaining possibility is that the specified level is below |
|
! lowest model level. If lowest model level value is missing, |
|
! set interpolated value to missing. |
|
|
|
IF (datain(i,j,1) .EQ. rmsg) THEN |
|
tempout(i,j) = rmsg |
|
CYCLE |
|
END IF |
|
|
|
! If the field comming in is not a pressure,temperature or height |
|
! field we can set the output array to the value at the lowest |
|
! model level. |
|
|
|
tempout(i,j) = datain(i,j,1) |
|
|
|
! For the special cases of pressure on height levels or height on |
|
! pressure levels, or temperature-related variables on pressure or |
|
! height levels, perform a special extrapolation based on |
|
! US Standard Atmosphere. Here we calcualate the surface pressure |
|
! with the altimeter equation. This is how RIP calculates the |
|
! surface pressure. |
|
IF (icase .GT. 0) THEN |
|
plhsl = pres(i,j,1) * 0.01D0 !pressure at lowest model level |
|
zlhsl = ght(i,j,1) !grid point height a lowest model level |
|
ezlhsl = EXP(-zlhsl/SCLHT) |
|
tlhsl = tk(i,j,1) !temperature in K at lowest model level |
|
zsurf = terrain(i,j) |
|
qvapor = MAX((qvp(i,j,1)*.001D0),1.e-15) |
|
! virtual temperature |
|
! vt = tlhsl * (eps + qvapor)/(eps*(1.0 + qvapor)) |
|
! psurf = plhsl * (vt/(vt+USSALR * (zlhsl-zsurf)))**rconst |
|
psurf = sfp(i,j) |
|
psurfsm = smsfp(i,j) |
|
ezsurf = EXP(-zsurf/SCLHT) |
|
|
|
! The if for checking above ground |
|
IF ((cvcord .EQ. 'z' .AND. vlev .LT. ezsurf) .OR. & |
|
(cvcord .EQ. 'p' .AND. vlev .LT. psurf)) THEN |
|
|
|
! We are below the lowest data level but above the ground. |
|
! Use linear interpolation (linear in prs and exp-height). |
|
|
|
IF (cvcord .EQ. 'p') THEN |
|
plev = vlev |
|
ezlev = ((plev - plhsl)*& |
|
ezsurf + (psurf - plev)*ezlhsl)/(psurf - plhsl) |
|
zlev = -SCLHT*LOG(ezlev) |
|
IF (icase .EQ. 2) THEN |
|
tempout(i,j) = zlev |
|
CYCLE |
|
END IF |
|
|
|
ELSE IF (cvcord .EQ. 'z') THEN |
|
ezlev = vlev |
|
zlev = -SCLHT*LOG(ezlev) |
|
plev = ((ezlev - ezlhsl)*& |
|
psurf + (ezsurf - ezlev)*plhsl)/(ezsurf - ezlhsl) |
|
IF (icase .EQ. 1) THEN |
|
tempout(i,j) = plev |
|
CYCLE |
|
END IF |
|
END IF |
|
|
|
ELSE !else for checking above ground |
|
ptarget = psurfsm - 150.D0 |
|
dpmin = 1.E4 |
|
DO k=1,nz |
|
ripk = nz-k+1 |
|
dp = ABS((pres(i,j,ripk) * 0.01D0) - ptarget) |
|
IF (dp .GT. dpmin) THEN |
|
EXIT |
|
END IF |
|
dpmin = MIN(dpmin, dp) |
|
END DO |
|
|
|
kupper = k-1 |
|
|
|
ripk = nz - kupper + 1 |
|
pbot = MAX(plhsl,psurf) |
|
zbot = MIN(zlhsl,zsurf) |
|
pratio = pbot/(pres(i,j,ripk) * 0.01D0) |
|
tbotextrap = tk(i,j,ripk)*(pratio)**EXPON |
|
! virtual temperature |
|
vt = tbotextrap * (EPS + qvapor)/(EPS*(1.0D0 + qvapor)) |
|
IF (cvcord .EQ. 'p') THEN |
|
plev = vlev |
|
zlev = zbot + vt/USSALR*(1. - (vlev/pbot)**EXPON) |
|
IF (icase .EQ. 2) THEN |
|
tempout(i,j) = zlev |
|
CYCLE |
|
END IF |
|
ELSE IF (cvcord .EQ. 'z') THEN |
|
zlev = -sclht*LOG(vlev) |
|
plev = pbot*(1. + USSALR/vt*(zbot - zlev))**EXPONI |
|
IF (icase .EQ. 1) THEN |
|
tempout(i,j) = plev |
|
CYCLE |
|
END IF |
|
END IF |
|
END IF !end if for checking above ground |
|
END IF !for icase gt 0 |
|
|
|
IF (icase .GT. 2) THEN !extrapolation for temperature |
|
tlev = tlhsl + (zlhsl - zlev)*USSALR |
|
qvapor = MAX(qvp(i,j,1), 1.e-15) |
|
gammam = GAMMA*(1. + GAMMAMD*qvapor) |
|
IF (icase .EQ. 3) THEN |
|
tempout(i,j) = tlev - CELKEL |
|
ELSE IF (icase .EQ. 4) THEN |
|
tempout(i,j) = tlev |
|
! Potential temperature - theta |
|
ELSE IF (icase .EQ. 5) THEN |
|
tempout(i,j) = tlev*(1000.D0/plev)**gammam |
|
! extraolation for equivalent potential temperature |
|
ELSE IF (icase .EQ. 6) THEN |
|
e = qvapor*plev/(EPS + qvapor) |
|
tlcl = TLCLC1/(LOG(tlev**TLCLC2/e) - TLCLC3) + TLCLC4 |
|
tempout(i,j)=tlev*(1000.D0/plev)**(gammam)*& |
|
EXP((THTECON1/tlcl - THTECON2)*& |
|
qvapor*(1. + THTECON3*qvapor)) |
|
END IF |
|
END IF |
|
|
|
!333 CONTINUE |
|
|
|
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 |
|
|
|
!$OMP PARALLEL DO COLLAPSE(2) |
|
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 |
|
|
|
RETURN |
|
|
|
END SUBROUTINE wrf_vintrp
|
|
|