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.
553 lines
21 KiB
553 lines
21 KiB
! The kind of code only a scientist could love. |
|
! TODO: The cape routine needs work to remove the GOTOs |
|
|
|
!====================================================================== |
|
! |
|
! !IROUTINE: TVIRTUAL -- 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) |
|
! |
|
|
|
REAL(KIND=8) FUNCTION tvirtual(temp,ratmix) |
|
IMPLICIT NONE |
|
REAL(KIND=8),INTENT(IN) :: temp,ratmix |
|
REAL(KIND=8),PARAMETER :: EPS = .622D0 |
|
|
|
tvirtual = temp*(EPS + ratmix)/(EPS*(1.D0 + ratmix)) |
|
RETURN |
|
END FUNCTION tvirtual |
|
|
|
REAL(KIND=8) FUNCTION tonpsadiabat(thte,prs,PSADITHTE,PSADIPRS,PSADITMK,GAMMA,& |
|
throw_exception) |
|
IMPLICIT NONE |
|
EXTERNAL throw_exception |
|
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 |
|
|
|
REAL(KIND=8) :: fracjt |
|
REAL(KIND=8) :: fracjt2 |
|
REAL(KIND=8) :: fracip |
|
REAL(KIND=8) :: fracip2 |
|
|
|
INTEGER :: ip, ipch, jt, jtch |
|
|
|
! 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 |
|
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 |
|
|
|
! JT = -1 |
|
!213 CONTINUE |
|
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 |
|
|
|
! IP = -1 |
|
!215 CONTINUE |
|
IF (jt.EQ.-1 .OR. ip.EQ.-1) THEN |
|
! Need an exception here |
|
CALL throw_exception('capecalc3d: ','Outside of lookup table bounds. prs,thte=',prs,thte) |
|
!STOP |
|
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 |
|
CALL throw_exception('capecalc3d: ','Tried to access missing temperature in lookup table.',& |
|
'Prs and Thte probably unreasonable. prs,thte=',prs,thte) |
|
!STOP |
|
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 |
|
|
|
|
|
! 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,miy,mjx,mkzh,ter_follow) |
|
|
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs |
|
REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: sfp |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(OUT) :: pf |
|
INTEGER,INTENT(IN) :: ter_follow,miy,mjx,mkzh |
|
|
|
INTEGER :: i,j,k |
|
|
|
! do j=1,mjx-1 Artifact of MM5 |
|
DO j = 1,mjx |
|
! do i=1,miy-1 staggered grid |
|
DO i = 1,miy |
|
DO k = 1,mkzh |
|
IF (k .EQ. mkzh) THEN |
|
! terrain-following data |
|
IF (ter_follow .EQ. 1) THEN |
|
pf(i,j,k) = sfp(i,j) |
|
! pressure-level data |
|
ELSE |
|
pf(i,j,k) = .5D0 * (3.D0*prs(i,j,k) - prs(i,j,k-1)) |
|
END IF |
|
ELSE |
|
pf(i,j,k) = .5D0 * (prs(i,j,k+1) + prs(i,j,k)) |
|
END IF |
|
END DO |
|
END DO |
|
END DO |
|
|
|
RETURN |
|
END SUBROUTINE dpfcalc |
|
|
|
!====================================================================== |
|
! |
|
! !IROUTINE: capecalc3d -- Calculate CAPE and CIN |
|
! |
|
! !DESCRIPTION: |
|
! |
|
! If i3dflag=1, 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). If i3dflag=0, then it |
|
! 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. |
|
! |
|
! In the case of i3dflag=0, |
|
! CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of |
|
! the cape and cin arrays. Also, if i3dflag=0, 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 |
|
! |
|
SUBROUTINE f_computecape(prs,tmk,qvp,ght,ter,sfp,cape,cin,& |
|
PSADITHTE,PSADIPRS,PSADITMK,cmsg,i3dflag,ter_follow,& |
|
throw_exception,miy,mjx,mkzh) |
|
|
|
IMPLICIT NONE |
|
EXTERNAL throw_exception |
|
INTEGER,INTENT(IN) :: miy,mjx,mkzh,i3dflag,ter_follow |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: prs |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: tmk |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: qvp |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(IN) :: ght |
|
REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) :: ter |
|
REAL(KIND=8),DIMENSION(miy,mjx),INTENT(IN) ::sfp |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(INOUT) :: cape |
|
REAL(KIND=8),DIMENSION(miy,mjx,mkzh),INTENT(INOUT) :: cin |
|
REAL(KIND=8),DIMENSION(150),INTENT(IN) :: PSADITHTE,PSADIPRS |
|
REAL(KIND=8),DIMENSION(150,150),INTENT(IN) :: PSADITMK |
|
REAL(KIND=8),INTENT(IN) :: cmsg |
|
|
|
! local variables |
|
INTEGER :: i,j,k,ilcl,kel,kk,klcl,klev,klfc,kmax,kpar,kpar1,kpar2 |
|
REAL(KIND=8) :: davg,ethmax,q,t,p,e,eth,tlcl,zlcl |
|
REAL(KIND=8) :: pavg,tvirtual,p1,p2,pp1,pp2,th,totthe,totqvp,totprs |
|
REAL(KIND=8) :: cpm,deltap,ethpari,gammam,ghtpari,qvppari,prspari,tmkpari |
|
REAL(KIND=8) :: facden,fac1,fac2,qvplift,tmklift,tvenv,tvlift,ghtlift |
|
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(miy,mjx,mkzh) :: prsf |
|
|
|
! constants |
|
INTEGER,PARAMETER :: IUP = 6 |
|
REAL(KIND=8),PARAMETER :: CELKEL = 273.15d0 |
|
REAL(KIND=8),PARAMETER :: GRAV = 9.81d0 |
|
! hpa |
|
REAL(KIND=8),PARAMETER :: EZERO = 6.112d0 |
|
REAL(KIND=8),PARAMETER :: ESLCON1 = 17.67d0 |
|
REAL(KIND=8),PARAMETER :: ESLCON2 = 29.65d0 |
|
REAL(KIND=8),PARAMETER :: EPS = 0.622d0 |
|
! j/k/kg |
|
REAL(KIND=8),PARAMETER :: RGAS = 287.04d0 |
|
! j/k/kg note: not using bolton's value of 1005.7 |
|
REAL(KIND=8),PARAMETER :: CP = 1004.d0 |
|
REAL(KIND=8),PARAMETER :: GAMMA = RGAS/CP |
|
! cp_moist=cp*(1.+cpmd*qvp) |
|
REAL(KIND=8),PARAMETER :: CPMD = .887d0 |
|
! rgas_moist=rgas*(1.+rgasmd*qvp) |
|
REAL(KIND=8),PARAMETER :: RGASMD = .608d0 |
|
! gamma_moist=gamma*(1.+gammamd*qvp) |
|
REAL(KIND=8),PARAMETER :: GAMMAMD = RGASMD - CPMD |
|
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 |
|
! k |
|
REAL(KIND=8),PARAMETER :: THTECON1 = 3376.d0 |
|
REAL(KIND=8),PARAMETER :: THTECON2 = 2.54d0 |
|
REAL(KIND=8),PARAMETER :: THTECON3 = .81d0 |
|
|
|
! To get rid of compiler warnings |
|
tmkpari = 0 |
|
qvppari = 0 |
|
klev = 0 |
|
klcl = 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 dpfcalc(prs,sfp,prsf,miy,mjx,mkzh,ter_follow) |
|
|
|
! before looping, set lookup table for getting temperature on |
|
! a pseudoadiabat. |
|
|
|
!call dlookup_table(psadithte,psadiprs,psaditmk,psafile) |
|
|
|
! do j=1,mjx-1 |
|
DO j = 1,mjx |
|
! do i=1,miy-1 |
|
DO i = 1,miy |
|
cape(i,j,1) = 0.d0 |
|
cin(i,j,1) = 0.d0 |
|
|
|
IF (i3dflag .EQ. 1) THEN |
|
kpar1 = 2 |
|
kpar2 = mkzh |
|
ELSE |
|
! find parcel with max theta-e in lowest 3 km agl. |
|
ethmax = -1.d0 |
|
DO k = mkzh,1,-1 |
|
IF (ght(i,j,k)-ter(i,j) .LT. 3000.d0) THEN |
|
q = MAX(qvp(i,j,k), 1.d-15) |
|
t = tmk(i,j,k) |
|
p = prs(i,j,k) |
|
e = q*p/(EPS + q) |
|
tlcl = TLCLC1/ (LOG(t**TLCLC2/e)-TLCLC3) + TLCLC4 |
|
eth = t * (1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q))*& |
|
EXP((THTECON1/tlcl - THTECON2)*q*(1.d0 + THTECON3*q)) |
|
IF (eth .GT. ethmax) THEN |
|
klev = k |
|
ethmax = eth |
|
END IF |
|
END IF |
|
END DO |
|
kpar1 = klev |
|
kpar2 = klev |
|
|
|
! Establish average properties of that parcel |
|
! (over depth of approximately davg meters) |
|
|
|
! davg=.1 |
|
davg = 500.d0 |
|
pavg = davg*prs(i,j,kpar1)*& |
|
GRAV/(RGAS*tvirtual(tmk(i,j,kpar1), qvp(i,j,kpar1))) |
|
p2 = MIN(prs(i,j,kpar1)+.5d0*pavg, prsf(i,j,mkzh)) |
|
p1 = p2 - pavg |
|
totthe = 0.d0 |
|
totqvp = 0.d0 |
|
totprs = 0.d0 |
|
DO k = mkzh,2,-1 |
|
IF (prsf(i,j,k) .LE. p1) GOTO 35 |
|
IF (prsf(i,j,k-1) .GE. p2) GOTO 34 |
|
p = prs(i,j,k) |
|
pup = prsf(i,j,k) |
|
pdn = prsf(i,j,k-1) |
|
q = MAX(qvp(i,j,k),1.d-15) |
|
th = tmk(i,j,k)*(1000.d0/p)**(GAMMA*(1.d0 + GAMMAMD*q)) |
|
pp1 = MAX(p1,pdn) |
|
pp2 = MIN(p2,pup) |
|
IF (pp2 .GT. pp1) THEN |
|
deltap = pp2 - pp1 |
|
totqvp = totqvp + q*deltap |
|
totthe = totthe + th*deltap |
|
totprs = totprs + deltap |
|
END IF |
|
34 CONTINUE |
|
END DO |
|
35 CONTINUE |
|
qvppari = totqvp/totprs |
|
tmkpari = (totthe/totprs)*& |
|
(prs(i,j,kpar1)/1000.d0)**(GAMMA*(1.d0+GAMMAMD*qvp(i,j,kpar1))) |
|
END IF |
|
|
|
DO kpar = kpar1,kpar2 |
|
|
|
! Calculate temperature and moisture properties of parcel |
|
! (note, qvppari and tmkpari already calculated above for 2d case.) |
|
|
|
IF (i3dflag .EQ. 1) THEN |
|
qvppari = qvp(i,j,kpar) |
|
tmkpari = tmk(i,j,kpar) |
|
END IF |
|
prspari = prs(i,j,kpar) |
|
ghtpari = ght(i,j,kpar) |
|
gammam = GAMMA * (1.d0 + GAMMAMD*qvppari) |
|
cpm = CP * (1.d0 + CPMD*qvppari) |
|
|
|
e = MAX(1.d-20,qvppari*prspari/(EPS + qvppari)) |
|
tlcl = TLCLC1/(LOG(tmkpari**TLCLC2/e) - TLCLC3) + TLCLC4 |
|
ethpari = tmkpari*(1000.d0/prspari)**(GAMMA*(1.d0 + GAMMAMD*qvppari))*& |
|
EXP((THTECON1/tlcl - THTECON2)*qvppari*(1.d0 + THTECON3*qvppari)) |
|
zlcl = ghtpari + (tmkpari - tlcl)/(GRAV/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 (ghtpari .GE. zlcl) THEN |
|
! Initial parcel already saturated or supersaturated. |
|
ilcl = 2 |
|
klcl = 1 |
|
END IF |
|
|
|
DO k = kpar,1,-1 |
|
! For arrays that go bottom to top |
|
33 kk = kk + 1 |
|
|
|
! Model level is below lcl |
|
IF (ght(i,j,k).lt.zlcl) THEN |
|
qvplift = qvppari |
|
tmklift = tmkpari - GRAV/cpm*(ght(i,j,k) - ghtpari) |
|
tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) |
|
tvlift = tvirtual(tmklift, qvplift) |
|
ghtlift = ght(i,j,k) |
|
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. |
|
tmklift = tlcl |
|
qvplift = qvppari |
|
facden = ght(i,j,k) - ght(i,j,k+1) |
|
fac1 = (zlcl-ght(i,j,k+1))/facden |
|
fac2 = (ght(i,j,k)-zlcl)/facden |
|
tmkenv = tmk(i,j,k+1)*fac2 + tmk(i,j,k)*fac1 |
|
qvpenv = qvp(i,j,k+1)*fac2 + qvp(i,j,k)*fac1 |
|
tvenv = tvirtual(tmkenv, qvpenv) |
|
tvlift = tvirtual(tmklift, qvplift) |
|
ghtlift = zlcl |
|
ilcl = 1 |
|
ELSE |
|
tmklift = tonpsadiabat(ethpari, prs(i,j,k), PSADITHTE, PSADIPRS,& |
|
PSADITMK, GAMMA, throw_exception) |
|
eslift = EZERO*EXP(ESLCON1*(tmklift - CELKEL)/(tmklift - ESLCON2)) |
|
qvplift = EPS*eslift/(prs(i,j,k) - eslift) |
|
tvenv = tvirtual(tmk(i,j,k), qvp(i,j,k)) |
|
tvlift = tvirtual(tmklift, qvplift) |
|
ghtlift = ght(i,j,k) |
|
END IF |
|
! Buoyancy |
|
buoy(kk) = GRAV*(tvlift - tvenv)/tvenv |
|
zrel(kk) = ghtlift - ghtpari |
|
|
|
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 |
|
GOTO 33 |
|
END IF |
|
END DO |
|
|
|
kmax = kk |
|
IF (kmax .GT. 150) THEN |
|
! Need an exception here |
|
CALL throw_exception('capecalc3d: kmax got too big. kmax=',kmax) |
|
!STOP |
|
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. |
|
|
|
DO k = kmax,klcl,-1 |
|
IF (buoy(k) .GE. 0.d0) THEN |
|
! k of equilibrium level |
|
kel = k |
|
GOTO 50 |
|
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. |
|
|
|
! cape(i,j,kpar) = -0.1d0 |
|
! cin(i,j,kpar) = -0.1d0 |
|
cape(i,j,kpar) = cmsg |
|
cin(i,j,kpar) = cmsg |
|
klfc = kmax |
|
|
|
GOTO 102 |
|
|
|
50 CONTINUE |
|
|
|
! 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) = -0.1d0 |
|
IF (cape(i,j,kpar) .LT. 100.d0) cin(i,j,kpar) = cmsg |
|
102 CONTINUE |
|
|
|
END DO |
|
|
|
IF (i3dflag .EQ. 0) THEN |
|
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) + ghtpari - ter(i,j) |
|
! meters agl |
|
cin(i,j,mkzh-2) = zrel(klfc) + ghtpari - ter(i,j) |
|
|
|
ENDIF |
|
END DO |
|
END DO |
|
|
|
RETURN |
|
END SUBROUTINE f_computecape
|
|
|