      SUBROUTINE BUOYANA(Z,P,PF,T,Q,PRESSG,ILEV,ILEVM,ILG,IL1,IL2,
     1                   MSG,RRL,GRAV,CPRES,EPS1,EPS2,A,B,RGAS,RGASV,
     2                   TVFA,RGOCP,MAXI,LCL,LNB,LPBL,TL,CAPE,CAPEP,
     3                   TC,QSTC,HMN,PBLT)
C
C     * NOV 24/2006 - M.LAZARE.    - CALLS NEW TINCLD3 INSTEAD OF TINCLD2.
C     *                            - SELECTIVE PROMOTION OF VARIABLES TO
C     *                              REAL*8 TO SUPPORT 32-BIT MODE.`
C     * JUN 15/2006 - M.LAZARE.    NEW VERSION FOR GCM15F:
C     *                            - NOW USES INTERNAL WORK ARRAYS
C     *                              INSTEAD OF PASSED WORKSPACE.    
C     *                            - "ZERO" DATA CONSTANT ADDED
C     *                              AND USED IN CALL TO INTRINSICS. 
C     * MAY 04/2006 - M.LAZARE/    PREVIOUS VERSION BUOYAN9 FOR GCM15E:
C     *               K.VONSALZEN. - CALLS NEW TINCLD2.
C     *                              THIS REQUIRES PASSING IN EPS2.
C     * MAY 05/2001 - K.VONSALZEN. PREVIOUS VERSION BUOYAN8 FOR GCM15D.
C                                                   
C     * CALCULATES PROPERTIES OF UNDILUTED AIR PARCEL ASCENDING
C     * FROM THE PBL TO ITS LEVEL OF NEUTRAL BUOYANCY.
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)

      PARAMETER (ITRIGG=0)
C
C     * I/O FIELDS:
C
      REAL, DIMENSION(ILG,ILEV)   :: Q,T,TC,QC,QSTC,Z,PF,P,HMN
      REAL, DIMENSION(ILG)        :: TL,CAPE,CAPEP,PRESSG,PBLT
      INTEGER,DIMENSION(ILG)      :: LCL,LNB,MAXI,LPBL
C
C     * INTERNAL WORK FIELDS (THE FIRST SET ARE PROMOTED TO
C     * INTEGRATE WITH THE NEW TINCLD3, IE FROM SHALLOW):
C
      REAL*8, DIMENSION(ILG,ILEV) :: TC8,TDRY,HMNC,QC8,QSTC8,Z8,P8
      REAL, DIMENSION(ILG,ILEV)   :: WC,BUOY,FACT
      REAL, DIMENSION(ILG)        :: QM,HMNM,MAXGR
C
      INTEGER,DIMENSION(ILG)      :: ISKIP,IDOC,LTM
C
      DATA ZERO    / 0.     /
C-----------------------------------------------------------------------
      DO 20 IL=IL1,IL2                      
         LPBL (IL)=NINT(PBLT(IL))                                
         MAXI (IL)=ILEV
         IDOC (IL)=1
         ISKIP(IL)=0
         MAXGR(IL)=9.E+20
         CAPE (IL)=0.
         CAPEP(IL)=0.
  20  CONTINUE
C
      DO 100 L=1,ILEV
      DO 100 IL=IL1,IL2  
         HMN(IL,L)=CPRES*T(IL,L)+GRAV*Z(IL,L)+RRL*Q(IL,L)
 100  CONTINUE    
C
      DO 120 IL=IL1,IL2                                                     
         QM  (IL)=Q   (IL,MAXI(IL))
         HMNM(IL)=HMN (IL,MAXI(IL))
 120  CONTINUE
C
C    * CALCULATE TEMPERATURE OF UNDILUTED AIR PARCEL.
C
      DO 140 L=MSG+1,ILEV 
      DO 140 IL=IL1,IL2 
         QC  (IL,L)=QM  (IL)
         HMNC(IL,L)=HMNM(IL)
 140  CONTINUE
C
      DO L=1,ILEV
      DO IL=IL1,IL2
        Z8   (IL,L)=Z   (IL,L)
        P8   (IL,L)=P   (IL,L)
        QC8  (IL,L)=QC  (IL,L)
      ENDDO
      ENDDO
C 
      CALL TINCLD3(TC8,ILG,ILEV,MSG-1,IL1,IL2,MAXI,TDRY,HMNC,QC8,QSTC8,
     1             Z8,P8,RRL,GRAV,CPRES,EPS1,EPS2,A,B,ISKIP)
C
      DO L=1,ILEV
      DO IL=IL1,IL2
        TC  (IL,L)=TC8  (IL,L)
        QSTC(IL,L)=QSTC8(IL,L)
      ENDDO
      ENDDO
C
C    * BUOYANCY OF UNDIULTED AIR PARCEL.
C
      DO 200 L=MSG+1,ILEV
      DO 200 IL=IL1,IL2 
         BUOY(IL,L)=-1.E-10
 200  CONTINUE
C
      L=ILEV
      DO 210 IL=IL1,IL2 
         IF ( L.LE.MAXI(IL) ) THEN
             TV =T(IL,L)*(1.+TVFA*Q(IL,L))
             TPV= TC(IL,L)*(1.+TVFA*QC(IL,L)
     1           -(1.+TVFA)*MAX(QC(IL,L)-QSTC(IL,L),ZERO))
             BUOY(IL,L)=(TPV-TV)/TV
             FACT(IL,L)=RGAS*BUOY(IL,L)*T(IL,L)
     1                 *LOG(0.01*PRESSG(IL)/PF(IL,L))
             BUOY(IL,L)=GRAV*BUOY(IL,L)
         ENDIF
 210  CONTINUE
      DO 220 L=MSG+1,ILEVM
      DO 220 IL=IL1,IL2 
         IF ( L.LE.MAXI(IL) ) THEN
             TV =T(IL,L)*(1.+TVFA*Q(IL,L))
             TPV= TC(IL,L)*(1.+TVFA*QC(IL,L)
     1           -(1.+TVFA)*MAX(QC(IL,L)-QSTC(IL,L),ZERO))
             BUOY(IL,L)=(TPV-TV)/TV
             FACT(IL,L)=RGAS*BUOY(IL,L)*T(IL,L)
     1                 *LOG(PF(IL,L+1)/PF(IL,L))
             BUOY(IL,L)=BUOY(IL,L)*GRAV
         ENDIF
 220  CONTINUE
C
C    * DETERMINE CLOUD BASE AS LIFTING CONDENSATION LEVEL (LCL).
C
      DO 240 IL=IL1,IL2
         LCL(IL)=MSG+1
         LTM(IL)=0
 240  CONTINUE
      DO 250 L=ILEVM,MSG+1,-1
      DO 250 IL=IL1,IL2
         IF (     L.LT.MAXI(IL) 
     1      .AND. QC(IL,L  ).GE.QSTC(IL,L  )   
     3      .AND. LTM(IL).EQ.0      ) THEN
            LTM(IL)=1
            LCL(IL)=L
         ENDIF
 250  CONTINUE
C
C     * DETERMINE CLOUD TOP AS LEVEL OF NEUTRAL BUOYANCY (LNB).
C
C
C     * VERTICAL VELOCITY OF PARCEL.  
C
      ADEF=-9.
      DO 295 L=MSG+1,ILEV                                                       
      DO 295 IL=IL1,IL2
         WC(IL,L)=ADEF
 295  CONTINUE
      DO 300 IL=IL1,IL2
         L=LCL(IL)
         TL(IL)=TC(IL,L)
         WC(IL,L)=0.4
 300  CONTINUE
      DO 310 L=ILEVM,MSG+1,-1 
      DO 310 IL=IL1,IL2 
         IF ( L.LT.LCL(IL) .AND. WC(IL,L+1).GE. 0.            ) THEN
            DZ=Z(IL,L)-Z(IL,L+1)
            DXDG=( BUOY(IL,L)-BUOY(IL,L+1) )/DZ
            ATMP=(BUOY(IL,L+1)-DXDG*Z(IL,L+1))*DZ
     1          +.5*DXDG*(Z(IL,L)**2-Z(IL,L+1)**2)
            ATMP=WC(IL,L+1)**2+2.*ATMP
            WC(IL,L)=SIGN(SQRT(ABS(ATMP)),ATMP)
         ENDIF
 310  CONTINUE
C
C     * ASSIGN TOP OF CLOUD LAYER TO LEVEL OF NEUTRAL BUOYANCY (LNB).
C     * THE LNB IS DEFINED AS THE FIRST LEVEL AT WHICH THE CLOUDY AIR
C     * IS NEUTRALLY BUOYANT BELOW THAT LEVEL AT WHICH THE MEAN VERTICAL
C     * VELOCITY IN THE CLOUD DROPS TO ZERO. 
C
      LDEF=ILEV
      DO 500 IL=IL1,IL2
         LNB(IL)=LDEF
         LTM(IL)=LDEF
 500  CONTINUE
      DO 520 L=ILEVM,MSG+1,-1
      DO 520 IL=IL1,IL2 
         IF ( L.LT.LCL(IL) .AND. LNB(IL).EQ.LDEF ) THEN
            IF ( WC(IL,L).LE.0. .AND.WC(IL,L+1).GT.0. ) THEN
               IF ( LTM(IL).NE.LDEF ) THEN
                  LNB(IL)=LTM(IL)
               ELSE
                  LNB(IL)=L
               ENDIF
            ELSE IF ( BUOY(IL,L).LE.0. .AND. BUOY(IL,L+1).GT.0. ) THEN
               LTM(IL)=L
            ENDIF
         ENDIF
 520  CONTINUE
      DO 540 IL=IL1,IL2 
         IF ( LNB(IL).EQ.LDEF ) THEN
            IF ( LTM(IL).NE.LDEF ) THEN 
               LNB(IL)=LTM(IL)
            ELSE
               LNB(IL)=ILEV+1
            ENDIF
         ENDIF
 540  CONTINUE
C                                                                            
C     * CALCULATE CONVECTIVE AVAILABLE POTENTIAL ENERGY (CAPE).             
C                    
      DO 600 L=MSG+1,ILEV                                                   
      DO 600 IL=IL1,IL2                                                       
        IF(IDOC(IL).NE.0 .AND. L.LT.MAXI(IL) .AND. L.GE.LNB(IL)) THEN         
           CAPE(IL)=CAPE(IL)+FACT(IL,L)
           IF ( L.LE.LCL(IL) ) CAPEP(IL)=CAPEP(IL)+FACT(IL,L)
        ENDIF                                                                
  600 CONTINUE
C
      RETURN
      END
