          SUBROUTINE KFPARA(NCUYES,ICUYES,J,LSB,HTOP,HBOT
     1,                     CNVTOP,CNVBOT,NSHALL)
C     ******************************************************************
C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C                .      .    .
C SUBPROGRAM:    KFPARA      CONVECTIVE PRECIPITATION PARAMETERIZATION
C   PRGRMMR: KAIN            ORG: W/NP2      DATE: 00-04-13
C
C ABSTRACT:
C     KFPARA IS THE CENTRAL SUBROUTINE IN THE K-F PARAMETERIZATION.
C     SHALLOW CONVECTION IS INCLUDED WITHOUT CAPE DEPENDENCE.
C
C
C PROGRAM HISTORY LOG:
C   ??-??-??  KAIN       - ORIGINATOR
C   00-04-13  BLACK      - INCORPORATED INTO THE ETA MODEL
C
C USAGE: CALL KFPARA FROM SUBROUTINE KFDRIVE
C
C   INPUT ARGUMENT LIST:
C     NCUYES - NUMBER OF GRID POINTS TO CHECK ON THIS J-SLICE
C     ICUYES - ARRAY CONTAINING THE I-VALUES FOR GRID POINTS TO CHECK
C     J      - THE ROW NUMBER
C     LSB    - ARRAY CONTAINING l VALUES FOR EACH NCUYES POINT AT WHICH
C              CHECKS FOR CONVECTIVE INITAITION SHOULD BEGIN
C
C   OUTPUT ARGUMENT LIST:
C     HTOP   - ARRAY USED TO STORE L VALUE FOR CLOUD TOP (FOR RADTN)
C     HBOT   - ARRAY USED TO STORE L VALUE FOR CLOUD BOTTOM (FOR RADTN)
C     CNVTOP - ARRAY USED TO STORE L VALUE FOR CLOUD TOP (FOR CHKOUT)
C     CNVBOT - ARRAY USED TO STORE L VALUE FOR CLOUD BOTTOM (FOR CHKOUT)
C     NSHALL - COUNTER FOR NUMBER OF SHALLOW CONVECTION POINTS ACTIVATED
C
C   OUTPUT FILES:
C     NONE
C
C   SUBPROGRAMS CALLED:
C
C     UNIQUE:
C        CONLOAD
C        DTFRZNEW
C        ENVIRTHT
C        PROF5
C
C     LIBRARY:
C        NONE
C
C   COMMON BLOCKS: CTLBLK
C                  LOOPS
C                  MASKS
C                  VRBLS
C                  DYNAM
C                  CNVCLD
C                  PVRBLS
C                  ACMCLH
C                  KFFDBK
C                  KFLUT
C                  INDX
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 90
C   MACHINE : IBM SP
C$$$
C**********************************************************************C
C                                                                      C
C ETA INPUT:   TEMPERATURE (T, K);    SPECIFIC HUMIDITY (Q, KG/KG)     C
C              PD, RES, AND PT USED TO CALCULATE PRESSURE (PASCAL)     C
C              HORIZONTAL WIND SPEED (U AND V,  M/S)                   C
C              HORIZONTAL GRID SPACING (DX,  M)                        C
C              MODEL TIME STEP (DT,  SECONDS)                          C
C              NUMBER OF TIME STEPS INTEGRATED (NTSD,  NO UNITS)       C
C              INSTANTANEOUS SURFACE SENSIBLE HEAT FLUX (TWBS, W/M^2)  C
C              INSTANTANEOUS SURFACE LATENT HEAT FLUX (QWBS, W/M^2)    C
C              ETA MODEL LEVELS (ETA,  NO UNITS)                       C
C              HALF ETA MODEL LEVELS (AETA, NO UNITS)                  C
C              MODEL LEVEL OF TOP OF PBL (PBLLEV)                      C
C              DIMENSIONS OF MODEL GRID (IMX,JMX,KMX FOR IM,JM,LM)     C
C              ARRAYS SPECIFYING WHERE H AND V ARRAYS ARE BELOW        C
C                                GROUND (LMH AND VMH,  NO UNITS)       C
C                                                                      C
C    OUTPUT:   CONVECTIVE TENDENCIES OF TEMPERATURE (DTDT), WATER      C
C              VAPOR (DQDT), CLOUD WATER (DQCDT), RAIN WATER (DQRDT)   C
C   !!!!! NOTE:  cloud water and rainwater arrays (DQCDT and DQRDT)    C
C                arrays removed for implementation in ETA model !!!!!!!C
C              AND PRECIPITATION RATE (RAINCV)                         C
C              HORIZ. LOCATION OF ACTIVE CONVECTION (NCA, NO UNITS)    C
C                                                                      C
C**********************************************************************C
C                                                                      C
C     REFERENCES:                                                      C
C                                                                      C
C         KAIN AND FRITSCH (1993):  "CONVECTIVE PARAMETERIZATION IN    C
C         MESOSCALE MODELS:  THE KAIN-FRITSCH SCHEME" IN THE REPRESEN- C
C         TATION OF CUMULUS CONVECTION IN NUMERICAL MODELS, A.M.S.     C
C         MONOGRAPH, K.A. EMANUEL AND D.J. RAYMOND, EDS., 165-170.     C
C                                                                      C
C         FRITSCH AND KAIN (1993):  "CONVECTIVE PARAMETERIZATION IN    C
C         MESOSCALE MODELS:  THE FRITSCH-CHAPPELL SCHEME" IN THE REP-  C
C         RESENTATION OF CUMULUS CONVECTION IN NUMERICAL MODELS, A.M.S.C
C         MONOGRAPH, K.A. EMANUEL AND D.J. RAYMOND, EDS., 165-170.     C
C                                                                      C
C         STENSRUD AND FRITSCH (1994), MON. WEA. REV., 2084-2104.      C
C                                                                      C
C         FRITSCH AND CHAPPELL (1980), J. ATMOS. SCI., 1722-1733.      C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C----------------------------------------------------------------------
      INCLUDE "parmeta"
      INCLUDE "mpp.h"
C----------------------------------------------------------------------
      PARAMETER(LP1=LM+1,JAM=6+2*(JM-10))
      PARAMETER(KX=LM,KXP1=KX+1,ILX=IM-1,JLX=JM-1)
C----------------------------------------------------------------------
      INCLUDE "ACMCLH.comm"
C----------------------------------------------------------------------
      INCLUDE "CTLBLK.comm"
C----------------------------------------------------------------------
      INCLUDE "DYNAM.comm"
C----------------------------------------------------------------------
      INCLUDE "MASKS.comm"
C----------------------------------------------------------------------
      INCLUDE "VRBLS.comm"
C----------------------------------------------------------------------
      INCLUDE "PVRBLS.comm"
C----------------------------------------------------------------------
      INCLUDE "LOOPS.comm"
C----------------------------------------------------------------------
      INCLUDE "CNVCLD.comm"
C----------------------------------------------------------------------
      INCLUDE "KFFDBK.comm"
C----------------------------------------------------------------------
      INCLUDE "KFLUT.comm"
C----------------------------------------------------------------------
      INCLUDE "INDX.comm"
C----------------------------------------------------------------------
                             R E A L
     1 HTOP  (IDIM1:IDIM2,JDIM1:JDIM2),HBOT  (IDIM1:IDIM2,JDIM1:JDIM2)
     2,CNVTOP(IDIM1:IDIM2,JDIM1:JDIM2),CNVBOT(IDIM1:IDIM2,JDIM1:JDIM2)
     3,QUER(KX),TKE(KX) 
     4,CLDHGT(KX),QSD(KX),DILFRC(KX),DDILFRC(KX)
     5,THTEEG(KX),TGU(KX),QGU(KX)
C
                             I N T E G E R
     1 ICUYES(IDIM1:IDIM2),LSB(IDIM1:IDIM2)
C----------------------------------------------------------------------
C
C***  DEFINE LOCAL VARIABLES
C     
                             R E A L
     1 P0(KXP1),Z00(KX),T0(KX),TV0(KX),Q0(KXP1),U0(KX),V0(KX)
     2,TU(KX),TVU(KX),QU(KX),TZ(KX),TVD(KX),QD(KX),QES(KX),THTES(KX)
     3,TG(KX),TVG(KX),QG(KX),WU(KX),WD(KX),EMS(KX),EMSD(KX),W0(KX)
     4,UMF(KX) ,UER(KX) ,UDR(KX) ,DMF(KX) ,DER(KX) ,DDR(KX) ,DZQ(KX)
     5,UMF2(KX),UER2(KX),UDR2(KX),DMF2(KX),DER2(KX),DDR2(KX),DZA(KX)
     6,THTA0(KXP1),THETEE(KX),THTAU(KX),THETEU(KX),THTAD(KX)
     7,THETED(KX)
     8,QLIQ(KX),QICE(KX),QLQOUT(KX),QICOUT(KX),PPTLIQ(KX),PPTICE(KX)
     9,DETLQ(KX),DETIC(KX),DETLQ2(KX),DETIC2(KX),RATIO(KX),RATIO2(KX)
C
                             R E A L
     1 DOMGDP(KX),EXN(KX),RHOE(KX),TVQU(KX),DP(KX),RH(KX)
     2,EQFRC(KX),WSPD(KX),QDT(KX),FXM(KX),THTAG(KX),THTESG(KX)
     3,THPA(KX),THFXIN(KX),THFXOUT(KX),QPA(KX),QFXIN(KX),QFXOUT(KX)
     4,QLPA(KX),QLFXIN(KX),QLFXOUT(KX),QIPA(KX),QIFXIN(KX),QIFXOUT(KX)
     5,QRPA(KX),QRFXIN(KX),QRFXOUT(KX),QSPA(KX),QSFXIN(KX),QSFXOUT(KX)
     6,QL0(KX),QLG(KX),QI0(KX),QIG(KX),QR0(KX),QRG(KX),QS0(KX),QSG(KX)
C
                             R E A L
     1 OMG(KXP1),RAINFB(KX),SNOWFB(KX)
C----------------------------------------------------------------------     
C
C...DEFINE CONSTANTS NEEDED FOR KFPARA SUBROUTINE
C
      DATA P00,T00,G,CP/1.E5,273.16,9.8,1004.6/
      DATA RLF/3.339E5/
      DATA RHIC,RHBC/1.,0.90/
      DATA PIE,TTFRZ,TBFRZ,C5/3.141592654,268.16,248.16,1.0723E-3/
      DATA RATE,RV/0.03,461.5/
      DATA ALIQ,BLIQ,CLIQ,DLIQ/613.3,17.502,4780.8,32.19/
      DATA AICE,BICE,CICE,DICE/613.2,22.452,6133.0,0.61/
      DATA XLV0,XLV1,XLS0,XLS1/3.147E6,2369.,2.905E6,259.532/
C----------------------------------------------------------------------     
C****************************************************************************
C----------------------------------------------------------------------     
C    
C***  OPTION TO FEED CONVECTIVELY GENERATED RAINWATER  
C***  INTO GRID-RESOLVED RAINWATER (OR SNOW/GRAUPEL)  
C***  FIELD.  'FBFRC' IS THE FRACTION OF AVAILABLE   
C***  PRECIPITATION TO BE FED BACK (0.0 - 1.0)
C
      FBFRC=0.0  
C
C----------------------------------------------------------------------     
      ROVG=R/G
      GDRY=-G/CP
      DT2 =2.*DT
C----------------------------------------------------------------------     
C
C***  SPECIFY DOWNDRAFT MASS FLUX AS FRACTION
C***  OF UPDRAFT MASS FLUX (DMFFRC)
C
      DMFFRC=0.9
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
C***
C***  LOOP OVER ALL POTENTIAL CONVECTIVE GRID POINTS
C***
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
      NCCNT=0
C
      DO 900 NC=1,NCUYES
C
C***  VARIABLES TO ALLOW SHALLOW CONVECTION
C
      NCHM=0
      ISHALL=0
      RAD=1500.
      FBFRC=0.  
      DPMIN=5.E3
C
      I=ICUYES(NC)
      XTIME=NTSD*DT/60.
      DXSQ=DX(I,J)*DX(I,J)
      P200=PD(I,J)+PT-2.E4
      P300=PD(I,J)+PT-3.E4
      P400=PD(I,J)+PT-4.E4
      ML=0
C
C***  DEFINE NUMBER OF LAYERS ABOVE GROUND LEVEL - CALL THIS KL
C
      KL=LMH(I,J)
      NLEV=KX-KL
      KLM=KL-1
C
C***  INPUT A VERTICAL SOUNDING
C
C***  THE FOLLOWING LOOP IS CRUCIAL.  THE KF SCHEME IS SET UP TO HAVE 
C***  LEVELS STARTING WITH 1 AT THE LOWEST MODEL LEVEL AND GOING UP, 
C***  WHICH IS INVERTED FROM WHAT SIGMA AND ETA COORDINATE MODELS 
C***  TEND TO USE.  THUS, THE FIRST STEP IS TO SWITCH THE ORDER OF THE
C***  VARIABLES USED WITHIN THE KF SCHEME.  NOTE THAT THE SCHEME NEEDS
C***  THE FOLLOWING VARIABLES
C
C     INPUT:   TEMPERTURE (T0, K) ;    SPECIFIC HUMIDITY (Q0, KG/KG) ; 
C              HORIZONTAL WIND SPEED (U0 AND V0, M/S) ;                
C              PRESSURE (P0, PASCAL) ;  HEIGHT (Z0, M);                
C              VERTICAL MOTION (W0, M/S).                              
C
      DO 15 K=1,KL
      NK=KX-K+1-NLEV
      P0(K)=PD(I,J)*RES(I,J)*AETA(NK)+PT
      T0(K)=T(I,J,NK)
      Q0(K)=Q(I,J,NK)
      TKE(K)=Q2(I,J,NK)
