          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"
C                   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"
C                                                         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),P0E(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),OMGV(KXP1)
C
C-------------variables for MOMENTUM FLUXES ------>>>>>>>>>>>>>>
C
                             R E A L
     1 UPA(KX),VPA(KX),UFXIN(KX),UFXOUT(KX),VFXIN(KX),VFXOUT(KX)
     2,UUP(KX),UDW(KX),VUP(KX),VDW(KX),UG(KX),VG(KX),DMFV(KX),DDRV(KX)
     3,UMFV(KX),UDRV(KX),DERV(KX),UERV(KX),DOMGVDP(KX),FXMV(KX)
C
C----------------------------------------------------<<<<<<<<<<<<<<<<
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)*VTM(I+IHE(J),J,NK)+U(I+IHW(J),J,NK)*
     1 VTM(I+IHW(J),J,NK)+U(I,J+1,NK)*VTM(I,J+1,NK)+
     2 U(I,J-1,NK)*VTM(I,J-1,NK))/SUMV
C
      V0(K)=(V(I+IHE(J),J,NK)*VTM(I+IHE(J),J,NK)+V(I+IHW(J),J,NK)*
     1 VTM(I+IHW(J),J,NK)+V(I,J+1,NK)*VTM(I,J+1,NK)+
     2 V(I,J-1,NK)*VTM(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
      IF(WKL.LT.0.0001)THEN
        DTLCL=0.
      ELSE
        DTLCL=4.64*WKL**0.33
      ENDIF
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
      DTTOT = DTLCL+DTRH
      IF(DTTOT.GT.1.E-4)THEN
        GDT=2.*G*DTTOT*500./TVEN
        WLCL=1.+0.5*SQRT(GDT)
        WLCL = AMIN1(WLCL,3.)
      ELSE
        WLCL=1.
      ENDIF
      PLCL=P0(K)+(P0(KLCL)-P0(K))*DLP
      WTW=WLCL*WLCL
      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=0.01*DXSQ
      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
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(WTW.LT.1.E-3)THEN
        GO TO 275
      ELSE
        WU(NK1)=SQRT(WTW)
      ENDIF
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
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
C
C
        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)
      TVDIFF = ABS(TU10-TVQU(NK1))
      IF(TVDIFF.LT.1.e-3)THEN
C
C
        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.
C
C
        GO TO 175
      ELSEIF(EQFRC(NK1).EQ.0.)THEN
        EE2=0.
        UD2=1.
C
C
        GO TO 175
      ELSE
C
C*** SUBROUTINE PROF5 INTEGRATES OVER THE GAUSSIAN DIST TO DETERMINE THE
C*** FRACTIONAL ENTRAINMENT AND DETRAINMENT RATES
C
C
        CALL PROF5(EQFRC(NK1),EE2,UD2)
C
      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.
C---------------------------------------------->>>>>>>>>>>>>>>>>>>>>>
        UMF(NK)=0.
C-----------------------------------------------<<<<<<<<<<<<<<<<<<<<<
        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.
C070114------------------------------->>>>>>>>>>>>>>>>>>>>
	UMFV(NK)=0.
	UERV(NK)=0.
C070114-------------------------------<<<<<<<<<<<<<<<<<<<<
      ENDIF
C
C070114------------------------------->>>>>>>>>>>>>>>>>>>>
      UDRV(NK)=0.
C070114-------------------------------<<<<<<<<<<<<<<<<<<<<
      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
C070114----------------------------------->>>>>>>>>>>>>>>>>>>>>>>>
        UMFV(NK)=0.
        UDRV(NK)=0.
        UERV(NK)=0.
C070114-----------------------------------<<<<<<<<<<<<<<<<<<<<<<<<
        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)
C----------ZERANDO ACIMA DO TOPO---------->>>>>>>>>>>>>>>>>>>
        UPA(NK)=U0(NK)
	VPA(NK)=V0(NK)
        UUP(NK)=0.
	VUP(NK)=0.
        UG(NK)=U0(NK)
	VG(NK)=V0(NK)
	FXM(NK)=0.
