      SUBROUTINE SCPRP8(MC,HMNC,QC,XC,HMNF,QF,XF,ZF,TG,HMNG,QG,XG,
     1                  ZG,PG,DG,ILEV,ILEVM,ILG,ILGG,MSG,NTRAC,ADEF,
     2                  NTRACA,ISO2,ISO4,IHPO,INH3,IHNO3,IO3,ICO2,
     3                  IDUA,IDUC,ISSA,ISSC,ISULF,ITRAC,LCL,LMAG,LTM,
     4                  LNB,FCT,QSTC,PF,COND,AUTO,RU,RUA,ITRCCH,
     5                  RRL,GRAV,CPRES,EPS1,EPS2,A,B,
     6                  BUOY,TDRY,TFE,WC,TSTAR,
     7                  TAUC,DZ,TC,DC,EC,PE,PD,ALPHA,ALPHM,HCT,QCT,QLCT,
     8                  XCT,NEQP,ACHPA,XCTA,XGA,SC,RHOA,AMHVAL,
     9                  RGAS,RGASV,TVFA,GAMMA,DADZF,WRK,AFMAX,ITRPHS,
     A                  ISKIP,PHSROW,ISHALL,ISVCONV,ISVCHEM,
     B                  SO4CROW,DCT,DZO,DTO,QLWC,SCDN,IPAM) 
C-----------------------------------------------------------------------
C     * FEB 10/2015 - M.LAZARE/    New routine for GCM18:
C     *               K.VONSALZEN. - GCROW removed (wasn't being used).
C     *                            - Scavenging now allowed for SO2.
C     *                            - Bugfix to possible accuracy 
C     *                              issue in ZINTE, especially at 32-bit.
C     * JUN 16/2013 - K.VONSALZEN. PREVIOUS VERSION SCPRP7 FOR GCM17:
C     *                            - ADD IN PLA OPTION (THROUGH "IPAM").
C     * APR 23/2010 - K.VONSALZEN. PREVIOUS VERSION SCPRP6 FOR GCM15I+:
C     *                            - MODIFIED DIAGNOSTIC FIELDS.
C     *                            - ASRSO4 INCREASED FROM 0.75 TO 0.90.
C     *                            - SLIGHTLY LESS INDIRECT EFFECT (PCDNC).
C     * FEB 16/2009 - K.VONSALZEN. PREVIOUS VERSION SCPRP5 FOR GCM15H:
C     *                            - REVISED CLOUD DROPLET DENSITY
C     *                              (PCDNC) CALCULATION, CONSISTENT
C     *                              WITH LARGESCALE.
C     * NOV 25/2006 - M.LAZARE.    PREVIOUS VERSION SCPRP4 FOR GCM15G.
C     *                            - SELECTED PROMOTION OF VARIABLES TO
C     *                              64BIT AND REVISED EXPRESSIONS,
C     *                              TO WORK IN 32-BIT MODE.
C     *                            - CALLS NEW OEQUIP2 AND TINCLD3.
C     * JUN 14/2006 - M.LAZARE.    PREVIOUS VERSION SCPRP3 FOR GCM15F:
C     *                            - DATA CONSTANTS ADDED AND USED
C     *                              IN CALL TO INTRINSICS. 
C     * MAY 04/2006 - M.LAZARE/    PREVIOUS VERSION SCPRP3 FOR GCM15E:
C     *               K.VONSALZEN. - CALLS NEW TINCLD2.
C     *                              THIS REQUIRES PASSING IN EPS2.
C     * DEC 13/2005 - K.VONSALZEN/ PREVIOUS VERSION SCPRP2 FOR GCM15D.
C
C     * CALCULATES CLOUD PROPERTIES AND FLUXES IN TRANSIENT SHALLOW 
C     * CUMULUS CLOUDS.
C-----------------------------------------------------------------------
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
C     
C     * WORK FIELDS (FIRST SET ARE PROMOTED FOR CONSISTENCY
C     * WITH REST OF SHALLOW). 
C
      REAL*8,  DIMENSION(ILG,ILEV)         :: BUOY,QF,PF,QC,QSTC,TC,
     1                                        HMNF,DZ,ZG,ZF,TDRY,
     2                                        HMNC,QCT,HCT,DG,PE,PD,
     3                                        ALPHA,WC,TSTAR
      REAL*8,  DIMENSION(ILG)              :: FCT,TAUC
      REAL*8   ATL,ATH,WTHRES,DHDF,DTDF,DQDF,DQSDT
C      
      REAL,   DIMENSION(ILG,ILEV)         :: AUTO,TFE,EC,TG,
     1                                       RHOA,MC,DADZF,WRK,
     2                                       QG,PG,ALPHM,HMNG,DC,
     3                                       COND,QLCT,AFMAX,DCT,
     4                                       DZO,QLWC,SCDN,PHSROW
C
      REAL,    DIMENSION(ILG,ILEV,NTRAC)   :: XC,XF,XCT,XG,SC,RU
C
C     * CHEMISTRY-RELATED FIELDS.
C
      REAL SO4CROW(ILG,ILEV)
      REAL,    DIMENSION(ILG,ILEV,NTRACA)  :: XGA,XCTA,RUA
C
      REAL ACHPA(ILG,ILEV,NEQP)
C     
C     * SINGLE-LEVEL WORK FIELDS.
C
      REAL      AMHVAL(ILG)
C
      INTEGER    LCL(ILG), LTM(ILG),  LNB(ILG), LMAG(ILG), ISKIP(ILG),
     1        ISHALL(ILG)
      INTEGER ITRPHS(NTRAC),ITRCCH(NTRAC)
C
C     * OTHER FIELDS AND PARAMETERS.
C
      LOGICAL KCTDET, KEQPROB, KCHEM, KRIT1, KRIT2
C
C     * LOCAL CONSTANTS.
C
      PARAMETER ( AMU    =4.5     )    ! determines lateral mixing
      PARAMETER ( APHI   =4.E-04  )    ! determines cloud-top mixing
      PARAMETER ( AEPS   =0.      ) 
      PARAMETER ( KCTDET =.FALSE. )    ! switch for cloud-top detrainmen
      PARAMETER ( KEQPROB=.FALSE. )    ! switch for version EQPROB
      PARAMETER ( YCOM3L =1.E+03       )
      PARAMETER ( YEPS   =1.E-20       )
      PARAMETER ( YSMALL =1.E-03       )
      PARAMETER ( YRHOWAT=1.E+03       )
      PARAMETER ( YSMASS =32.06E-03    )
      PARAMETER ( YNA    =6.022045E+23 )
      PARAMETER ( BETA   =.5           )
      PARAMETER ( ZCRAUT = 1.          )
      PARAMETER ( YAUTEX = 4.7         )
      PARAMETER ( YAPN   = 0.02        )
      PARAMETER ( ASRSO4 = 0.90        )                                
      PARAMETER ( ASRCO2 = 0.          )
      PARAMETER ( ASRPHOB= 0.          )
      PARAMETER ( YCUTL  = -3.05E-02   )  ! Accuracy 3% for exponential
      PARAMETER ( YCUTH  =  2.96E-02   )  ! Accuracy 3% for exponential
      PARAMETER ( ZCDNMIN= 1.E6        )
C
      REAL*8 ZINTE, ARG1, ARG2, YEFPA, APA, AEXPF, APP, AG, APDF0, 
     1       ATEMPO, AETMP, SMALLZINTE
      REAL*8 YFPA, YFPB, YFPC 
C
      DATA YFPA / -2.1066  / 
      DATA YFPB / 1.0443   / 
      DATA YFPC / -0.87296 / 
C
      DATA ZERO   / 0.    /
      DATA ONE    / 1.    /
      DATA OPEM1  /1.E-1  /
      DATA OPEM2  /1.E-2  /
      DATA OPEM10 /1.E-10 /
      DATA TFHUN  /2400.  /
C
C     * INTERNAL FUNCTIONS.
C
      ZINTE(ARG1,ARG2)=(1./ARG2)*(ARG1 - 5.*(1./ARG2)*(ARG1 -
     1     4.*(1./ARG2)*(ARG1 - 3.*(1./ARG2)*(ARG1 -
     2     2.*(1./ARG2)*(ARG1 + (1./ARG2)*(1-ARG1))))))
C
C-----------------------------------------------------------------------
C
C     * INITIALIZATIONS.
C
      YEFPA=EXP(YFPA) 
      APA0=YFPB*ZINTE(YEFPA,YFPA)-YFPC/6.
      TVFAP1=TVFA+1.
      DO 10 IL=1,ILGG
         TAUC(IL)=ADEF
 10   CONTINUE
      DO 20 L=ILEV,MSG+2,-1 
      DO 20 IL=1,ILGG 
         AFMAX(IL,L)=0.
         ALPHA(IL,L)=0.
         ALPHM(IL,L)=0.
         DADZF(IL,L)=0.
         PE   (IL,L)=0.
         PD   (IL,L)=0.
         DG   (IL,L)=1.
         WRK  (IL,L)=0.
         COND (IL,L)=0.
         AUTO (IL,L)=0.
         QLCT (IL,L)=0.
         QLWC (IL,L)=0.
 20   CONTINUE
      IF ( ITRAC.NE.0 ) THEN
         DO 30 N=1,NTRAC
         DO 30 L=MSG+1,ILEV
         DO 30 IL=1,ILGG
            SC(IL,L,N)=0.
            RU(IL,L,N)=0.
 30      CONTINUE
         IF ( ISULF.NE.0 ) THEN
            DO 40 N=1,NTRACA
            DO 40 L=MSG+1,ILEV
            DO 40 IL=1,ILGG
               RUA(IL,L,N)=0.
 40         CONTINUE
         END IF
      END IF
C
C     * FRACTIONAL ENTRAINMENT RATE E=AMU*DBU/DZ FOR LATERAL 
C     * ENTRAINMENT AND PD FOR LATERAL DETRAINMENT. ALL CLOUD 
C     * PROPERTIES REFER TO THE UNDILUTED PARCEL ASCENDING FROM 
C     * CLOUD BASE.
C
      DO 50 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LCL(IL)+1
            TV=TFE(IL,L)*(1.+TVFA*QF(IL,L))
            TPV=TC(IL,L)*(1.-QC(IL,L)+TVFAP1*QSTC(IL,L))
            IF ( ZF(IL,L).LE.ZG(IL,L-1) ) THEN
               DTEDZ=(TG(IL,L-1)-TG(IL,L))/(ZG(IL,L-1)-ZG(IL,L))
               DREDZ=(QG(IL,L-1)-QG(IL,L))/(ZG(IL,L-1)-ZG(IL,L))
               DPDZ =(PG(IL,L-1)-PG(IL,L))/(ZG(IL,L-1)-ZG(IL,L))
            ELSE
               DTEDZ=(TG(IL,L-2)-TG(IL,L-1))/(ZG(IL,L-2)-ZG(IL,L-1))
               DREDZ=(QG(IL,L-2)-QG(IL,L-1))/(ZG(IL,L-2)-ZG(IL,L-1))
               DPDZ =(PG(IL,L-2)-PG(IL,L-1))/(ZG(IL,L-2)-ZG(IL,L-1))
            ENDIF
            DHDZ=0.
            DRTDZ=0.
            DQSDP=-QSTC(IL,L)/(PF(IL,L)-EXP(A-B/TC(IL,L)))
            DQSDT=B*QSTC(IL,L)*(QSTC(IL,L)+EPS1)/(EPS1*TC(IL,L)**2)
            DTPDZ=-(GRAV+RRL*DQSDP*DPDZ-DHDZ)/(CPRES+RRL*DQSDT)
            DRPDZ=DQSDP*DPDZ+DQSDT*DTPDZ
            DBCDZ=(TPV/TV)*((TVFAP1*DRPDZ-DRTDZ)
     1               /(1.-QC(IL,L)+TVFAP1*QSTC(IL,L))
     2               +DTPDZ/TC(IL,L)-(TVFA*DREDZ/(1.+TVFA*QF(IL,L))
     3               +DTEDZ/TFE(IL,L)))*GRAV
C
C         *** PUT BOUNDS ON PE AND PD.
C
            PE(IL,L)=MIN(MAX( AMU*DBCDZ,ZERO),OPEM2)
            PD(IL,L)=MIN(MAX(-AMU*DBCDZ,ZERO),OPEM2)
C
C         *** SET LATERAL DETRAINMENT TO ZERO.
C
C            PD(IL,L)=0.
         ENDIF
 50   CONTINUE
      DO 70 L=ILEV,MSG+2,-1
      DO 70 IL=1,ILGG
         IF (       L.LT.(LCL(IL)+1) .AND. L.GE.LTM(IL) 
     1        .AND. ISKIP(IL).EQ.0 ) THEN
            TV=TFE(IL,L)*(1.+TVFA*QF(IL,L))
            TPV=TC(IL,L)*(1.-QC(IL,L)+TVFAP1*QSTC(IL,L))
            DTEDZ=( TG(IL,L-1)-TFE(IL,L) )/( ZG(IL,L-1)-ZF(IL,L) )
            DREDZ=( QG(IL,L-1)- QF(IL,L) )/( ZG(IL,L-1)-ZF(IL,L) )
            DPDZ =( PG(IL,L-1)- PF(IL,L) )/( ZG(IL,L-1)-ZF(IL,L) )
            DHDZ=0.
            DRTDZ=0.
            DQSDP=-QSTC(IL,L)/(PF(IL,L)-EXP(A-B/TC(IL,L)))
            DQSDT=B*QSTC(IL,L)*(QSTC(IL,L)+EPS1)/(EPS1*TC(IL,L)**2)
            DTPDZ=-(GRAV+RRL*DQSDP*DPDZ-DHDZ)/(CPRES+RRL*DQSDT)
            DRPDZ=DQSDP*DPDZ+DQSDT*DTPDZ
            DBCDZ=(TPV/TV)*((TVFAP1*DRPDZ-DRTDZ)
     1               /(1.-QC(IL,L)+TVFAP1*QSTC(IL,L))
     2               +DTPDZ/TC(IL,L)-(TVFA*DREDZ/(1.+TVFA*QF(IL,L))
     3               +DTEDZ/TFE(IL,L)))*GRAV
C
C         *** PUT BOUNDS ON PE AND PD.
C
            PE(IL,L)=MIN(MAX( AMU*DBCDZ,ZERO),OPEM2)
            PD(IL,L)=MIN(MAX(-AMU*DBCDZ,ZERO),OPEM2)
C
C         *** SET LATERAL DETRAINMENT TO ZERO.
C
C            PD(IL,L)=0.
         ENDIF
 70   CONTINUE
C
C     * PROBABILITY DENSITY FUNCTION AND ALPHA FOR CLOUD
C     * INHOMOGENEITIES, EXPRESSED IN TERMS OF UNDILUTE PARCEL 
C     * PROPERTIES.
C
      DO 80 L=ILEV,MSG+2,-1
      DO 80 IL=1,ILGG
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0 ) THEN
C
C         *** MIXING FRACTION AT THE EDGE OF THE CLOUD.
C
            DHDF=HMNF(IL,L)-HMNC(IL,L)
            DQDF=QF(IL,L)-QC(IL,L)
