rip_cape.f

c======================================================================
c
c !IROUTINE: capecalc3d -- Calculate CAPE and CIN
c
c !DESCRIPTION:
c
c   If i3dflag=1, this routine calculates CAPE and CIN (in m**2/s**2,
c   or J/kg) for every grid point in the entire 3D domain (treating
c   each grid point as a parcel).  If i3dflag=0, then it
c   calculates CAPE and CIN only for the parcel with max theta-e in
c   the column, (i.e. something akin to Colman's MCAPE).  By "parcel",
c   we mean a 500-m deep parcel, with actual temperature and moisture
c   averaged over that depth.
c
c   In the case of i3dflag=0,
c   CAPE and CIN are 2D fields that are placed in the k=mkzh slabs of
c   the cape and cin arrays.  Also, if i3dflag=0, LCL and LFC heights
c   are put in the k=mkzh-1 and k=mkzh-2 slabs of the cin array.
c
c ASSUMPTIONS:
c
c !REVISION HISTORY:
c     2005-May-15 - Mark T. Stoelinga - oringinal version from RIP4
c     2005-Nov-28 - J. Schramm - modified to run outside of RIP4 with
c                                NCL
c
c !INTERFACE:
c ------------------------------------------------------------------
C NCLFORTSTART
      SUBROUTINE DCAPECALC3D(PRS,TMK,QVP,GHT,TER,SFP,CAPE,CIN,MIY,MJX,
     +                       MKZH,I3DFLAG,TER_FOLLOW,PSAFILE)
c
      IMPLICIT NONE
      INTEGER MIY,MJX,MKZH,I3DFLAG,TER_FOLLOW
      DOUBLE PRECISION PRS(MIY,MJX,MKZH)
      DOUBLE PRECISION TMK(MIY,MJX,MKZH)
      DOUBLE PRECISION QVP(MIY,MJX,MKZH)
      DOUBLE PRECISION GHT(MIY,MJX,MKZH)
      DOUBLE PRECISION TER(MIY,MJX)
      DOUBLE PRECISION SFP(MIY,MJX)
      DOUBLE PRECISION CAPE(MIY,MJX,MKZH)
      DOUBLE PRECISION CIN(MIY,MJX,MKZH)
      CHARACTER*(*) PSAFILE
C NCLEND
c Local variables
      INTEGER I,J,K,ILCL,IUP,KEL,KK,KLCL,KLEV,KLFC,KMAX,KPAR,KPAR1,KPAR2
      INTEGER MM,NN
      DOUBLE PRECISION DAVG,ETHMAX,Q,T,P,E,ETH,TLCL,ZLCL
      DOUBLE PRECISION CP,EPS,GAMMA,GAMMAMD,RGAS,RGASMD,TLCLC1,TLCLC2,
     +                 TLCLC3,TLCLC4
      DOUBLE PRECISION CPMD,THTECON1,THTECON2,THTECON3
      DOUBLE PRECISION CELKEL,EZERO,ESLCON1,ESLCON2,GRAV,THECON1,
     +                 THECON2,THECON3
      DOUBLE PRECISION PAVG,VIRTUAL,P1,P2,PP1,PP2,TH,TOTTHE,TOTQVP,
     +                 TOTPRS
      DOUBLE PRECISION CPM,DELTAP,ETHPARI,GAMMAM,GHTPARI,QVPPARI,
     +                 PRSPARI,TMKPARI
      DOUBLE PRECISION FACDEN,FAC1,FAC2,QVPLIFT,TMKLIFT,TVENV,TVLIFT,
     +                 GHTLIFT
      DOUBLE PRECISION ESLIFT,TMKENV,QVPENV,TONPSADIABAT
      DOUBLE PRECISION BENAMIN,DZ,Z,PUP,PDN
      DOUBLE PRECISION BUOY(150),ZREL(150),BENACCUM(150),
     +                 PRSF(MIY,MJX,MKZH)
      DOUBLE PRECISION PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
c
C The comments were taken from a Mark Stoelinga email, 23 Apr 2007,
C in response to a user getting the "Outside of lookup table bounds"
C error message. 
C
C TMKPARI  - Initial temperature of parcel, K
C    Values of 300 okay. (Not sure how much from this you can stray.)
C
C PRSPARI - Initial pressure of parcel, hPa
C    Values of 980 okay. (Not sure how much from this you can stray.)
C
C THTECON1, THTECON2, THTECON3
C     These are all constants, the first in K and the other two have
C     no units.  Values of 3376, 2.54, and 0.81 were stated as being
C     okay.
C
C TLCL - The temperature at the parcel's lifted condensation level, K
C        should be a reasonable atmospheric temperature around 250-300 K
C        (398 is "way too high")
C
C QVPPARI - The initial water vapor mixing ratio of the parcel,
C           kg/kg (should range from 0.000 to 0.025)
C

c Constants
      IUP = 6
      CELKEL = 273.15D0
      GRAV = 9.81D0
C hPa
      EZERO = 6.112D0
      ESLCON1 = 17.67D0
      ESLCON2 = 29.65D0
      EPS = 0.622D0
C J/K/kg
      RGAS = 287.04D0
C  J/K/kg  Note: not using Bolton's value of 1005.7
      CP = 1004.D0
      GAMMA = RGAS/CP
C  cp_moist=cp*(1.+cpmd*qvp)
      CPMD = .887D0
C  rgas_moist=rgas*(1.+rgasmd*qvp)
      RGASMD = .608D0
C  gamma_moist=gamma*(1.+gammamd*qvp)
      GAMMAMD = RGASMD - CPMD
      TLCLC1 = 2840.D0
      TLCLC2 = 3.5D0
      TLCLC3 = 4.805D0
      TLCLC4 = 55.D0
C  K
      THTECON1 = 3376.D0
      THTECON2 = 2.54D0
      THTECON3 = .81D0
c
c  Calculated the pressure at full sigma levels (a set of pressure
c  levels that bound the layers represented by the vertical grid points)

      CALL DPFCALC(PRS,SFP,PRSF,MIY,MJX,MKZH,TER_FOLLOW)
c
c  Before looping, set lookup table for getting temperature on
c  a pseudoadiabat.
c
      CALL DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,PSAFILE)