C----------------------------------------<<<<<<<<<<<<<<<<<<<
        QLG(NK)=0.
        QIG(NK)=0.
        QRG(NK)=0.
        QSG(NK)=0.
        OMG(NK)=0.
      ENDDO
C
      OMG(KXP1)=0.
      OMGV(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
C------DEF UDW VDW NO LFS------------->>>>>>>>>>>>>>>>>>>>>>>
C
      UDW(LFS)=U0(LFS)
      VDW(LFS)=V0(LFS)
C
C-------------------------------------<<<<<<<<<<<<<<<<<<<<<<<
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------DEF UDW VDW NO KSTART-------------->>>>>>>>>>>>>>>>>>>>>>>>>
C
      UDW(ND)=U0(KSTART)
      VDW(ND)=V0(KSTART)
C
C-----------------------------------------<<<<<<<<<<<<<<<<<<<<<<<<<
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.
C
C-------DEF UDW VDW DA SUP AO TOPO-----NO-DOWNDRAFT->>>>>>>>>>>>>>>>>
C
          UDW(NDK)=0.
          VDW(NDK)=0.
C
C-------------------------------------<<<<<<<<<<<<<<<<<<<<<<<<<
        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
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.
C
C------DEF UDW VDW DA SUP AO LDB-------->>>>>>>>>>>>>>>>>>>>>
          UDW(NK)=0.
          VDW(NK)=0.
C---------------------------------------<<<<<<<<<<<<<<<<<<<<<	  
        ENDDO
      ENDIF
Cc
      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.
C
C-------DEF UDW VDW ACIMA DO LFS------->>>>>>>>>>>>>>>>>>>>>>>>
C
        UDW(NK)=0.
        VDW(NK)=0.
C
C--------------------------------------<<<<<<<<<<<<<<<<<<<<<<<<
      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.1.e-3)THEN
          AINCM1=EMS(NK)/((UER(NK)-DER(NK))*TIMEC)
          AINCMX=AMIN1(AINCMX,AINCM1)
        ENDIF
      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
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)
          ABSOMG = ABS(OMG(NK))
          ABSOMGTC = ABSOMG*TIMEC
          FRDP = 0.75*DP(NK-1)
          IF(ABSOMGTC.GT.FRDP)THEN
            DTT1 = FRDP/ABSOMG
            DTT=AMIN1(DTT,DTT1)
          ENDIF
        ENDIF
      ENDDO
C
C=====================================================>>>>>>>>>>>>>>>>
C
C070114---zero momentum fluxes--->>>>>>>>>>
C
       DO NK=1,KX
          UMFV(NK)=0.
          UERV(NK)=0.
          UDRV(NK)=0.
          DMFV(NK)=0.
          DERV(NK)=0.
          DDRV(NK)=0.
       ENDDO
C070114-------------------------------------<<<<<<<<<
C
C*********************************************************************
C                                                                    *
C         CALCULATE UDW & VDW within cloud layer --> LFS-1,KSTART,-1      *
C                                                                    *
C*********************************************************************
C
C IF NO FLUX AT KSTART(CLOUD BASE) ---> NO DRAWDRAFT - UDW VDW -> ZERO
C
      IF(ABS(DMF(KSTART)).GT.1E-3)THEN
C
       DO ND=LFS-1,KSTART,-1
          ND1=ND+1
C070115--------------------------------->>>>>>>>>>>
        IF(ND.EQ.(LFS-1))THEN
	  DMFV(ND1)=DMF(LFS)
          DERV(ND1)=DER(LFS)
          DDRV(ND1)=DDR(LFS)
	ENDIF
	  DMFV(ND)=DMF(ND)
          DERV(ND)=DER(ND)
          DDRV(ND)=DDR(ND)