C
            DQSDT=B*QSTC(IL,L)*(QSTC(IL,L)+EPS1)/(EPS1*TC(IL,L)**2)
            DTDF=DHDF/(CPRES+RRL*DQSDT)
            AFMAX(IL,L)=(QC(IL,L)-QSTC(IL,L))/(DQSDT*DTDF-DQDF)
            AFMAX(IL,L)=MAX(MIN(AFMAX(IL,L),ONE),ZERO)
C
C         *** CLOUD-TOP ENTAINMENT PARAMETER AND HORIZONTALLY AVERAGED
C         *** BUOYANCY IN THE CLOUD FOR AN EXPONENTIALLY DECAYING PDF.
C
            AFC=0.
            AF2C=0.
            IF ( .NOT.KEQPROB ) THEN
               IF ( APHI.GT.0. ) THEN
                  IF ( TSTAR(IL,L).NE.ADEF .AND. TSTAR(IL,L).GT.YSMALL
     1                 .AND. AFMAX(IL,L).GT.0. ) THEN
                     AG=EXP(-APHI*TSTAR(IL,L))
                     APP=AG/(1.-AG)
                     AEXPF=EXP(-APP*AFMAX(IL,L))
                     IF ( AEXPF.EQ.0. ) THEN
                        AFC=0.
                        AF2C=0.
                     ELSE IF ( AEXPF.EQ.1. ) THEN
                        AFC=.5*AFMAX(IL,L)
                        AF2C=(AFMAX(IL,L)**2)/3.
                     ELSE
                        AFC=1./APP-AFMAX(IL,L)*AEXPF/(1.-AEXPF)
                        AFC=MAX(MIN(AFC,.5*AFMAX(IL,L)),ZERO)
                        AF2C=2./APP**2-(2.*AFMAX(IL,L)/APP
     1                      +AFMAX(IL,L)**2)*AEXPF/(1.-AEXPF)
                        AF2C=MAX(MIN(AF2C,(AFMAX(IL,L)**2)/3.),ZERO)
                     END IF
                  ELSE IF ( TSTAR(IL,L).EQ.ADEF ) THEN
                     AFC=.5*AFMAX(IL,L)
                     AF2C=(AFMAX(IL,L)**2)/3.
                  ELSE
                     AFC=0.
                     AF2C=0.
                  END IF
               END IF 
            ELSE
               AFC=.5*AFMAX(IL,L)
               AF2C=(AFMAX(IL,L)**2)/3.
            END IF
            IF ( AFC.NE.1. ) THEN
               WRK(IL,L)=AF2C
               ALPHA(IL,L)=AFC/(1.-AFC)
            ELSE
               AFC=0.999999
               WRK(IL,L)=0.999999**2
               ALPHA(IL,L)=AFC/(1.-AFC)
            END IF
         END IF
 80   CONTINUE
C
C     * CLOUD CORE PROPERTIES OF TOTAL WATER AND MOIST STATIC ENERGY.
C
      DO 140 L=ILEVM,MSG+2,-1
      DO 140 IL=1,ILGG
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0 ) THEN
            ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                     *PD(IL,L+1)
            ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                     *PD(IL,L  )
            IF ( ATL*ATH.GT.0. ) THEN
C
C            *** MAXIMUM VALUE OF PD. THIS IS BASED ON THE ASSUMPTION
C            *** THAT LATERAL DETRAINMENT CANNOT INCREASE THE VALUES
C            *** OF THE PROPERTIES IN THE CLOUD CORE TO VALUES ABOVE
C            *** THOSE OF THE INIITIATING PARCEL OF AIR OR DECREASE
C            *** THE VALUES BELOW ZERO.
C
               ACTMP=QCT(IL,L+1)-QG(IL,L)
               IF ( ALPHA(IL,L).NE.0. .AND. PD(IL,L).NE.0. 
     1             .AND. ACTMP.LT.0. ) THEN
                  ARAT=-ACTMP/(QG(IL,L)*0.99)
                  APDMAX=(ATL-2.*LOG(ARAT)/DZ(IL,L)+PE(IL,L))
     1                  *((ALPHA(IL,L)+1.)/ALPHA(IL,L))
                  IF ( PD(IL,L).GT.APDMAX ) THEN
                     PD(IL,L)=APDMAX
                     ATH=PE(IL,L)-(ALPHA(IL,L)/(1.+ALPHA(IL,L)))
     1                            *PD(IL,L)
                  END IF
               END IF
               IF ( ATL*ATH.GT.0. ) THEN
                  DG(IL,L)=EXP(-.5*DZ(IL,L)*(ATH+ATL))
                  QCT(IL,L)=QG  (IL,L)+(QCT(IL,L+1)-  QG(IL,L))
     1                     *DG(IL,L)
                  HCT(IL,L)=HMNG(IL,L)+(HCT(IL,L+1)-HMNG(IL,L))
     1                     *DG(IL,L)
               ELSE
                  QCT(IL,L)=QCT(IL,L+1)
                  HCT(IL,L)=HCT(IL,L+1)
               END IF
            ELSE
               QCT(IL,L)=QCT(IL,L+1)
               HCT(IL,L)=HCT(IL,L+1)
            END IF
         ENDIF
 140  CONTINUE