C
C***  SATURATION VAPOR PRESSURE (ES) IS CALCULATED USING FORMULA GIVEN BY
C***  BUCK (1981).  IF Q0 IS ABOVE SATURATION VALUE USING THIS METHOD,
C***  REDUCE IT TO SATURATION LEVEL.
C
      ES=ALIQ*EXP((BLIQ*T0(K)-CLIQ)/(T0(K)-DLIQ))
      QES(K)=0.622*ES/(P0(K)-ES)
      Q0(K)=AMIN1(QES(K),Q0(K))
C
      IF(Q0(K).LE.0.1E-8)THEN
        Q0(K)=0.1E-8
      ENDIF
C
      RH(K)=Q0(K)/QES(K)  
C
C***  SET HYDROMETEOR CONCENTRATIONS TO ZERO INITIALLY.  ALTHOUGH SOME
C***  HYDROMETEORS MAY ACTUALLY BE PRESENT, THIS IS DONE TO DISALLOW
C***  ENTRAINMENT OF THESE INTO CONVECTIVE UPDRAFTS AND DOWNDRAFTS.
C***  NECESSARY BECAUSE CONVECTIVE SCHEME OPERATES OVER MULTIPLE TIME
C***  STEPS ASSUMING STEADY-STATE ENVIRONMENTAL CONDITIONS AND THERE
C***  IS NO WAY OF GUARANTEEING THAT LIQUID WATER WILL CONTINUE TO
C***  REMAIN AVAILABLE FOR ENTRAINMENT EVEN IF IT IS THERE INITIALLY.
C
      QL0(K)=0.
      QI0(K)=0.
      QR0(K)=0.
      QS0(K)=0.
      DILFRC(K)=1.
C
C***  CALCULATE WIND AT H-POINT AS AVERAGE OF 4 SURROUNDING V-POINT VALUES
C
      SUMV=VTM(I+IHE(J),J,NK)+VTM(I+IHW(J),J,NK)+VTM(I,J+1,NK)+
     1     VTM(I,J-1,NK)

Cmp - correct for possible HTM/VTM discrepancies
	IF (SUMV .gt. 0) THEN

      U0(K)=(U(I+IHE(J),J,NK)+U(I+IHW(J),J,NK)+U(I,J+1,NK)+
     1       U(I,J-1,NK))/SUMV
      V0(K)=(V(I+IHE(J),J,NK)+V(I+IHW(J),J,NK)+V(I,J+1,NK)+
     1       V(I,J-1,NK))/SUMV

	ELSE

      U0(K)=0.
      V0(K)=0.

	ENDIF
C
      TV0(K)=T0(K)*(1.+0.608*Q0(K))
      RHOE(K)=P0(K)/(R*TV0(K))
C
      DZQ(K)=ROVG*TV0(K)*ALOG((PD(I,J)*RES(I,J)*ETA(NK+1)+PT)/
     1      (PD(I,J)*RES(I,J)*ETA(NK)+PT))
      DP(K)=(ETA(NK+1)-ETA(NK))*PD(I,J)*RES(I,J)
      W0(K)=W0AVG(I,J,NK)
C
C***  DZQ IS DZ BETWEEN ETA SURFACES. DZA IS DZ BETWEEN MODEL HALF LEVEL.
C***  DP IS THE PRESSURE INTERVAL BETWEEN FULL ETA LEVELS.
C
      IF(P0(K).GE.500E2)L5=K
      IF(P0(K).GE.400E2)L4=K
      IF(P0(K).GE.P300)LLFC=K 
      IF(T0(K).GT.T00)ML=K
      CLDHGT(K)=0.
C
   15 CONTINUE
C----------------------------------------------------------------------     
C
C***  FILL REMAINING LAYERS WITH ZEROES
C
      DO K=KL+1,KX
        P0(K)=0.
        T0(K)=0.
        Q0(K)=0.
        QES(K)=0.
        U0(K)=0.
        V0(K)=0.
        QL0(K)=0.
        QI0(K)=0.
        QR0(K)=0.
        QS0(K)=0.
        TV0(K)=0.
        RHOE(K)=0.
        DZQ(K)=0.
        DP(K)=0.
        Z00(K)=0.
        W0(K)=0.
        DZA(K)=0.
        DILFRC(K)=0.
      ENDDO
C
      Z00(1)=.5*DZQ(1)
C
      DO K=2,KL
        Z00(K)=Z00(K-1)+0.5*(DZQ(K)+DZQ(K-1))
        DZA(K-1)=Z00(K)-Z00(K-1)
      ENDDO
C
      DZA(KL)=0.
      NLOOP=0.
C
      KMIX=LSB(I)
C----------------------------------------------------------------------     
   25 LOW=KMIX
      NLOOP=NLOOP+1
C
      IF(NLOOP.GT.50)THEN
        PRINT*,'I, J, LOW, ISHALL, NCHM =',I,J,LOW,ISHALL,NCHM
        IF(NLOOP.GT.100)THEN
          GO TO 600
          STOP 'NLOOP'
        ENDIF
      ENDIF
C
C***  IF PARCEL ORIGINATES FROM 300 MB ABOVE SURFACE OR HIGHER THEN SCHEME
C***  IS NOT USED - CHECK NEXT GRID POINT
C
C***  VARIABLES FOR SHALLOW CONVECTION
C      
      IF(LOW.GT.LLFC)THEN 
        IF(ISHALL.EQ.1)THEN
          CHMAX=0.
          NCHM=0
C
          DO NK=1,LLFC
            IF(CLDHGT(NK).GT.CHMAX)THEN
              NCHM=NK
              CHMAX=CLDHGT(NK)
            ENDIF
          ENDDO
C
          KMIX=NCHM
          GO TO 25
        ENDIF
        GO TO 900   
      ENDIF      
C----------------------------------------------------------------------     
C***
      LC=LOW
      IF(ISHALL.EQ.1.AND.LC.EQ.NCHM)THEN
        FBFRC=1.
C
C***  COMMENT OUT STATEMENT BELOW SO THAT SHALLOW CLOUD DRAWS
C***  FROM SAME LAYER AS DEEP CLOUD TO PREVENT SHALLOW CLOUD
C***  FROM BECOMING TOO DEEP!
C
C       DPMIN=2.5E3
      ENDIF
C----------------------------------------------------------------------     
C
C***  ASSUME THAT IN ORDER TO SUPPORT A DEEP UPDRAFT YOU NEED A LAYER OF
C***  UNSTABLE AIR 50 TO 100 MB DEEP.  TO APPROXIMATE THIS, ISOLATE A
C***  GROUP OF ADJACENT INDIVIDUAL MODEL LAYERS, WITH THE BASE AT LEVEL
C***  LC, SUCH THAT THE COMBINED DEPTH OF THESE LAYERS IS AT LEAST 60 MB.
C
C***  INTRODUCE A MOISTURE PERTURBATION IN UPDRAFT SOURCE LAYERS
C***  IF RH > 75%.  NOTE THAT THIS APPROACH SHOULD PERHAPS BE
C***  MODIFIED IF THE GRID-SCALE CONDENSATION THRESHOLD IS CHANGED.
C
      NLAYRS=0
      DPTHMX=0.
C
      DO NK=LC,KX
        DPTHMX=DPTHMX+DP(NK)
        NLAYRS=NLAYRS+1
        QEF=0.
        QUER(NK)=Q0(NK)*(1.+QEF)  
        IF(DPTHMX.GT.DPMIN)GO TO 50
      ENDDO
C
      GO TO 900
   50 KPBL=LC+NLAYRS-1   
C
C***  DETERMINE WHAT LEVEL TO START WITH FOR THE NEXT
C***  MIXTURE IN CASE THE CURRENT MIXTURE, WITH BASE AT
C***  LEVEL LC, IS NOT BUOYANT.
C***  INSTEAD OF CHECKING MIXTURES USING EVERY SINGLE LAYER,
C***  MOVE UP IN INCREMENTS OF AT LEAST 20 MB.
C
      PM15=P0(LC)-15.E2
C
      DO NK=LC+1,KL
        IF(P0(NK).LT.PM15)THEN
          KMIX=NK
          GO TO 75
        ENDIF
      ENDDO
C
      GO TO 900
   75 CONTINUE      
C----------------------------------------------------------------------     
C***
C***  FOR COMPUTATIONAL SIMPLICITY WITHOUT MUCH LOSS IN ACCURACY,
C***  MIX TEMPERATURE INSTEAD OF THETA FOR EVALUATING CONVECTIVE
C***  INITIATION (TRIGGERING) POTENTIAL.
C
      TMIX=0.
      QMIX=0.
      ZMIX=0.
      PMIX=0.
C
C***  FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
C***  MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
C***  LAYERS.
C
      DO NK=LC,KPBL
        TMIX=TMIX+DP(NK)*T0(NK)
        QMIX=QMIX+DP(NK)*QUER(NK)  
        ZMIX=ZMIX+DP(NK)*Z00(NK)
        PMIX=PMIX+DP(NK)*P0(NK)
      ENDDO
C
      TMIX=TMIX/DPTHMX
      QMIX=QMIX/DPTHMX
      ZMIX=ZMIX/DPTHMX
      PMIX=PMIX/DPTHMX
      EMIX=QMIX*PMIX/(0.622+QMIX)
C
C***  FIND THE TEMPERATURE OF THE MIXTURE AT ITS LCL
C
      ASTRT=1.E-3
      AINC=0.075
      A1=EMIX/ALIQ
      TP=(A1-ASTRT)/AINC
      INDLU=INT(TP)+1
      VALUE=(INDLU-1)*AINC+ASTRT
      AINTRP=(A1-VALUE)/AINC
      TLOG=AINTRP*ALU(INDLU+1)+(1-AINTRP)*ALU(INDLU)
      TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
      TLCL=TDPT-(0.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))
     1         *(TMIX-TDPT)
      TLCL=AMIN1(TLCL,TMIX)
      TVLCL=TLCL*(1.+0.608*QMIX)
Celim      TVQU(K)=TVLCL
      ZLCL=ZMIX+(TLCL-TMIX)/GDRY
C
      DO NK=LC,KL
        KLCL=NK
        IF(ZLCL.LE.Z00(NK))GO TO 100
      ENDDO 
C
      GO TO 900
  100 K=KLCL-1
C----------------------------------------------------------------------     
C
C***  CALCULATE DLP USING Z INSTEAD OF LOG(P)
C
      DLP=(ZLCL-Z00(K))/(Z00(KLCL)-Z00(K))
C     
C***  ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL.
C     
      TENV=T0(K)+(T0(KLCL)-T0(K))*DLP
      QENV=Q0(K)+(Q0(KLCL)-Q0(K))*DLP
      TVEN=TENV*(1.+0.608*QENV)
C     
C***  CHECK TO SEE IF CLOUD IS BUOYANT USING FRITSCH-CHAPPELL TRIGGER
C***  FUNCTION DESCRIBED IN KAIN AND FRITSCH (1992).  W0 IS AN
C***  APROXIMATE VALUE FOR THE RUNNING-MEAN GRID-SCALE VERTICAL
C***  VELOCITY, WHICH GIVES SMOOTHER FIELDS OF CONVECTIVE INITIATION
C***  THAN THE INSTANTANEOUS VALUE.  FORMULA RELATING TEMPERATURE
C***  PERTURBATION TO VERTICAL VELOCITY HAS BEEN USED WITH THE MOST
C***  SUCCESS AT GRID LENGTHS NEAR 25 KM.  FOR DIFFERENT GRID-LENGTHS,
C***  ADJUST VERTICAL VELOCITY TO EQUIVALENT VALUE FOR 25 KM GRID
C***  LENGTH, ASSUMING LINEAR DEPENDENCE OF W ON GRID LENGTH.
C     
      IF(ZLCL.LT.2.E3)THEN
        WKLCL=0.02*ZLCL/2.E3
      ELSE
        WKLCL=0.02
      ENDIF
C
      WKL=(W0(K)+(W0(KLCL)-W0(K))*DLP)*DX(I,J)
     1      /25.E3-WKLCL
      WABS=ABS(WKL)
C
      IF(WABS.EQ.0.)THEN
        WSIGNE=1.
        DTLCL=0.
        GO TO 120
      ENDIF
C
      WSIGNE=WKL/WABS
      DTLCL=4.64*WSIGNE*WABS**0.33
      DTLCL=AMAX1(DTLCL,0.)
  120 CONTINUE
C----------------------------------------------------------------------     
C
C***  GIVE PARCEL AN EXTRA TEMPERATURE PERTURBATION BASED
C***  THE THRESHOLD RH FOR CONDENSATION.
C***  FOR NOW, JUST ASSUME U00=0.75.
C
      U00=0.75
C
      IF(U00.LT.1.)THEN
        QSLCL=QES(K)+(QES(KLCL)-QES(K))*DLP
        RHLCL=QENV/QSLCL
        DQSDT=QMIX*(CLIQ-BLIQ*DLIQ)/((TLCL-DLIQ)*(TLCL-DLIQ))
        IF(RHLCL.GE.0.75.AND.RHLCL.LE.0.95)THEN
          DTRH=0.25*(RHLCL-0.75)*QMIX/DQSDT
        ELSEIF(RHLCL.GT.0.95)THEN
          DTRH=(1./RHLCL-1.)*QMIX/DQSDT
        ELSE
          DTRH=0.
        ENDIF
      ENDIF   