C
C        IF(ABS(DMF(ND)).GT.1E-3)THEN
C        UDW(ND)=(UDW(ND1)*DMF(ND1)+0.5*(U0(ND)+U0(ND1))*
C     1          DER(ND))/DMF(ND)
CC
C        VDW(ND)=(VDW(ND1)*DMF(ND1)+0.5*(V0(ND)+V0(ND1))*
C     1          DER(ND))/DMF(ND)
        IF(ABS(DMFV(ND)).GT.1E-3)THEN
        UDW(ND)=(UDW(ND1)*DMFV(ND1)+0.5*(U0(ND)+U0(ND1))*
     1          DERV(ND))/DMFV(ND)
C
        VDW(ND)=(VDW(ND1)*DMFV(ND1)+0.5*(V0(ND)+V0(ND1))*
     1          DERV(ND))/DMFV(ND)
C070115---------------------------------<<<<<<<<<<<
        ELSE
        UDW(ND)=UDW(ND1)
	VDW(ND)=VDW(ND1)
        ENDIF
       ENDDO
C
C*********************************************************************
C                                                                    *
C           CALCULO DE UDW E VDW ABAIXO DA NUVEM --> LDT,LDB,-1      *
C                                                                    *
C*********************************************************************
C---------------PRINT-------------
C       IF(MYPE.EQ.0)THEN
C       IF(NTSD.LT.100)THEN
C       IF(J.EQ.209)THEN
C      DMFf=-10.
C      ALFA=ALOG((DMF(KSTART)-DMFf)/DMFf)/(P0(KSTART)-P0(LDB))
C      WRITE(96,*)'ND(LDT-->LDB)  DMFV   DERV   DDRV'
C      WRITE(96,*)KSTART,DMF(KSTART),DER(KSTART),DDR(KSTART)
C      DO ND=LDT,LDB,-1
C         ND1=ND+1
C      IF (ND.EQ.LDT)DMFV(ND1)=DMF(KSTART)
C        DMFV(ND)=DMFf*EXP(ALFA*(P0(ND)-P0(LDB)))+DMFf
C	 DDRV(ND)=DMFV(ND1)-DMFV(ND)
C	 DERV(ND)=0.
C      IF (ND.EQ.LDB)THEN
C        DMFV(ND)=0.
C        DDRV(ND)=DDRV(ND)+2*DMFf
C      ENDIF
C      WRITE(96,*)ND,DMFV(ND),DERV(ND),DDRV(ND)
C      ENDDO
C       ENDIF
C       ENDIF
C       ENDIF
C
C-------CALCULO DMFV - DDRV  -   DERV  ---------->>>>>>>>>>>>>>> 
C
C070120=======================================>>>>>>>>>>>>>>>>>>>>>>
C
C VALOR 10 SIMILAR A SITUACAO DO FIM DA CORRENTE ASCENDENTE
C
C070120      DMFf=-10.
C070120      IF(ABS(DMF(KSTART)).GT.ABS(DMFf))THEN
C070120         ALFA=ALOG((DMF(KSTART)-DMFf)/DMFf)/(P0(KSTART)-P0(LDB))
C
C070120       DO ND=LDT,LDB,-1
C070120          ND1=ND+1
C070120         IF (ND.EQ.LDT)THEN
C070120            DMFV(ND1)=DMF(KSTART)
C070120            DDRV(ND1)=DDR(KSTART)
C070120            DERV(ND1)=DER(KSTART)
C070120         ENDIF
C070119     DMFV(ND)=DMFf*EXP(ALFA*(P0(ND)-P0(LDB)))+DMFf
C070120          DMFV(ND)=DMFf*EXP(ALFA*(P0(ND)-P0(LDB)))-DMFf
C070114	    DDRV(ND)=DMFV(ND1)-DMFV(ND)
C070120          DDRV(ND)=DMFV(ND)-DMFV(ND1)
C070120          DERV(ND)=0.
C070120       ENDDO
C
C------------AJUSTA OS VALORES PARA LDB------------
C
C070120            DMFV(LDB)=0.
C070115----------------------------------->>>>>>>>>>
C            DDRV(LDB)=DDRV(LDB)+2*DMFf
C070119            DDRV(LDB)=0.
C070119            DDRV(LDB+1)=-DMFV(LDB+1)
C070115-----------------------------------<<<<<<<<<<
C
C-------SE EXISTE CORRECAO: LDT --> LDB
C
C070120       DO ND=LDT,LDB,-1
C070120         ND1=ND+1
C
C------CALCULO UDW VDW COM CORRECAO----------------------->>>>>>>>>
C
C070114        IF(ABS(2*DMFV(ND)+DDRV(ND)).GE.1E-3)THEN
C070114        UDW(ND)=UDW(ND1)*((2*DMFV(ND1)-DDRV(ND))/
C070114    1          (2*DMFV(ND)+DDRV(ND)))
CC
C070114        VDW(ND)=VDW(ND1)*((2*DMFV(ND1)-DDRV(ND))/
C070114     1          (2*DMFV(ND)+DDRV(ND)))
C070120        IF(ABS(2*DMFV(ND)-DDRV(ND)).GE.1E-3)THEN
C070120        UDW(ND)=UDW(ND1)*((2*DMFV(ND1)+DDRV(ND))/
C070120     1          (2*DMFV(ND)-DDRV(ND)))
C
C070120        VDW(ND)=VDW(ND1)*((2*DMFV(ND1)+DDRV(ND))/
C070120     1          (2*DMFV(ND)-DDRV(ND)))
C070120        ELSE
C070120	UDW(ND)=UDW(ND1)
C070120	VDW(ND)=VDW(ND1)
C070120        ENDIF
C
C070120       ENDDO
C
C*********************************************************<<<<<<<<	 
C
C070120      ELSE
C
C-------SE NAO EXISTE A CORRECAO:  LDT--> LDB
C     
C-----------------------------------------------------<<<<<<<<<<<<<<
C070115------------------------->>>>>>>>>>
       DO ND=LDT,LDB,-1
          ND1=ND+1