C
C     * OBTAIN TEMPERATURE OF CLOUD CORE FROM MOIST STATIC ENERGY.
C
      DO 160 IL=1,ILGG 
         LTM(IL)=LCL(IL)+1
 160  CONTINUE
      CALL TINCLD3(TC,ILG,ILEV,MSG,1,ILGG,LTM,TDRY,HCT,QCT,QSTC,ZF,
     1             PF,RRL,GRAV,CPRES,EPS1,EPS2,A,B,ISKIP)
C
C     * CLOUD CORE DRY STATIC ENERGY.
C
      DO 180 L=ILEVM,MSG+2,-1
      DO 180 IL=1,ILGG
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0 ) THEN
            DCT(IL,L)=HCT(IL,L)-RRL*QSTC(IL,L)
         END IF
 180  CONTINUE
C
C     * DERIVED CLOUD PROPERTIES. 
C
      DO 220 L=ILEV,MSG+2,-1
      DO 220 IL=1,ILGG
         IF ( L.LE.(LCL(IL)+1) .AND. ISKIP(IL).EQ.0 ) THEN
C
C         *** MIXING FRACTION AT THE EDGE OF THE CLOUD.
C
            DHDF=HMNF(IL,L)-HCT(IL,L)
            DQDF=QF(IL,L)-QCT(IL,L)
            IF ( TC(IL,L).GT.0. ) THEN
               DQSDT=B*QSTC(IL,L)*(QSTC(IL,L)+EPS1)/(EPS1*TC(IL,L)**2)
            ELSE
               DQSDT=0.
            END IF
            DTDF=DHDF/(CPRES+RRL*DQSDT)
            AFMAX(IL,L)=(QCT(IL,L)-QSTC(IL,L))/(DQSDT*DTDF-DQDF)
            AFMAX(IL,L)=MAX(MIN(AFMAX(IL,L),ONE),ZERO)
C
C         *** PUT AN UPPER BOUND ON ALPHA IF FC EXCEEDS AFMAX.
C
            ALPHAR=ALPHA(IL,L)/(1.+ALPHA(IL,L))
            AFC=MAX(MIN(ALPHAR,AFMAX(IL,L)),ZERO)
            AF2C=MAX(MIN(WRK(IL,L),AFMAX(IL,L)**2),ZERO)
            ALPHA(IL,L)=AFC/(1.-AFC)
C
C         *** MEAN BUOYANCY IN THE CLOUD.
C
            ATC=TC(IL,L)+DTDF*AFC
            ATQC =TC(IL,L)*QCT (IL,L)+(DTDF*QCT (IL,L)
     1           +DQDF      *TC(IL,L))*AFC+DTDF   *DQDF *AF2C
            ATQSC=TC(IL,L)*QSTC(IL,L)+(DTDF*QSTC(IL,L)
     1           +DQSDT*DTDF*TC(IL,L))*AFC+DTDF**2*DQSDT*AF2C
            TPV=ATC-ATQC+TVFAP1*ATQSC
            TV=TFE(IL,L)*(1.+TVFA*QF(IL,L))
            BUOY(IL,L)=GRAV*( TPV-TV )/( TV*(1.+GAMMA) )
C
C         *** LIQUID WATER CONTENT (WITHOUT PRECIPITATION).
C
            QLCT(IL,L)=MAX(QCT(IL,L)-QSTC(IL,L),0.*QCT(IL,L))
         ENDIF
 220  CONTINUE
C
C     * MEAN AND VERTICAL GRADIENT OF CLOUD-TOP ENTRAINMENT PARAMETER.
C
      DO 280 L=ILEVM,MSG+2,-1
      DO 280 IL=1,ILGG
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0  ) THEN
            ALPHM(IL,L)=.5*( ALPHA(IL,L)+ALPHA(IL,L+1) )
         ENDIF 
 280  CONTINUE
      DO 300 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LCL(IL)+1
            DADZF(IL,L)=( ALPHA(IL,L-1)-ALPHA(IL,L) )
     1                 /(    ZF(IL,L-1)-   ZF(IL,L) )
         ENDIF
 300  CONTINUE
      DO 320 L=ILEV,MSG+3,-1
      DO 320 IL=1,ILGG
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0  ) THEN
            DADZF(IL,L)=2.*( ALPHM(IL,L-1)-ALPHM(IL,L) )
     1                    /( ZF(IL,L-1)-ZF(IL,L+1) )
         ENDIF
 320  CONTINUE
      DO 330 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=MSG+2
            DADZF(IL,L)=( ALPHA(IL,L)-ALPHA(IL,L+1) )
     1                 /(    ZF(IL,L)-   ZF(IL,L+1) )
         ENDIF
 330  CONTINUE
C
C     * HORIZONTALLY AVERAGED VERTICAL VELOCITY.
C
      DO 340 L=ILEVM,MSG+2,-1
      DO 340 IL=1,ILGG
         IF ( L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0  ) THEN
            ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                     *PD(IL,L+1)
            ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                     *PD(IL,L  )
            ATH=ATH+DADZF(IL,L  )/(1.+ALPHA(IL,L  ))
            ATL=ATL+DADZF(IL,L+1)/(1.+ALPHA(IL,L+1))
            IF ( WC(IL,L+1).LT. 0. ) THEN
               ATMP=ADEF
            ELSE IF ( ATL*ATH.LE.0. ) THEN
               ATMP=.5*DZ(IL,L)*(BUOY(IL,L)+BUOY(IL,L+1))
               ATMP=WC(IL,L+1)**2+2.*ATMP/(1.+ALPHM(IL,L))
            ELSE 
               DGT=DZ(IL,L)*(ATH+ATL)
               PLYH=BUOY(IL,L  )/(ATH*(1.+ALPHA(IL,L  )))
               PLYL=BUOY(IL,L+1)/(ATL*(1.+ALPHA(IL,L+1)))
               DPLDG=( PLYH-PLYL )/DGT
               ATMP=PLYH-DPLDG+(DPLDG+WC(IL,L+1)**2-PLYL)/EXP(DGT)
            ENDIF
            WC(IL,L)=SIGN(SQRT(ABS(ATMP)),ATMP)
         ENDIF
 340  CONTINUE
C
C     * SKIP PATHOLOGICAL CASE.
C
      DO 350 IL=1,ILGG
         L=LCL(IL)
         WTHRES=DZO(IL,L-1)/MAX(DTO,TFHUN)
         IF ( WC(IL,L).LE.WTHRES .AND. ISKIP(IL).EQ.0 ) ISKIP(IL)=2
 350  CONTINUE
C
C     * ASSIGN TOP OF CLOUD LAYER TO LEVEL OF NEUTRAL BUOYANCY (LNB).
C     * THE LNB IS DEFINED AS THE FIRST LEVEL AT WHICH THE CLOUDY AIR
C     * IS NEUTRALLY BUOYANT BELOW THAT LEVEL AT WHICH THE MEAN VERTICAL
C     * VELOCITY IN THE CLOUD DROPS TO ZERO. 
C
      LDEF=ILEV
      DO 360 IL=1,ILGG
         LNB(IL)=LDEF
         LTM(IL)=LDEF
         FCT(IL)=0.
 360  CONTINUE
      DO 370 L=ILEVM-1,MSG+2,-1
      DO 370 IL=1,ILGG 
         IF (       L.LT.LCL(IL) .AND. LTM(IL).EQ.LDEF 
     1        .AND. ISKIP(IL).EQ.0  ) THEN
            WTHRES=DZO(IL,L-1)/MAX(DTO,TFHUN)	
            IF ( WC(IL,L).LE.WTHRES .AND. WC(IL,L+1).GT.WTHRES ) THEN
               LTM(IL)=L+1
               FCT(IL)=(WTHRES-WC(IL,L))/(WC(IL,L+1)-WC(IL,L))
            ELSE IF ( WC(IL,L+1).LE.WTHRES .AND. (L+1).LT.LCL(IL)) THEN
               LTM(IL)=L+2
               FCT(IL)=0.
            ENDIF
         ENDIF
 370  CONTINUE
      DO 380 IL=1,ILGG 
         IF ( LNB(IL).EQ.LDEF .AND. ISKIP(IL).EQ.0 ) THEN
            IF ( LTM(IL).NE.LDEF ) THEN 
               LNB(IL)=LTM(IL)
            ELSE
               ISKIP(IL)=4
            ENDIF
         ENDIF
 380  CONTINUE