C
      IF(TLCL+DTLCL+DTRH.GT.TENV)GO TO 150
      IF(ISHALL.EQ.1.AND.LC.EQ.NCHM)GO TO 900
      IF(KMIX.LE.LLFC)GO TO 25
C
C***  SHALLOW CONVECTION CHANGES
C
      IF(ISHALL.EQ.1)THEN
        CHMAX=0.
        NCHM=0
C
        DO NK=1,LLFC
          IF(CLDHGT(NK).GT.CHMAX)THEN
            NCHM=NK
            CHMAX=CLDHGT(NK)
          ENDIF
        ENDDO
C
        KMIX=NCHM
        GO TO 25
      ENDIF       
C***
      GO TO 900
  150 CONTINUE
C----------------------------------------------------------------------     
C     
C***  CONVECTIVE TRIGGERING CRITERIA HAS BEEN SATISFIED.
C***  COMPUTE EQUIVALENT POTENTIAL TEMPERATURE.
C***  (THETEU) AND VERTICAL VELOCITY OF THE RISING PARCEL AT THE LCL.
C     
      THMIX=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))
      THETEU(K)=THMIX*
     1          EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
      ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
      TVAVG=0.5*(TV0(KLCL)+TENV*(1.+0.608*QENV))
C
C***  MODIFY CALCULATION OF INITIAL PARCEL VERTICAL VELOCITY
C
      GDT=2.*G*(DTLCL+DTRH)*500./TVEN
      WLCL=1.+0.5*SQRT(ABS(GDT)+1.E-10)
      WLCL=AMIN1(WLCL,3.)
      PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
      WTW=WLCL*WLCL
      IF(WLCL.LT.0.)GO TO 25
      TVLCL=TLCL*(1.+0.608*QMIX)
      RHOLCL=PLCL/(R*TVLCL)
      PLCL0=PLCL
C     
      LCL=KLCL
      LET=LCL
C
      IF(WKL.LT.0.)THEN
        RAD=1000.
      ELSEIF(WKL.GT.0.1)THEN
        RAD=2000.
      ELSE
        RAD=1000.+1.E4*WKL
      ENDIF
C     
C*******************************************************************
C                                                                  *
C                 COMPUTE UPDRAFT PROPERTIES                       *
C                                                                  *
C*******************************************************************
C     
C***
C***  ESTIMATE INITIAL UPDRAFT MASS FLUX (UMF(K))
C***     
      WU(K)=WLCL
      AU0=PIE*RAD*RAD
      UMF(K)=RHOLCL*AU0
      VMFLCL=UMF(K)
      UPOLD=VMFLCL
      UPNEW=UPOLD
C     
C***  RATIO2 IS THE DEGREE OF GLACIATION IN THE CLOUD (0 TO 1),
C***  UER IS THE ENVIR ENTRAINMENT RATE, ABE IS AVAILABLE
C***  BUOYANT ENERGY, TRPPT IS THE TOTAL RATE OF PRECIPITATION
C***  PRODUCTION.
C     
      RATIO2(K)=0.
      UER(K)=0.
      ABE=0.
      TRPPT=0.
      TU(K)=TLCL
      TVU(K)=TVLCL
      QU(K)=QMIX
      EQFRC(K)=1.
      QLIQ(K)=0.
      QICE(K)=0.
      QLQOUT(K)=0.
      QICOUT(K)=0.
      DETLQ(K)=0.
      DETIC(K)=0.
      PPTLIQ(K)=0.
      PPTICE(K)=0.
      IFLAG=0
      KFRZ=LC
C     
C***  THE AMOUNT OF CONV AVAIL POT ENERGY (CAPE) IS CALCULATED WITH
C***  RESPECT TO UNDILUTE PARCEL ASCENT; EQ POT TEMP OF UNDILUTE
C***  PARCEL IS THTUDL, UNDILUTE TEMPERATURE IS GIVEN BY TUDL.
C     
      THTUDL=THETEU(K)
      TUDL=TLCL
C     
C***  TTEMP IS USED DURING CALCULATION OF THE LINEAR GLACIATION
C***  PROCESS; IT IS INITIALLY SET TO THE TEMPERATURE AT WHICH
C***  FREEZING IS SPECIFIED TO BEGIN.  WITHIN THE GLACIATION
C***  INTERVAL, IT IS SET EQUAL TO THE UPDRAFT TEMP AT THE
C***  PREVIOUS MODEL LEVEL.
C     
      TTEMP=TTFRZ
C     
C***  ENTER THE LOOP FOR UPDRAFT CALCULATIONS.  CALCULATE UPDRAFT TEMP,
C***  MIXING RATIO, VERTICAL MASS FLUX, LATERAL DETRAINMENT OF MASS AND
C***  MOISTURE, PRECIPITATION RATES AT EACH MODEL LEVEL.
C     
C     
C***  DIAGNOSTIC STUFF
C
      PLCL0=PLCL
      LFC=0
      PMIX0=PMIX
      REI=0.
      DILBE=0.
      UDLBE=0.
C
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
      DO 250 NK=K,KLM
      NK1=NK+1
      RATIO2(NK1)=RATIO2(NK)
      FRC1=0.
      TU(NK1)=T0(NK1)
      THETEU(NK1)=THETEU(NK)
      QU(NK1)=QU(NK)
      QLIQ(NK1)=QLIQ(NK)
      QICE(NK1)=QICE(NK)
C
      CALL TPMIX2(P0(NK1),THETEU(NK1),TU(NK1),QU(NK1),QLIQ(NK1),
     1            QICE(NK1),QNEWLQ,QNEWIC,RATIO2(NK1),
     2            XLV1,XLV0)
C----------------------------------------------------------------------     
C
C***  CHECK TO SEE IF UPDRAFT TEMP IS ABOVE THE TEMPERATURE AT WHICH
C***  GLACIATION IS ASSUMED TO INITIATE; IF IT IS, CALCULATE THE
C***  FRACTION OF REMAINING LIQUID WATER TO FREEZE.  FRZ IS THE
C***  TEMP AT WHICH FREEZING BEGINS, TBFRZ THE TEMP BELOW WHICH ALL
C***  LIQUID WATER IS FROZEN AT EACH LEVEL.
C
      IF(TU(NK1).LE.TTFRZ)THEN
C
        IF(TU(NK1).GT.TBFRZ)THEN
          IF(TTEMP.GT.TTFRZ)TTEMP=TTFRZ
          FRC1=(TTEMP-TU(NK1))/(TTEMP-TBFRZ)
        ELSE
          FRC1=1.
          IFLAG=1
        ENDIF
C
        TTEMP=TU(NK1)
C
C***  DETERMINE THE EFFECTS OF LIQUID WATER FREEZING WHEN TEMPERATURE
C***  IS BELOW TTFRZ
C
        QFRZ=(QLIQ(NK1)+QNEWLQ)*FRC1
        QNEWIC=QNEWIC+QNEWLQ*FRC1
        QNEWLQ=QNEWLQ-QNEWLQ*FRC1
        QICE(NK1)=QICE(NK1)+QLIQ(NK1)*FRC1
        QLIQ(NK1)=QLIQ(NK1)-QLIQ(NK1)*FRC1
        CALL DTFRZNEW(TU(NK1),P0(NK1),THETEU(NK1),QU(NK1),QFRZ,
     1                QICE(NK1),ALIQ,BLIQ,CLIQ,DLIQ)
      ENDIF
C----------------------------------------------------------------------     
C
      TVU(NK1)=TU(NK1)*(1.+0.608*QU(NK1))
C
C***  CALCULATE UPDRAFT VERTICAL VELOCITY AND PRECIPITATION FALLOUT...
C
      IF(NK.EQ.K)THEN
        BE=(TVLCL+TVU(NK1))/(TVEN+TV0(NK1))-1.
        BOTERM=2.*(Z00(NK1)-ZLCL)*G*BE/1.5
        DZZ=Z00(NK1)-ZLCL
      ELSE
        BE=(TVU(NK)+TVU(NK1))/(TV0(NK)+TV0(NK1))-1.
        BOTERM=2.*DZA(NK)*G*BE/1.5
        DZZ=DZA(NK)
      ENDIF
C
      ENTERM=2.*REI*WTW/UPOLD
C
C***  DIAGNOSTICS
C
      IF(TVU(NK1).GT.TV0(NK1))THEN
        IF(TVU(NK).LT.TV0(NK).OR.NK1.EQ.KLCL)LFC=NK1
      ENDIF
C
      WSQ=WTW
C
C----------------------------------------------------------------------     
      CALL CONDLOAD(QLIQ(NK1),QICE(NK1),WTW,DZZ,BOTERM,ENTERM,
     1              RATE,QNEWLQ,QNEWIC,QLQOUT(NK1),QICOUT(NK1))
C----------------------------------------------------------------------     
C
      WABS=SQRT(ABS(WTW))
      WU(NK1)=WTW/WABS
C
C***  IF VERT VELOCITY IS LESS THAN ZERO, EXIT THE UPDRAFT LOOP AND,
C***  IF CLOUD IS TALL ENOUGH, FINALIZE UPDRAFT CALCULATIONS.
C
      IF(WU(NK1).LT.0.)GO TO 275
C
C***  CALL SUBROUTINE TO CALCULATE ENVIRONMENTAL EQUIVALENT POTENTIAL
C***  TEMP...WITHIN GLACIATION INTERVAL, THETAE MUST BE CALCULATED
C***  WITH RESPECT TO THE SAME DEGREE OF GLACIATION FOR ALL ENTRAINING
C***  AIR.
C
C***  FOR LOOKUP TABLE VERSION, CALCULATE THETAE WITH RESPECT TO
C***  LIQUID WATER AT ALL LEVELS.
C
      IF(NK1.LE.KPBL)THEN
        CALL ENVIRTHT(P0(NK1),T0(NK1),QUER(NK1),THETEE(NK1),0.,
     1                RL,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
      ELSE
        CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),0.,
     1                RL,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
      ENDIF  
C
C***  REI IS THE RATE OF ENVIRONMENTAL INFLOW
C
      REI=VMFLCL*DP(NK1)*0.03/RAD
      TVQU(NK1)=TU(NK1)*(1.+0.608*QU(NK1)-QLIQ(NK1)-QICE(NK1))
C
      IF(NK.EQ.K)THEN
        DILBE=((TVLCL+TVQU(NK1))/(TVEN+TV0(NK1))-1.)*DZZ
      ELSE
        DILBE=((TVQU(NK)+TVQU(NK1))/(TV0(NK)+TV0(NK1))-1.)*DZZ
      ENDIF
C
      IF(DILBE.GT.0.)ABE=ABE+DILBE*G
C----------------------------------------------------------------------     
C
C***  IF CLOUD PARCELS ARE VIRTUALLY COLDER THAN THE ENVIRONMENT, 
C***  NO ENTRAINMENT IS ALLOWED AT THIS LEVEL.
C
      IF(TVQU(NK1).LE.TV0(NK1))THEN
C
C***  USE A MINIMUM ENTRAINMENT RATE (UER) OF 0.5*REI
C
        UER(NK1)=0.5*REI
        UDR(NK1)=1.5*REI
        EE2=0.5
        UD2=1.
        EQFRC(NK1)=0.
        GO TO 200
      ENDIF
C
      LET=NK1
      TTMP=TVQU(NK1)
C
C***  DETERMINE THE CRITICAL MIXED FRACTION OF UPDRAFT 
C***  AND ENVIRONMENTAL AIR
C
      F1=0.95
      F2=1.-F1
      THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
      QTMP=F1*Q0(NK1)+F2*QU(NK1)
      TMPLIQ=F2*QLIQ(NK1)
      TMPICE=F2*QICE(NK1)

      CALL TPMIX2(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,
     1            QNEWLQ,QNEWIC,RATIO2(NK1),
     2            XLV1,XLV0)
C
      TU95=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
C
      IF(TU95.GT.TV0(NK1))THEN
        EE2=1.
        UD2=0.
        EQFRC(NK1)=1.
        GO TO 175
      ENDIF
C
      F1=0.10
      F2=1.-F1
      THTTMP=F1*THETEE(NK1)+F2*THETEU(NK1)
      QTMP=F1*Q0(NK1)+F2*QU(NK1)
      TMPLIQ=F2*QLIQ(NK1)
      TMPICE=F2*QICE(NK1)
C
      CALL TPMIX2(P0(NK1),THTTMP,TTMP,QTMP,TMPLIQ,TMPICE,
     1            QNEWLQ,QNEWIC,RATIO2(NK1),
     2            XLV1,XLV0)
C
      TU10=TTMP*(1.+0.608*QTMP-TMPLIQ-TMPICE)
C
      IF(TU10.EQ.TVQU(NK1))THEN
        EE2=1.
        UD2=0.
        EQFRC(NK1)=1.0
        GO TO 175
      ENDIF
C
      EQFRC(NK1)=(TV0(NK1)-TVQU(NK1))*F1/(TU10-TVQU(NK1))
      EQFRC(NK1)=AMAX1(0.0,EQFRC(NK1))
      EQFRC(NK1)=AMIN1(1.0,EQFRC(NK1))
C
      IF(EQFRC(NK1).EQ.1)THEN
        EE2=1.
        UD2=0.
        GO TO 175
      ELSEIF(EQFRC(NK1).EQ.0.)THEN
        EE2=0.
        UD2=1.
        GO TO 175
      ELSE
C
C***  SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
C***  FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES
C
        CALL PROF5(EQFRC(NK1),EE2,UD2)
      ENDIF