C
         IF (ND.EQ.LDT)THEN
            DMFV(ND1)=DMF(KSTART)
            DDRV(ND1)=DDR(KSTART)
            DERV(ND1)=DER(KSTART)
         ENDIF
C  
          DMFV(ND)=DMF(ND)
          DDRV(ND)=DDR(ND)
          DERV(ND)=DER(ND)
C
	 IF(ND.EQ.LDB)THEN
	  DDRV(ND)=DMFV(ND)-DMFV(ND1)
	  DMFV(ND)=0.
         ENDIF
       ENDDO
       DO ND=LDT,LDB,-1
          ND1=ND+1
C070115-------------------------<<<<<<<<<<<
C
C------------CALCULO UDW VDW SEM CORRECAO------------------>>>>>>
C
C070114       IF(ABS(2*DMFV(ND)+DDRV(ND)).GE.1E-3)THEN
C070114        UDW(ND)=UDW(ND1)*((2*DMFV(ND1)-DDRV(ND))/
C070114     1          (2*DMFV(ND)+DDRV(ND)))
C070114
C070114        VDW(ND)=VDW(ND1)*((2*DMFV(ND1)-DDRV(ND))/
C070114     1          (2*DMFV(ND)+DDRV(ND)))
       IF(ABS(2*DMFV(ND)-DDRV(ND)).GE.1E-3)THEN
        UDW(ND)=UDW(ND1)*((2*DMFV(ND1)+DDRV(ND))/
     1          (2*DMFV(ND)-DDRV(ND)))
C
        VDW(ND)=VDW(ND1)*((2*DMFV(ND1)+DDRV(ND))/
     1          (2*DMFV(ND)-DDRV(ND)))
       ELSE
	UDW(ND)=UDW(ND1)
	VDW(ND)=VDW(ND1)
C
       ENDIF
       ENDDO