C
C     * CALCULATE PROPERTIES AT TOP OF THE CLOUD.
C
      DO 385 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LNB(IL)-1
            BUOY (IL,L)=FCT(IL)*BUOY (IL,L+1)+(1.-FCT(IL))*BUOY (IL,L)
            PF   (IL,L)=FCT(IL)*PF   (IL,L+1)+(1.-FCT(IL))*PF   (IL,L)
            ZF   (IL,L)=FCT(IL)*ZF   (IL,L+1)+(1.-FCT(IL))*ZF   (IL,L)
            DZ   (IL,L)=ZF(IL,L)-ZF(IL,L+1)
            WC   (IL,L)=FCT(IL)*WC   (IL,L+1)+(1.-FCT(IL))*WC   (IL,L)
            WC   (IL,L)=MAX(WC(IL,L),0.*WC(IL,L))
            QSTC (IL,L)=FCT(IL)*QSTC (IL,L+1)+(1.-FCT(IL))*QSTC (IL,L)
            QLCT (IL,L)=FCT(IL)*QLCT (IL,L+1)+(1.-FCT(IL))*QLCT (IL,L)
            QF   (IL,L)=FCT(IL)*QF   (IL,L+1)+(1.-FCT(IL))*QF   (IL,L)
            HMNF (IL,L)=FCT(IL)*HMNF (IL,L+1)+(1.-FCT(IL))*HMNF (IL,L)
            TC   (IL,L)=FCT(IL)*TC   (IL,L+1)+(1.-FCT(IL))*TC   (IL,L)
            TFE  (IL,L)=FCT(IL)*TFE  (IL,L+1)+(1.-FCT(IL))*TFE  (IL,L)
            RHOA (IL,L)=FCT(IL)*RHOA (IL,L+1)+(1.-FCT(IL))*RHOA (IL,L)
            ALPHA(IL,L)=FCT(IL)*ALPHA(IL,L+1)+(1.-FCT(IL))*ALPHA(IL,L)
            PE   (IL,L)=FCT(IL)*PE   (IL,L+1)+(1.-FCT(IL))*PE   (IL,L)
            PD   (IL,L)=FCT(IL)*PD   (IL,L+1)+(1.-FCT(IL))*PD   (IL,L)
            ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                     *PD(IL,L+1)
            ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                     *PD(IL,L  )
            IF ( ATL*ATH.GT.0. ) THEN
               DG(IL,L)=EXP(-.5*DZ(IL,L)*(ATH+ATL))
            ELSE
               DG(IL,L)=1.
            END IF
         ENDIF
 385  CONTINUE
      IF ( ITRAC.NE.0 ) THEN
         DO 390 N=1,NTRAC
         DO 390 IL=1,ILGG
            IF ( ISKIP(IL).EQ.0 ) THEN
               L=LNB(IL)-1
               XF (IL,L,N)=    FCT(IL) *XF (IL,L+1,N)
     1                    +(1.-FCT(IL))*XF (IL,L,N  )
            END IF
 390     CONTINUE
      END IF
C
C     * CONDENSATION RATE.
C
      DO 395 L=ILEVM,MSG+2,-1
      DO 395 IL=1,ILGG
         IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1        .AND. ISKIP(IL).EQ.0                              ) THEN
            ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                     *PD(IL,L+1)
            ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                     *PD(IL,L  )
            IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN
               APA=.5*(WC(IL,L+1)*ATL+WC(IL,L)*ATH)
               COND(IL,L)=-APA*((QSTC(IL,L)-DG(IL,L)*QSTC(IL,L+1))
     1                        /(1.-DG(IL,L))-QG(IL,L))
            ELSE
               DT=2.*DZ(IL,L)/(WC(IL,L+1)+WC(IL,L))
               COND(IL,L)=-(QSTC(IL,L)-QSTC(IL,L+1))/DT
            END IF
         END IF
 395  CONTINUE
C
C     * PRELIMINARY CONCENTRATIONS FOR CALCULATION OF
C     * DROPLET NUMBER CONCENTRATION IN AUTOCONVERSION.
C
      IF ( ISULF.EQ.1 ) THEN
         DO 396 N=1,NTRAC 
           IF ( ((N.EQ.ISO4).OR.(N.EQ.ISSA).OR.(N.EQ.ISSC))
     1                                       .AND.ITRPHS(N).GT.0 ) THEN
             DO 397 L=ILEVM,MSG+2,-1
             DO 397 IL=1,ILGG
                IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1               .AND. ISKIP(IL).EQ.0                     ) THEN
                   ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                            *PD(IL,L+1)
                   ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                            *PD(IL,L  )
                   IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN
                      APA=.5*(WC(IL,L+1)*ATL+WC(IL,L)*ATH)
                      APB=1.-DG(IL,L)
                      XCT(IL,L,N)=APA*(XG(IL,L,N)+(XCT(IL,L+1,N)
     1                                -XG(IL,L,N))*DG(IL,L))
     2                                 /(APA+APB*RU(IL,L,N))
                      XCT(IL,L,N)=MAX(XCT(IL,L,N),ZERO)
                   ELSE
                      DT=2.*DZ(IL,L)/(WC(IL,L+1)+WC(IL,L))
                      XCT(IL,L,N)=XCT(IL,L+1,N)/(1.+DT*RU(IL,L,N))
                   ENDIF
                ENDIF
 397         CONTINUE
           ENDIF
 396     CONTINUE
      ELSE
         DO 398 L=ILEVM,MSG+2,-1
         DO 398 IL=1,ILGG
            IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1           .AND. ISKIP(IL).EQ.0                     ) THEN
               XCT(IL,L,N)=SO4CROW(ISHALL(IL),L)
            ENDIF
 398     CONTINUE
      ENDIF
C
C     * MASS SCALING FACTOR. ASSUME THAT ALL SULPHATE AEROSOL
C     * IS IN THE FORM OF SO4.
C     
      SO4MASS=96.058/32.064
      SSAMASS=22.990/58.443
      SSCMASS=22.990/58.443
C
C
      DO 400 L=ILEVM,MSG+2,-1
      DO 400 IL=1,ILGG
         IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1        .AND. QLCT(IL,L).GT.0. .AND. ISKIP(IL).EQ.0       ) THEN
            ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                     *PD(IL,L+1)
            ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                     *PD(IL,L  )
C
            IF ( IPAM.EQ.0 ) THEN
              ZRHO0=100.*PF(IL,L)/(RGAS*TFE(IL,L))
              CSO4=XCT(IL,L,ISO4)*SO4MASS*ZRHO0*1.E+09
              CSS =XCT(IL,L,ISSA)*SSAMASS*ZRHO0*1.E+09
     1            +XCT(IL,L,ISSC)*SSCMASS*ZRHO0*1.E+09
              SURAER=CSO4+0.27*CSS
              PCDNC=100.*CSO4**0.2
              PCDNC=PCDNC*0.6
              PCDNC=MAX(1.E+06*PCDNC,ZCDNMIN)
            ELSE
              PCDNC=1.E+08
            ENDIF
            SCDN(ISHALL(IL),L)=PCDNC
            AUTCF=(ZCRAUT*6.E+28*YAPN*(PCDNC
     1           *1.E-06)**(-3.3)*(1.E-03*RHOA(IL,L))**YAUTEX)
     2           /RHOA(IL,L)
C
C           * ACCOUNT FOR CLOUD INHOMOGENEITY. THIS IS BASED
C           * ON BEHENG'S FORMULATION OF THE AUTOCONVERSION RATE.

            IF ( .NOT.KEQPROB ) THEN
               APA=1.
               IF ( APHI.GT.0. ) THEN
                  IF ( TSTAR(IL,L).NE.ADEF .AND. TSTAR(IL,L).GT.YSMALL
     1                 .AND. AFMAX(IL,L).GT.0. ) THEN
                     AG=EXP(-APHI*TSTAR(IL,L))
                     IF ( AG.GT. 23.1/(23.1+AFMAX(IL,L)) ) THEN
                        AEXPF=1.E-10
                     ELSE
                        APP=AG/(1.-AG)
                        AEXPF=EXP(-APP*AFMAX(IL,L))
                     ENDIF
                     IF ( ABS(AEXPF).LE.1.E-10 ) THEN
                        APA=1.
                     ELSE IF ( AEXPF.EQ.1. ) THEN
                        APA=APA0
                     ELSE
                        APDF0=APP/(1.-AEXPF)
                        ATEMPO=YFPA+APP*AFMAX(IL,L)
                        AETMP =YEFPA/AEXPF
                        IF ( ATEMPO.LT.YCUTH .AND. ATEMPO.GT.YCUTL 
     1                                                          ) THEN
                           APA=YFPB/6.
                        ELSE
                           APA=YFPB*ZINTE(AETMP,ATEMPO)
                        END IF
                        ATEMPO=APP*AFMAX(IL,L)
                        AETMP =1./AEXPF
                        IF ( ATEMPO.LT.YCUTH .AND. ATEMPO.GT.YCUTL 
     1                                                          ) THEN
                           APA=APA-YFPC/6.
                        ELSE
                           APA=APA-YFPC*ZINTE(AETMP,ATEMPO)
                        END IF
                        APA=APA*AEXPF*AFMAX(IL,L)*APDF0
                     END IF
                  ELSE IF ( TSTAR(IL,L).EQ.ADEF ) THEN
                     APA=APA0
                  ELSE
                     APA=1.
                  END IF
               END IF 
            ELSE
               APA=APA0
            END IF
               if(apa.lt.0.) then
                 print *, '0apa = ',apa
                 call abort
               endif
            AUTCF=AUTCF*APA
C
            IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN
               APA=.5*(WC(IL,L+1)*ATL+WC(IL,L)*ATH)
               APB=(1.-DG(IL,L))/APA
               APC=-(1.-DG(IL,L))*COND(IL,L)/APA-DG(IL,L)*QLCT(IL,L+1)
               QLCT(IL,L)=(1.-DG(IL,L))*COND(IL,L)/APA
     1                   +DG(IL,L)*QLCT(IL,L+1)
               QLCT(IL,L)=MIN(QLCT(IL,L)
     1                       ,(MAX(-APC/(APB*AUTCF),ZERO))**(1/YAUTEX))
C
C            *** FIRST ITERATION.
C
               AUTO(IL,L)=SIGN(AUTCF*ABS(QLCT(IL,L))**YAUTEX,QLCT(IL,L))
               AFUN=QLCT(IL,L)+APC+APB*AUTO(IL,L)
               IF ( QLCT(IL,L).GT.(1./(APB*AUTCF))**(1./(YAUTEX-1.)) 
     1                                                            ) THEN 
                  ADFUN=AUTCF/(YAUTEX*AUTO(IL,L)/QLCT(IL,L))+APB*AUTCF
                  QLCT(IL,L)=(AUTO(IL,L)/AUTCF-AFUN/ADFUN)**(1./YAUTEX)
               ELSE IF ( QLCT(IL,L).NE.0. ) THEN
                  ADFUN=1.+YAUTEX*APB*AUTO(IL,L)/QLCT(IL,L)
                  QLCT(IL,L)=QLCT(IL,L)-AFUN/ADFUN
               ELSE
                  QLCT(IL,L)=QLCT(IL,L)-AFUN
               END IF
               QLCT(IL,L)=MIN(QLCT(IL,L)
     1                       ,(MAX(-APC/(APB*AUTCF),ZERO))**(1/YAUTEX))
