      SUBROUTINE PARPRP4(BUOY,TC,QSTC,HMNC,HMNM,QC,QM,QF,ILEV,ILEVM,
     1                   ILG,ILGG,MSG,LMAG,LCL,FCB,TDRY,ZF,TSTAR,
     2                   WC,DZ,PF,TFE,ISKIP,RRL,GRAV,CPRES,EPS1,EPS2,
     3                   A,B,RGAS,RGASV,TVFA,ADEF,LTM,MB,PG,FACT)
C-------------------------------------------------------------------------
C     * FEB 02/2017 - M.LAZARE.    REMOVE (UNUSED) CAPEN.
C     * MAY 01/2015 - K.VONSALZEN: NEW VERSION FOR GCM18+:
C     *                            - Loops 250/260 modified to handle
C     *                              pathalogical case where LCL is
C     *                              exactly at the same height as a
C     *                              model level base (just happened 
C     *                              for the first time in almost 10 years).
C     *                              The modifications give bit-for-bit
C     *                              otherwise.
C     * NOV 25/2006 - M.LAZARE.    PREVIOUS VERSION PARPRP3 FOR GCM15F->GCM17.
C                                  - CALLS NEW TINCLD3 INSTEAD OF TINCLD2.
C     *                            - SELECTIVE PROMOTION OF VARIABLES TO
C     *                              REAL*8 TO SUPPORT 32-BIT MODE.
C     * JUN 14/2006 - M.LAZARE.    NEW VERSION FOR GCM15F:
C     *                            - "ZERO","P999999" DATA CONSTANTS
C     *                              ADDED AND USED IN CALL TO 
C     *                              INTRINSICS. 
C     * MAY 04/2006 - M.LAZARE/    PREVIOUS VERSION PARPRP2 FOR GCM15E:
C     *               K.VONSALZEN. - CALLS NEW TINCLD2.
C     *                              THIS REQUIRES PASSING IN EPS2.
C     * MAY 05/2001 - K.VONSALZEN/ PREVIOUS VERSION PARPRP FOR GCM15D.
C     *               N.MCFARLANE.
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)
C     
C     * WORK FIELDS (FIRST SET ARE PROMOTED FOR CONSISTENCY
C     * WITH REST OF SHALLOW). 
C
      REAL*8,  DIMENSION(ILG,ILEV)         :: BUOY,QF,PF,QC,QSTC,TC,
     1                                        DZ,ZF,TDRY,HMNC,WC,
     2                                        TSTAR
      REAL*8,  DIMENSION(ILG)              :: FCB
      REAL*8                                  ZERO,P999999
C
      REAL,    DIMENSION(ILG,ILEV)         :: TFE,PG,FACT
      REAL,    DIMENSION(ILG)              :: HMNM,QM,MB
      INTEGER LMAG(ILG), LCL(ILG), ISKIP(ILG), LTM(ILG)
C
      DATA ZERO    / 0.     /
      DATA P999999 /0.999999/
C-----------------------------------------------------------------------
C
C    * CALCULATE TEMPERATURE OF UNDILUTED AIR PARCEL.
C
      DO 100 IL=1,ILGG 
         ISKIP(IL)=0
 100  CONTINUE
      DO 140 L=MSG+2,ILEV 
      DO 140 IL=1,ILGG 
         QC  (IL,L)=QM(IL)
         HMNC(IL,L)=HMNM(IL)
 140  CONTINUE
      CALL TINCLD3(TC,ILG,ILEV,MSG,1,ILGG,LMAG,TDRY,HMNC,QC,QSTC,ZF,
     1             PF,RRL,GRAV,CPRES,EPS1,EPS2,A,B,ISKIP)