C070120      ENDIF
C     
C      DO ND=LDT,LDB,-1
C
C      IF(mype.eq.0)then
C      IF(NTSD.EQ.10)THEN      
C      IF(NTSD.LT.100)THEN
C      IF(J.EQ.209)THEN
C      WRITE(96,*)'LDT=',LDT,' LDB=',LDB,' LTOP=',LTOP,
C     1                ' LET=',LET,'  I=',I,' J=',J
C      ENDIF
C      ENDIF
C      ENDIF
C
C      ENDDO
C070120=================================================<<<<<<<<<<<<<<
C
      ELSE
C
C-----ZERA TODA CAMADA LFS --> LDB - IF NO FLUX AT KSTART(CLOUD BASE)
C
	  DO ND=LFS,LDB,-1
	     UDW(ND)=0.
	     VDW(ND)=0.
	  ENDDO
      ENDIF
C=====================================================<<<<<<<<<<<<<<<<
C
      DO NK=1,LTOP
C-------------------------------->>>>>>>>>>>>>>>>>>>>>
C********************************************************************
C                                                                   *
C               AJUSTES INICIAIS PARA UUP E VUP                     * 
C                                                                   * 
C********************************************************************
        UPA(NK)=U0(NK)
	VPA(NK)=V0(NK)
C
C**INICIO ASSUMI DENTRO DA NUVEM PARCELA NAO MISTURA, ABAIXO=AMBIENTE
C
        IF(NK.LT.KSTART)THEN
          UUP(NK)=U0(NK)
          VUP(NK)=V0(NK)
C070115---------------------------->>>>>>>>>>>>>>>>>
          UMFV(NK)=UMF(NK)
          UERV(NK)=UER(NK)
          UDRV(NK)=UDR(NK)
C070115----------------------------<<<<<<<<<<<<<<<<<
        ELSE
          UUP(NK)=U0(KSTART)
          VUP(NK)=V0(KSTART)
        ENDIF
C
C--------------------------------<<<<<<<<<<<<<<<<<<<<<<
        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************ ZERO FLUXES ABOVE CLOUD TOP************************
C
      DO NK=LTOP+1,KX
         FXM(NK)=0.
         UFXIN(NK)=0.
	 UFXOUT(NK)=0.
	 VFXIN(NK)=0.
	 VFXOUT(NK)=0.
	 UPA(NK)=U0(NK)
	 VPA(NK)=V0(NK) 
      ENDDO
C----------------------------------<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
C=====================================================>>>>>>>>>>
C------------------------------------------>>>>>>>>>>>>>>>>>>>
C*************************************************************
C                                                            *
C    CALCULO DA VELOC HORIZONTAL NA CORR UP KSTART --> LTOP  *
C                                                            *
C*************************************************************
C
C IF NO FLUX AT KSTART(CLOUD BASE) ---> NO UPDRAFT - UUP VUP -> ZERO
C
      IF (ABS(UMF(KSTART)).GT.1E-3)THEN
C
C070115---------------------------------------->>>>>>>>>>>>>
       DO NK=KSTART,LET-1
          NK1=NK+1
          UMFV(NK)=UMF(NK)
          UERV(NK)=UER(NK)
          UDRV(NK)=UDR(NK)
	IF(NK.EQ.(LET-1))THEN
	  UMFV(NK1)=UMF(LET)
          UERV(NK1)=UER(LET)
          UDRV(NK1)=UDR(LET)
	ENDIF
       ENDDO
       DO NK=KSTART,LET-1
          NK1=NK+1       
C        IF(ABS(2*UMF(NK1)+UDR(NK1)).LT.1E-3)THEN
C        UUP(NK1)=UUP(NK)
C        VUP(NK1)=VUP(NK)
C        ELSE
C        UUP(NK1)=(UUP(NK)*(2*UMF(NK)-UDR(NK1))+UER(NK1)*
C     1          (U0(NK)+U0(NK1)))/(2*UMF(NK1)+UDR(NK1))
CC
C        VUP(NK1)=(VUP(NK)*(2*UMF(NK)-UDR(NK1))+UER(NK1)*
C     1          (V0(NK)+V0(NK1)))/(2*UMF(NK1)+UDR(NK1))
        IF(ABS(2*UMFV(NK1)+UDRV(NK1)).LT.1E-3)THEN
        UUP(NK1)=UUP(NK)
        VUP(NK1)=VUP(NK)
        ELSE
        UUP(NK1)=(UUP(NK)*(2*UMFV(NK)-UDRV(NK1))+UERV(NK1)*
     1          (U0(NK)+U0(NK1)))/(2*UMFV(NK1)+UDRV(NK1))