C
C            *** SECOND ITERATION.
C
               AUTO(IL,L)=SIGN(AUTCF*ABS(QLCT(IL,L))**YAUTEX,QLCT(IL,L))
               AFUN=QLCT(IL,L)+APC+APB*AUTO(IL,L)
               IF ( QLCT(IL,L).GT.(1./(APB*AUTCF))**(1./(YAUTEX-1.)) 
     1                                                            ) THEN 
                  ADFUN=AUTCF/(YAUTEX*AUTO(IL,L)/QLCT(IL,L))+APB*AUTCF
                  QLCT(IL,L)=(AUTO(IL,L)/AUTCF-AFUN/ADFUN)**(1./YAUTEX)
               ELSE IF ( QLCT(IL,L).NE.0. ) THEN
                  ADFUN=1.+YAUTEX*APB*AUTO(IL,L)/QLCT(IL,L)
                  QLCT(IL,L)=QLCT(IL,L)-AFUN/ADFUN
               ELSE
                  QLCT(IL,L)=QLCT(IL,L)-AFUN
               END IF
               QLCT(IL,L)=MIN(QLCT(IL,L)
     1                       ,(MAX(-APC/(APB*AUTCF),ZERO))**(1/YAUTEX))
C
C            *** THIRD ITERATION.
C
               AUTO(IL,L)=SIGN(AUTCF*ABS(QLCT(IL,L))**YAUTEX,QLCT(IL,L))
               AFUN=QLCT(IL,L)+APC+APB*AUTO(IL,L)
               IF ( QLCT(IL,L).GT.(1./(APB*AUTCF))**(1./(YAUTEX-1.)) 
     1                                                            ) THEN 
                  ADFUN=AUTCF/(YAUTEX*AUTO(IL,L)/QLCT(IL,L))+APB*AUTCF
                  QLCT(IL,L)=(AUTO(IL,L)/AUTCF-AFUN/ADFUN)**(1./YAUTEX)
               ELSE IF ( QLCT(IL,L).NE.0. ) THEN
                  ADFUN=1.+YAUTEX*APB*AUTO(IL,L)/QLCT(IL,L)
                  QLCT(IL,L)=QLCT(IL,L)-AFUN/ADFUN
               ELSE
                  QLCT(IL,L)=QLCT(IL,L)-AFUN
               END IF
               QLCT(IL,L)=MIN(QLCT(IL,L)
     1                       ,(MAX(-APC/(APB*AUTCF),ZERO))**(1/YAUTEX))
            ELSE
               DT=2.*DZ(IL,L)/(WC(IL,L)+WC(IL,L+1))
               QLCT(IL,L)=QLCT(IL,L+1)+DT*COND(IL,L)
               QLCT(IL,L)=MIN(QLCT(IL,L),(MAX((QLCT(IL,L+1)
     1                   +DT*COND(IL,L))/(DT*AUTCF),ZERO))**(1/YAUTEX))
C
C            *** FIRST ITERATION.
C
               AUTO(IL,L)=SIGN(AUTCF*ABS(QLCT(IL,L))**YAUTEX,QLCT(IL,L))
               AFUN=QLCT(IL,L)-QLCT(IL,L+1)+DT*(AUTO(IL,L)-COND(IL,L))
               IF ( QLCT(IL,L).GT.(1./(DT*AUTCF))**(1./(YAUTEX-1.))
     1                                                            ) THEN 
                  ADFUN=AUTCF/(YAUTEX*AUTO(IL,L)/QLCT(IL,L))+DT*AUTCF
                  QLCT(IL,L)=(AUTO(IL,L)/AUTCF-AFUN/ADFUN)**(1./YAUTEX)
               ELSE IF ( QLCT(IL,L).NE.0. ) THEN
                  ADFUN=1.+YAUTEX*DT*AUTO(IL,L)/QLCT(IL,L)
                  QLCT(IL,L)=QLCT(IL,L)-AFUN/ADFUN
               ELSE
                  QLCT(IL,L)=QLCT(IL,L)-AFUN
               END IF
               QLCT(IL,L)=MIN(QLCT(IL,L),(MAX((QLCT(IL,L+1)
     1                   +DT*COND(IL,L))/(DT*AUTCF),ZERO))**(1/YAUTEX))
C
C            *** SECOND ITERATION.
C
               AUTO(IL,L)=SIGN(AUTCF*ABS(QLCT(IL,L))**YAUTEX,QLCT(IL,L))
               AFUN=QLCT(IL,L)-QLCT(IL,L+1)+DT*(AUTO(IL,L)-COND(IL,L))
               IF ( QLCT(IL,L).GT.(1./(DT*AUTCF))**(1./(YAUTEX-1.))
     1                                                            ) THEN 
                  ADFUN=AUTCF/(YAUTEX*AUTO(IL,L)/QLCT(IL,L))+DT*AUTCF
                  QLCT(IL,L)=(AUTO(IL,L)/AUTCF-AFUN/ADFUN)**(1./YAUTEX)
               ELSE IF ( QLCT(IL,L).NE.0. ) THEN
                  ADFUN=1.+YAUTEX*DT*AUTO(IL,L)/QLCT(IL,L)
                  QLCT(IL,L)=QLCT(IL,L)-AFUN/ADFUN
               ELSE
                  QLCT(IL,L)=QLCT(IL,L)-AFUN
               END IF
               QLCT(IL,L)=MIN(QLCT(IL,L),(MAX((QLCT(IL,L+1)
     1                   +DT*COND(IL,L))/(DT*AUTCF),ZERO))**(1/YAUTEX))
C
C            *** THIRD ITERATION.
C
               AUTO(IL,L)=SIGN(AUTCF*ABS(QLCT(IL,L))**YAUTEX,QLCT(IL,L))
               AFUN=QLCT(IL,L)-QLCT(IL,L+1)+DT*(AUTO(IL,L)-COND(IL,L))
               IF ( QLCT(IL,L).GT.(1./(DT*AUTCF))**(1./(YAUTEX-1.))
     1                                                            ) THEN 
                  ADFUN=AUTCF/(YAUTEX*AUTO(IL,L)/QLCT(IL,L))+DT*AUTCF
                  QLCT(IL,L)=(AUTO(IL,L)/AUTCF-AFUN/ADFUN)**(1./YAUTEX)
               ELSE IF ( QLCT(IL,L).NE.0. ) THEN
                  ADFUN=1.+YAUTEX*DT*AUTO(IL,L)/QLCT(IL,L)
                  QLCT(IL,L)=QLCT(IL,L)-AFUN/ADFUN
               ELSE
                  QLCT(IL,L)=QLCT(IL,L)-AFUN
               END IF
               QLCT(IL,L)=MIN(QLCT(IL,L),(MAX((QLCT(IL,L+1)
     1                   +DT*COND(IL,L))/(DT*AUTCF),ZERO))**(1/YAUTEX))
            END IF
            QLCT(IL,L)=MAX(QLCT(IL,L),ZERO)
            AUTO(IL,L)=MAX(AUTO(IL,L),ZERO)
         ENDIF
 400  CONTINUE
C
C     * ALLOW ONLY NON-ZERO LIQUID WATER BELOW TOP OF THE CLOUD.
C
      DO 405 L=ILEVM,MSG+2,-1
      DO 405 IL=1,ILGG 
         IF (      L.LT.LCL(IL) .AND. L.GE.LNB(IL) .AND. ISKIP(IL).EQ.0 
     1       .AND. QLCT(IL,L).EQ.0.  ) THEN
            LNB(IL)=L+1
            FCT(IL)=0.
         ENDIF
 405  CONTINUE
C
C     * TOTAL CLOUD WATER CONTENT AFTER PRECIPITATION HAS BEEN 
C     * ACCOUNTED FOR. 
C
      DO 410 L=ILEVM,MSG+2,-1
      DO 410 IL=1,ILGG
         IF (      L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1       .AND. ISKIP(IL).EQ.0                               ) THEN
            QCT(IL,L)=QSTC(IL,L)+QLCT(IL,L)
         END IF 
 410  CONTINUE
C
C     * HORIZONALLY AVERAGED CLOUD PROPERTIES.
C
      DO L=ILEV,MSG+2,-1
      DO IL=1,ILGG
         IF (       L.LE.(LCL(IL)+1) .AND. L.GE.(LNB(IL)-1) 
     1        .AND. ISKIP(IL).EQ.0                              ) THEN
            DHDF=HMNF(IL,L)-HCT(IL,L)
            DQDF=QF(IL,L)-QCT(IL,L)
            IF ( TC(IL,L).GT.0. ) THEN
               DQSDT=B*QSTC(IL,L)*(QSTC(IL,L)+EPS1)/(EPS1*TC(IL,L)**2)
            ELSE
               DQSDT=0.
            END IF
            DTDF=DHDF/(CPRES+RRL*DQSDT)
            AFC=ALPHA(IL,L)/(1.+ALPHA(IL,L))
            QLWC(IL,L)=MAX(QCT(IL,L)-QSTC(IL,L)+(DQDF-DQSDT*DTDF)*AFC
     1                   ,0.*QCT(IL,L))
         ENDIF
      ENDDO
      ENDDO
C
      DO 415 L=ILEV,MSG+2,-1
      DO 415 IL=1,ILGG
         IF (       L.LE.(LCL(IL)+1) .AND. L.GE.(LNB(IL)-1) 
     1        .AND. ISKIP(IL).EQ.0                              ) THEN
            QC  (IL,L)=( QCT(IL,L)+ALPHA(IL,L)*QF(IL,L) )
     1                /( 1.+ALPHA(IL,L) )
            HMNC(IL,L)=( HCT(IL,L)+ALPHA(IL,L)*HMNF(IL,L) )
     1                /( 1.+ALPHA(IL,L) )
         ENDIF
 415  CONTINUE
