! 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