C
        VUP(NK1)=(VUP(NK)*(2*UMFV(NK)-UDRV(NK1))+UERV(NK1)*
     1          (V0(NK)+V0(NK1)))/(2*UMFV(NK1)+UDRV(NK1))
C070115---------------------------------------<<<<<<<<<<<<<<
        ENDIF
       ENDDO
C
C070120============================================>>>>>>>>>>>
C------------ DECAIMENTO EXPONENCIAL  UMF  LET --> LTOP ------
C
C070120      UMFf=10.
C070120       IF(ABS(UMF(LET)).GT.ABS(UMFf))THEN
C070120         ALFAUP=ALOG((UMF(LET)-UMFf)/UMFf)/(P0(LTOP)-P0(LET))
C
C070120        DO NU=LET,LTOP-1
C070120         NU1=NU+1
C070120        IF(NU.EQ.LET)THEN
C070120         UMFV(NU)=UMF(LET)
C070120	 UDRV(NU)=UDR(LET)
C070120	 UERV(NU)=UER(LET)
C070120        ENDIF
C070119         UMFV(NU1)=UMFf*EXP(ALFAUP*(P0(LTOP)-P0(NU1)))
C070120         UMFV(NU1)=UMFf*EXP(ALFAUP*(P0(LTOP)-P0(NU1)))-UMFf
C070120	 UDRV(NU1)=UMFV(NU)-UMFV(NU1)
C070120	 UERV(NU1)=0.
C
C-----------AJUSTA OS VALORES PARA LTOP------------
C
C070119        IF(NU1.EQ.LTOP)THEN
C070119           UMFV(NU1)=0.
C070119           UDRV(NU1)=UDRV(NU1)+UMFf
C070119        ENDIF
C070120        ENDDO
C070119-------AJUSTA OS VALORES PARA LTOP----------->>>>
C
C070120         UMFV(LTOP)=0.
C
C070119--------------------------------------------<<<<<
C----------SE EXISTE CORRECAO ---  LET ' LTOP -----------
C
C070120        DO NU= LET,LTOP-1
C070120           NU1=NU+1
C
C-----------CALCULO DE UUP VUP COM CORRECAO----------------
C
C070120         IF(ABS(2*UMFV(NU1)+UDRV(NU1)).LT.1E-3)THEN
C070120          UUP(NU1)=UUP(NU)
C070120          VUP(NU1)=VUP(NU)
C070120         ELSE
C070120          UUP(NU1)=(UUP(NU)*(2*UMFV(NU)-UDRV(NU1))+UERV(NU1)*
C070120     1          (U0(NU)+U0(NU1)))/(2*UMFV(NU1)+UDRV(NU1))
C
C070120          VUP(NU1)=(VUP(NU)*(2*UMFV(NU)-UDRV(NU1))+UERV(NU1)*
C070120     1          (V0(NU)+V0(NU1)))/(2*UMFV(NU1)+UDRV(NU1))
C070120         ENDIF
C070120        ENDDO
C
C070120       ELSE
C
C------SE NAO EXISTE CORRECAO: LET -> LTOP -----------------
C
        DO NU= LET,LTOP-1
               NU1=NU+1
C
        UMFV(NU)=UMF(NU)
        UDRV(NU)=UDR(NU)
        UERV(NU)=UER(NU)