C
C     * IN-CLOUD MIXING RATIOS OF CHEMICAL TRACERS.
C
      KCHEM=.FALSE.
      IF ( ITRAC.NE.0 ) THEN
         KCHEM=.TRUE.
C
C      *** NON-REACTIVE TRACERS.
C
         DO 420 N=1,NTRAC 
            IF ( IPAM.EQ.0 ) THEN
              KRIT1=ITRPHS(N).GT.0 
     1          .AND. ( (N.NE.ISO2 .AND. N.NE.ISO4 .AND. N.NE.IHPO) 
     2                  .OR. .NOT.KCHEM )
            ELSE
              KRIT1=ITRPHS(N).GT.0
            ENDIF
            IF ( KRIT1 ) THEN
               DO 430 L=ILEVM,MSG+2,-1
               DO 430 IL=1,ILGG
                  IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1                 .AND. ISKIP(IL).EQ.0                     ) THEN
                     ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                              *PD(IL,L+1)
                     ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                              *PD(IL,L  )
C
C                    *** PRELIMINARY SCAVENGING RATES IN 1/S.
C
                     ALWCVMR=(     BETA *QLCT(IL,L+1)*RHOA(IL,L+1)
     1                        +(1.-BETA)*QLCT(IL,L  )*RHOA(IL,L  ))
     2                          /YRHOWAT
                     IF (       L.LT.LCL(IL) .AND. WC(IL,L).GT.0.
     1                 .AND. ALWCVMR.GT.YEPS  ) THEN
                        APA=BETA*QLCT(IL,L+1)+(1.-BETA)*QLCT(IL,L)
                        APB=BETA*AUTO(IL,L+1)+(1.-BETA)*AUTO(IL,L)
C
                        IF ( IPAM.EQ.0 ) THEN
                          KRIT2=ITRCCH(N).EQ.-1
                        ELSE
                          KRIT2=ITRCCH(N).EQ.-1 .OR. N.EQ.IHPO
                        ENDIF
                        IF(KRIT2) THEN
                         RUFAC=ASRPHOB
                        ELSE IF ( ITRCCH(N).EQ.1 ) THEN
                         RUFAC=ASRSO4
                        ELSE
                         RUFAC=0.
                        ENDIF
                        RU (IL,L,N) = RUFAC*APB/APA
                     ENDIF
                     IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN
                        APA=.5*(WC(IL,L+1)*ATL+WC(IL,L)*ATH)
                        APB=1.-DG(IL,L)
                        XCT(IL,L,N)=APA*(XG(IL,L,N)+(XCT(IL,L+1,N)
     1                             -XG(IL,L,N))*DG(IL,L))
     2                             /(APA+APB*RU(IL,L,N))
                        XCT(IL,L,N)=MAX(XCT(IL,L,N),ZERO)
                     ELSE
                        DT=2.*DZ(IL,L)/(WC(IL,L+1)+WC(IL,L))
                        XCT(IL,L,N)=XCT(IL,L+1,N)/(1.+DT*RU(IL,L,N))
                     ENDIF
                     RU (IL,L,N) = RU(IL,L,N)*XCT(IL,L,N)
                  ENDIF
 430           CONTINUE
            ENDIF
 420     CONTINUE
         IF ( IPAM.EQ.0 ) THEN
C
C        * S(IV) OXIDATION. CHEMICAL SOURCES/SINKS ARE CONSIDERED 
C        * FOR SULPHUR DIOXIDE, HYDROGEN PEROXIDE, OZONE. THEY
C        * ARE ONLY APPLIED IF ALL OF THESE THREE TRACERS PLUS
C        * SULPHATE ARE AVAILABLE. THE CORRESPONDING SULPHATE 
C        * PRODUCTION RATE IS GIVEN BY AF1*[O3]*[SO2]+AF2*[H2O2]*[SO2]
C        * (IN KG-SULPHUR/KG-AIR/SEC), WHERE THE CONCENTRATIONS
C        * REFER TO THE TOTAL (GAS PHASE + DISSOLVED) SPECIES (IN
C        * KG-SULPHUR/KG-AIR).
C
         IF ( KCHEM ) THEN
C
C         *** HENRY'S LAW CONSTANTS AT AMBIENT TEMPERATURES.
C
            CALL OEQUIP2(ILG,ILEV,NEQP,MSG+2,1,ILGG,ACHPA,TFE)
C
C         *** MIXING RATIOS OF CHEMICAL SPECIES AT NEXT LEVEL ABOVE.
C
            DO 450 L=ILEVM,MSG+2,-1
              DO 440 IL=1,ILGG
C
C              * INITIALIZE WORK FIELDS.
C
               AMHVAL(IL)=0.
C
               IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1              .AND. ISKIP(IL).EQ.0                        ) THEN
                  ATL=PE(IL,L+1)-(ALPHA(IL,L+1)/(1.+ALPHA(IL,L+1)))
     1                           *PD(IL,L+1)
                  ATH=PE(IL,L  )-(ALPHA(IL,L  )/(1.+ALPHA(IL,L  )))
     1                           *PD(IL,L  )
                  AF1=0.
                  AF2=0.
                  DGN=DG(IL,L)
                  SO2T=MAX(XG(IL,L,ISO2)
     1                +(XCT(IL,L+1,ISO2)-XG(IL,L,ISO2))*DGN,ZERO)
                  HPOT=MAX(XG(IL,L,IHPO)
     1                +(XCT(IL,L+1,IHPO)-XG(IL,L,IHPO))*DGN,ZERO)
                  SO4T=MAX(XG(IL,L,ISO4)
     1                +(XCT(IL,L+1,ISO4)-XG(IL,L,ISO4))*DGN,ZERO)
C
C               *** VOLUME FRACTION OF CLOUD LIQUID WATER (IN M^3/M^3).
C
                  ALWCVMR=(     BETA *QLCT(IL,L+1)*RHOA(IL,L+1)
     1                     +(1.-BETA)*QLCT(IL,L  )*RHOA(IL,L  ))
     2                        /YRHOWAT
                  IF (       L.LT.LCL(IL) .AND. WC(IL,L).GT.0.
     1                 .AND. ALWCVMR.GT.YEPS ) THEN
C
C                  *** REACTION RATES AF1 AND AF2 (IN LITRE/MOL/SEC).
C                  *** THE PH AND CONCENTRATIONS RELEVANT TO AF1 AND AF2
C                  *** ARE OBTAINED FROM LINEAR INERPOLATION OF THE 
C                  *** RESULTS AT THE LEVEL BELOW AND AT THE ACTUAL
C                  *** LEVEL UNDER THE ASSUMPTION THAT IN-CLOUD OXIDATIO
C                  *** CAN BE OMITTED IN THE CALCULATION OF THE VALUES A
C                  *** THE ACTUAL LEVEL. 
C    
C                    Total (gas + aq. phase) concentrations. 
C                    Conversion kg-S/kg -> mol/l.
C
                     AMR2M =(BETA*RHOA(IL,L+1)+(1.-BETA)*RHOA(IL,L))
     1                     /(YSMASS*YCOM3L)
                     AGTSO2=AMR2M*(    BETA*XCT(IL,L+1,ISO2)
     1                           +(1.-BETA)*SO2T )
                     AGTSO4 =AMR2M*(  BETA*XCT (IL,L+1,ISO4 )
     1                           +(1.-BETA)*SO4T  )
                     AHNO3=MAX(ZERO,XGA(IL,L,IHNO3)
     1                   +(XCTA(IL,L+1,IHNO3)-XGA(IL,L,IHNO3))*DGN)
                     AGTHNO3=AMR2M*(  BETA*XCTA(IL,L+1,IHNO3)
     1                           +(1.-BETA)*AHNO3 )
                     ANH3 =MAX(ZERO,XGA(IL,L,INH3 )
     1                   +(XCTA(IL,L+1,INH3 )-XGA(IL,L,INH3 ))*DGN)
                     AGTNH3 =AMR2M*(  BETA*XCTA(IL,L+1,INH3 )
     1                           +(1.-BETA)*ANH3  )
                     ACO2 =MAX(ZERO,XGA(IL,L,ICO2 )
     1                   +(XCTA(IL,L+1,ICO2 )-XGA(IL,L,ICO2 ))*DGN)
                     AGTCO2 =AMR2M*(  BETA*XCTA(IL,L+1,ICO2 )
     1                           +(1.-BETA)*ACO2  )
C
C                    Sodium molarity (mol/l-water).
C
                     AMNA=0./ALWCVMR     
C
C                    Chlorine molarity (mol/l-water).
C
                     AMCL  =MAX(AMNA-2.E-06,ZERO) 
                     AGTHCL=AMCL*ALWCVMR 
C
C                  *** INITIAL GUESS FOR PH.
C
C                    Initial molarities in cloud water.
C
                     AMSO3=0.                  ! SO3(2-)
                     AMNH4=AGTNH3 /ALWCVMR      ! NH4(+)
                     AMNO3=AGTHNO3/ALWCVMR      ! NO3(-)
                     AMSO4=AGTSO4 /ALWCVMR      ! SO4(2-)
C
C                    Initial guess H+ molarity from initial ion 
C                    molarities.
C
                     ADELTA=AMNO3+2.*(AMSO4+AMSO3)+AMCL-AMNH4-AMNA
                     AFOH=1.E-14
                     AMH=0.5*( ADELTA+SQRT(ADELTA**2+4.*(AFOH
     1                  +ACHPA(IL,L,2)*AGTSO2+ACHPA(IL,L,5)*AGTCO2)) )
                     AMH=MAX(MIN(AMH,OPEM1),OPEM10)  
C
C                    Equilibrium parameters for soluble species.
C
                     ATVAL=ACHPA(IL,L,1)+ACHPA(IL,L,2)/AMH
     1                    +ACHPA(IL,L,3)/AMH**2
                     AFRAS=1./( 1.+ALWCVMR*ATVAL )
                     AFNH4=1./( 1.+AGTNH3/(1./ACHPA(IL,L,9)
     1                    +ALWCVMR*AMH) )
                     AFNO3=AGTHNO3/( 1./ACHPA(IL,L,8)+ALWCVMR/AMH )