C
  175 CONTINUE
C----------------------------------------------------------------------     
C
      IF(NK.EQ.K)THEN
        EE1=1.
        UD1=0.
      ENDIF
C
C***  NET ENTRAINMENT AND DETRAINMENT RATES ARE GIVEN BY THE AVERAGE 
C***  FRACTIONAL VALUES IN THE LAYER
C
C***  USE A MINIMUM ENTRAINMENT RATE (UER) OF 0.5*REI
C
      EE2=AMAX1(EE2,0.5)
      UD2=1.5*UD2
      UER(NK1)=0.5*REI*(EE1+EE2)
      UDR(NK1)=0.5*REI*(UD1+UD2)
C----------------------------------------------------------------------     
C
C***  IF THE CALCULATED UPDRAFT DETRAINMENT RATE IS GREATER THAN THE TOTAL
C***  UPDRAFT MASS FLUX, ALL CLOUD MASS DETRAINS, EXIT UPDRAFT CALCULATIONS.
C
  200 IF(UMF(NK)-UDR(NK1).LT.10.)THEN
C
C***  IF THE CALCULATED DETRAINED MASS FLUX IS GREATER THAN THE TOTAL 
C***  UPWARD MASS FLUX, IMPOSE TOTAL DETRAINMENT OF UPDRAFT MASS AT 
C***  THE PREVIOUS MODEL LVL
C
        IF(DILBE.GT.0.)ABE=ABE-DILBE*G
        LET=NK
        GO TO 275
      ENDIF
C
      EE1=EE2
      UD1=UD2
      UPOLD=UMF(NK)-UDR(NK1)
      UPNEW=UPOLD+UER(NK1)
      UMF(NK1)=UPNEW
      DILFRC(NK1)=UPNEW/UPOLD
C
C***  DETLQ AND DETIC ARE THE RATES OF DETRAINMENT OF LIQUID AND
C***  ICE IN THE DETRAINING UPDRAFT MASS
C
      DETLQ(NK1)=QLIQ(NK1)*UDR(NK1)
      DETIC(NK1)=QICE(NK1)*UDR(NK1)
      QDT(NK1)=QU(NK1)
C
      IF(NK1.LE.KPBL)THEN
        QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*QUER(NK1))/UPNEW
      ELSE
        QU(NK1)=(UPOLD*QU(NK1)+UER(NK1)*Q0(NK1))/UPNEW
      ENDIF   
C
      THETEU(NK1)=(THETEU(NK1)*UPOLD+THETEE(NK1)*UER(NK1))/UPNEW
      QLIQ(NK1)=QLIQ(NK1)*UPOLD/UPNEW
      QICE(NK1)=QICE(NK1)*UPOLD/UPNEW
C
C***  KFRZ IS THE HIGHEST MODEL LEVEL AT WHICH LIQUID CONDENSATE
C***  IS GENERATED...PPTLIQ IS THE RATE OF GENERATION (FALLOUT) OF
C***  LIQUID PRECIP AT A GIVEN MODEL LVL, PPTICE THE SAME FOR ICE,
C***  TRPPT IS THE TOTAL RATE OF PRODUCTION OF PRECIP UP TO THE
C***  CURRENT MODEL LEVEL...
C
      IF(ABS(RATIO2(NK1)-1.).GT.1.E-6)KFRZ=NK1
C
C***  REVERSE THE MOD THAT ALLOWS FEEDBACK OF RAIN/SNOW THAT ORIGINATES 
C***  IN DETRAINED AIR
C
      PPTLIQ(NK1)=QLQOUT(NK1)*UMF(NK)
      PPTICE(NK1)=QICOUT(NK1)*UMF(NK)
C
      TRPPT=TRPPT+PPTLIQ(NK1)+PPTICE(NK1)
      IF(NK1.LE.KPBL)UER(NK1)=UER(NK1)+VMFLCL*DP(NK1)/DPTHMX
  250 CONTINUE
C----------------------------------------------------------------------     
C
C***  CHECK CLOUD DEPTH.  IF CLOUD IS TALL ENOUGH, ESTIMATE THE 
C***  EQUILIBRIUM TEMPERATURE LEVEL (LET) AND ADJUST MASS FLUX PROFILE 
C***  AT CLOUD TOP SO THAT MASS FLUX DECREASES TO ZERO AS A LINEAR 
C***  FUNCTION OF PRESSURE BETWEEN THE LET AND CLOUD TOP.
C     
C***  LTOP IS THE MODEL LEVEL JUST BELOW THE LEVEL AT WHICH VERTICAL 
C***  VELOCITY FIRST BECOMES NEGATIVE
C     
  275 LTOP=NK
      CLDHGT(LC)=Z00(LTOP)-ZLCL 
C     
C***  IF CLOUD TOP HEIGHT IS LESS THAN THE SPECIFIED MINIMUM FOR DEEP 
C***  CONVECTION, SAVE VALUE TO CONSIDER THIS LEVEL AS SOURCE FOR 
C***  SHALLOW CONVECTION, GO BACK UP TO CHECK NEXT LEVEL.
C     
C***  TRY SPECIFYING MINIMUM CLOUD DEPTH AS A FUNCTION OF TLCL
C
      IF(TLCL.GT.293.)THEN
        CHMIN=4.E3
      ELSEIF(TLCL.LE.293. .and. TLCL.GE.273)THEN
        CHMIN=2.E3 + 100.*(TLCL-273.)
      ELSEIF(TLCL.LT.273.)THEN
        CHMIN=2.E3
      ENDIF
C
      KSTART=MAX0(KPBL,KLCL)
C
C----------------------------------------------------------------------     
C***  DO NOT ALLOW ANY CLOUD FROM THIS LAYER IF:
C
C   1.) IF THERE IS NO CAPE, OR 
C   2.) CLOUD TOP IS AT MODEL LEVEL JUST ABOVE LCL, OR
C   3.) CLOUD TOP IS WITHIN UPDRAFT SOURCE LAYER, OR
C   4.) CLOUD-TOP DETRAINMENT LAYER BEGINS WITHIN 
C       UPDRAFT SOURCE LAYER
C
      IF(LTOP.LE.KLCL.OR. 
     1   LTOP.LE.KPBL.OR.LET+1.LE.KPBL)THEN
        CLDHGT(LC)=0.
C
C***  IF THIS IS SELECTED SHALLOW SOURCE AND STILL DOES NOT MAKE A 
C***  SIGNIFICANT CLOUD, GO ON TO NEXT GRID POINT.
C
        IF(LC.EQ.NCHM)GO TO 900
C
C***  IF ALL LAYERS IN THE SPECIFIED PORTION OF LOWER ATMOSPHERE
C***  HAVE BEEN CHECKED, THEN...
C
        IF(KMIX.GT.LLFC)THEN
C
C***  IF NO POSSIBLE SHALLOW CLOUD LAYERS WERE FOUND, GO TO NEXT 
C***  GRID POINT
C
          IF(ISHALL.EQ.0)THEN
            GO TO 900
          ELSE
C
C***  IF SOME POTENTIAL SHALLOW CLOUD SOURCE LAYERS WERE FOUND, 
C***  FIND THE ONE THAT GIVES THE TALLEST CLOUD.
C
            GO TO 305
          ENDIF
        ELSE
C
C***  IF THERE ARE MORE LAYERS TO CHECK, RESET CLOUD CHARACTERISTICS,
C***  GO TO NEXT LEVEL.
          GO TO 310
        ENDIF 
C
C***  IF THIS LAYER HAS BEEN SELECTED AS A SHALLOW CONVECTIVE SOURCE,
C***  ALLOW SHALLOW CONVECTION.  (EVEN IF THE SHALLOW-CLOUD PARAMETERS
C***  ALLOW CLDHGT TO EXCEED CHMIN)
C
      ELSEIF(LC.EQ.NCHM)THEN
        GO TO 315
C
C***  IF CLOUD DEPTH IS GREATER THAN MINIMUM DEPTH CRITERION, AND
C***  THIS LAYER HAS NOT BEEN MARKED AS A SHALLOW CONVECTIVE SOURCE LAYER,
C***  ALLOW DEEP CONVECTION.
C
      ELSEIF(CLDHGT(LC).GT.CHMIN.AND.ABE.GT.1)THEN
        ISHALL=0
        GO TO 315
      ENDIF
C
C***  AT THIS POINT, WE HAVE A PARCEL THAT IS ABLE TO MAINTAIN UPWARD 
C***  MOMENTUM FOR AT LEAST A SHORT DISTANCE, BUT THE CLOUD FROM THIS 
C***  SOURCE LAYER (LC) IS NOT DEEP ENOUGH FOR "DEEP" (PRECIPITATING)
C***  CONVECTION.  SET ISHALL=1 TO SAVE LC AS A POSSIBLE FOR SOURCE 
C***  LAYER FOR SHALLOW CONVECTION.
C
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
C
C***  TO DISALLOW SHALLOW CONVECTION, COMMENT OUT NEXT LINE !!!!!!!!
C
      ISHALL=1
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
      IF(KMIX.GT.LLFC)THEN
        GO TO 305
      ELSE
        GO TO 310
      ENDIF
C
  305 CONTINUE
C
      IF(ISHALL.EQ.0)THEN
        GO TO 900
      ELSE
        CHMAX=0.
        NCHM=0
        DO NK=1,LLFC
          IF(CLDHGT(NK).GT.CHMAX)THEN
            NCHM=NK
            CHMAX=CLDHGT(NK)
          ENDIF
        ENDDO
        KMIX=NCHM
      ENDIF
C
  310 CONTINUE
C
      DO NK=K,LTOP
        UMF(NK)=0.
        UDR(NK)=0.
        UER(NK)=0.
        DETLQ(NK)=0.
        DETIC(NK)=0.
        PPTLIQ(NK)=0.
        PPTICE(NK)=0.
      ENDDO
C
      GO TO 25        
C
  315 CONTINUE   
      IF(ISHALL.EQ.1)THEN
        KSTART=MAX0(KPBL,KLCL)
        LET=KSTART
      ENDIF
C----------------------------------------------------------------------     
C     
C***  IF THE LET AND LTOP ARE THE SAME, DETRAIN ALL OF THE UPDRAFT MASS 
C***  FLUX AT THIS LEVEL.
C     
      IF(LET.EQ.LTOP)THEN
        UDR(LTOP)=UMF(LTOP)+UDR(LTOP)-UER(LTOP)
        DETLQ(LTOP)=QLIQ(LTOP)*UDR(LTOP)*UPNEW/UPOLD
        DETIC(LTOP)=QICE(LTOP)*UDR(LTOP)*UPNEW/UPOLD
C
C***  REVERSE THE MOD THAT ALLOWS FEEDBACK OF RAIN/SNOW THAT ORIGINATES 
C***  IN DETRAINED AIR.
C
        UER(LTOP)=0.
        UMF(LTOP)=0.
        GO TO 350
      ENDIF
C     
C***  BEGIN TOTAL DETRAINMENT AT THE LEVEL ABOVE THE LET
C     
      DPTT=0.
C
      DO NJ=LET+1,LTOP
        DPTT=DPTT+DP(NJ)
      ENDDO
C
      DUMFDP=UMF(LET)/DPTT
C----------------------------------------------------------------------     
C     
C***  ADJUST MASS FLUX PROFILES, DETRAINMENT RATES, AND PRECIPITATION FALL
C***  RATES TO REFLECT THE LINEAR DECREASE IN MASS FLX BETWEEN THE LET AND
C     
      DO 325 NK=LET+1,LTOP
C
C***  ENTRAINMENT IS ALLOWED AT EVERY LEVEL EXCEPT FOR LTOP, SO DISALLOW
C***  ENTRAINMENT AT LTOP AND ADJUST ENTRAINMENT RATES BETWEEN LET AND LTOP
C***  SO THE THE DILUTION FACTOR DUE TO ENTRAINMENT IS NOT CHANGED BUT 
C***  THE ACTUAL ENTRAINMENT RATE WILL CHANGE DUE DUE FORCED TOTAL 
C***  DETRAINMENT IN THIS LAYER.
C
      IF(NK.EQ.LTOP)THEN
        UDR(NK)=UMF(NK-1)
        UER(NK)=0.
        DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
        DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
      ELSE
        UMF(NK)=UMF(NK-1)-DP(NK)*DUMFDP
        UER(NK)=UMF(NK)*(1.-1./DILFRC(NK))
        UDR(NK)=UMF(NK-1)-UMF(NK)+UER(NK)
        DETLQ(NK)=UDR(NK)*QLIQ(NK)*DILFRC(NK)
        DETIC(NK)=UDR(NK)*QICE(NK)*DILFRC(NK)
      ENDIF
C
C***  REVERSE THE MOD THAT ALLOWS FEEDBACK OF RAIN/SNOW THAT ORIGINATES 
C***  IN DETRAINED AIR.
C
      IF(NK.GE.LET+2)THEN
        TRPPT=TRPPT-PPTLIQ(NK)-PPTICE(NK)
        PPTLIQ(NK)=UMF(NK-1)*QLQOUT(NK)
        PPTICE(NK)=UMF(NK-1)*QICOUT(NK)
        TRPPT=TRPPT+PPTLIQ(NK)+PPTICE(NK)
      ENDIF
C
  325 CONTINUE
C
  350 CONTINUE
