w3fa06.f

Go to the documentation of this file.

C> @file
C> @brief Calculation of the lifted index.
C> @author James Howcroft @date 1978-07-01
 
C> Given the pressure,temperature and relative humidity of
C> an air parcel at some point in the atmosphere, calculate the
C> lifted index of the parcel. Lifted index is defined as the
C> temperature difference between the observed 500mb temperature and
C> the supposed temperature that the parcel would obtain if it were
C> lifted dry-adiabatically to saturation and then moved moist
C> adiabatically to the 500mb level.
C>
C> Program history log:
C> - James Howcroft 1978-07-01
C> - Ralph Jones 1989-01-24 Change to microsoft fortran 4.10.
C> - Ralph Jones 1990-06-08 Change to sun fortran 1.3.
C> - Ralph Jones 1991-03-29 Convert to silicongraphics fortran.
C> - Ralph Jones 1993-03-29 Add save statement.
C> - Ralph Jones 1995-09-25 Put in w3 library on cray.
C>
C> @param[in] P Parcel pressure in millibars.
C> @param[in] T Parcel temperataure in degrees celsius.
C> @param[in] RH Parcel relative humidity in percent.
C> @param[in] T5 Temperature at the 500mb level in deg. celsius.
C> @param[out] TLI Lifted index in degrees celsius
C> tli = 9.9999  iteration diverges; return to user program.
C>
C> @author James Howcroft @date 1978-07-01
      SUBROUTINE w3fa06 (P,T,RH,T5,TLI)
C
       SAVE
C
       DATA            eps   /0.5/
       DATA            kout  /  6/
C
  300  FORMAT ('  *** ITERATION NOT CONVERGING IN W3FA06 ***')
  350  FORMAT ('      INPUT PARAMS ARE:',4f15.8,/
     1         '      CALCULATIONS ARE',7e15.8)
C
       potemp(t,p) = (t+273.16)*((1000./p)**0.2857)
C
       eep(t,p,es) = exp((596.73-0.601*t)*((0.622*es)/(p-es))
     1                / (0.24*(t+273.16)))
C
       unpot(te,p) = (((p/1000.)**0.2857)*te)-273.16
C
       vapres(t)   = 6.11*exp(17.2694*t/(t+237.3))
C
       CALL w3fa01 (p,t,rh,td,plcl,tlcl)
       IF (plcl .GT. 500.)  GO TO 30
       IF (plcl .LT. 500.)  GO TO 20
         tli = t5 - tlcl
         GO TO 80
   20  CONTINUE
C        LCL IS ABOVE THE 500MB LVL
         tli = t5 - unpot((potemp(tlcl,plcl)),500.)
         GO TO 80
   30  CONTINUE
C        USE STACKPOLE ALGORITHM (JAM VOL 6/1967 PP 464-7) TO FIND TGES
C        SO THAT (TGES,500) IS ON SAME MOIST ADIABAT AS (TLCL,PLCL).
         es     = vapres(tlcl)
         thd    = potemp(tlcl,(plcl-es))
         theta  = thd * eep(tlcl,plcl,es)
C          THETA IS THE PSEUDO-EQUIV POTENTIAL TEMP THRU (PLCL,TLCL).
C      NOW FIND TEMP WHERE THETA INTERSECTS 500MB SFC.
C          INITIALIZE FOR STACKPOLIAN ITERATION
         tges = t5
         dtt  = 10.
         piin = 1./(0.5**0.2857)
         a    = 0.
         istp = 0
C          START ITERATION.
   40  CONTINUE
         istp = istp + 1
         IF (istp .GT. 200)  GO TO 50
         sva = vapres(tges)
         ax  = a
         a   = (tges+273.16)*piin * eep(tges,500.,sva) - theta
         IF (abs(a) .LT. eps)  GO TO 70
         dtt = dtt * 0.5
         IF (a*ax.LT.0.0)  dtt = -dtt
         tp   = tges + dtt
         sva  = vapres(tp)
         ap   = (tp+273.16)*piin * eep(tp,500.,sva) - theta
         IF (abs(ap) .LT. eps)  GO TO 60
C          FIND NEXT ESTIMATE, DTT IS ADJUSTMENT FROM OLD TO NEW TGES.
           dtt  = a*dtt/(a-ap)
           IF (abs(dtt).LT.0.01)  dtt = sign(0.01,dtt)
           tges = tges + dtt
           IF (tges .GT. 50)  tges = 50.
           GO TO 40
C
   50  CONTINUE
C        DISASTER SECTION
         WRITE (kout,300)
         WRITE (kout,350) p,t,rh,t5,theta,ax,a,ap,tges,tp,sva
           tli = 9.9999
           GO TO 80
   60  CONTINUE
         tges = tp
   70  CONTINUE
         tli  = t5 - tges
   80  CONTINUE
       RETURN
       END