C
C                    First iteration H+ molarity calculation.
C
                     AFCL =AGTHCL/( 1./ACHPA(IL,L,10)+ALWCVMR/AMH )
                     AFHSO=AGTSO2*AFRAS*ACHPA(IL,L,2)
                     AFCO3=AGTCO2*ACHPA(IL,L,5) 
                     AMSO3=( AGTSO2*AFRAS*ACHPA(IL,L,3) )/AMH**2
                     ADELTA=AFNH4*( 2.*(AMSO4+AMSO3)-AMNA )
                     AMH=0.5*( ADELTA+SQRT(ADELTA**2+4.*AFNH4*(AFOH
     1                                    +AFCO3+AFHSO+AFNO3+AFCL)) )
                     AMH=MAX(MIN(AMH,OPEM1),OPEM10)
                     ATVAL=ACHPA(IL,L,1)+ACHPA(IL,L,2)/AMH
     1                    +ACHPA(IL,L,3)/AMH**2
                     AFRAS=1./( 1.+ALWCVMR*ATVAL )
                     AFNH4=1./( 1.+AGTNH3/(1./ACHPA(IL,L,9)
     1                    +ALWCVMR*AMH) )
                     AFNO3=AGTHNO3/( 1./ACHPA(IL,L,8)+ALWCVMR/AMH )
C
C                    Second iteration H+ molarity calculation.
C
                     AFCL =AGTHCL/( 1./ACHPA(IL,L,10)+ALWCVMR/AMH )
                     AFHSO=AGTSO2*AFRAS*ACHPA(IL,L,2)
                     AFCO3=AGTCO2*ACHPA(IL,L,5) 
                     AMSO3=( AGTSO2*AFRAS*ACHPA(IL,L,3) )/AMH**2
                     ADELTA=AFNH4*( 2.*(AMSO4+AMSO3)-AMNA )
                     AMH=0.5*( ADELTA+SQRT(ADELTA**2+4.*AFNH4*(AFOH
     1                                    +AFCO3+AFHSO+AFNO3+AFCL)) )
                     AMH=MAX(MIN(AMH,OPEM1),OPEM10)
                     ATVAL=ACHPA(IL,L,1)+ACHPA(IL,L,2)/AMH
     1                    +ACHPA(IL,L,3)/AMH**2
                     AFRAS=1./( 1.+ALWCVMR*ATVAL )
C
                     AMHVAL(IL)=AMH
C
C                    Ozone oxidation parameter (in litre-air/mol/sec).
C
                     ATVO3=1./( 1.+ALWCVMR*ACHPA(IL,L,7) )
                     AF1=( ACHPA(IL,L,11)+ACHPA(IL,L,12)/AMH )*AFRAS
     1                   *ATVAL*ATVO3*ACHPA(IL,L,7)*ALWCVMR
C
C                    Hydrogen peroxide oxidation parameter 
C                    (in litre-air/mol/sec).
C
                     AFRAH=1./( 1.+ALWCVMR*ACHPA(IL,L,6) )
                     AF2=( ACHPA(IL,L,13)/(0.1+AMH) )*AFRAS
     1                    *ACHPA(IL,L,1)*AFRAH*ACHPA(IL,L,6)*ALWCVMR
C
C                  *** CONVERSION FROM LITRE/MOL/SEC TO 
C                  *** KG-AIR/KG-SULPHUR/SEC.
C
                     AF1=AF1*AMR2M
                     AF2=AF2*AMR2M
C
C                  *** SCAVENGING RATIOS.
C
                     ASRSO2=MAX( 1.-AFRAS, ZERO )
                     ASRHPO=MAX( 1.-AFRAH, ZERO )
                     ASRO3 =MAX( 1.-ATVO3, ZERO )
                     ASRNH3=MAX( ALWCVMR*AMH*(1./AFNH4 - 1.)  
     1                     /MAX(AGTNH3,YEPS), ZERO )                      
                     ASRHNO=MAX( ALWCVMR*AFNO3 
     1                   /( MAX(AGTHNO3,YEPS)*AMH ), ZERO )               
C
C                  *** PRELIMINARY SCAVENGING RATES IN 1/S.
C
                     APA=BETA*QLCT(IL,L+1)+(1.-BETA)*QLCT(IL,L)
                     APB=BETA*AUTO(IL,L+1)+(1.-BETA)*AUTO(IL,L)
C
                     RU (IL,L,ISO4) = ASRSO4*APB/APA
                     RU (IL,L,ISO2) = ASRSO2*APB/APA
                     RU (IL,L,IHPO) = ASRHPO*APB/APA
                     RUA(IL,L,INH3) = ASRNH3*APB/APA
                     RUA(IL,L,IHNO3)= ASRHNO*APB/APA
                     RUA(IL,L,ICO2) = ASRCO2*APB/APA
                     RUA(IL,L,IO3)  = ASRO3 *APB/APA
                  ENDIF
                  APA=.5*(WC(IL,L+1)*ATL+WC(IL,L)*ATH)
                  APB=1.-DG(IL,L)
                  DT=2.*DZ(IL,L)/(WC(IL,L+1)+WC(IL,L))
                  IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN
                     XCTA(IL,L,INH3 )=APA*MAX(ZERO,XGA(IL,L,INH3 )
     1                   +(XCTA(IL,L+1,INH3 )-XGA(IL,L,INH3 ))*DGN)
     2                   /(APA+APB*RUA(IL,L,INH3 ))
                     XCTA(IL,L,IHNO3)=APA*MAX(ZERO,XGA(IL,L,IHNO3)
     1                   +(XCTA(IL,L+1,IHNO3)-XGA(IL,L,IHNO3))*DGN)
     2                   /(APA+APB*RUA(IL,L,IHNO3))
                     XCTA(IL,L,ICO2 )=APA*MAX(ZERO,XGA(IL,L,ICO2 )
     1                   +(XCTA(IL,L+1,ICO2 )-XGA(IL,L,ICO2 ))*DGN)
     2                   /(APA+APB*RUA(IL,L,ICO2 ))
                     XCTA(IL,L,IO3  )=APA*MAX(ZERO,XGA(IL,L,IO3  )
     1                   +(XCTA(IL,L+1,IO3  )-XGA(IL,L,IO3  ))*DGN)
     2                   /(APA+APB*RUA(IL,L,IO3  ))
                  ELSE
                     XCTA(IL,L,INH3 )=XCTA(IL,L+1,INH3 )
     1                               /(1.+DT*RUA(IL,L,INH3 ))
                     XCTA(IL,L,IHNO3)=XCTA(IL,L+1,IHNO3)
     1                               /(1.+DT*RUA(IL,L,IHNO3))
                     XCTA(IL,L,ICO2 )=XCTA(IL,L+1,ICO2 )
     1                               /(1.+DT*RUA(IL,L,ICO2 ))
                     XCTA(IL,L,IO3  )=XCTA(IL,L+1,IO3  )
     1                               /(1.+DT*RUA(IL,L,IO3  ))
                  END IF
C
C               *** SEMI-IMPLICIT SOLUTION OF FULL IN-CLOUD BUDGET FOR
C               *** GIVEN AF1 AND AF2. IT IS ASSUMED THAT OZONE 
C               *** CONCENTRATIONS DO NOT CHANGE DUE TO CHEMICAL 
C               *** REACTIONS OWING TO THE RELATIVELY LOW MAGNITUDE OF 
C               *** IN-CLOUD S(IV) OXIDATION. 
C
                  IF ( (AF1*XCTA(IL,L,IO3)*SO2T+AF2*HPOT*SO2T).GT.0.
     1               ) THEN
                     IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN
                        AP1=APB*(AF1*XCTA(IL,L,IO3)+RU(IL,L,ISO2))+APA
                        AP2=APB*AF2
                        AP3=AP1+AP2*(HPOT-SO2T)
                        AP4=4.*SO2T*AP1*AP2
                        AP5=APA/(2.*AP1*AP2)
                        AP6=(APA+APB*(AF1*XCTA(IL,L,IO3)+RU(IL,L,ISO2)
     1                     +AF2*(HPOT-SO2T)))*AP5
                        XCT(IL,L,ISO2)=-AP6+ABS(AP5)*SQRT(AP3**2+AP4)
                        XCT(IL,L,ISO2)=MAX(XCT(IL,L,ISO2),ZERO)
                        XCT(IL,L,IHPO)=APA*HPOT/(APA
     1                      +APB*(AF2*XCT(IL,L,ISO2)+RU(IL,L,IHPO)))
                        XCT(IL,L,IHPO)=MAX(XCT(IL,L,IHPO),ZERO)
                        XCT(IL,L,ISO4)=( APB*( XG(IL,L,ISO4)
     1                     +(AF1*XCTA(IL,L,IO3 )+AF2*XCT(IL,L,IHPO)
     2                      )*XCT(IL,L,ISO2)/APA )
     3                      +DG(IL,L)*XCT(IL,L+1,ISO4) )
     4                      /( 1.+APB*RU(IL,L,ISO4)/APA )
                        XCT(IL,L,ISO4)=MAX(XCT(IL,L,ISO4),ZERO)
                     ELSE
                        AP1=0.5*( 1.+DT*(RU(IL,L,ISO2)
     1                     +AF1*XCTA(IL,L,IO3 )+AF2*(XCT(IL,L+1,IHPO)
     2                      -XCT(IL,L+1,ISO2))) )
     3                      /(DT*AF2*(DT*(RU(IL,L,ISO2)
     4                       +AF1*XCTA(IL,L,IO3))+1.))
                        AP2=-XCT(IL,L+1,ISO2)
     1                     /(DT*AF2*(DT*(RU(IL,L,ISO2)
     2                       +AF1*XCTA(IL,L,IO3))+1.))
                        XCT(IL,L,ISO2)=-AP1+SQRT(AP1**2-AP2)
                        XCT(IL,L,ISO2)=MAX(XCT(IL,L,ISO2),ZERO)
                        XCT(IL,L,IHPO)=XCT(IL,L+1,IHPO)
     1                                /(1.+DT*(RU(IL,L,IHPO)
     2                                +AF2*XCT(IL,L,ISO2)))
                        XCT(IL,L,IHPO)=MAX(XCT(IL,L,IHPO),ZERO)
                        XCT(IL,L,ISO4)=( XCT(IL,L+1,ISO4)
     1                      +DT*((AF1*XCTA(IL,L,IO3)+AF2*XCT(IL,L,IHPO))
     2                         *XCT(IL,L,ISO2)) )
     3                      /( 1.+DT*RU(IL,L,ISO4) )
                        XCT(IL,L,ISO4)=MAX(XCT(IL,L,ISO4),ZERO)
                     ENDIF