C
C    * BUOYANCY OF UNDIULTED AIR PARCEL.
C
      DO 200 L=MSG+1,ILEV
      DO 200 IL=1,ILGG 
         BUOY(IL,L)=-1.E-10
 200  CONTINUE
      DO 210 L=MSG+2,ILEV
      DO 210 IL=1,ILGG 
         IF ( L.LE.LMAG(IL) ) THEN
             TV =TFE(IL,L)*(1.+TVFA*QF(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)=GRAV*( TPV-TV )/TV
             FACT(IL,L)=RGAS*BUOY(IL,L)*TFE(IL,L)
     1                 *LOG(PG(IL,L)/PG(IL,L-1))/GRAV
         ENDIF
 210  CONTINUE
C
C    * DETERMINE CLOUD BASE AS LIFTING CONDENSATION LEVEL (LCL).
C
      DO 240 IL=1,ILGG
         LCL(IL)=MSG+2
         LTM(IL)=0
         FCB(IL)=0.
 240  CONTINUE
      DO 250 L=ILEVM,MSG+2,-1
      DO 250 IL=1,ILGG
         IF (     L.LT.LMAG(IL) .AND. LTM(IL).EQ.0      ) THEN
            IF (       QC(IL,L  ).GT.QSTC(IL,L  )  
     1           .AND. QC(IL,L+1).LE.QSTC(IL,L+1)       ) THEN
               LTM(IL)=1
               LCL(IL)=L
               FCB(IL)=(QC(IL,L  )-QSTC(IL,L  ))
     1              /( (QC(IL,L  )-QSTC(IL,L  ))
     2                -(QC(IL,L+1)-QSTC(IL,L+1)) )
               FCB(IL)=MAX(MIN(FCB(IL),P999999),ZERO)
            ELSE IF (   QC(IL,L  ).EQ.QSTC(IL,L  )      ) THEN
               LTM(IL)=1
               LCL(IL)=L+1
               FCB(IL)=1.
            ENDIF
         ENDIF
 250  CONTINUE
      DO 260 IL=1,ILGG
         IF ( LTM(IL).EQ.0 .OR. LCL(IL) .GE. ILEV ) THEN
            ISKIP(IL)=1
         ENDIF
 260  CONTINUE
      DO 270 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LCL(IL)
            BUOY(IL,L+1)=FCB(IL)*BUOY(IL,L+1)+(1.-FCB(IL))*BUOY(IL,L)
            ZF  (IL,L+1)=FCB(IL)*ZF  (IL,L+1)+(1.-FCB(IL))*ZF  (IL,L)
            DZ  (IL,L  )=ZF(IL,L)-ZF(IL,L+1)
         ENDIF
 270  CONTINUE
C
C     * VERTICAL VELOCITY OF PARCEL.  
C
      DO 295 L=MSG+1,ILEV                                                       
      DO 295 IL=1,ILGG
         WC(IL,L)=ADEF
 295  CONTINUE
      DO 300 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LCL(IL)
C
C           *** SET CLOUD-BASE VELOCITY TO 0.25 M/S.
C
            WC(IL,L+1)=0.25
         ENDIF
 300  CONTINUE
      DO 310 L=ILEVM,MSG+2,-1 
      DO 310 IL=1,ILGG 
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0 
     1       .AND. WC(IL,L+1).GE. 0.            ) THEN
            DXDG=( BUOY(IL,L)-BUOY(IL,L+1) )/DZ(IL,L)
            ATMP=(BUOY(IL,L+1)-DXDG*ZF(IL,L+1))*DZ(IL,L)
     1          +.5*DXDG*(ZF(IL,L)**2-ZF(IL,L+1)**2)
            ATMP=WC(IL,L+1)**2+2.*ATMP
            WC(IL,L)=SIGN(SQRT(ABS(ATMP)),ATMP)
         ENDIF
 310  CONTINUE
      DO 330 IL=1,ILGG 
         LTM(IL)=MSG+2
 330  CONTINUE
      DO 350 L=ILEVM,MSG+2,-1 
      DO 350 IL=1,ILGG 
         IF (       L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0 
     1        .AND. WC(IL,L).NE.ADEF) THEN
            LTM(IL)=L
         ENDIF
 350  CONTINUE
C
C     * TIME AT WHICH PARCEL REACHES A CERTAIN LEVEL.
C
      DO 390 L=MSG+1,ILEV
      DO 390 IL=1,ILGG
         TSTAR(IL,L)=ADEF
 390  CONTINUE
      DO 400 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LCL(IL)+1
            TSTAR(IL,L)=0.
         END IF
 400  CONTINUE
      DO 410 L=ILEVM,MSG+2,-1 
      DO 410 IL=1,ILGG 
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0
     1       .AND. L.GE.(LTM(IL)+1) ) THEN
            TSTAR(IL,L)=TSTAR(IL,L+1)+2.*DZ(IL,L)
     1                               /(WC(IL,L)+WC(IL,L+1))
         ENDIF
 410  CONTINUE
C
      RETURN
      END