C070114--------------------------------->>>>>>>>>
	IF(NU.EQ.LTOP-1)THEN
	UDRV(NU1)=UMFV(NU)-UMFV(NU1)
	UMFV(NU1)=0.
	ENDIF
C070114--------------------------------->>>>>>>>>
C
C-------- CALCULO DE UUP VUP SEM CORRECAO----------------
C
         IF(ABS(2*UMFV(NU1)+UDRV(NU1)).LT.1E-3)THEN
          UUP(NU1)=UUP(NU)
          VUP(NU1)=VUP(NU)
         ELSE
          UUP(NU1)=(UUP(NU)*(2*UMFV(NU)-UDRV(NU1))+UERV(NU1)*
     1          (U0(NU)+U0(NU1)))/(2*UMFV(NU1)+UDRV(NU1))
C
          VUP(NU1)=(VUP(NU)*(2*UMFV(NU)-UDRV(NU1))+UERV(NU1)*
     1          (V0(NU)+V0(NU1)))/(2*UMFV(NU1)+UDRV(NU1))
         ENDIF
        ENDDO
C
C070121       ENDIF
C
      ELSE
C
C----------ZERA TODA CAMADA SE NENHUM FLUXO EM KSTART
C                  NAO HA CORRENTE ASCENDENTE
C
C070114        DO NK=KSTART,LTOP
        DO NK=1,LTOP
         UUP(NK)=0.
         VUP(NK)=0.
        ENDDO
      ENDIF
C
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
C---------------------------------------------->>>>>>>>>>>>>>
C
C********ZERANDO OS FLUXO NAS INTERFACES ABAIXO DO TOPO************
C
        UFXIN(NK)=0.
	UFXOUT(NK)=0.
	VFXIN(NK)=0.
	VFXOUT(NK)=0.
C----------DEF UDW VDW PARA ANULAR O EFEITO DOWNDRAFT----------
CC     	UDW(NK)=0.
CC        VDW(NK)=0.
C
C----------------------------------------------<<<<<<<<<<<<<<
        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==================================================>>>>>>>>>>>>>>>>>>>
C************************************************************
C                                                           *
C       CALCULO SEPARADO PARA VENTO DO FLUXO NA CAMADA      *
C                                                           *
C************************************************************
C-------------------------------------------------
C        CALCULO DO OMG PARA O VENTO
C-------------------------------------------------
      DO NK=1,LTOP
C-------CORRECAO PARA O FLUXO ACIMA KSTART, ABAIXO LET
        IF((NK.GE.KSTART).AND.(NK.LE.LFS))THEN
            DDRV(NK)=DDR(NK)
            DERV(NK)=DER(NK)
        ENDIF
        IF(NK.LE.LET)THEN
            UDRV(NK)=UDR(NK)
            UERV(NK)=UER(NK)
        ENDIF
C--------------------------------------------
        DOMGVDP(NK)=-(UERV(NK)-DERV(NK)-UDRV(NK)-DDRV(NK))*EMSD(NK)
          IF(NK.GT.1)THEN
           OMGV(NK)=OMGV(NK-1)-DP(NK-1)*DOMGVDP(NK-1)
          ELSE
	   OMGV(NK)=0.
	  ENDIF
        FXMV(NK)=OMGV(NK)*DXSQ/G
      ENDDO
      DO NK=LTOP+1,KX
         FXMV(NK)=0.
      ENDDO
C----------------------------------------------------
C
      DO NK=2,LTOP
        IF(OMGV(NK).LE.0.)THEN
C       NK1=NK+1
C
C-------------------------CAMADA 2 DOWN------------>>>>>>>>>>
C
	  UFXIN(NK)=-FXMV(NK)*UPA(NK-1)
	  VFXIN(NK)=-FXMV(NK)*VPA(NK-1)
CCC	  UFXIN(NK1)=-FXM(NK1)*UPA(NK1)
CCC	  VFXIN(NK1)=-FXM(NK1)*VPA(NK1)
C
C--------------------------------------------------<<<<<<<<<<
C
C-------------------------CAMADA 1 DOWN------------>>>>>>>>>>
C
	  UFXOUT(NK-1)=UFXOUT(NK-1)+UFXIN(NK)
	  VFXOUT(NK-1)=VFXOUT(NK-1)+VFXIN(NK)