c
C   do j=1,mjx-1
      DO J = 1,MJX
C   do i=1,miy-1
          DO I = 1,MIY
              CAPE(I,J,1) = 0.D0
              CIN(I,J,1) = 0.D0
c
              IF (I3DFLAG.EQ.1) THEN
                  KPAR1 = 2
                  KPAR2 = MKZH
              ELSE
c
c      Find parcel with max theta-e in lowest 3 km AGL.
c
                  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
c
c      Establish average properties of that parcel
c         (over depth of approximately davg meters)
c
c         davg=.1
                  DAVG = 500.D0
                  PAVG = DAVG*PRS(I,J,KPAR1)*GRAV/
     +                   (RGAS*VIRTUAL(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) GO TO 35
                      IF (PRSF(I,J,K-1).GE.P2) GO TO 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
c
              DO KPAR = KPAR1,KPAR2
c
c   Calculate temperature and moisture properties of parcel
c     (Note, qvppari and tmkpari already calculated above for 2D case.)
c
                  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)
c
                  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)
c
c   Calculate buoyancy and relative height of lifted parcel at
c   all levels, and store in bottom up arrays.  Add a level at the LCL,
c   and at all points where buoyancy is zero.
c
C  for arrays that go bottom to top
                  KK = 0
                  ILCL = 0
                  IF (GHTPARI.GE.ZLCL) THEN
c
c      initial parcel already saturated or supersaturated.
c
                      ILCL = 2
                      KLCL = 1
                  END IF
                  DO K = KPAR,1,-1
C  for arrays that go bottom to top
   33                 KK = KK + 1
C  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 = VIRTUAL(TMK(I,J,K),QVP(I,J,K))
                          TVLIFT = VIRTUAL(TMKLIFT,QVPLIFT)
                          GHTLIFT = GHT(I,J,K)
                      ELSE IF (GHT(I,J,K).GE.ZLCL .AND. ILCL.EQ.0) THEN
c
c     This model level and previous model level straddle the LCL,
c     so first create a new level in the bottom-up array, at the LCL.
c
                          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 = VIRTUAL(TMKENV,QVPENV)
                          TVLIFT = VIRTUAL(TMKLIFT,QVPLIFT)
                          GHTLIFT = ZLCL
                          ILCL = 1
                      ELSE
                          TMKLIFT = TONPSADIABAT(ETHPARI,PRS(I,J,K),
     +                              PSADITHTE,PSADIPRS,PSADITMK,GAMMA)
                          ESLIFT = EZERO*EXP(ESLCON1* (TMKLIFT-CELKEL)/
     +                             (TMKLIFT-ESLCON2))
                          QVPLIFT = EPS*ESLIFT/ (PRS(I,J,K)-ESLIFT)
                          TVENV = VIRTUAL(TMK(I,J,K),QVP(I,J,K))
                          TVLIFT = VIRTUAL(TMKLIFT,QVPLIFT)
                          GHTLIFT = GHT(I,J,K)
                      END IF