C
C                  *** FINAL MEAN CHEMCIAL SINKS AND SOURCES IN THE CLOU
C                  *** (IN KG/KG/SEC).
C
                     SC(IL,L,ISO2)=-( AF1*XCTA(IL,L,IO3 )
     1                               +AF2*XCT (IL,L,IHPO) )
     2                              *XCT(IL,L,ISO2)
                     SC(IL,L,IHPO)=-AF2*XCT(IL,L,IHPO)*XCT(IL,L,ISO2)
                     SC(IL,L,ISO4)=-SC(IL,L,ISO2)
                  ELSE
                     IF ( ATL*ATH.GT.0. .AND. DG(IL,L).NE.1.) THEN 
                        XCT(IL,L,ISO2)=APA*SO2T/(APA+APB*RU(IL,L,ISO2))
                        XCT(IL,L,IHPO)=APA*HPOT/(APA+APB*RU(IL,L,IHPO))
                        XCT(IL,L,ISO4)=APA*SO4T/(APA+APB*RU(IL,L,ISO4))
                     ELSE
                        XCT(IL,L,ISO2)=XCT(IL,L+1,ISO2)
     1                                /(1.+DT*RU(IL,L,ISO2))
                        XCT(IL,L,IHPO)=XCT(IL,L+1,IHPO)
     1                                /(1.+DT*RU(IL,L,IHPO))
                        XCT(IL,L,ISO4)=XCT(IL,L+1,ISO4)
     1                                /(1.+DT*RU(IL,L,ISO4))
                     ENDIF
                  ENDIF
                  RU (IL,L,ISO2) = RU(IL,L,ISO2)*XCT(IL,L,ISO2)
                  RU (IL,L,IHPO) = RU(IL,L,IHPO)*XCT(IL,L,IHPO)
                  RU (IL,L,ISO4) = RU(IL,L,ISO4)*XCT(IL,L,ISO4)
               ENDIF
  440        CONTINUE
C
             IF(ISVCHEM.NE.0)                       THEN  
              DO 445 IL=1,ILGG
                IF ( AMHVAL(IL) > 1.E-20 ) THEN
                  PHSROW(ISHALL(IL),L)=-LOG10(AMHVAL(IL))
                ELSE
                  PHSROW(ISHALL(IL),L)=0.
                ENDIF
  445         CONTINUE
             ENDIF
  450      CONTINUE
         ENDIF
         ENDIF
         DO 465 N=1,NTRAC 
            IF ( ITRPHS(N).GT.0 ) THEN
               DO 460 L=ILEV,MSG+2,-1 
               DO 460 IL=1,ILGG 
                  IF (       L.LE.(LCL(IL)+1) .AND. L.GE.(LNB(IL)-1)
     1                 .AND. ISKIP(IL).EQ.0  ) THEN
                     XC(IL,L,N)=( XCT(IL,L,N)+ALPHA(IL,L)*XF(IL,L,N) )
     1                         /( 1.+ALPHA(IL,L) )
                  ENDIF
 460           CONTINUE
            ENDIF
 465     CONTINUE
      ENDIF
C
C     * RESET IF THE TOP OF THE CLOUD IS ABOVE THE FREEZING LEVEL.
C
      DO 600 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 ) THEN
           IF ( TFE(IL,LNB(IL)-1).LT.273.  )     ISKIP(IL)=5
        ENDIF
 600  CONTINUE
C
C     * TIME AT WHICH THE CLOUD TOP REACHES A CERTAIN LEVEL.
C
      DO 620 L=MSG+1,ILEV
      DO 620 IL=1,ILGG
         TSTAR(IL,L)=0.
 620  CONTINUE
      IF ( KCTDET ) THEN
         DO 640 L=ILEVM,MSG+2,-1 
         DO 640 IL=1,ILGG 
            IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1           .AND. ISKIP(IL).EQ.0  ) THEN
               TSTAR(IL,L)=TSTAR(IL,L+1)+2.*DZ(IL,L)
     1                             /( WC(IL,L)+WC(IL,L+1) )
            ENDIF
 640     CONTINUE
      ELSE
         DO 660 L=ILEVM,MSG+2,-1 
         DO 660 IL=1,ILGG 
            IF (       L.LE.LCL(IL) .AND. L.GE.(LNB(IL)-1)
     1           .AND. ISKIP(IL).EQ.0  ) THEN
               TSTAR(IL,L)=TSTAR(IL,L+1)+2.*DZ(IL,L)
     1                          /( WC(IL,L  )*(1.+ALPHA(IL,L  ))
     2                            +WC(IL,L+1)*(1.+ALPHA(IL,L+1)) )
            ENDIF
 660     CONTINUE
      END IF
C
C     * CLOUD LIFE TIME.
C
      DO 700 IL=1,ILGG 
         IF ( ISKIP(IL).EQ.0 ) THEN
            L=LNB(IL)-1
            TAUC(IL)=TSTAR(IL,L)
         ENDIF
 700  CONTINUE
C
C     * DIMENSIONLESS UPDRAFT MASS FLUX.
C
      DO 710 L=MSG+1,ILEV
      DO 710 IL=1,ILGG
         MC (IL,L)=0.
         WRK(IL,L)=0.
 710  CONTINUE
      DO 720 L=MSG+2,ILEV
      DO 720 IL=1,ILGG
         IF ( ISKIP(IL).EQ.0 .AND. L.LE.LMAG(IL) 
     1                       .AND. L.GT.LCL(IL)                 ) THEN
            MC(IL,L)=1.
            WRK(IL,L)=1.
         END IF 
 720  CONTINUE
      DO 740 L=ILEVM,MSG+2,-1
      DO 740 IL=1,ILGG 
         IF (       L.LE.LCL(IL) .AND. L.GE.LNB(IL) 
     1        .AND. ISKIP(IL).EQ.0 .AND. WC(IL,L).GT.0. ) THEN
            ATL=PE(IL,L+1)-PD(IL,L+1)
            ATH=PE(IL,L  )-PD(IL,L  )
            DL=ATH-ATL
            DGT=( ATL*DZ(IL,L)-DL*ZF(IL,L+1) )
     1          +.5*DL*( ZF(IL,L)**2-ZF(IL,L+1)**2 )/DZ(IL,L)
            WRK(IL,L)=WRK(IL,L+1)*EXP(DGT)
            MC(IL,L)=WRK(IL,L)*( 1.-TSTAR(IL,L)/TAUC(IL) )
         ENDIF
 740  CONTINUE
      DO 760 IL=1,ILGG 
         L=LNB(IL)-1
         IF (       L.LE.LCL(IL) .AND. ISKIP(IL).EQ.0 
     1        .AND. WC(IL,L).GT.0. ) THEN
            ATL=PE(IL,L+1)-PD(IL,L+1)
            ATH=PE(IL,L  )-PD(IL,L  )
            DL=ATH-ATL
            DGT=( ATL*DZ(IL,L)-DL*ZF(IL,L+1) )
     1          +.5*DL*( ZF(IL,L)**2-ZF(IL,L+1)**2 )/DZ(IL,L)
            WRK(IL,L)=WRK(IL,L+1)*EXP(DGT)
         END IF
 760  CONTINUE
C
C     * DETRAINMENT AND ENTRAINMENT RATES (FOR DIAGNOSTIC PURPOSES).
C
      DO 780 L=MSG+1,ILEV
      DO 780 IL=1,ILGG
         EC(IL,L)=0.
         DC(IL,L)=0.
 780  CONTINUE
      IF ( KCTDET ) THEN
         DO 800 L=MSG+2,ILEV
         DO 800 IL=1,ILGG  
            IF (       L.LE.(LCL(IL)+1) .AND. L.GE.(LNB(IL)-1)
     1           .AND. ISKIP(IL).EQ.0 .AND. WC(IL,L).GT.0. ) THEN
               DC(IL,L)=WRK(IL,L)*( 1.+ALPHA(IL,L) )
     1                 /( WC(IL,L)*TAUC(IL) )  +MC(IL,L)*PD(IL,L)
               EC(IL,L)=WRK(IL,L)*( ALPHA(IL,L)
     1                 /( WC(IL,L)*TAUC(IL) ) )+MC(IL,L)*PE(IL,L)
            ENDIF
 800     CONTINUE
      ELSE
         DO 820 L=MSG+2,ILEV
         DO 820 IL=1,ILGG  
            IF (       L.LE.(LCL(IL)+1) .AND. L.GE.(LNB(IL)-1)
     1           .AND. ISKIP(IL).EQ.0 .AND. WC(IL,L).GT.0. ) THEN
               DC(IL,L)=WRK(IL,L)/( WC(IL,L)*TAUC(IL) )
     1                 +MC(IL,L)*PD(IL,L) 
               EC(IL,L)=WRK(IL,L)*( ALPHA(IL,L)/( WC(IL,L)*TAUC(IL) 
     1                 *(1.+ALPHA(IL,L)) ) )+MC(IL,L)*PE(IL,L)
            ENDIF
 820     CONTINUE
      END IF
C
      RETURN 
      END