CCC	  UFXOUT(NK)=-FXM(NK)*UPA(NK)
CCC	  VFXOUT(NK)=-FXM(NK)*VPA(NK)
C--------------------------------------------------<<<<<<<<<<
        ELSE
C
C-------------------------CAMADA 2 UP-------------->>>>>>>>>>
C
	  UFXOUT(NK)=FXMV(NK)*UPA(NK)
	  VFXOUT(NK)=FXMV(NK)*VPA(NK)
C
C--------------------------------------------------<<<<<<<<<<
C
C-------------------------CAMADA 1 UP-------------->>>>>>>>>>
C
	  UFXIN(NK-1)=UFXIN(NK-1)+UFXOUT(NK)
	  VFXIN(NK-1)=VFXIN(NK-1)+VFXOUT(NK)
CCC	  UFXIN(NK-1)=FXM(NK-1)*UPA(NK-1)
CCC	  VFXIN(NK-1)=FXM(NK-1)*VPA(NK-1)
C--------------------------------------------------<<<<<<<<<<
        ENDIF
      ENDDO     
C========================================================<<<<<<<<<<<<<<
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
CCCCCCCCCCCCCCCCC------------------------------->>>>>>>>>>>>>>>>>>
C****************************************************************
C                                                               * 
C              CALCULO SEPARADO DO UPDATE U E V                 *
C                                                               *
C****************************************************************
C
      DO NK=1,LTOP
C
C
C----------------UPDATE U E V---------------------->>>>>>>>>>>>>
C
	UPA(NK)=UPA(NK)+(UFXIN(NK)+UDRV(NK)*UUP(NK)+DDRV(NK)*UDW(NK)
     1          -UFXOUT(NK)-(UERV(NK)-DERV(NK))*U0(NK))*
     2          DTIME*EMSD(NK)
C
	VPA(NK)=VPA(NK)+(VFXIN(NK)+UDRV(NK)*VUP(NK)+DDRV(NK)*VDW(NK)
     1          -VFXOUT(NK)-(UERV(NK)-DERV(NK))*V0(NK))*
     2          DTIME*EMSD(NK)
C
      ENDDO
CCCCCCCCCCCCCCCCC-------------------------------<<<<<<<<<<<<<<<<<
  525 CONTINUE
C----------------------------------------------------------------------
      DO NK=1,LTOP
C------------------------------------------->>>>>>>>>>>>>>>>>>
C
        UG(NK)=UPA(NK)
        VG(NK)=VPA(NK)
C
C-------------------------------------------<<<<<<<<<<<<<<<<<<<
        THTAG(NK)=THPA(NK)
        QG(NK)=QPA(NK)
      ENDDO
C----------------------------------------------------------------------
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
        IF(ABS(AINC-AINCOLD).LT.0.0001)THEN
          NOITR=1
          AINC=AINCOLD
          GO TO 575
        ENDIF
        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
        IF(DABE.LT.1.e-4)THEN
          NOITR=1
          AINC=AINCOLD
!          CYCLE iter
        ELSE
          AINC=AINC*STAB*ABE/DABE
        ENDIF
      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
      IF(PPTFLX.GT.0.)THEN
        RELERR=ERR2*QINIT/(PPTFLX*TIMEC)
      ELSE
        RELERR=0.
      ENDIF
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
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--------------------------------------->>>>>>>>>>>>>>>>>>>>>
C*********************************************************************
C                                                                    *
C                     CALCULO DE DUDT E DVDT                         *
C                                                                    *
C*********************************************************************
C
      DUDT(I,J,NK)=(UG(K)-U0(K))/TIMEC
      DVDT(I,J,NK)=(VG(K)-V0(K))/TIMEC
C
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
      close(95)
      close(96)
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