C  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
c
c   Parcel ascent curve crosses sounding curve, so create a new level
c   in the bottom-up array at the crossing.
c
                          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
                          GO TO 33
                      END IF
                  END DO
                  KMAX = KK
                  IF (KMAX.GT.150) THEN
                      print *,
     +                  'capecalc3d: kmax got too big. kmax=',KMAX
                      STOP
                  END IF
c
c   If no LCL was found, set klcl to kmax.  It is probably not really
c   at kmax, but this will make the rest of the routine behave
c   properly.
c
                  IF (ILCL.EQ.0) KLCL=KMAX
c
c   Get the accumulated buoyant energy from the parcel's starting
c   point, at all levels up to the top level.
c
                  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
c
c     Determine equilibrium level (EL), which we define as the highest
c     level of non-negative buoyancy above the LCL. Note, this may be
c     the top level if the parcel is still buoyant there.
c
                  DO K = KMAX,KLCL,-1
                      IF (BUOY(K).GE.0.D0) THEN
C  k of equilibrium level
                          KEL = K
                          GO TO 50
                      END IF
                  END DO
c
c   If we got through that loop, then there is no non-negative
c   buoyancy above the LCL in the sounding.  In these situations,
c   both CAPE and CIN will be set to -0.1 J/kg.  Also, where CAPE is
c   non-zero, CAPE and CIN will be set to a minimum of +0.1 J/kg, so
c   that the zero contour in either the CIN or CAPE fields will
c   circumscribe regions of non-zero CAPE.
c
                  CAPE(I,J,KPAR) = -0.1D0
                  CIN(I,J,KPAR) = -0.1D0
                  KLFC = KMAX
c
                  GO TO 102
c
   50             CONTINUE
c
c   If there is an equilibrium level, then CAPE is positive.  We'll
c   define the level of free convection (LFC) as the point below the
c   EL, but at or above the LCL, where accumulated buoyant energy is a
c   minimum.  The net positive area (accumulated buoyant energy) from
c   the LFC up to the EL will be defined as the CAPE, and the net
c   negative area (negative of accumulated buoyant energy) from the
c   parcel starting point to the LFC will be defined as the convective
c   inhibition (CIN).
c
c   First get the LFC according to the above definition.
c
                  BENAMIN = 9D9
                  KLFC = KMAX
                  DO K = KLCL,KEL
                      IF (BENACCUM(K).LT.BENAMIN) THEN
                          BENAMIN = BENACCUM(K)
                          KLFC = K
                      END IF
                  END DO
c
c   Now we can assign values to cape and cin
c
                  CAPE(I,J,KPAR) = MAX(BENACCUM(KEL)-BENAMIN,0.1D0)
                  CIN(I,J,KPAR) = MAX(-BENAMIN,0.1D0)
c
c   CIN is uninteresting when CAPE is small (< 100 J/kg), so set
c   CIN to -.1 in that case.
c
                  IF (CAPE(I,J,KPAR).LT.100.D0) CIN(I,J,KPAR) = -0.1D0
  102             CONTINUE
c
              END DO
c
              IF (I3DFLAG.EQ.0) THEN
                  CAPE(I,J,MKZH) = CAPE(I,J,KPAR1)
                  CIN(I,J,MKZH) = CIN(I,J,KPAR1)
C  meters AGL
                  CIN(I,J,MKZH-1) = ZREL(KLCL) + GHTPARI - TER(I,J)
C  meters AGL
                  CIN(I,J,MKZH-2) = ZREL(KLFC) + GHTPARI - TER(I,J)
              END IF
c
          END DO
      END DO
c
      RETURN
      END
c                                                                     c
c*********************************************************************c
c                                                                     c
C NCLFORTSTART
      FUNCTION TONPSADIABAT(THTE,PRS,PSADITHTE,PSADIPRS,PSADITMK,GAMMA)
      DOUBLE PRECISION TONPSADIABAT
      DOUBLE PRECISION THTE
      DOUBLE PRECISION PRS
      DOUBLE PRECISION PSADITHTE
      DOUBLE PRECISION PSADIPRS
      DOUBLE PRECISION PSADITMK
      DOUBLE PRECISION GAMMA
      DOUBLE PRECISION FRACJT
      DOUBLE PRECISION FRACJT2
      DOUBLE PRECISION FRACIP
      DOUBLE PRECISION FRACIP2
      DIMENSION PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)
C NCLEND
c                                                                     c
c   This function gives the temperature (in K) on a moist adiabat
c   (specified by thte in K) given pressure in hPa.  It uses a
c   lookup table, with data that was generated by the Bolton (1980)
c   formula for theta_e.
c
c     First check if pressure is less than min pressure in lookup table.
c     If it is, assume parcel is so dry that the given theta-e value can
c     be interpretted as theta, and get temperature from the simple dry
c     theta formula.
c
      IF (PRS.LE.PSADIPRS(150)) THEN
          TONPSADIABAT = THTE* (PRS/1000.D0)**GAMMA
          RETURN
      END IF