C----------------------------------------------------------------------     
C     
C***  SEND UPDRAFT CHARACTERISTICS TO OUTPUT FILES
C     
C***  WHEN MOIST PERTURBATION IS ADDED TO UPDRAFT AND KLCL.LE.KPBL,
C***  RESET THETEE SO THAT MOISTURE PERTURBATION ONLY AFFECTS UPDRAFTS.
C
      IF(KLCL.LE.KPBL)THEN
        DO NK1=KLCL,KPBl
          CALL ENVIRTHT(P0(NK1),T0(NK1),Q0(NK1),THETEE(NK1),0.,
     1                  RL,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
        ENDDO
      ENDIF  
C
      XTIME=NTSD*DT/60.
C     
C***  EXTEND THE UPDRAFT MASS FLUX PROFILE DOWN TO THE SOURCE LAYER FOR
C***  THE UPDRAFT AIR.  ALSO, DEFINE THETAE FOR LEVELS BELOW THE LCL.
C     
      DO 360 NK=1,K
C
      IF(NK.GE.LC)THEN
        IF(NK.EQ.LC)THEN
          UMF(NK)=VMFLCL*DP(NK)/DPTHMX
          UER(NK)=VMFLCL*DP(NK)/DPTHMX
        ELSEIF(NK.LE.KPBL)THEN
          UER(NK)=VMFLCL*DP(NK)/DPTHMX
          UMF(NK)=UMF(NK-1)+UER(NK)
        ELSE
          UMF(NK)=VMFLCL
          UER(NK)=0.
        ENDIF
        TU(NK)=TMIX+(Z00(NK)-ZMIX)*GDRY
        QU(NK)=QMIX
        WU(NK)=WLCL
      ELSE
        TU(NK)=0.
        QU(NK)=0.
        UMF(NK)=0.
        WU(NK)=0.
        UER(NK)=0.
      ENDIF
C
      UDR(NK)=0.
      QDT(NK)=0.
      QLIQ(NK)=0.
      QICE(NK)=0.
      QLQOUT(NK)=0.
      QICOUT(NK)=0.
      PPTLIQ(NK)=0.
      PPTICE(NK)=0.
      DETLQ(NK)=0.
      DETIC(NK)=0.
      RATIO2(NK)=0.
      EE=Q0(NK)*P0(NK)/(0.622+Q0(NK))
      TLOG=ALOG(EE/ALIQ)
      TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
      TSAT=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(T0(NK)-T00))*(
     1     T0(NK)-TDPT)
      THTA=T0(NK)*(1.E5/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
      THETEE(NK)=THTA*
     1           EXP((3374.6525/TSAT-2.5403)*Q0(NK)*(1.+0.81*Q0(NK)))
      THTES(NK)=THTA*
     1          EXP((3374.6525/T0(NK)-2.5403)*QES(NK)*(1.+0.81*
     2               QES(NK)))
      EQFRC(NK)=1.0
  360 CONTINUE
C----------------------------------------------------------------------     
C
C***  FIND HEIGHT OF LFC, CALCULATE CIN.
C
      IF(LFC.GT.0)THEN
        CIN1=0.
        IF(LFC.EQ.KLCL)THEN
          IF(TVLCL.GT.TVEN)THEN
            CIN1=0.
          ELSE
            Z00LFC=Z00(LFC)
            IF(Z00LFC.NE.ZLCL)THEN
              DTUDZ=(TVU(LFC)-TVLCL)/(Z00LFC-ZLCL)
              DTEDZ=(TV0(LFC)-TVEN )/(Z00LFC-ZLCL)
              ZLFC=ZLCL+(TVEN-TVLCL)/(DTUDZ-DTEDZ)
              TVLFC=TVEN+(ZLFC-ZLCL)*DTEDZ
            ELSE
              ZLFC=ZLCL
              TVLFC=TVEN
            ENDIF
            CIN1=G*(ZLFC-ZLCL)*(TVLCL-TVEN)/(TVEN+TVLFC)
          ENDIF
        ELSE
          CIN1=CIN1+
     1         G*(Z00(KLCL)-ZLCL)*((TVU(KLCL)+TVLCL)/
     2         (TV0(KLCL)+TVEN)-1.)     
          DO NK=KLCL+1,LFC
            IF(NK.EQ.LFC)THEN
              DTUDZ=(TVU(LFC)-TVU(LFC-1))/(Z00(LFC)-Z00(LFC-1))
              DTEDZ=(TV0(LFC)-TV0(LFC-1))/(Z00(LFC)-Z00(LFC-1))
              ZLFC=Z00(LFC-1)+(TV0(LFC-1)-TV0(LFC-1))/(DTUDZ-DTEDZ)
              TVEN=TV0(LFC-1)+(ZLFC-Z00(LFC-1))*DTEDZ
              CIN1=CIN1+G*(ZLFC-Z00(LFC-1))*(TVU(LFC-1)-TV0(LFC-1))/
     1                (TVEN+TV0(LFC-1))
            ELSE
              CIN1=CIN1+G*(Z00(NK)-Z00(NK-1))*
     1              ((TVU(NK)+TVU(NK-1))/(TV0(NK)+TV0(NK-1))-1.)
            ENDIF
          ENDDO
        ENDIF
C
      ENDIF
C----------------------------------------------------------------------     
C     
      LTOP1=LTOP+1
      LTOPM1=LTOP-1
C     
C***  DEFINE VARIABLES ABOVE CLOUD TOP
C     
      DO NK=LTOP1,KX
        UMF(NK)=0.
        UDR(NK)=0.
        UER(NK)=0.
        QDT(NK)=0.
        QLIQ(NK)=0.
        QICE(NK)=0.
        QLQOUT(NK)=0.
        QICOUT(NK)=0.
        DETLQ(NK)=0.
        DETIC(NK)=0.
        PPTLIQ(NK)=0.
        PPTICE(NK)=0.
C
        IF(NK.GT.LTOP1)THEN
          TU(NK)=0.
          QU(NK)=0.
          WU(NK)=0.
        ENDIF
C
        THTA0(NK)=0.
        THTAU(NK)=0.
        EMS(NK)=0.
        EMSD(NK)=0.
        TG(NK)=T0(NK)
        QG(NK)=Q0(NK)
        QLG(NK)=0.
        QIG(NK)=0.
        QRG(NK)=0.
        QSG(NK)=0.
        OMG(NK)=0.
      ENDDO
C
      OMG(KXP1)=0.
      P150=P0(KLCL)-1.50E4
C----------------------------------------------------------------------     
      DO 375 NK=1,LTOP
      EMS(NK)=DP(NK)*DXSQ/G
      EMSD(NK)=1./EMS(NK)
C     
C***  INITIALIZE SOME VARIABLES TO BE USED LATER IN THE
C***  VERTICAL ADVECTION SCHEME
C     
      EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QDT(NK)))
      THTAU(NK)=TU(NK)*EXN(NK)
      EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*Q0(NK)))
      THTA0(NK)=T0(NK)*EXN(NK)
C     
C***  LVF IS THE LEVEL AT WHICH MOISTURE FLUX IS ESTIMATED AS THE 
C***  BASIS FOR PRECIPITATION EFFICIENCY CALCULATIONS.
C     
      IF(P0(NK).GT.P150)LVF=NK
      DDILFRC(NK)=1./DILFRC(NK)
      OMG(NK)=0.
  375 CONTINUE
C----------------------------------------------------------------------     
      LVF=MIN0(LVF,LET)     
      USR=UMF(LVF+1)*(QU(LVF+1)+QLIQ(LVF+1)+QICE(LVF+1))
      USR=AMIN1(USR,TRPPT)
C     
C***  COMPUTE CONVECTIVE TIME SCALE(TIMEC). THE MEAN WIND AT THE LCL
C***  AND MIDTROPOSPHERE IS USED.
C     
      WSPD(KLCL)=SQRT(U0(KLCL)*U0(KLCL)+V0(KLCL)*V0(KLCL))
      WSPD(L5)=SQRT(U0(L5)*U0(L5)+V0(L5)*V0(L5))
      WSPD(LTOP)=SQRT(U0(LTOP)*U0(LTOP)+V0(LTOP)*V0(LTOP))
      VCONV=0.5*(WSPD(KLCL)+WSPD(L5))
      TIMEC=DX(I,J)/VCONV
      TADVEC=TIMEC
      TIMEC=AMAX1(1800.,TIMEC)
      TIMEC=AMIN1(3600.,TIMEC)
      IF(ISHALL.EQ.1)TIMEC=2400.
      NIC=NINT(TIMEC/(0.5*DT2))
      TIMEC=FLOAT(NIC)*0.5*DT2
C   
C***  COMPUTE WIND SHEAR AND PRECIPITATION EFFICIENCY.
C     
      IF(WSPD(LTOP).GT.WSPD(KLCL))THEN
        SHSIGN=1.
      ELSE
        SHSIGN=-1.
      ENDIF
C
      VWS=(U0(LTOP)-U0(KLCL))*(U0(LTOP)-U0(KLCL))+(V0(LTOP)-V0(KLCL))*
     1    (V0(LTOP)-V0(KLCL))
      VWS=1.E3*SHSIGN*SQRT(VWS)/(Z00(LTOP)-Z00(LCL))
      PEF=1.591+VWS*(-.639+VWS*(9.53E-2-VWS*4.96E-3))
      PEF=AMAX1(PEF,.2)
      PEF=AMIN1(PEF,.9)
C     
C***  PRECIPITATION EFFICIENCY IS A FUNCTION OF THE HEIGHT OF CLOUD BASE
C     
      CBH=(ZLCL-Z00(1))*3.281E-3
C
      IF(CBH.LT.3.)THEN
        RCBH=.02
      ELSE
        RCBH=0.96729352+CBH*(-.70034167+CBH*(.162179896+CBH*
     1         (-1.2569798E-2+CBH*(4.2772E-4-CBH*5.44E-6))))
      ENDIF
C
      IF(CBH.GT.25)RCBH=2.4
      PEFCBH=1./(1.+RCBH)
      PEFCBH=AMIN1(PEFCBH,.9)
C     
C***  MEAN PRECIP EFFICIENCY IS USED TO COMPUTE RAINFALL.
C     
      PEFF=.5*(PEF+PEFCBH)
      PEFF2=PEFF   
C
C*****************************************************************
C                                                                *
C                  COMPUTE DOWNDRAFT PROPERTIES                  *
C                                                                *
C*****************************************************************
C     
      TDER=0.
C
      IF(ISHALL.EQ.1)THEN
        LFS=1
        GO TO 450
      ENDIF   
C
C***  START DOWNDRAFT ABOUT 150 MB ABOVE CLOUD BASE
C
      KSTART=KPBL+1     
C
      DO NK=KSTART+1,KL
        DPPP=P0(KSTART)-P0(NK)
        IF(DPPP.GT.150.E2)THEN
          KLFS=NK
          GO TO 405
        ENDIF
      ENDDO
C
  405 CONTINUE
C
      KLFS=MIN0(KLFS,LET-1)
      LFS=KLFS
C
C***  IF LFS IS NOT AT LEAST 50 MB ABOVE CLOUD BASE (IMPLYING THAT THE 
C***  LEVEL OF EQUIL TEMP, LET, IS JUST ABOVE CLOUD BASE) DO NOT ALLOW 
C***  A DOWNDRAFT.
C
      IF((P0(KSTART)-P0(LFS)).LT.50.E2)THEN
        TDER=0.
        GO TO 450
      ENDIF
C
      THETED(LFS)=THETEE(LFS)
      QD(LFS)=Q0(LFS)
C
C***  CALL TPMIX2DD TO FIND WET-BULB TEMP, QV
C
      CALL TPMIX2DD(P0(LFS),THETED(LFS),TZ(LFS),QSS)
C
      THTAD(LFS)=TZ(LFS)*(P00/P0(LFS))**(0.2854*(1.-0.28*QSS))
C     
C***  TAKE A FIRST GUESS AT THE INITIAL DOWNDRAFT MASS FLUX...
C     
      TVD(LFS)=TZ(LFS)*(1.+0.608*QSS)
      RDD=P0(LFS)/(R*TVD(LFS))
      A1=(1.-PEFF)*AU0
      DMF(LFS)=-A1*RDD
      DER(LFS)=DMF(LFS)
      DDR(LFS)=0.
      RHBAR=RH(LFS)*DP(LFS)
      DPTT=DP(LFS)
C
      DO ND=LFS-1,KSTART,-1
        ND1=ND+1
        DER(ND)=DER(LFS)*EMS(ND)/EMS(LFS)
        DDR(ND)=0.
        DMF(ND)=DMF(ND1)+DER(ND)
        THETED(ND)=(THETED(ND1)*DMF(ND1)+THETEE(ND)*DER(ND))/DMF(ND)
        QD(ND)=(QD(ND1)*DMF(ND1)+Q0(ND)*DER(ND))/DMF(ND)    
        DPTT=DPTT+DP(ND)
        RHBAR=RHBAR+RH(ND)*DP(ND)
      ENDDO
C
      RHBAR=RHBAR/DPTT
      DMFFRC=2.*(1.-RHBAR)
      DPDD=0.
C----------------------------------------------------------------------     
C
C***  CALCULATE MELTING EFFECT.
C***  FIRST, COMPUTE TOTAL FROZEN PRECIPITATION GENERATED.
C
      PPTMLT=0.
C
      DO NK=KLCL,LTOP
        PPTMLT=PPTMLT+PPTICE(NK)
      ENDDO
C
C----------------------------------------------------------------------     
      IF(LC.LT.ML)THEN
