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.
983 lines
37 KiB
983 lines
37 KiB
!====================================================================== |
|
! |
|
! !IROUTINE: VIRTUAL -- Calculate virtual temperature (K) |
|
! |
|
! !DESCRIPTION: |
|
! |
|
! This function returns a single value of virtual temperature in |
|
! K, given temperature in K and mixing ratio in kg/kg. For an |
|
! array of virtual temperatures, use subroutine VIRTUAL_TEMP. |
|
! |
|
! !INPUT: |
|
! RATMIX - water vapor mixing ratio (kg/kg) |
|
! TEMP - temperature (K) |
|
! |
|
! !OUTPUT: |
|
! TV - Virtual temperature (K) |
|
! |
|
|
|
! NCLFORTSTART |
|
REAL(KIND=8) FUNCTION TVIRTUAL(temp, ratmix) |
|
USE wrf_constants, ONLY : EPS |
|
|
|
!f2py threadsafe |
|
|
|
IMPLICIT NONE |
|
REAL(KIND=8), INTENT(IN) :: temp, ratmix |
|
|
|
! NCLEND |
|
|
|
TVIRTUAL = temp*(EPS + ratmix)/(EPS*(1.D0 + ratmix)) |
|
|
|
RETURN |
|
|
|
END FUNCTION TVIRTUAL |
|
|
|
! NCLFORTSTART |
|
REAL(KIND=8) FUNCTION TONPSADIABAT(thte, prs, psadithte, psadiprs, psaditmk, gamma,& |
|
errstat, errmsg) |
|
USE wrf_constants, ONLY : ALGERR |
|
!!$OMP DECLARE SIMD (TONPSADIABAT) |
|
!!uniform(thte,prs,psadithte,psadiprs,psaditmk) |
|
!f2py threadsafe |
|
!f2py intent(in,out) :: cape, cin |
|
|
|
IMPLICIT NONE |
|
REAL(KIND=8), INTENT(IN) :: thte |
|
REAL(KIND=8), INTENT(IN) :: prs |
|
REAL(KIND=8), DIMENSION(150), INTENT(IN) :: psadithte |
|
REAL(KIND=8), DIMENSION(150), INTENT(IN) :: psadiprs |
|
REAL(KIND=8), DIMENSION(150,150), INTENT(IN) :: psaditmk |
|
REAL(KIND=8), INTENT(IN) :: gamma |
|
INTEGER, INTENT(INOUT) :: errstat |
|
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
|
|
|
! NCLEND |
|
|
|
REAL(KIND=8) :: fracjt |
|
REAL(KIND=8) :: fracjt2 |
|
REAL(KIND=8) :: fracip |
|
REAL(KIND=8) :: fracip2 |
|
|
|
INTEGER :: l1, h1, mid1, rang1, l2, h2, mid2, rang2 |
|
INTEGER :: ip, jt |
|
|
|
! This function gives the temperature (in K) on a moist adiabat |
|
! (specified by thte in K) given pressure in hPa. It uses a |
|
! lookup table, with data that was generated by the Bolton (1980) |
|
! formula for theta_e. |
|
|
|
! First check if pressure is less than min pressure in lookup table. |
|
! If it is, assume parcel is so dry that the given theta-e value can |
|
! be interpretted as theta, and get temperature from the simple dry |
|
! theta formula. |
|
|
|
IF (prs .LE. psadiprs(150)) THEN |
|
TONPSADIABAT = thte * (prs/1000.D0)**gamma |
|
RETURN |
|
END IF |
|
|
|
! Otherwise, look for the given thte/prs point in the lookup table. |
|
|
|
jt = -1 |
|
l1 = 1 |
|
h1 = 149 |
|
rang1 = h1 - l1 |
|
mid1 = (h1 + l1) / 2 |
|
DO WHILE(rang1 .GT. 1) |
|
if(thte .GE. psadithte(mid1)) then |
|
l1 = mid1 |
|
else |
|
h1 = mid1 |
|
end if |
|
rang1 = h1 - l1 |
|
mid1 = (h1 + l1) / 2 |
|
END DO |
|
jt = l1 |
|
|
|
! DO jtch = 1, 150-1 |
|
! IF (thte .GE. psadithte(jtch) .AND. thte .LT. psadithte(jtch+1)) THEN |
|
! jt = jtch |
|
! EXIT |
|
! !GO TO 213 |
|
! END IF |
|
! END DO |
|
|
|
ip = -1 |
|
l2 = 1 |
|
h2 = 149 |
|
rang2 = h2 - l2 |
|
mid2 = (h2 + l2) / 2 |
|
DO WHILE(rang2 .GT. 1) |
|
if(prs .LE. psadiprs(mid2)) then |
|
l2 = mid2 |
|
else |
|
h2 = mid2 |
|
end if |
|
rang2 = h2 - l2 |
|
mid2 = (h2 + l2) / 2 |
|
END DO |
|
ip = l2 |
|
|
|
! ip = -1 |
|
! DO ipch = 1, 150-1 |
|
! IF (prs .LE. psadiprs(ipch) .AND. prs .GT. psadiprs(ipch+1)) THEN |
|
! ip = ipch |
|
! EXIT |
|
! !GO TO 215 |
|
! END IF |
|
! END DO |
|
|
|
IF (jt .EQ. -1 .OR. ip .EQ. -1) THEN |
|
! Set the error and return |
|
TONPSADIABAT = -1 |
|
errstat = ALGERR |
|
WRITE(errmsg, *) "capecalc3d: Outside of lookup table bounds. prs,thte=", prs, thte |
|
RETURN |
|
END IF |
|
|
|
fracjt = (thte-psadithte(jt)) / (psadithte(jt+1)-psadithte(jt)) |
|
fracjt2 = 1.D0 - fracjt |
|
fracip = (psadiprs(ip)-prs) / (psadiprs(ip)-psadiprs(ip+1)) |
|
fracip2 = 1.D0 - fracip |
|
|
|
IF (psaditmk(ip,jt) .GT. 1D9 .OR. psaditmk(ip+1,jt) .GT. 1D9 .OR. & |
|
psaditmk(ip,jt+1) .GT. 1D9 .OR. psaditmk(ip+1,jt+1) .GT. 1D9) THEN |
|
! Set the error and return |
|
TONPSADIABAT = -1 |
|
errstat = ALGERR |
|
WRITE(errmsg, *) "capecalc3d: Tried to access missing temperature in lookup table. ", & |
|
"Prs and Thte probably unreasonable. prs,thte=", prs, thte |
|
RETURN |
|
END IF |
|
|
|
TONPSADIABAT = fracip2*fracjt2*psaditmk(ip,jt) + fracip*fracjt2*psaditmk(ip+1,jt) + & |
|
fracip2*fracjt*psaditmk(ip,jt+1) + fracip*fracjt*psaditmk(ip+1,jt+1) |
|
|
|
RETURN |
|
|
|
END FUNCTION TONPSADIABAT |
|
|
|
!NCLFORTSTART |
|
SUBROUTINE DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, fname, errstat, errmsg) |
|
USE wrf_constants, ONLY : ALGERR |
|
|
|
!f2py threadsafe |
|
|
|
REAL(KIND=8), DIMENSION(150), INTENT(INOUT) :: psadithte, psadiprs |
|
REAL(KIND=8), DIMENSION(150,150), INTENT(INOUT) :: psaditmk |
|
CHARACTER(LEN=*), INTENT(IN) :: fname |
|
INTEGER, INTENT(INOUT) :: errstat |
|
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
|
|
|
!NCLEND |
|
|
|
! Locals |
|
INTEGER :: iustnlist, i, nthte, nprs, ip, jt |
|
|
|
! FNAME = 'psadilookup.dat' |
|
iustnlist = 33 |
|
OPEN (UNIT=iustnlist, FILE=fname, FORM='formatted', STATUS='old') |
|
|
|
DO i = 1,14 |
|
READ (iustnlist, FMT=*) |
|
END DO |
|
|
|
READ (iustnlist, FMT=*) nthte, nprs |
|
|
|
IF (nthte .NE. 150 .OR. nprs .NE. 150) THEN |
|
errstat = ALGERR |
|
errmsg = "Number of pressure or theta_e levels in lookup table file not 150" |
|
RETURN |
|
END IF |
|
|
|
READ (iustnlist, FMT="(5D15.7)") (psadithte(jt),jt=1,nthte) |
|
READ (iustnlist, FMT="(5D15.7)") (psadiprs(ip),ip=1,nprs) |
|
READ (iustnlist, FMT="(5D15.7)") ((psaditmk(ip,jt),ip=1,nprs),jt=1,nthte) |
|
|
|
CLOSE (iustnlist) |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DLOOKUP_TABLE |
|
|
|
|
|
! Historically, this routine calculated the pressure at full sigma |
|
! levels when RIP was specifically designed for MM4/MM5 output. |
|
! With the new generalized RIP (Feb '02), this routine is still |
|
! intended to calculate a set of pressure levels that bound the |
|
! layers represented by the vertical grid points, although no such |
|
! layer boundaries are assumed to be defined. The routine simply |
|
! uses the midpoint between the pressures of the vertical grid |
|
! points as the bounding levels. The array only contains mkzh |
|
! levels, so the pressure of the top of the uppermost layer is |
|
! actually excluded. The kth value of pf is the lower bounding |
|
! pressure for the layer represented by kth data level. At the |
|
! lower bounding level of the lowest model layer, it uses the |
|
! surface pressure, unless the data set is pressure-level data, in |
|
! which case it assumes the lower bounding pressure level is as far |
|
! below the lowest vertical level as the upper bounding pressure |
|
! level is above. |
|
SUBROUTINE DPFCALC(prs, sfp, pf, mix, mjy, mkzh, ter_follow) |
|
|
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(IN) :: prs |
|
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) :: sfp |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy), INTENT(OUT) :: pf |
|
INTEGER, INTENT(IN) :: ter_follow,mix,mjy,mkzh |
|
|
|
INTEGER :: i,j,k |
|
|
|
DO j = 1,mjy |
|
DO i = 1,mix |
|
DO k = 1,mkzh |
|
IF (k .EQ. mkzh) THEN |
|
! terrain-following data |
|
IF (ter_follow .EQ. 1) THEN |
|
pf(k,i,j) = sfp(i,j) |
|
! pressure-level data |
|
ELSE |
|
pf(k,i,j) = .5D0 * (3.D0*prs(k,i,j) - prs(k-1,i,j)) |
|
END IF |
|
ELSE |
|
pf(k,i,j) = .5D0 * (prs(k+1,i,j) + prs(k,i,j)) |
|
END IF |
|
END DO |
|
END DO |
|
END DO |
|
|
|
RETURN |
|
|
|
END SUBROUTINE DPFCALC |
|
|
|
!====================================================================== |
|
! |
|
! !IROUTINE: capecalc3d -- Calculate CAPE and CIN |
|
! |
|
! !DESCRIPTION: |
|
! |
|
! This routine calculates CAPE and CIN (in m**2/s**2, |
|
! or J/kg) for every grid point in the entire 3D domain (treating |
|
! each grid point as a parcel). |
|
! |
|
|
|
|
|
! Important! The z-indexes must be arranged so that mkzh (max z-index) is the |
|
! surface pressure. So, pressure must be ordered in ascending order before |
|
! calling this routine. Other variables must be ordered the same (p,tk,q,z). |
|
|
|
! Also, be advised that missing data values are not checked during the computation. |
|
! Also also, Pressure must be hPa |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DCAPECALC3D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& |
|
cmsg,mix,mjy,mkzh,ter_follow,& |
|
psafile, errstat, errmsg) |
|
USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & |
|
EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & |
|
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 |
|
|
|
!USE omp_lib |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: cape, cin |
|
|
|
INTEGER, INTENT(IN) :: mix, mjy, mkzh, ter_follow |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: prs |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: tmk |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: qvp |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: ght |
|
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) :: ter |
|
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin |
|
REAL(KIND=8), INTENT(IN) :: cmsg |
|
CHARACTER(LEN=*), INTENT(IN) :: psafile |
|
INTEGER, INTENT(INOUT) :: errstat |
|
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
|
|
|
! NCLFORTEND |
|
|
|
! local variables |
|
INTEGER :: i, j, k, ilcl, kel, kk, klcl, klev, klfc, kmax, kpar |
|
REAL(KIND=8) :: tlcl, zlcl |
|
REAL(KIND=8) :: ethpari, qvppari, tmkpari |
|
REAL(KIND=8) :: facden, qvplift, tmklift, tvenv, tvlift, ghtlift |
|
REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat |
|
REAL(KIND=8) :: benamin, dz |
|
REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf |
|
REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs |
|
REAL(KIND=8), DIMENSION(150,150) :: psaditmk |
|
LOGICAL :: elfound |
|
REAL :: t1,t2 |
|
|
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new |
|
|
|
! To remove compiler warnings |
|
tmkpari = 0 |
|
qvppari = 0 |
|
klev = 0 |
|
klcl = 0 |
|
kel = 0 |
|
|
|
|
|
! the comments were taken from a mark stoelinga email, 23 apr 2007, |
|
! in response to a user getting the "outside of lookup table bounds" |
|
! error message. |
|
|
|
! tmkpari - initial temperature of parcel, k |
|
! values of 300 okay. (not sure how much from this you can stray.) |
|
|
|
! prspari - initial pressure of parcel, hpa |
|
! values of 980 okay. (not sure how much from this you can stray.) |
|
|
|
! thtecon1, thtecon2, thtecon3 |
|
! these are all constants, the first in k and the other two have |
|
! no units. values of 3376, 2.54, and 0.81 were stated as being |
|
! okay. |
|
|
|
! tlcl - the temperature at the parcel's lifted condensation level, k |
|
! should be a reasonable atmospheric temperature around 250-300 k |
|
! (398 is "way too high") |
|
|
|
! qvppari - the initial water vapor mixing ratio of the parcel, |
|
! kg/kg (should range from 0.000 to 0.025) |
|
! |
|
|
|
! calculated the pressure at full sigma levels (a set of pressure |
|
! levels that bound the layers represented by the vertical grid points) |
|
|
|
!CALL cpu_time(t1) |
|
!CALL OMP_SET_NUM_THREADS(16) |
|
!$OMP PARALLEL DO |
|
DO j = 1,mjy |
|
DO i = 1,mix |
|
DO k = 1,mkzh |
|
prs_new(k,i,j) = prs(i,j,k) |
|
tmk_new(k,i,j) = tmk(i,j,k) |
|
qvp_new(k,i,j) = qvp(i,j,k) |
|
ght_new(k,i,j) = ght(i,j,k) |
|
END DO |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
CALL DPFCALC(prs_new, sfp, prsf, mix, mjy, mkzh, ter_follow) |
|
|
|
! before looping, set lookup table for getting temperature on |
|
! a pseudoadiabat. |
|
|
|
CALL DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, psafile, errstat, errmsg) |
|
|
|
IF (errstat .NE. 0) THEN |
|
RETURN |
|
END IF |
|
|
|
!$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, & |
|
!$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 |
|
|
|
!!$OMP SIMD |
|
DO kpar = 2, mkzh |
|
|
|
! Calculate temperature and moisture properties of parcel |
|
! (note, qvppari and tmkpari already calculated above for 2d case.) |
|
|
|
tlcl = TLCLC1/(LOG(tmk_new(kpar,i,j)**TLCLC2/(MAX(1.D-20,qvp_new(kpar,i,j)*prs_new(kpar,i,j)/ & |
|
(EPS + qvp_new(kpar,i,j))))) - TLCLC3) + TLCLC4 |
|
|
|
ethpari = tmk_new(kpar,i,j)*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvp_new(kpar,i,j)))* & |
|
EXP((THTECON1/tlcl - THTECON2)*qvp_new(kpar,i,j)*(1.D0 + THTECON3*qvp_new(kpar,i,j))) |
|
|
|
zlcl = ght_new(kpar,i,j) + (tmk_new(kpar,i,j) - tlcl)/(G/CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) |
|
|
|
! DO k = kpar,1,-1 |
|
! tmklift_new(k) = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& |
|
! psaditmk, GAMMA, errstat, errmsg) |
|
! END DO |
|
! Calculate buoyancy and relative height of lifted parcel at |
|
! all levels, and store in bottom up arrays. add a level at the lcl, |
|
! and at all points where buoyancy is zero. |
|
! |
|
! For arrays that go bottom to top |
|
kk = 0 |
|
ilcl = 0 |
|
|
|
IF (ght_new(kpar,i,j) .GE. zlcl) THEN |
|
! Initial parcel already saturated or supersaturated. |
|
ilcl = 2 |
|
klcl = 1 |
|
END IF |
|
|
|
!!$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 |
|
|
|
! Model level is below lcl |
|
IF (ght_new(k,i,j) .LT. zlcl) THEN |
|
tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) * (ght_new(k,i,j) - ght_new(kpar,i,j)) |
|
tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j))) |
|
tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j))) |
|
ghtlift = ght_new(k,i,j) |
|
ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN |
|
! This model level and previous model level straddle the lcl, |
|
! so first create a new level in the bottom-up array, at the lcl. |
|
facden = 1.0/(ght_new(k,i,j) - ght_new(k+1,i,j)) |
|
tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) |
|
qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) |
|
tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv)) |
|
tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j))) |
|
ghtlift = zlcl |
|
ilcl = 1 |
|
ELSE |
|
tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& |
|
psaditmk, GAMMA, errstat, errmsg) |
|
eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) |
|
qvplift = EPS*eslift/(prs_new(k,i,j) - eslift) |
|
tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j))) |
|
tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift)) |
|
ghtlift = ght_new(k,i,j) |
|
END IF |
|
! Buoyancy |
|
buoy(kk) = G*(tvlift - tvenv)/tvenv |
|
zrel(kk) = ghtlift - ght_new(kpar,i,j) |
|
IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN |
|
! Parcel ascent curve crosses sounding curve, so create a new level |
|
! in the bottom-up array at the crossing. |
|
kk = kk + 1 |
|
buoy(kk) = buoy(kk-1) |
|
zrel(kk) = zrel(kk-1) |
|
buoy(kk-1) = 0.D0 |
|
zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& |
|
(buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) |
|
END IF |
|
IF (ilcl .EQ. 1) THEN |
|
klcl = kk |
|
ilcl = 2 |
|
CYCLE |
|
END IF |
|
|
|
END DO |
|
|
|
kmax = kk |
|
! IF (kmax .GT. 150) THEN |
|
! print *,'kmax got too big' |
|
! errstat = ALGERR |
|
! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax |
|
! RETURN |
|
! END IF |
|
|
|
! If no lcl was found, set klcl to kmax. it is probably not really |
|
! at kmax, but this will make the rest of the routine behave |
|
! properly. |
|
IF (ilcl .EQ. 0) klcl=kmax |
|
|
|
! Get the accumulated buoyant energy from the parcel's starting |
|
! point, at all levels up to the top level. |
|
benaccum(1) = 0.0D0 |
|
benamin = 9d9 |
|
DO k = 2,kmax |
|
dz = zrel(k) - zrel(k-1) |
|
benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k)) |
|
IF (benaccum(k) .LT. benamin) THEN |
|
benamin = benaccum(k) |
|
END IF |
|
END DO |
|
! Determine equilibrium level (el), which we define as the highest |
|
! level of non-negative buoyancy above the lcl. note, this may be |
|
! the top level if the parcel is still buoyant there. |
|
|
|
elfound = .FALSE. |
|
DO k = kmax,klcl,-1 |
|
IF (buoy(k) .GE. 0.D0) THEN |
|
! k of equilibrium level |
|
kel = k |
|
elfound = .TRUE. |
|
EXIT |
|
END IF |
|
END DO |
|
|
|
! If we got through that loop, then there is no non-negative |
|
! buoyancy above the lcl in the sounding. in these situations, |
|
! both cape and cin will be set to -0.1 j/kg. (see below about |
|
! missing values in v6.1.0). also, where cape is |
|
! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so |
|
! that the zero contour in either the cin or cape fields will |
|
! circumscribe regions of non-zero cape. |
|
|
|
! In v6.1.0 of ncl, we added a _fillvalue attribute to the return |
|
! value of this function. at that time we decided to change -0.1 |
|
! to a more appropriate missing value, which is passed into this |
|
! routine as cmsg. |
|
|
|
IF (.NOT. elfound) THEN |
|
!print *,'el not found' |
|
cape(i,j,kpar) = cmsg |
|
cin(i,j,kpar) = cmsg |
|
klfc = kmax |
|
CYCLE |
|
END IF |
|
|
|
! If there is an equilibrium level, then cape is positive. we'll |
|
! define the level of free convection (lfc) as the point below the |
|
! el, but at or above the lcl, where accumulated buoyant energy is a |
|
! minimum. the net positive area (accumulated buoyant energy) from |
|
! the lfc up to the el will be defined as the cape, and the net |
|
! negative area (negative of accumulated buoyant energy) from the |
|
! parcel starting point to the lfc will be defined as the convective |
|
! inhibition (cin). |
|
|
|
! First get the lfc according to the above definition. |
|
benamin = 9D9 |
|
klfc = kmax |
|
DO k = klcl,kel |
|
IF (benaccum(k) .LT. benamin) THEN |
|
benamin = benaccum(k) |
|
klfc = k |
|
END IF |
|
END DO |
|
|
|
! Now we can assign values to cape and cin |
|
|
|
cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0) |
|
cin(i,j,kpar) = MAX(-benamin, 0.1D0) |
|
|
|
! cin is uninteresting when cape is small (< 100 j/kg), so set |
|
! cin to -0.1 (see note about missing values in v6.1.0) in |
|
! that case. |
|
|
|
! In v6.1.0 of ncl, we added a _fillvalue attribute to the return |
|
! value of this function. at that time we decided to change -0.1 |
|
! to a more appropriate missing value, which is passed into this |
|
! routine as cmsg. |
|
|
|
IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg |
|
|
|
END DO |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
!CALL cpu_time(t2) |
|
!print *,'Time taken in seconds ',(t2-t1) |
|
RETURN |
|
END SUBROUTINE DCAPECALC3D |
|
|
|
!====================================================================== |
|
! |
|
! !IROUTINE: capecalc2d -- Calculate CAPE and CIN |
|
! |
|
! !DESCRIPTION: |
|
! |
|
! Calculates CAPE and CIN only for the parcel with max theta-e in |
|
! the column, (i.e. something akin to Colman's MCAPE). By "parcel", |
|
! we mean a 500-m deep parcel, with actual temperature and moisture |
|
! averaged over that depth. |
|
! |
|
! CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of |
|
! the cape and cin arrays. Also, LCL and LFC heights |
|
! are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array. |
|
! |
|
|
|
|
|
! Important! The z-indexes must be arranged so that mkzh (max z-index) is the |
|
! surface pressure. So, pressure must be ordered in ascending order before |
|
! calling this routine. Other variables must be ordered the same (p,tk,q,z). |
|
|
|
! Also, be advised that missing data values are not checked during the |
|
! computation. |
|
! Also also, Pressure must be hPa |
|
|
|
! NCLFORTSTART |
|
SUBROUTINE DCAPECALC2D(prs,tmk,qvp,ght,ter,sfp,cape,cin,& |
|
cmsg,mix,mjy,mkzh,ter_follow,& |
|
psafile, errstat, errmsg) |
|
USE wrf_constants, ONLY : ALGERR, CELKEL, G, EZERO, ESLCON1, ESLCON2, & |
|
EPS, RD, CP, GAMMA, CPMD, RGASMD, GAMMAMD, TLCLC1, & |
|
TLCLC2, TLCLC3, TLCLC4, THTECON1, THTECON2, THTECON3 |
|
|
|
!USE omp_lib |
|
IMPLICIT NONE |
|
|
|
!f2py threadsafe |
|
!f2py intent(in,out) :: cape, cin |
|
|
|
INTEGER, INTENT(IN) :: mix, mjy, mkzh, ter_follow |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: prs |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: tmk |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: qvp |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(IN) :: ght |
|
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) :: ter |
|
REAL(KIND=8), DIMENSION(mix,mjy), INTENT(IN) ::sfp |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cape |
|
REAL(KIND=8), DIMENSION(mix,mjy,mkzh), INTENT(OUT) :: cin |
|
REAL(KIND=8), INTENT(IN) :: cmsg |
|
CHARACTER(LEN=*), INTENT(IN) :: psafile |
|
INTEGER, INTENT(INOUT) :: errstat |
|
CHARACTER(LEN=*), INTENT(INOUT) :: errmsg |
|
|
|
! NCLFORTEND |
|
|
|
|
|
! local variables |
|
INTEGER :: i, j, k, ilcl, kel, kk, klcl, klev, klfc, kmax, kpar, kpar1, kpar2 |
|
REAL(KIND=8) :: ethmax, q, p, e, tlcl, zlcl |
|
REAL(KIND=8) :: pavg, tvirtual, p1, p2, pp1, pp2, th, totthe, totqvp, totprs |
|
REAL(KIND=8) :: cpm, deltap, ethpari, gammam, qvppari, tmkpari |
|
REAL(KIND=8) :: facden, qvplift, tmklift, tvenv, tvlift, ghtlift, fac1, fac2 |
|
REAL(KIND=8) :: eslift, tmkenv, qvpenv, tonpsadiabat |
|
REAL(KIND=8) :: benamin, dz, pup, pdn |
|
REAL(KIND=8), DIMENSION(150) :: buoy, zrel, benaccum |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prsf |
|
REAL(KIND=8), DIMENSION(150) :: psadithte, psadiprs |
|
REAL(KIND=8), DIMENSION(150,150) :: psaditmk |
|
LOGICAL :: elfound |
|
INTEGER :: nthreads |
|
REAL(KIND=8), DIMENSION(mkzh) :: eth_temp |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: prs_new |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: tmk_new |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: qvp_new |
|
REAL(KIND=8), DIMENSION(mkzh,mix,mjy) :: ght_new |
|
|
|
! To remove compiler warnings |
|
errstat = 0 |
|
tmkpari = 0 |
|
qvppari = 0 |
|
klev = 0 |
|
klcl = 0 |
|
kel = 0 |
|
deltap = 0 |
|
|
|
|
|
! the comments were taken from a mark stoelinga email, 23 apr 2007, |
|
! in response to a user getting the "outside of lookup table bounds" |
|
! error message. |
|
|
|
! tmkpari - initial temperature of parcel, k |
|
! values of 300 okay. (not sure how much from this you can stray.) |
|
|
|
! prspari - initial pressure of parcel, hpa |
|
! values of 980 okay. (not sure how much from this you can stray.) |
|
|
|
! thtecon1, thtecon2, thtecon3 |
|
! these are all constants, the first in k and the other two have |
|
! no units. values of 3376, 2.54, and 0.81 were stated as being |
|
! okay. |
|
|
|
! tlcl - the temperature at the parcel's lifted condensation level, k |
|
! should be a reasonable atmospheric temperature around 250-300 k |
|
! (398 is "way too high") |
|
|
|
! qvppari - the initial water vapor mixing ratio of the parcel, |
|
! kg/kg (should range from 0.000 to 0.025) |
|
! |
|
|
|
!$OMP PARALLEL DO COLLAPSE(3) |
|
DO j = 1,mjy |
|
DO i = 1,mix |
|
DO k = 1,mkzh |
|
prs_new(k,i,j) = prs(i,j,k) |
|
tmk_new(k,i,j) = tmk(i,j,k) |
|
qvp_new(k,i,j) = qvp(i,j,k) |
|
ght_new(k,i,j) = ght(i,j,k) |
|
END DO |
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
|
|
|
|
! calculated the pressure at full sigma levels (a set of pressure |
|
! levels that bound the layers represented by the vertical grid points) |
|
CALL DPFCALC(prs_new, sfp, prsf, mix, mjy, mkzh, ter_follow) |
|
|
|
! before looping, set lookup table for getting temperature on |
|
! a pseudoadiabat. |
|
|
|
CALL DLOOKUP_TABLE(psadithte, psadiprs, psaditmk, psafile, errstat, errmsg) |
|
|
|
IF (errstat .NE. 0) THEN |
|
RETURN |
|
END IF |
|
|
|
!CALL OMP_SET_NUM_THREADS(16) |
|
!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) |
|
DO j = 1,mjy |
|
DO i = 1,mix |
|
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 |
|
DO k = 1, mkzh |
|
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))))*& |
|
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)))) |
|
END IF |
|
END DO |
|
klev = mkzh |
|
DO k = 1,mkzh |
|
IF (eth_temp(k) .GT. ethmax) THEN |
|
klev = k |
|
ethmax = eth_temp(k) |
|
END IF |
|
END DO |
|
|
|
kpar1 = klev |
|
kpar2 = klev |
|
|
|
|
|
! Establish average properties of that parcel |
|
! (over depth of approximately davg meters) |
|
|
|
!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 |
|
totthe = 0.D0 |
|
totqvp = 0.D0 |
|
totprs = 0.D0 |
|
DO k = mkzh,2,-1 |
|
IF (prsf(k,i,j) .LE. p1) EXIT !GOTO 35 |
|
IF (prsf(k-1,i,j) .GE. p2) CYCLE !GOTO 34 |
|
p = prs_new(k,i,j) |
|
pup = prsf(k,i,j) |
|
pdn = prsf(k-1,i,j) |
|
!q = MAX(qvp_new(k,i,j),1.D-15) |
|
th = 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))) |
|
pp1 = MAX(p1,pdn) |
|
pp2 = MIN(p2,pup) |
|
IF (pp2 .GT. pp1) THEN |
|
! deltap = pp2 - pp1 |
|
totqvp = totqvp + MAX(qvp_new(k,i,j),1.D-15)*(pp2 - pp1) |
|
totthe = totthe + th*(pp2 - pp1) |
|
totprs = totprs + (pp2 - pp1) |
|
END IF |
|
END DO |
|
qvppari = totqvp/totprs |
|
tmkpari = (totthe/totprs)*& |
|
(prs_new(kpar1,i,j)/1000.D0)**(GAMMA*(1.D0+GAMMAMD*qvp_new(kpar1,i,j))) |
|
|
|
!CALL CPU_TIME(t3) |
|
DO kpar = kpar1, kpar2 |
|
|
|
! Calculate temperature and moisture properties of parcel |
|
! (note, qvppari and tmkpari already calculated above for 2d |
|
! case.) |
|
|
|
!prspari = prs_new(kpar,i,j) |
|
!ghtpari = ght_new(kpar,i,j) |
|
gammam = GAMMA * (1.D0 + GAMMAMD*qvppari) |
|
cpm = CP * (1.D0 + CPMD*qvppari) |
|
|
|
e = MAX(1.D-20,qvppari*prs_new(kpar,i,j)/(EPS + qvppari)) |
|
tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4 |
|
ethpari = tmkpari*(1000.D0/prs_new(kpar,i,j))**(GAMMA*(1.D0 + GAMMAMD*qvppari))*& |
|
EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.D0 + THTECON3*qvppari)) |
|
zlcl = ght_new(kpar,i,j) + (tmkpari - tlcl)/(G/cpm) |
|
|
|
! Calculate buoyancy and relative height of lifted parcel at |
|
! all levels, and store in bottom up arrays. add a level at the |
|
! lcl, |
|
! and at all points where buoyancy is zero. |
|
! |
|
! |
|
! For arrays that go bottom to top |
|
kk = 0 |
|
ilcl = 0 |
|
|
|
IF (ght_new(kpar,i,j) .GE. zlcl) THEN |
|
! Initial parcel already saturated or supersaturated. |
|
ilcl = 2 |
|
klcl = 1 |
|
END IF |
|
|
|
k = kpar |
|
DO k = kpar,1,-1 |
|
! For arrays that go bottom to top |
|
kk = kk + 1 |
|
|
|
! Model level is below lcl |
|
IF (ght_new(k,i,j) .LT. zlcl) THEN |
|
tmklift = tmk_new(kpar,i,j) - G/(CP * (1.D0 + CPMD*qvp_new(kpar,i,j))) * (ght_new(k,i,j) - ght_new(kpar,i,j)) |
|
tvenv = tmk_new(k,i,j)*(EPS + qvp_new(k,i,j))/(EPS*(1.D0 + qvp_new(k,i,j))) |
|
tvlift = tmklift*(EPS + qvp_new(kpar,i,j))/(EPS*(1.D0 + qvp_new(kpar,i,j))) |
|
ghtlift = ght_new(k,i,j) |
|
ELSE IF (ght(i,j,k) .GE. zlcl .AND. ilcl .EQ. 0) THEN |
|
! This model level and previous model level straddle the lcl, |
|
! so first create a new level in the bottom-up array, at the lcl. |
|
facden = 1/(ght_new(k,i,j) - ght_new(k+1,i,j)) |
|
tmkenv = tmk_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + tmk_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) |
|
qvpenv = qvp_new(k+1,i,j)*((ght_new(k,i,j)-zlcl)*facden) + qvp_new(k,i,j)*((zlcl-ght_new(k+1,i,j))*facden) |
|
tvenv = tmkenv* (EPS + qvpenv) / (EPS * (1.D0 + qvpenv)) |
|
tvlift = tlcl* (EPS + qvp_new(kpar,i,j)) / (EPS *(1.D0 + qvp_new(kpar,i,j))) |
|
ghtlift = zlcl |
|
ilcl = 1 |
|
ELSE |
|
tmklift = TONPSADIABAT(ethpari, prs_new(k,i,j), psadithte, psadiprs,& |
|
psaditmk, GAMMA, errstat, errmsg) |
|
eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) |
|
qvplift = EPS*eslift/(prs_new(k,i,j) - eslift) |
|
tvenv = tmk_new(k,i,j) * (EPS + qvp_new(k,i,j)) / (EPS * (1.D0 + qvp_new(k,i,j))) |
|
tvlift = tmklift*(EPS + qvplift) / (EPS * (1.D0 + qvplift)) |
|
ghtlift = ght_new(k,i,j) |
|
END IF |
|
! Buoyancy |
|
buoy(kk) = G*(tvlift - tvenv)/tvenv |
|
zrel(kk) = ghtlift - ght_new(kpar,i,j) |
|
IF ((kk .GT. 1) .AND. (buoy(kk)*buoy(kk-1) .LT. 0.0D0)) THEN |
|
! Parcel ascent curve crosses sounding curve, so create a new level |
|
! in the bottom-up array at the crossing. |
|
kk = kk + 1 |
|
buoy(kk) = buoy(kk-1) |
|
zrel(kk) = zrel(kk-1) |
|
buoy(kk-1) = 0.D0 |
|
zrel(kk-1) = zrel(kk-2) + buoy(kk-2)/& |
|
(buoy(kk-2) - buoy(kk))*(zrel(kk) - zrel(kk-2)) |
|
END IF |
|
IF (ilcl .EQ. 1) THEN |
|
klcl = kk |
|
ilcl = 2 |
|
CYCLE |
|
END IF |
|
|
|
END DO |
|
|
|
kmax = kk |
|
! IF (kmax .GT. 150) THEN |
|
! errstat = ALGERR |
|
! WRITE(errmsg, *) 'capecalc3d: kmax got too big. kmax=',kmax |
|
! RETURN |
|
! END IF |
|
|
|
! If no lcl was found, set klcl to kmax. it is probably not |
|
! really |
|
! at kmax, but this will make the rest of the routine behave |
|
! properly. |
|
IF (ilcl .EQ. 0) klcl=kmax |
|
|
|
! Get the accumulated buoyant energy from the parcel's starting |
|
! point, at all levels up to the top level. |
|
benaccum(1) = 0.0D0 |
|
benamin = 9d9 |
|
DO k = 2,kmax |
|
dz = zrel(k) - zrel(k-1) |
|
benaccum(k) = benaccum(k-1) + .5D0*dz*(buoy(k-1) + buoy(k)) |
|
IF (benaccum(k) .LT. benamin) THEN |
|
benamin = benaccum(k) |
|
END IF |
|
END DO |
|
! Determine equilibrium level (el), which we define as the highest |
|
! level of non-negative buoyancy above the lcl. note, this may be |
|
! the top level if the parcel is still buoyant there. |
|
|
|
elfound = .FALSE. |
|
DO k = kmax,klcl,-1 |
|
IF (buoy(k) .GE. 0.D0) THEN |
|
! k of equilibrium level |
|
kel = k |
|
elfound = .TRUE. |
|
EXIT |
|
END IF |
|
END DO |
|
|
|
! If we got through that loop, then there is no non-negative |
|
! buoyancy above the lcl in the sounding. in these situations, |
|
! both cape and cin will be set to -0.1 j/kg. (see below about |
|
! missing values in v6.1.0). also, where cape is |
|
! non-zero, cape and cin will be set to a minimum of +0.1 j/kg, so |
|
! that the zero contour in either the cin or cape fields will |
|
! circumscribe regions of non-zero cape. |
|
! In v6.1.0 of ncl, we added a _fillvalue attribute to the return |
|
! value of this function. at that time we decided to change -0.1 |
|
! to a more appropriate missing value, which is passed into this |
|
! routine as cmsg. |
|
|
|
IF (.NOT. elfound) THEN |
|
cape(i,j,kpar) = cmsg |
|
cin(i,j,kpar) = cmsg |
|
klfc = kmax |
|
CYCLE |
|
END IF |
|
|
|
|
|
! If there is an equilibrium level, then cape is positive. |
|
! we'll |
|
! define the level of free convection (lfc) as the point below |
|
! the |
|
! el, but at or above the lcl, where accumulated buoyant energy |
|
! is a |
|
! minimum. the net positive area (accumulated buoyant energy) |
|
! from |
|
! the lfc up to the el will be defined as the cape, and the net |
|
! negative area (negative of accumulated buoyant energy) from |
|
! the |
|
! parcel starting point to the lfc will be defined as the |
|
! convective |
|
! inhibition (cin). |
|
|
|
! First get the lfc according to the above definition. |
|
benamin = 9D9 |
|
klfc = kmax |
|
DO k = klcl,kel |
|
IF (benaccum(k) .LT. benamin) THEN |
|
benamin = benaccum(k) |
|
klfc = k |
|
END IF |
|
END DO |
|
|
|
! Now we can assign values to cape and cin |
|
|
|
cape(i,j,kpar) = MAX(benaccum(kel)-benamin, 0.1D0) |
|
cin(i,j,kpar) = MAX(-benamin, 0.1D0) |
|
|
|
! cin is uninteresting when cape is small (< 100 j/kg), so set |
|
! cin to -0.1 (see note about missing values in v6.1.0) in |
|
! that case. |
|
|
|
! In v6.1.0 of ncl, we added a _fillvalue attribute to the return |
|
! value of this function. at that time we decided to change -0.1 |
|
! to a more appropriate missing value, which is passed into this |
|
! routine as cmsg. |
|
|
|
IF (cape(i,j,kpar) .LT. 100.D0) cin(i,j,kpar) = cmsg |
|
|
|
END DO |
|
|
|
cape(i,j,mkzh) = cape(i,j,kpar1) |
|
cin(i,j,mkzh) = cin(i,j,kpar1) |
|
! meters agl |
|
cin(i,j,mkzh-1) = zrel(klcl) + ght_new(kpar,i,j) - ter(i,j) |
|
! meters agl |
|
cin(i,j,mkzh-2) = zrel(klfc) + ght_new(kpar,i,j) - ter(i,j) |
|
|
|
END DO |
|
END DO |
|
!$OMP END PARALLEL DO |
|
RETURN |
|
END SUBROUTINE DCAPECALC2D
|
|
|