c
c   Otherwise, look for the given thte/prs point in the lookup table.
c
      DO JTCH = 1,150 - 1
          IF (THTE.GE.PSADITHTE(JTCH) .AND.
     +        THTE.LT.PSADITHTE(JTCH+1)) THEN
              JT = JTCH
              GO TO 213
          END IF
      END DO
      JT = -1
  213 CONTINUE
      DO IPCH = 1,150 - 1
          IF (PRS.LE.PSADIPRS(IPCH) .AND. PRS.GT.PSADIPRS(IPCH+1)) THEN
              IP = IPCH
              GO TO 215
          END IF
      END DO
      IP = -1
  215 CONTINUE
      IF (JT.EQ.-1 .OR. IP.EQ.-1) THEN
         print *,'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
          print *,'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)
c
      RETURN
      END
c                                                                     c
c*********************************************************************c
      SUBROUTINE DLOOKUP_TABLE(PSADITHTE,PSADIPRS,PSADITMK,FNAME)
      DOUBLE PRECISION PSADITHTE
      DOUBLE PRECISION PSADIPRS
      DOUBLE PRECISION PSADITMK
c   Set up lookup table for getting temperature on a pseudoadiabat.
c   (Borrow the unit number for the stationlist, just for the moment.)
c
C      CHARACTER*15 FNAME
      CHARACTER*(*) FNAME
      DIMENSION PSADITHTE(150),PSADIPRS(150),PSADITMK(150,150)

C      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
          WRITE (IUP,FMT=*)
     +      'Number of pressure or theta_e levels in lookup table'
          WRITE (IUP,FMT=*) 'file not = 150.  Check lookup table file.'
          STOP
      END IF
      READ (IUSTNLIST,FMT=173) (PSADITHTE(JT),JT=1,NTHTE)
      READ (IUSTNLIST,FMT=173) (PSADIPRS(IP),IP=1,NPRS)
      READ (IUSTNLIST,FMT=173) ((PSADITMK(IP,JT),IP=1,NPRS),JT=1,NTHTE)
  173 FORMAT (5D15.7)
      CLOSE (IUSTNLIST)

      RETURN
      END
c                                                                     c
c*********************************************************************c
c                                                                     c
      SUBROUTINE DPFCALC(PRS,SFP,PF,MIY,MJX,MKZH,TER_FOLLOW)
      DOUBLE PRECISION PRS
      DOUBLE PRECISION SFP
      DOUBLE PRECISION PF
c
c     Historically, this routine calculated the pressure at full sigma
c     levels when RIP was specifically designed for MM4/MM5 output.
c     With the new generalized RIP (Feb '02), this routine is still
c     intended to calculate a set of pressure levels that bound the
c     layers represented by the vertical grid points, although no such
c     layer boundaries are assumed to be defined.  The routine simply
c     uses the midpoint between the pressures of the vertical grid
c     points as the bounding levels.  The array only contains mkzh
c     levels, so the pressure of the top of the uppermost layer is
c     actually excluded.  The kth value of pf is the lower bounding
c     pressure for the layer represented by kth data level.  At the
c     lower bounding level of the lowest model layer, it uses the
c     surface pressure, unless the data set is pressure-level data, in
c     which case it assumes the lower bounding pressure level is as far
c     below the lowest vertical level as the upper bounding pressure
c     level is above.
c
      DIMENSION PRS(MIY,MJX,MKZH),SFP(MIY,MJX),PF(MIY,MJX,MKZH)
      INTEGER TER_FOLLOW
c
C  do j=1,mjx-1  Artifact of MM5
      DO J = 1,MJX
C  do i=1,miy-1  staggered grid
          DO I = 1,MIY
              DO K = 1,MKZH
                  IF (K.EQ.MKZH) THEN
C  terrain-following data
                      IF (TER_FOLLOW.EQ.1) THEN
                          PF(I,J,K) = SFP(I,J)
C  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
c
      RETURN
      END
c*********************************************************************c
c
C NCLFORTSTART
      FUNCTION VIRTUAL(TEMP,RATMIX)
      DOUBLE PRECISION VIRTUAL
      DOUBLE PRECISION EPS
c
c   This function returns virtual temperature in K, given temperature
c      in K and mixing ratio in kg/kg.
c
      DOUBLE PRECISION TEMP,RATMIX
C NCLEND
      EPS = 0.622D0
      VIRTUAL = TEMP* (EPS+RATMIX)/ (EPS* (1.D0+RATMIX))
      RETURN
      END