C
C***  FOR NOW, CALCULATE MELTING EFFECT AS IF DMF=-UMF AT KLCL, I.E., 
C***  AS IF DMFFRC=1.  OTHERWISE, FOR SMALL DMFFRC, DTMELT GETS TOO LARGE.
C
        DTMELT=RLF*PPTMLT/(CP*UMF(KLCL))
      ELSE
        DTMELT=0.
      ENDIF
C
      LDT=MIN0(LFS-1,KSTART-1)
      CALL TPMIX2DD(P0(KSTART),THETED(KSTART),TZ(KSTART),QSS)
      TZ(KSTART)=TZ(KSTART)-DTMELT
      ES=ALIQ*EXP((BLIQ*TZ(KSTART)-CLIQ)/(TZ(KSTART)-DLIQ))
      QSS=0.622*ES/(P0(KSTART)-ES)
      THETED(KSTART)=TZ(KSTART)*(1.E5/P0(KSTART))**
     1              (0.2854*(1.-0.28*QSS))*
     2         EXP((3374.6525/TZ(KSTART)-2.5403)*QSS*(1.+0.81*QSS))
C
      LDT=MIN0(LFS-1,KSTART-1)
C
C----------------------------------------------------------------------     
      DO 425 ND=LDT,1,-1
      DPDD=DPDD+DP(ND)
      THETED(ND)=THETED(KSTART)
      QD(ND)=QD(KSTART)       
C
C***  CALL TPMIX2DD TO FIND WET BULB TEMP, SATURATION MIXING RATIO
C
      CALL TPMIX2DD(P0(ND),THETED(ND),TZ(ND),QSS)
      QSD(ND)=QSS
C
C***  SPECIFY RH DECREASE OF 10%/KM IN DOWNDRAFT
C
      RHH=1.-0.2/1000.*(Z00(KSTART)-Z00(ND))
C
C***  ADJUST DOWNDRAFT TEMP, Q TO SPECIFIED RH
C
      IF(RHH.LT.1.)THEN
        DSSDT=(CLIQ-BLIQ*DLIQ)/((TZ(ND)-DLIQ)*(TZ(ND)-DLIQ))
        RL=XLV0-XLV1*TZ(ND)
        DTMP=RL*QSS*(1.-RHH)/(CP+RL*RHH*QSS*DSSDT)
        T1RH=TZ(ND)+DTMP
        ES=RHH*ALIQ*EXP((BLIQ*T1RH-CLIQ)/(T1RH-DLIQ))
        QSRH=0.622*ES/(P0(ND)-ES)
C
C***  CHECK TO SEE IF MIXING RATIO AT SPECIFIED RH IS LESS THAN 
C***  ACTUAL MIXING RATIO.  IF SO, ADJUST TO GIVE ZERO EVAPORATION.
C
        IF(QSRH.LT.QD(ND))THEN
          QSRH=QD(ND)
          T1RH=TZ(ND)+(QSS-QSRH)*RL/CP
        ENDIF
C
        TZ(ND)=T1RH
        QSS=QSRH
        QSD(ND)=QSS
      ENDIF         
C
      TVD(ND)=TZ(ND)*(1.+0.608*QSD(ND))
C
      IF(TVD(ND).GT.TV0(ND).OR.ND.EQ.1)THEN
        LDB=ND
        GO TO 430
      ENDIF
C
  425 CONTINUE
C----------------------------------------------------------------------     
  430 CONTINUE
C
      IF((P0(LDB)-P0(LFS)).LT.50.E2)THEN   ! No Downdraft allowed! 
        TDER=0.
        GO TO 450
      ENDIF
C
      TDER=0.
C     
C***  CALCULATE AN EVAPORATION RATE FOR GIVEN MASS FLUX
C     
      DO ND=LDT,LDB,-1
        ND1=ND+1
        DDR(ND)=-DMF(KSTART)*DP(ND)/DPDD
        DER(ND)=0.
        DMF(ND)=DMF(ND1)+DDR(ND)
        TDER=TDER+(QSD(nd)-QD(ND))*DDR(ND)
        QD(ND)=QSD(ND)
        THTAD(ND)=TZ(ND)*(P00/P0(ND))**(0.2854*(1.-0.28*QD(ND)))
      ENDDO
C
  450 CONTINUE
C----------------------------------------------------------------------     
C     
C***  IF DOWNDRAFT DOES NOT EVAPORATE ANY WATER FOR SPECIFIED RELATIVE
C***  HUMIDITY, NO DOWNDRAFT IS ALLOWED.
C     
      IF(TDER.LT.1.)THEN
        PPTFLX=TRPPT
        CPR=TRPPT
        TDER=0.
        CNDTNF=0.
        UPDINC=1.
        LDB=LFS
C
        DO NDK=1,LTOP
          DMF(NDK)=0.
          DER(NDK)=0.
          DDR(NDK)=0.
          THTAD(NDK)=0.
          WD(NDK)=0.
          TZ(NDK)=0.
          QD(NDK)=0.
        ENDDO
C
        AINCM2=100.
        GO TO 475
      ENDIF
C     
C***  ADJUST DOWNDRAFT MASS FLUX SO THAT EVAPORATION RATE IN DOWNDRAFT
C***  IS CONSISTENT WITH PRECIPITATION EFFICIENCY RELATIONSHIP.
C     
C***  DDINC IS THE FACTOR BY WHICH TO INCREASE THE FIRST-GUESS DOWNDRAFT 
C***  MOMENTUM FLUX TO SATISFY THE PRECIP EFFICIENCY RELATIONSHIP.
C***  UPDINC IS THE FACTOR BY WHICH TO INCREASE THE UPDRAFT MASS FLUX BELOW 
C***  THE LFS TO ACCOUNT FOR TRANSFER OF MASS FROM UPDRAFT TO DOWNDRAFT.
C     
      DDINC=-DMFFRC*UMF(KLCL)/DMF(KSTART)
      UPDINC=1.
C
      IF(TDER*DDINC.GT.TRPPT)THEN
        DDINC=TRPPT/TDER
      ENDIF
C
      TDER=TDER*DDINC
C
      DO NK=LDB,LFS
        DMF(NK)=DMF(NK)*DDINC
        DER(NK)=DER(NK)*DDINC
        DDR(NK)=DDR(NK)*DDINC
      ENDDO
C     
      CPR=TRPPT
      PPTFLX=TRPPT-TDER
C     
C***  ADJUST UPDRAFT MASS FLUX, MASS DETRAINMENT RATE, AND LIQUID WATER 
C***  AND DETRAINMENT RATES TO BE CONSISTENT WITH THE TRANSFER OF THE 
C***  ESTIMATE FROM THE UPDRAFT TO THE DOWNDRAFT AT THE LFS.
C     
      DO NK=LC,LFS
        UMF(NK)=UMF(NK)*UPDINC
        UDR(NK)=UDR(NK)*UPDINC
        UER(NK)=UER(NK)*UPDINC
        PPTLIQ(NK)=PPTLIQ(NK)*UPDINC
        PPTICE(NK)=PPTICE(NK)*UPDINC
        DETLQ(NK)=DETLQ(NK)*UPDINC
      ENDDO
C     
C***  ZERO OUT THE ARRAYS FOR DOWNDRAFT DATA AT LEVELS ABOVE AND BELOW 
C***  THE DOWNDRAFT
C     
      IF(LDB.GT.1)THEN
        DO NK=1,LDB-1
          DMF(NK)=0.
          DER(NK)=0.
          DDR(NK)=0.
          WD(NK)=0.
          TZ(NK)=0.
          QD(NK)=0.
          THTAD(NK)=0.
        ENDDO
      ENDIF
C
      DO NK=LFS+1,KX
        DMF(NK)=0.
        DER(NK)=0.
        DDR(NK)=0.
        WD(NK)=0.
        TZ(NK)=0.
        QD(NK)=0.
        THTAD(NK)=0.
      ENDDO
C
      DO NK=LDT+1,LFS-1
        TZ(NK)=0.
        QD(NK)=0.
        THTAD(NK)=0. 
      ENDDO
C----------------------------------------------------------------------     
  475 CONTINUE
C----------------------------------------------------------------------     
C
C***  SET LIMITS ON THE UPDRAFT AND DOWNDRAFT MASS FLUXES SO THAT THE 
C***  INFLOW INTO CONVECTIVE DRAFTS FROM A GIVEN LAYER IS NO MORE THAN 
C***  IS AVAILABLE IN THAT LAYER INITIALLY.
C     
      AINCMX=1000.
      LMAX=MAX0(KLCL,LFS)
C
      DO NK=LC,LMAX
        IF((UER(NK)-DER(NK)).GT.0.)AINCM1=EMS(NK)/((UER(NK)-DER(NK))*
     1                                    TIMEC)
        AINCMX=AMIN1(AINCMX,AINCM1)
      ENDDO
C
      AINC=1.
      IF(AINCMX.LT.AINC)AINC=AINCMX
C     
C***  SAVE THE RELEVENT VARIABLES FOR A UNIT UPDRFT AND DOWNDRFT.
C***  THEY WILL BE ADJUSTED ITERATIVELY BY THE FACTOR AINC TO SATISFY 
C***  THE STABILIZATION CLOSURE.
C     
      NCOUNT=0
      TDER2=TDER
      PPTFL2=PPTFLX
C
      DO NK=1,LTOP
        DETLQ2(NK)=DETLQ(NK)
        DETIC2(NK)=DETIC(NK)
        UDR2(NK)=UDR(NK)
        UER2(NK)=UER(NK)
        DDR2(NK)=DDR(NK)
        DER2(NK)=DER(NK)
        UMF2(NK)=UMF(NK)
        DMF2(NK)=DMF(NK)
      ENDDO
C
      FABE=1.
      STAB=0.95
      NOITR=0
      ISTOP=0
C
      IF(AINC/AINCMX.GT.0.999)THEN
        NCOUNT=0
        GO TO 575
      ENDIF
C----------------------------------------------------------------------     
C
C***  SHALLOW CONVECTION EFFECTS
C
      IF(ISHALL.EQ.1)THEN
C
C***  FIND THE MAXIMUM tke VALUE BETWEEN LC AND KLCL
C
        TKEMAX=0.
        DO K=LC,KLCL
          NK=KX-K+1-NLEV
          TKEMAX=AMAX1(TKEMAX,Q2(I,J,NK))
        ENDDO
C
        TKEMAX=AMIN1(TKEMAX,10.)
        TKEMAX=AMAX1(TKEMAX,5.)
C
C***  DPMIN WAS CHANGED FOR SHALLOW CONVECTION SO THAT IT IS THE
C***  THE SAME AS FOR DEEP CONVECTION (5.E3).  sINCE THIS DOUBLES
C***  (ROUGHLY) THE VALUE OF DPTHMX, ADD A FACTOR OF 0.5 TO THE
C***  CALCULATION OF EVAC.
C
        EVAC=0.5*TKEMAX*0.1
        AINC=EVAC*DPTHMX*DX(I,J)*DX(I,J)/(VMFLCL*G*TIMEC)
        GO TO 575
      ENDIF
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
C     
C*****************************************************************
C                                                                *
C           COMPUTE PROPERTIES FOR COMPENSATIONAL SUBSIDENCE     *
C                                                                *
C*****************************************************************
C     
  500 NCOUNT=NCOUNT+1
C     
C***  DETERMINE OMEGA VALUE NECESSARY AT TOP AND BOTTOM OF EACH LAYER 
C***  TO SATISFY MASS CONTINUITY.
C     
      DTT=TIMEC
C
      DO NK=1,LTOP
        DOMGDP(NK)=-(UER(NK)-DER(NK)-UDR(NK)-DDR(NK))*EMSD(NK)
        IF(NK.GT.1)THEN
          OMG(NK)=OMG(NK-1)-DP(NK-1)*DOMGDP(NK-1)
          DTT1=0.75*DP(NK-1)/(ABS(OMG(NK))+1.E-10)
          DTT=AMIN1(DTT,DTT1)
        ENDIF
      ENDDO
C
      DO NK=1,LTOP
        THPA(NK)=THTA0(NK)
        QPA(NK)=Q0(NK)
        NSTEP=NINT(TIMEC/DTT+1)
        DTIME=TIMEC/FLOAT(NSTEP)
        FXM(NK)=OMG(NK)*DXSQ/G
      ENDDO
C----------------------------------------------------------------------     
C     
C***  DO AN UPSTREAM/FORWARD-IN-TIME ADVECTION OF THETA, QV
C     
      DO 525 NTC=1,NSTEP
C     
C***  ASSIGN THETA AND Q VALUES AT THE TOP AND BOTTOM OF EACH LAYER BASED
C***  SIGN OF OMEGA.
C     
      DO NK=1,LTOP
        THFXIN(NK)=0.
        THFXOUT(NK)=0.
        QFXIN(NK)=0.
        QFXOUT(NK)=0.
      ENDDO
C
      DO NK=2,LTOP
        IF(OMG(NK).LE.0.)THEN
          THFXIN(NK)=-FXM(NK)*THPA(NK-1)
          QFXIN(NK)=-FXM(NK)*QPA(NK-1)
          THFXOUT(NK-1)=THFXOUT(NK-1)+THFXIN(NK)
          QFXOUT(NK-1)=QFXOUT(NK-1)+QFXIN(NK)
        ELSE
          THFXOUT(NK)=FXM(NK)*THPA(NK)
          QFXOUT(NK)=FXM(NK)*QPA(NK)
          THFXIN(NK-1)=THFXIN(NK-1)+THFXOUT(NK)
          QFXIN(NK-1)=QFXIN(NK-1)+QFXOUT(NK)
        ENDIF
      ENDDO
C     
C***  UPDATE THE THETA AND QV VALUES AT EACH LEVEL
C     
      DO NK=1,LTOP
        THPA(NK)=THPA(NK)+(THFXIN(NK)+UDR(NK)*THTAU(NK)+DDR(NK)*
     1           THTAD(NK)-THFXOUT(NK)-(UER(NK)-DER(NK))*THTA0(NK))*
     2           DTIME*EMSD(NK)
        IF(NK.GE.LC.AND.NK.LE.KPBL)THEN
          QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)-
     1            QFXOUT(NK)-UER(NK)*QUER(NK)+DER(NK)*Q0(NK))*
     2            DTIME*EMSD(NK)
        ELSE
          QPA(NK)=QPA(NK)+(QFXIN(NK)+UDR(NK)*QDT(NK)+DDR(NK)*QD(NK)-
     1            QFXOUT(NK)-(UER(NK)-DER(NK))*Q0(NK))*DTIME*EMSD(NK)
        ENDIF  
      ENDDO
C
  525 CONTINUE
C----------------------------------------------------------------------     
      DO NK=1,LTOP
        THTAG(NK)=THPA(NK)
        QG(NK)=QPA(NK)
      ENDDO
C----------------------------------------------------------------------     
C     
C***  CHECK TO SEE IF MIXING RATIO DIPS BELOW ZERO ANYWHERE.  IF SO, 
C***  BORROW MOISTURE FROM ADJACENT LAYERS TO BRING IT BACK UP ABOVE 
C***  ZERO.
C     
      DO 530 NK=1,LTOP
C
      IF(QG(NK).LT.0.)THEN
C
        IF(NK.EQ.1)THEN                            
          PRINT *,'!!!!! PROBLEM WITH KF SCHEME:  '
          PRINT *,'QG = 0 AT THE SURFACE!!!!!!!'   
          STOP 'QG'                                
        ENDIF                                      
C
        NK1=NK+1
        IF(NK.EQ.LTOP)NK1=KLCL
        TMA=QG(NK1)*EMS(NK1)
        TMB=QG(NK-1)*EMS(NK-1)
        TMM=(QG(NK)-1.E-9)*EMS(NK  )
        BCOEFF=-TMM/((TMA*TMA)/TMB+TMB)
        ACOEFF=BCOEFF*TMA/TMB
        TMB=TMB*(1.-BCOEFF)
        TMA=TMA*(1.-ACOEFF)
C
        IF(NK.EQ.LTOP)THEN
          QVDIFF=(QG(NK1)-TMA*EMSD(NK1))*100./QG(NK1)
C
          IF(ABS(QVDIFF).GT.1.)THEN
            PRINT*,'!!!WARNING!!! CLOUD BASE WATER VAPOR CHANGES BY ',
     1             QVDIFF,
     2             '% WHEN MOISTURE IS BORROWED TO PREVENT NEGATIVE ',
     3             'VALUES IN KAIN-FRITSCH'
          ENDIF
C
        ENDIF
C
        QG(NK)=1.E-9
        QG(NK1)=TMA*EMSD(NK1)
        QG(NK-1)=TMB*EMSD(NK-1)
      ENDIF
C
  530 CONTINUE
C----------------------------------------------------------------------     
      TOPOMG=(UDR(LTOP)-UER(LTOP))*DP(LTOP)*EMSD(LTOP)
C
      IF(ABS(TOPOMG-OMG(LTOP)).GT.1.E-3)THEN
        WRITE(99,*)'ERR:  MASS NOT BALANCED IN KF SCHEME; TOPOMG, OMG ='
     1   ,TOPOMG,OMG(LTOP)
        WRITE(6,*)'I, J, NTSD =',I,J,NTSD
        WRITE(6,*)'ERR:  MASS NOT BALANCED IN KF SCHEME; TOPOMG, OMG ='
     1   ,TOPOMG,OMG(LTOP)
        ISTOP=1
        GO TO 600
      ENDIF
C     
C***  CONVERT THETA TO T
C     
      DO NK=1,LTOP
        EXN(NK)=(P00/P0(NK))**(0.2854*(1.-0.28*QG(NK)))
        TG(NK)=THTAG(NK)/EXN(NK)
        TVG(NK)=TG(NK)*(1.+0.608*QG(NK))
      ENDDO
C
C***  SHALLOW 
C
      IF(ISHALL.EQ.1)THEN
        GO TO 600
      ENDIF
C     
C*******************************************************************
C                                                                  *
C     COMPUTE NEW CLOUD AND CHANGE IN AVAILABLE BUOYANT ENERGY.    *
C                                                                  *
C*******************************************************************
C     
C***  THE FOLLOWING COMPUTATIONS ARE SIMILAR TO THAT FOR UPDRAFT
C     
      TMIX=0.
      QMIX=0.
C
C***  FIND THE THERMODYNAMIC CHARACTERISTICS OF THE LAYER BY
C***  MASS-WEIGHTING THE CHARACTERISTICS OF THE INDIVIDUAL MODEL
C***  LAYERS.
C
      DO NK=LC,KPBL
        TMIX=TMIX+DP(NK)*TG(NK)
        QMIX=QMIX+DP(NK)*QG(NK)  
      ENDDO
C
      TMIX=TMIX/DPTHMX
      QMIX=QMIX/DPTHMX
      ES=ALIQ*EXP((TMIX*BLIQ-CLIQ)/(TMIX-DLIQ))
      QSS=0.622*ES/(PMIX-ES)
C     
C***  REMOVE SUPERSATURATION FOR DIAGNOSTIC PURPOSES, IF NECESSARY
C     
      IF(QMIX.GT.QSS)THEN
        RL=XLV0-XLV1*TMIX
        CPM=CP*(1.+0.887*QMIX)
        DSSDT=QSS*(CLIQ-BLIQ*DLIQ)/((TMIX-DLIQ)*(TMIX-DLIQ))
        DQ=(QMIX-QSS)/(1.+RL*DSSDT/CPM)
        TMIX=TMIX+RL/CP*DQ
        QMIX=QMIX-DQ
        TLCL=TMIX
      ELSE
        QMIX=AMAX1(QMIX,0.)
        EMIX=QMIX*PMIX/(0.622+QMIX)
        ASTRT=1.E-3
        BINC=0.075
        A1=EMIX/ALIQ
        TP=(A1-ASTRT)/BINC
        INDLU=INT(TP)+1
        VALUE=(INDLU-1)*BINC+ASTRT
        AINTRP=(A1-VALUE)/BINC
        TLOG=AINTRP*ALU(INDLU+1)+(1-AINTRP)*ALU(INDLU)
        TDPT=(CLIQ-DLIQ*TLOG)/(BLIQ-TLOG)
        TLCL=TDPT-(.212+1.571E-3*(TDPT-T00)-4.36E-4*(TMIX-T00))*(TMIX-
     1       TDPT)
        TLCL=AMIN1(TLCL,TMIX)
      ENDIF
C
      TVLCL=TLCL*(1.+0.608*QMIX)
      ZLCL=ZMIX+(TLCL-TMIX)/GDRY
C
      DO NK=LC,KL
        KLCL=NK
        IF(ZLCL.LE.Z00(NK))GO TO 550
      ENDDO
C
      PRINT*,'PROBLEM...LCL cannot be found after conv adjustment!!'
      STOP 'KF BAD ADJUST'
C
  550 CONTINUE
      K=KLCL-1
      DLP=(ZLCL-Z00(K))/(Z00(KLCL)-Z00(K))
C     
C***  ESTIMATE ENVIRONMENTAL TEMPERATURE AND MIXING RATIO AT THE LCL
C     
      TENV=TG(K)+(TG(KLCL)-TG(K))*DLP
      QENV=QG(K)+(QG(KLCL)-QG(K))*DLP
      TVEN=TENV*(1.+0.608*QENV)
      PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
      THETEU(K)=TMIX*(1.E5/PMIX)**(0.2854*(1.-0.28*QMIX))*
     1          EXP((3374.6525/TLCL-2.5403)*QMIX*(1.+0.81*QMIX))
      ES=ALIQ*EXP((TENV*BLIQ-CLIQ)/(TENV-DLIQ))
      QESE=0.622*ES/(PLCL-ES)
      THTESG(K)=TENV*(1.E5/PLCL)**(0.2854*(1.-0.28*QESE))*
     1          EXP((3374.6525/TENV-2.5403)*QESE*(1.+0.81*QESE))
C----------------------------------------------------------------------     
C     
C***  COMPUTE ADJUSTED ABE(ABEG)
C     
      ABEG=0.
      THTUDL=THETEU(K)
      THEU_ADJ=THTUDL
C
      DO 560 NK=K,LTOPM1
      NK1=NK+1
      THETEU(NK1)=THETEU(NK)
C
      CALL TP_CAPE(P0(NK1),THETEU(NK1),TGU(NK1),QGU(NK1))
C
      TVQU(NK1)=TGU(NK1)*(1.+0.608*QGU(NK1)-QLIQ(NK1)-QICE(NK1))
C
      IF(NK.EQ.K)THEN
        DZZ=Z00(KLCL)-ZLCL
        DILBE=((TVLCL+TVQU(NK1))/(TVEN+TVG(NK1))-1.)*DZZ
      ELSE
        DZZ=DZA(NK)
        DILBE=((TVQU(NK)+TVQU(NK1))/(TVG(NK)+TVG(NK1))-1.)*DZZ
      ENDIF
C
      IF(DILBE.GT.0.)ABEG=ABEG+DILBE*G
C
C*** DILUTE BY ENTRAINMENT BY THE RATE AS ORIGINAL UPDRAFT.
C
      CALL ENVIRTHT(P0(NK1),TG(NK1),QG(NK1),THTEEG(NK1),0.,
     1              RL,ALIQ,BLIQ,CLIQ,DLIQ,AICE,BICE,CICE,DICE)
      THETEU(NK1)=
     1     THETEU(NK1)*DDILFRC(NK1)+THTEEG(NK1)*(1.-DDILFRC(NK1))
  560 CONTINUE
C----------------------------------------------------------------------     
C     
C***  ASSUME AT LEAST 90% OF CAPE (ABE) IS REMOVED BY CONVECTION DURING
C***  THE PERIOD TIMEC
C     
      IF(NOITR.EQ.1)THEN
        GO TO 600
      ENDIF
C
      DABE=AMAX1(ABE-ABEG,0.1*ABE)
      FABE=ABEG/(ABE+1.E-8)
C
      IF(FABE.GT.1..AND.ISHALL.EQ.0)THEN
c         WRITE(98,*)'UPDRAFT/DOWNDRAFT COUPLET INCREASES CAPE AT THIS'
c         WRITE(98,*)'GRID POINT; NO CONVECTION ALLOWED!'
          GO TO 900
        ENDIF
C
      IF(NCOUNT.NE.1)THEN
        DFDA=(FABE-FABEOLD)/(AINC-AINCOLD)
        IF(DFDA.GT.0.)THEN
          NOITR=1
          AINC=AINCOLD
          GO TO 575
        ENDIF
      ENDIF
C
      AINCOLD=AINC
      FABEOLD=FABE
C
      IF(AINC/AINCMX.GT.0.999.AND.FABE.GT.1.05-STAB)THEN
        GO TO 600
      ENDIF
C
      IF(FABE.LE.1.05-STAB.AND.FABE.GE.0.95-STAB)GO TO 600
C
      IF(NCOUNT.GT.10)THEN
        GO TO 600
      ENDIF
C----------------------------------------------------------------------     
C     
C***  IF MORE THAN 10% OF THE ORIGINAL CAPE REMAINS, INCREASE THE 
C***  CONVECTIVE MASS FLUX BY THE FACTOR AINC
C     
      IF(FABE.EQ.0.)THEN
        AINC=AINC*0.5
      ELSE
        AINC=AINC*STAB*ABE/(DABE+1.E-8)
      ENDIF
C
  575 AINC=AMIN1(AINCMX,AINC)
C
C***  IF AINC BECOMES VERY SMALL, EFFECTS OF CONVECTION 
C***  WILL BE MINIMAL SO JUST IGNORE IT.
C
      IF(AINC.LT.0.05)THEN
        GO TO 900         
      ENDIF
C
C     AINC=AMAX1(AINC,0.05)   
      TDER=TDER2*AINC
      PPTFLX=PPTFL2*AINC
C
      DO NK=1,LTOP
        UMF(NK)=UMF2(NK)*AINC
        DMF(NK)=DMF2(NK)*AINC
        DETLQ(NK)=DETLQ2(NK)*AINC
        DETIC(NK)=DETIC2(NK)*AINC
        UDR(NK)=UDR2(NK)*AINC
        UER(NK)=UER2(NK)*AINC
        DER(NK)=DER2(NK)*AINC
        DDR(NK)=DDR2(NK)*AINC
      ENDDO
C***  
C***  GO BACK UP FOR ANOTHER ITERATION...
C***  
      GO TO 500
  600 CONTINUE
C----------------------------------------------------------------------     
C     
C***  COMPUTE HYDROMETEOR TENDENCIES AS IS DONE FOR T, QV
C     
C***  FRC2 IS THE FRACTION OF TOTAL CONDENSATE  
C***  GENERATED THAT GOES INTO PRECIPITIATION  
C
      IF(CPR.GT.0.)THEN 
        FRC2=PPTFLX/(CPR*AINC)    
      ELSE
        FRC2=0.
      ENDIF
C
      DO NK=1,LTOP
        QLPA(NK)=QL0(NK)
        QIPA(NK)=QI0(NK)
        QRPA(NK)=QR0(NK)
        QSPA(NK)=QS0(NK)
        RAINFB(NK)=PPTLIQ(NK)*AINC*FBFRC*FRC2   
        SNOWFB(NK)=PPTICE(NK)*AINC*FBFRC*FRC2  
      ENDDO
C
C----------------------------------------------------------------------     
      DO 625 NTC=1,NSTEP
C     
C***  ASSIGN HYDROMETEORS CONCENTRATIONS AT THE TOP AND BOTTOM OF EACH 
C***  LAYER BASED ON THE SIGN OF OMEGA
C     
      DO NK=1,LTOP
        QLFXIN(NK)=0.
        QLFXOUT(NK)=0.
        QIFXIN(NK)=0.
        QIFXOUT(NK)=0.
        QRFXIN(NK)=0.
        QRFXOUT(NK)=0.
        QSFXIN(NK)=0.
        QSFXOUT(NK)=0.
      ENDDO
C
      DO NK=2,LTOP
        IF(OMG(NK).LE.0.)THEN
          QLFXIN(NK)=-FXM(NK)*QLPA(NK-1)
          QIFXIN(NK)=-FXM(NK)*QIPA(NK-1)
          QRFXIN(NK)=-FXM(NK)*QRPA(NK-1)
          QSFXIN(NK)=-FXM(NK)*QSPA(NK-1)
          QLFXOUT(NK-1)=QLFXOUT(NK-1)+QLFXIN(NK)
          QIFXOUT(NK-1)=QIFXOUT(NK-1)+QIFXIN(NK)
          QRFXOUT(NK-1)=QRFXOUT(NK-1)+QRFXIN(NK)
          QSFXOUT(NK-1)=QSFXOUT(NK-1)+QSFXIN(NK)
        ELSE
          QLFXOUT(NK)=FXM(NK)*QLPA(NK)
          QIFXOUT(NK)=FXM(NK)*QIPA(NK)
          QRFXOUT(NK)=FXM(NK)*QRPA(NK)
          QSFXOUT(NK)=FXM(NK)*QSPA(NK)
          QLFXIN(NK-1)=QLFXIN(NK-1)+QLFXOUT(NK)
          QIFXIN(NK-1)=QIFXIN(NK-1)+QIFXOUT(NK)
          QRFXIN(NK-1)=QRFXIN(NK-1)+QRFXOUT(NK)
          QSFXIN(NK-1)=QSFXIN(NK-1)+QSFXOUT(NK)
        ENDIF
      ENDDO
C----------------------------------------------------------------------     
C     
C***  UPDATE THE HYDROMETEOR CONCENTRATION VALUES AT EACH LEVEL
C     
      DO 620 NK=1,LTOP
        QLPA(NK)=QLPA(NK)+(QLFXIN(NK)+DETLQ(NK)-QLFXOUT(NK))*DTIME*
     1           EMSD(NK)
        QIPA(NK)=QIPA(NK)+(QIFXIN(NK)+DETIC(NK)-QIFXOUT(NK))*DTIME*
     1           EMSD(NK)
C
C***  REVERSE THE MOD THAT ALLOWS FEEDBACK OF RAIN/SNOW THAT ORIGINATES 
C***  IN DETRAINED AIR.
C
c     QRPA(NK)=QRPA(NK)+(QRFXIN(NK)+QLQOUT(NK)*UDR(NK)-QRFXOUT(NK)
      QRPA(NK)=QRPA(NK)+(QRFXIN(NK)-QRFXOUT(NK)
     1        +RAINFB(NK))*DTIME*EMSD(NK)
C
C***  REVERSE THE MOD THAT ALLOWS FEEDBACK OF RAIN/SNOW THAT ORIGINATES 
C***  IN DETRAINED AIR.
C
c     QSPA(NK)=QSPA(NK)+(QSFXIN(NK)+QICOUT(NK)*UDR(NK)-QSFXOUT(NK)
      QSPA(NK)=QSPA(NK)+(QSFXIN(NK)-QSFXOUT(NK)
     1        +SNOWFB(NK))*DTIME*EMSD(NK)
C
  620 CONTINUE
  625 CONTINUE
C----------------------------------------------------------------------     
C----------------------------------------------------------------------     
      DO NK=1,LTOP
        QLG(NK)=QLPA(NK)
        QIG(NK)=QIPA(NK)
        QRG(NK)=QRPA(NK)
        QSG(NK)=QSPA(NK)
      ENDDO
C----------------------------------------------------------------------     
C
C***  CLEAN THINGS UP, CALCULATE CONVECTIVE FEEDBACK TENDENCIES FOR 
C***  THIS GRID POINT.
C     
C***  SEND FINAL PARAMETERIZED VALUES TO OUTPUT FILES.
C     
c 640 IF(IPRNT)THEN 
c       IF(ISTOP.EQ.1)THEN
c         WRITE(99,*)
c         WRITE(99,*)'At T(H), I, J =',FLOAT(NTSD)*72./3600.,I,J
c         WRITE(99,*)'P(LC), DTP, WKL, WKLCL =',P0(LC)/100.,
c    1                TLCL+DTLCL+DTRH-TENV,WKL,WKLCL
c         WRITE(99,*)'TLCL, DTLCL, DTRH, TENV =',TLCL,DTLCL,
c    1                DTRH,TENV   
c         WRITE(99,1025)KLCL,ZLCL,DTLCL,LTOP,P0(LTOP),IFLAG,
c    1    TMIX-T00,PMIX,QMIX,ABE
c         WRITE(99,1030)P0(LET)/100.,P0(LTOP)/100.,VMFLCL,PLCL/100.,
c    1    WLCL,CLDHGT(LC)
c         WRITE(99,1035)PEF,PEFCBH,LC,LET,WKL,VWS 
c         WRITE(99,*)'PRECIP EFFICIENCY =',PEFF 
c         WRITE(99,1080)LFS,LDB,LDT,TIMEC,TADVEC,NSTEP,NCOUNT,FABE,AINC
c       ENDIF
C     
C***  IF CALCULATIONS ABOVE SHOW AN ERROR IN THE MASS BUDGET, PRINT OUT 
C***  FOR USE LATER FOR DIAGNOSTIC PURPOSES, THEN ABORT RUN.
C     
c       IF(ISHALL.EQ.0.OR.ISTOP.EQ.1)THEN
c         WRITE(98)I,J,IYR,IMO,IDY,IHR,IMN,KL
c         WRITE(98) P0,T0,Q0,U0,V0,W0,DP,TKE
c         IF(ISTOP.EQ.1)THEN
c           STOP 'KAIN-FRITSCH'
c         ENDIF
c       ENDIF
c     ENDIF
C----------------------------------------------------------------------     
C
      CNDTNF=(1.-EQFRC(LFS))*(QLIQ(LFS)+QICE(LFS))*DMF(LFS)
C***
C***  EVALUATE MOISTURE BUDGET
C***
      QINIT=0.
      QFNL=0.
      DPT=0.
C
      DO NK=1,LTOP
        DPT=DPT+DP(NK)
        QINIT=QINIT+Q0(NK)*EMS(NK)
        QFNL=QFNL+QG(NK)*EMS(NK)
        QFNL=QFNL+(QLG(NK)+QIG(NK)+QRG(NK)+QSG(NK))*EMS(NK)
      ENDDO
C
      QFNL=QFNL+PPTFLX*TIMEC*(1.-FBFRC) 
c     QFNL=QFNL+PPTFLX*TIMEC          
C
      ERR2=(QFNL-QINIT)*100./QINIT
C
      IF(ABS(ERR2).GT.0.05.AND.ISTOP.EQ.0)THEN 
        WRITE(99,*)'!!!!!!!! MOISTURE BUDGET ERROR IN KFPARA !!!'
        WRITE(99,666)QINIT,QFNL,ERR2
  666   FORMAT(' QINIT=',E12.5,' QFNL=',E12.5,' ERR2=',E12.5)
c       IPRNT=.TRUE.
c       ISTOP=1
c       GO TO 640
        STOP 'QVERR'
      ENDIF
C
      RELERR=ERR2*QINIT/(PPTFLX*TIMEC+1.E-10)
C----------------------------------------------------------------------     
C     
C***  FEEDBACK TO RESOLVABLE SCALE TENDENCIES.
C     
C***  IF THE ADVECTIVE TIME PERIOD (TADVEC) IS LESS THAN SPECIFIED MINIMUM
C***  TIMEC, ALLOW FEEDBACK TO OCCUR ONLY DURING TADVEC.
C     
      IF(TADVEC.LT.TIMEC)NIC=NINT(TADVEC/(0.5*DT2))
      NCA(I,J)=NIC
C
      IF(ISHALL.EQ.1)THEN
        TIMEC=2400.
        NCA(I,J)=NCLDCK
        NSHALL=NSHALL+1
      ENDIF 
C
      NKFP=NKFP+1
C
C----------------------------------------------------------------------     
      DO 675 K=1,KX-NLEV
      NK=KX-K+1-NLEV  
C
C***  EVAPORATE OR SUBLIMATE ALL HYDROMETEORS SO THAT THE
C***  FEEDBACKS ARE TEMPERATURE AND WATER VAPOR TENDENCIES.  THIS MAY
C***  CREATE SUPERSATURATED VALUES OF TG, BUT THESE WILL BE REMOVED BY
C***  NORMAL SUPERSATURATION-REMOVAL MECHANISMS.
C
      QLG(K)=QLG(K)+QIG(K)+QRG(K)+QSG(K)
C
      IF(ISHALL.EQ.1)THEN
        RL=XLV0-XLV1*TG(K)            
        QG(K)=QG(K)+QLG(K)
        TG(K)=TG(K)-RL*QLG(K)/CP
        QLG(K)=0.
      ENDIF
C
      DTDT(I,J,NK)=(TG(K)-T0(K))/TIMEC
      DQDT(I,J,NK)=(QG(K)-Q0(K))/TIMEC
      DQCDT(I,J,NK)=(QLG(K)-QL0(K))/TIMEC  
C
      IF(ISHALL.EQ.0)THEN
        NCAD(I,J)=NCA(I,J)
        PSRC(I,J)=PMIX0/100.
        PCLB(I,J)=PLCL0/100.
        EMST=DPTHMX*DXSQ/G
        UMFB(I,J)=100.*VMFLCL*TIMEC*AINC/EMST
        CIN(I,J)=CIN1
      ENDIF
C
  675 CONTINUE
C
C--------------SAVE CLOUD TOP AND BOTTOM FOR RADIATION------------------
C
      LTOP=KL-LTOP+1
      LBOT=KL-LCL+1
      HTOP(I,J)=MIN(FLOAT(LTOP),HTOP(I,J))
      HBOT(I,J)=MAX(FLOAT(LBOT),HBOT(I,J))
      CNVTOP(I,J)=MIN(FLOAT(LTOP),CNVTOP(I,J))
      CNVBOT(I,J)=MAX(FLOAT(LBOT),CNVBOT(I,J))
C
      RAINCV(I,J)=.1*.5*DT2*PPTFLX*(1.-FBFRC)/DXSQ    
      RNC=RAINCV(I,J)*NIC
      NCCNT=NCCNT+1
  900 CONTINUE
C----------------------------------------------------------------------     
C***  END OF MARCH THROUGH THE POINTS
C----------------------------------------------------------------------     
c      WRITE(*)'AT NTSD =',NTSD,',NO. OF KF POINTS ACTIVATED =',
c    1             NCCNT
C----------------------------------------------------------------------     
      RETURN
      END
C**********************************************************************
      SUBROUTINE TP_CAPE(P,THES,TU,QU)
C
      PARAMETER(KFNT=250,KFNP=220)
C
      COMMON/KFLUT/ TTAB(KFNT,KFNP),QSTAB(KFNT,KFNP),THE0K(KFNP),
     1              ALU(200),RDPR,RDTHK,PTOP
C----------------------------------------------------------------------     
C***  SCALING PRESSURE & TT TABLE INDEX
C----------------------------------------------------------------------     
      TP=(P-PTOP)*RDPR
      QQ=TP-AINT(TP)
      IPTB=INT(TP)+1
C
C----------------------------------------------------------------------     
C***  BASE AND SCALING FACTOR FOR THE
C----------------------------------------------------------------------     
C
C----------------------------------------------------------------------     
C***  SCALING THE & TT TABLE INDEX
C----------------------------------------------------------------------     
      BTH=(THE0K(IPTB+1)-THE0K(IPTB))*QQ+THE0K(IPTB)
      TTH=(THES-BTH)*RDTHK
      PP =TTH-AINT(TTH)
      ITHTB=INT(TTH)+1
C
      T00=TTAB(ITHTB  ,IPTB  )
      T10=TTAB(ITHTB+1,IPTB  )
      T01=TTAB(ITHTB  ,IPTB+1)
      T11=TTAB(ITHTB+1,IPTB+1)
C
      Q00=QSTAB(ITHTB  ,IPTB  )
      Q10=QSTAB(ITHTB+1,IPTB  )
      Q01=QSTAB(ITHTB  ,IPTB+1)
      Q11=QSTAB(ITHTB+1,IPTB+1)
C
C----------------------------------------------------------------------     
C***  PARCEL TEMPERATURE
C----------------------------------------------------------------------     
C
      TU=(T00+(T10-T00)*PP+(T01-T00)*QQ
     1       +(T00-T10-T01+T11)*PP*QQ)
C
      QU=(Q00+(Q10-Q00)*PP+(Q01-Q00)*QQ
     1       +(Q00-Q10-Q01+Q11)*PP*QQ)
C----------------------------------------------------------------------     
      RETURN
      END

