      SUBROUTINE CONTRA9 (X,XU,XD,DXDT,DUDT,DVDT,XGA,SUMDX,CLDCV,WRK,
     1                    IDEEP,DUO,EUO,EDO,MU,MD,MDO,QL,Q,T,P,PFG,DZ,
     2                    DP,PRESSG,MB,EPS0,JB,JT,MX,JD,ITRPHS,ITRCCH,
     3                    XSRC,XSNK,PHDROW,DSR,NTRAC,
     4                    NTRACA,IU,IV,INH3,IHNO3,
     6                    IO3,ICO2,ISO2,ISO4,IHPO,ISSA,ISSC,IDUA,IDUC, 
     7                    C0FAC,EPS1,RGAS,GRAV,ROG,ILG,ILEV,IL1G,IL2G,
     8                    MSG,DT,NEQP,ISVCHEM,IPAM)
C-----------------------------------------------------------------------
C     * FEB 01/2018 - M.LAZARE.    BUGFIXES TO HANDLE CASES WHERE EUO~DUO. 
C     * FEB 10/2015 - K.VONSALZEN. New routine for GCM18:
C     *                            - Scavenging now allowed for SO2.
C     * JUN 20/2013 - K.VONSALZEN/ PREVIOUS VERSION CONTRA8 FOR GCM17:
C     *               M.LAZARE.    - ADD SUPPORT FOR PLA.
C     *                            - COSMETIC CLEANUP OF LOOP NUMBERS.
C     *                            - REMOVE "ISULF", INCLUDING CALL.
C     * APR 30/2012 - M.LAZARE.    PREVIOUS VERSION CONTRA7 FOR GCM16:
C     *                            - CCU,CCD DECREASED FROM 0.55 TO ZERO.
C     *                            - CALLS NEW XTEND5.
C     * APR 23/2010 - M.LAZARE/    PREVIOUS VERSION CONTRA6 FOR GCM15I:
C     *               K.VONSALZEN. - CALLS NEW XTEND4.
C     *                            - EXPLICIT SOURCE/SINK TERMS.
C     *                            - REVISED DIAGNOSTIC FIELDS.
C     * JUL 15/2009 - M.LAZARE/    BUGFIX TO TRACER DIMENSION FOR
C     *               K.VONSALZEN. INPUT ARRAY "X" (NTRAC INSTEAD OF 
C     *                            NTRACA). COSMETIC ONLY FOR NON-SCM
C     *                            BECAUSE WE PASS BY ADDRESS NOT VALUE.
C     * FEB 16/2009 - M.LAZARE/    PREVIOUS VERSION CONTRA5 FOR GCM15H:
C     *               K.VONSALZEN. ACCURACY RESTRICTION ON CALCULATION
C     *                            OF FACV TO HANDLE 0/0 UNUSUAL CASE. 
C     * DEC 18/2007 - M.LAZARE/    PREVIOUS VERSION CONTRA4 FOR GCM15G:
C     *               K.VONSALZEN. ADD CALCULATIONS FOR DIAGNOSTIC FIELD
C     *                            WDD4ROW.
C     * NOV 23/2006 - M.LAZARE/    PREVIOUS VERSION CONTRA3 FOR GCM15F:
C     *               K.VONSALZEN. - VARIOUS CHANGES TO AVOID PATHALOGICAL
C     *                              CASES FOR SMALL VALUES IN 32-BIT 
C     *                              MODE, PARTICULARILY IN DEFINITION OF
C     *                              "ATMP". 
C     * JUL 30/2006 - M.LAZARE.    - C0FAC NOW PASSED IN FROM CONVECTION
C     *                              DRIVER, TO ENSURE CONSISTENT USEAGE
C     *                              IN CLDPRP AND CONTRA.
C     * JUN 15/2006 - M.LAZARE.    - "ZERO","ONE","OPEM1","OPEM10" 
C     *                              DATA CONSTANTS ADDED AND USED IN 
C     *                              CALLS TO INTRINSICS. 
C     *                            - NOW USES INTERNAL WORK ARRAYS
C     *                              INSTEAD OF PASSED WORKSPACE.    
C     *                            - CALLS REVISED XTEND3.
C     * DEC 13/2005 - K. VON SALZEN. PREVIOUS VERSION CONTRA2 FOR GCM15D/E:
C     *                              - ADDS CODE FOR MOMENTUM, INCLUDING
C     *                                CALLING REVISED XTEND2 AND PASSING
C     *                                IN EDO FROM NEW CLDPRP9. 
C     * MAY 14/2004 - M.LAZARE/      PREVIOUS VERSION CONTRA FOR GCM15B/C.
C     *               K. VON SALZEN.    
C     * SEP 18/2000 - K. VON SALZEN. ORIGINAL VERSION BASED ON NARCM CODE.
C-----------------------------------------------------------------------
      IMPLICIT REAL   (A-H,O-Z),
     +INTEGER (I-N)
C   
C     * OUTPUT ARRAYS.
C
      REAL  , DIMENSION(ILG,ILEV,NTRAC)  :: XU,XD,DXDT
      REAL  , DIMENSION(ILG,NTRAC)       :: SUMDX,XSRC,XSNK
      REAL  , DIMENSION(ILG,ILEV)        :: DUDT,DVDT,PHDROW
      REAL  , DIMENSION(ILG)             :: WRK
C
C     * INPUT TRACER FIELDS.
C
      REAL  , DIMENSION(ILG,ILEV,NTRAC ) :: X
      REAL  , DIMENSION(ILG,ILEV,NTRACA) :: XGA 
C
C     * OTHER INPUT FIELDS.
C 
      REAL  , DIMENSION(ILG,ILEV)        :: DUO,EUO,DSR,EDO,PFG,MDO,QL,
     1                                      Q,P,DZ,DP,CLDCV,MD,T,MU
      REAL  , DIMENSION(ILG)             :: MB,EPS0,PRESSG
      INTEGER,DIMENSION(ILG)             :: JB,JT,MX,JD,IDEEP
      INTEGER,DIMENSION(NTRAC)           :: ITRPHS,ITRCCH
C
C     * INTERNAL WORK ARRAYS.
C
      REAL  , DIMENSION(ILG,ILEV,NTRAC)  :: SC,RU,XF,WRKX
      REAL  , DIMENSION(ILG,ILEV,NTRACA) :: XUA,XDA,RUA,XFA,SCA,WRKXA
      REAL  , DIMENSION(ILG,ILEV,NEQP)   :: ACHPA
      REAL  , DIMENSION(ILG,ILEV)        :: ALWCVMR,WRKS  
      REAL  , DIMENSION(ILG)             :: AMHVAL,FACV,ATMPV,
     1                                      C0V,MBC,AMBCOR
C
      LOGICAL KCHEM,KRIT1,KRIT2
C
C     * PARAMETERS.
C
      PARAMETER (YCOM3L = 1.E+03  )
      PARAMETER (YEPS   = 1.E-20  )
      PARAMETER (YRHOW  = 1.E+03  )
      PARAMETER (YSMALL = 1.E-33  )
      PARAMETER (YSMASS =32.06E-03)
      PARAMETER (BETA   =.5       )
      PARAMETER (CCU=0.)
      PARAMETER (CCD=0.)
C
C     * COMMON DECLARATIONS.
C
      COMMON /HTCP  / TFREEZ,T2S,AI,BI,AW,BW,SLP
C
      DATA ZERO    / 0.     /
      DATA ONE     / 1.     /
      DATA OPEM1   / 1.E-1  /
      DATA OPEM10  / 1.E-10 /
C-----------------------------------------------------------------------
      KCHEM=.TRUE.
      ASRSO4 = 0.9                                      
      ASRPHOB= 0.
C
C     * INVOKE CFL TO LIMIT CLOUD BASE MASS FLUX. THIS ENSURES
C     * STABILITY OF THE EXPLICIT CALCULATIONS OF TRACER TENDENCIES
C     * IN SUBROUTINE XTEND3.
C
      DO 10 IL=IL1G,IL2G                                                
        MBC(IL)=MB(IL)
        AMBCOR(IL)=0.
   10 CONTINUE
      DO 15 L=MSG+2,ILEV                                               
      DO 15 IL=IL1G,IL2G                                                
        IF((MU(IL,L)+MD(IL,L)).NE.0.)                        THEN    
          MBC(IL)=MIN( MBC(IL) , MB(IL)*MIN(DP(IL,L  )/DSR(IL,L  ),
     1                                      DP(IL,L-1)/DSR(IL,L-1))
     1                         /(DT*ABS(MU(IL,L)+MD(IL,L))) )                
        ENDIF                                                           
   15 CONTINUE                                        
      DO 20 IL=IL1G,IL2G                                                
        IF ( MB(IL).GT.0. ) THEN
             AMBCOR(IL)=MBC(IL)/MB(IL)
        ENDIF
   20 CONTINUE
C
C     * INITIALIZATION OF WORK FIELDS.
C     * INITIALIZE UPDRAFT AND DOWNDRAFT PROPERTIES.
C
      DO 25 N=1,NTRAC
      DO 25 L=MSG+1,ILEV                                       
      DO 25 IL=IL1G,IL2G                                       
         RU(IL,L,N) = 0.                                      
         SC(IL,L,N) = 0.                                      
         XU(IL,L,N) = X(IL,L,N)
         XD(IL,L,N) = X(IL,L,N)
   25 CONTINUE   
C
C     * SET CLOUD BASE UPDRAFT PROPERTIES.
C
      DO 30 N=1,NTRAC
      DO 30 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.) THEN
          XU(IL,JB(IL),N) = X(IL,MX(IL),N)
        ENDIF
   30 CONTINUE
C
C     * CALCULATE UPDRAFT PROPERTIES FOR PASSIVE TRACERS.
C
      DO N=1,NTRAC
       IF(ITRPHS(N).GT.0 .AND. ITRCCH(N).EQ.0) THEN    
         DO J=ILEV-1,MSG+1,-1 
         DO IL=IL1G,IL2G
          IF(EPS0(IL).GT.0. .AND. (J.GT.JT(IL) .AND. J.LT.JB(IL)) .AND.
     1       MU(IL,J+1).NE.0. )THEN 
             IF ( EUO(IL,J).LE.10.*ABS(EUO(IL,J)-DUO(IL,J))) THEN
                AEXP=EUO(IL,J)/(DUO(IL,J)-EUO(IL,J))
                ATMP=(ABS(MU(IL,J)/MU(IL,J+1)))**AEXP
             ELSE
                ATMP=EXP(-MB(IL)*EUO(IL,J)*DP(IL,J)/MU(IL,J+1))
             END IF
             XU(IL,J,N)=X(IL,J,N)+(XU(IL,J+1,N)-X(IL,J,N))*ATMP
          ENDIF
         ENDDO
         ENDDO
       ENDIF
      ENDDO
C
      IF ( IPAM.EQ.0 ) THEN
        DO L=MSG+1,ILEV                                       
        DO IL=IL1G,IL2G                                       
          ALWCVMR(IL,L) = 0.                                             
        ENDDO
        ENDDO
      ENDIF
C
      DO 40 N=1,NTRACA                                         
      DO 40 L=MSG+1,ILEV                                       
      DO 40 IL=IL1G,IL2G                                       
        RUA(IL,L,N) = 0.                                      
        XUA(IL,L,N) = XGA(IL,L,N)                             
        XDA(IL,L,N) = XGA(IL,L,N)                             
        SCA(IL,L,N) = 0.
   40 CONTINUE                                                 
      DO 160 N=1,NTRACA
      DO 160 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.) THEN
           XUA(IL,JB(IL),N) = XGA(IL,MX(IL),N)
        ENDIF
  160 CONTINUE
      DO 161 N=IU,IV
      DO 161 J=ILEV-1,MSG+1,-1 
      DO 161 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0. .AND. (J.GT.JT(IL) .AND. J.LT.JB(IL)) .AND.
     1     MU(IL,J+1).NE.0. )THEN 
C
C          * UNPERTURBED WIND IN UPDRAFT.
C
            IF ( EUO(IL,J).LE.10.*ABS(EUO(IL,J)-DUO(IL,J))) THEN
              AEXP=EUO(IL,J)/(DUO(IL,J)-EUO(IL,J))
              ATMP=(ABS(MU(IL,J)/MU(IL,J+1)))**AEXP
           ELSE
              ATMP=EXP(-MB(IL)*EUO(IL,J)*DP(IL,J)/MU(IL,J+1))
           END IF
           XUA(IL,J,N)=XGA(IL,J,N)+(XUA(IL,J+1,N)-XGA(IL,J,N))*ATMP
        ENDIF
  161 CONTINUE
      DO 162 N=IU,IV
      DO 162 J=ILEV-1,MSG+1,-1 
      DO 162 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0. .AND. (J.GT.JT(IL) .AND. J.LT.JB(IL)) .AND.
     1     MU(IL,J+1).NE.0. )THEN 
            IF ( EUO(IL,J).LE.10.*ABS(EUO(IL,J)-DUO(IL,J))) THEN
              AEXP=EUO(IL,J)/(DUO(IL,J)-EUO(IL,J))
              ATMP=(ABS(MU(IL,J)/MU(IL,J+1)))**AEXP
           ELSE
              ATMP=EXP(-MB(IL)*EUO(IL,J)*DP(IL,J)/MU(IL,J+1))
           END IF
C
C         * WIND IN UPDRAFT FROM PRESSURE GRADIENT FORCE PERTURBATION.
C
          XUA(IL,J,N)=CCU*XGA(IL,J,N)+(1.-CCU)*XUA(IL,J,N)
C  
C         * DIAGNOSE PRESSURE GRADIENT FORCE TERM UNDER THE ASSUMPTION
C         * OF A CONSTANT SOURCE/SINK IN EACH GRID CELL.
C
          IF ( ABS(1.-ATMP).GT.1.E-02 ) THEN
            SCAT=EUO(IL,J)*(XUA(IL,J,N)-XGA(IL,J,N)
     1          -ATMP*(XUA(IL,J+1,N)-XGA(IL,J,N)))/(1.-ATMP)
          ELSE
            AMUO=MU(IL,J)/MB(IL)
            SCAT=AMUO*(XUA(IL,J,N)-XUA(IL,J+1,N))/DZ(IL,J)
          ENDIF
          SCA(IL,J,N)=SCAT*MBC(IL)*DZ(IL,J)*DSR(IL,J)/DP(IL,J)
        ENDIF
  162 CONTINUE
C
      IF ( KCHEM.AND.IPAM.EQ.0 ) THEN
C
C       *** COMPUTE HENRY'S LAW CONSTANTS.
C
        CALL OEQUIP2(ILG,ILEV,NEQP,MSG+1,IL1G,IL2G,ACHPA,T)               
C                                                                       
C       *** LIQUID WATER VOLUME MIXING RATIOS.
C
        DO 170 L=ILEV,MSG+2,-1                                            
        DO 170 IL=IL1G,IL2G                                               
         ZRHOL=100.*P(IL,L  )/(RGAS*T(IL,L  ))
         ZRHOH=100.*P(IL,L-1)/(RGAS*T(IL,L-1))
         IF(EPS0(IL).GT.0. .AND. (L.GT.JT(IL)                           
     1     .AND. L.LT.JB(IL)) .AND. MU(IL,L).GT.0.                      
     2     .AND. (QL(IL,L) .GT. 1.E-09) 
     3     .AND. (CLDCV(IL,L).GT.1.E-06) ) THEN                         
            ALWCVMR(IL,L) = QL(IL,L)*.5*(ZRHOL+ZRHOH)/YRHOW           
         END IF                                                         
  170   CONTINUE   
      ENDIF
C
C     * UPDRAFT MIXING RATIOS OF CHEMICAL SPECIES.
C
      DO 260 L=ILEV-1,MSG+2,-1                                         
        DO 240 IL=IL1G,IL2G
          AMHVAL(IL)=0.
          ATMPV(IL)=0.
          FACV(IL)=0.
          C0V(IL)=0. 
          IF(      L.GT.JT(IL) .AND. L.LT.JB(IL) .AND. EPS0(IL).GT.0.  
     1      .AND. MU(IL,L).GT.0. ) THEN
            IF ( EUO(IL,L).LE.10.*ABS(EUO(IL,L)-DUO(IL,L))) THEN
               AEXP=EUO(IL,L)/(DUO(IL,L)-EUO(IL,L))
               ATMP=(ABS(MU(IL,L)/MU(IL,L+1)))**AEXP
            ELSE
               ATMP=EXP(-MB(IL)*EUO(IL,L)*DP(IL,L)/MU(IL,L+1))
            END IF
            AFAKT=0.
            ATMPV(IL)=ATMP
            IF ( EUO(IL,L) .NE. 0. ) THEN 
              FACV(IL)=(1.-ATMP)/EUO(IL,L)
            ELSE
              FACV(IL)=0.
            ENDIF
            C0V(IL)=0.
            IF ( IPAM.EQ.0 )                                        THEN
             SO2T=MAX(X(IL,L,ISO2)+(XU(IL,L+1,ISO2)-X(IL,L,ISO2))
     1                                                      *ATMP,ZERO)
             SO4T=MAX(X(IL,L,ISO4)+(XU(IL,L+1,ISO4)-X(IL,L,ISO4))
     1                                                      *ATMP,ZERO)
             HPOT=MAX(X(IL,L,IHPO)+(XU(IL,L+1,IHPO)-X(IL,L,IHPO))
     1                                                      *ATMP,ZERO)
             AFAKT=CLDCV(IL,L)*DP(IL,L)/(DSR(IL,L)*MAX(MBC(IL),YSMALL)
     1                                  *DZ(IL,L))
             AF1=0.
             AF2=0.
             IF ( ALWCVMR(IL,L).GT.YEPS .AND. KCHEM ) THEN
C
C              *** REACTION RATES AF1 AND AF2 (IN LITRE/MOL/SEC) AND
C              *** SCAVENGING RATES. CLOUD WATER PH AND CONCENTRATIONS 
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 OXIDATION
C              *** CAN BE OMITTED IN THE CALCULATION OF THE VALUES AT
C              *** THE ACTUAL LEVEL IN A FIRST APPROXIMATION. 
C    
C              Total (gas + aq. phase) concentrations. 
C              Conversion kg-S/kg -> mol/l.
C  
               ZRHOH=100.*P(IL,L  )/(RGAS*T(IL,L  ))
               ZRHOL=100.*P(IL,L+1)/(RGAS*T(IL,L+1))
               AMR2M =(BETA*ZRHOL+(1.-BETA)*ZRHOH)
     1               /(YSMASS*YCOM3L)
               AGTSO2 =AMR2M*( BETA*XU (IL,L+1,ISO2 )+(1.-BETA)*SO2T  )
               AGTSO4 =AMR2M*( BETA*XU (IL,L+1,ISO4 )+(1.-BETA)*SO4T  )
               AHNO3=XGA(IL,L,IHNO3)+(XUA(IL,L+1,IHNO3)-XGA(IL,L,IHNO3)
     1                               )*ATMP
               AGTHNO3=AMR2M*( BETA*XUA(IL,L+1,IHNO3)+(1.-BETA)*AHNO3 )
               ANH3 =XGA(IL,L,INH3 )+(XUA(IL,L+1,INH3 )-XGA(IL,L,INH3 )
     1                               )*ATMP
               AGTNH3 =AMR2M*( BETA*XUA(IL,L+1,INH3 )+(1.-BETA)*ANH3  )
               ACO2 =XGA(IL,L,ICO2 )+(XUA(IL,L+1,ICO2 )-XGA(IL,L,ICO2 )
     1                               )*ATMP
               AGTCO2 =AMR2M*( BETA*XUA(IL,L+1,ICO2 )+(1.-BETA)*ACO2  )
C
C              Sea salt species.
C
               ASSMR=XU(IL,L+1,ISSA)+XU(IL,L+1,ISSC)
               ACONVF=(BETA*ZRHOL+(1.-BETA)*ZRHOH)
     1               /(ALWCVMR(IL,L)*YCOM3L)
C
C            *** INITIAL GUESS FOR PH.
C
C              Initial molarities in cloud water.
C   
               AMSO3=0.                         ! SO3(2-)
               AMNH4=AGTNH3 /ALWCVMR(IL,L)      ! NH4(+)
               AMNO3=AGTHNO3/ALWCVMR(IL,L)      ! NO3(-)
               AMSO4=AGTSO4 /ALWCVMR(IL,L)      ! SO4(2-)
C
C              Calculate molarities of sea salt species assuming a typical
C              composition of sea water. No anions other than sulphate and
C              chloride are considered. All cations, except H+, are lumped
C              together. The pH of the sea water is assumed to be 8.2.
C              H+ is accounted for in the ion balance and in the calculation
C              of hydrogen chlorine in the gas phase.
C
               AMSSSO4=0.077*ACONVF/96.058E-03*ASSMR
               AMSSCL =0.552*ACONVF/35.453E-03*ASSMR
               AMSSANI=0.
               AMSSCAT=(AMSSCL+2.*AMSSSO4+AMSSANI)*(1.-1.048E-08)
               AMSO4=AMSO4+AMSSSO4
               AMCL=AMSSCL
               AMHCL=AMSSCL*6.31E-09/ACHPA(IL,L,10)
               AGTHCL=(AMSSCL+AMHCL)*ALWCVMR(IL,L)
C              AGTHCL=AMSSCL*ALWCVMR(IL,L)   
C
C              Initial guess H+ molarity from initial ion molarities.
C   
               ADELTA=AMNO3+2.*(AMSO4+AMSO3)+AMCL+AMSSANI-AMNH4
     1               -AMSSCAT
               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+ACHPA(IL,L,3)
     1              /AMH**2
               AFRAS=1./( 1.+ALWCVMR(IL,L)*ATVAL )
               AFNH4=1./( 1.+AGTNH3/(1./ACHPA(IL,L,9)
     1                      +ALWCVMR(IL,L)*AMH) )
               AFNO3=AGTHNO3/( 1./ACHPA(IL,L,8)+ALWCVMR(IL,L)/AMH )
C
C              First iteration H+ molarity calculation.
C
               AFCL =AGTHCL/( 1./ACHPA(IL,L,10)+ALWCVMR(IL,L)/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)+AMSSANI-AMSSCAT)
               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+ACHPA(IL,L,3)
     1              /AMH**2
               AFRAS=1./( 1.+ALWCVMR(IL,L)*ATVAL )
               AFNH4=1./( 1.+AGTNH3/(1./ACHPA(IL,L,9)
     1                      +ALWCVMR(IL,L)*AMH) )
               AFNO3=AGTHNO3/( 1./ACHPA(IL,L,8)+ALWCVMR(IL,L)/AMH )
C 
C              Second iteration H+ molarity calculation.
C
               AFCL =AGTHCL/( 1./ACHPA(IL,L,10)+ALWCVMR(IL,L)/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)+AMSSANI-AMSSCAT)
               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+ACHPA(IL,L,3)
     1              /AMH**2
               AFRAS=1./( 1.+ALWCVMR(IL,L)*ATVAL )
               AMHVAL(IL)=AMH
C
C              Ozone oxidation parameter (in litre-air/mol/sec).
C
               ATVO3=1./( 1.+ALWCVMR(IL,L)*ACHPA(IL,L,7) )
               AF1=( ACHPA(IL,L,11)+ACHPA(IL,L,12)/AMH )*AFRAS
     1               *ATVAL*ATVO3*ACHPA(IL,L,7)*ALWCVMR(IL,L)
C
C              Hydrogen peroxide oxidation parameter 
C              (in litre-air/mol/sec).
C
               AFRAH=1./( 1.+ALWCVMR(IL,L)*ACHPA(IL,L,6) )
               AF2=( ACHPA(IL,L,13)/(0.1+AMH) )*AFRAS
     1               *ACHPA(IL,L,1)*AFRAH*ACHPA(IL,L,6)*ALWCVMR(IL,L)
C
C            *** CONVERSION FROM LITRE/MOL/SEC TO 
C            *** KG-AIR/KG-SULPHUR/SEC.
C
               AF1=AF1*AMR2M
               AF2=AF2*AMR2M
C
C            *** CONVERSION TO HPA/(SEC*M*CLOUD_BASE_FLX)*KG-AIR/KG-S,
C            *** HORIZONTALLY AVERAGED IN GRID CELL.
C
               AF1=AF1*AFAKT
               AF2=AF2*AFAKT
C
C            *** SCAVENGING RATIOS, I.E. FRACTION OF TRACER THAT
C            *** IS DISSOLVED.
C
               ASRSO2=MAX( 1.-AFRAS, ZERO )
               ASRHPO=MAX( 1.-AFRAH, ZERO )
               ASRO3 =MAX( 1.-ATVO3, ZERO )
               ASRNH3=MAX(ALWCVMR(IL,L)*AMH/(1./ACHPA(IL,L,9)
     1               +ALWCVMR(IL,L)*AMH),ZERO)
               ASRHNO=MAX(ALWCVMR(IL,L)/(AMH/ACHPA(IL,L,8)
     1               +ALWCVMR(IL,L)),ZERO)
               ASRCO2=MAX(ALWCVMR(IL,L)*ACHPA(IL,L,5)/AMH,ZERO)
C                                                                       
C            *** CALCULATE TRACER SCAVENGING PROPORTIONAL TO CONVERSION
C            *** OF CLOUD WATER INTO PRECIPITATION.
C
               FAC=MAX(T(IL,JB(IL))-.5*(T(IL,L+1)+T(IL,L)),ZERO)          
     1            /MAX(T(IL,JB(IL))-TFREEZ, ONE)                         
               C0=0.
               IF ( MBC(IL).GT.0. ) 
     1            C0=C0FAC*MIN(FAC**2, ONE)*MU(IL,L)/MBC(IL)
               C0V(IL)=C0
C
C            *** SCAVENGING RATES.
C
               RU (IL,L,ISO4) = ASRSO4 *C0                              
               RU (IL,L,ISO2) = ASRSO2 *C0                              
               RU (IL,L,IHPO) = ASRHPO *C0                              
               RUA(IL,L,INH3) = ASRNH3 *C0                            
               RUA(IL,L,IHNO3)= ASRHNO *C0                            
               RUA(IL,L,ICO2) = ASRCO2 *C0                            
               RUA(IL,L,IO3)  = ASRO3  *C0                        
             ENDIF
C
C          *** UPDRAFT TRACER CONCENTRATIONS FOR ADDITIONAL SPECIES.
C
             XUA(IL,L,INH3) =(XGA(IL,L,INH3)+(XUA(IL,L+1,INH3)-
     1                        XGA(IL,L,INH3))*ATMP)/
     2                       (1.+RUA(IL,L,INH3)*FACV(IL))
             XUA(IL,L,IHNO3)=(XGA(IL,L,IHNO3)+(XUA(IL,L+1,IHNO3)-
     1                        XGA(IL,L,IHNO3))*ATMP)/
     2                       (1.+RUA(IL,L,IHNO3)*FACV(IL))
             XUA(IL,L,ICO2) =(XGA(IL,L,ICO2)+(XUA(IL,L+1,ICO2)-
     1                        XGA(IL,L,ICO2))*ATMP)/
     2                       (1.+RUA(IL,L,ICO2)*FACV(IL))
             XUA(IL,L,IO3)  =(XGA(IL,L,IO3)+(XUA(IL,L+1,IO3)-
     1                        XGA(IL,L,IO3))*ATMP)/
     2                       (1.+RUA(IL,L,IO3)*FACV(IL))
C
C          *** UPDRAFT TRACER CONCENTRATIONS FOR REACTIVE SPECIES.
C          *** NO CHEMICAL SOURCES/SINKS ARE ASSUMED FOR THE 
C          *** PATHOLOGICAL CASE THAT ENTRAINMENT EQUALS DETRAINMENT.
C
             IF ( ((AF1*XUA(IL,L,IO3)+AF2*HPOT)*SO2T).GT.0.
     1            .AND.(.5*(T(IL,L)+T(IL,L+1))).GT.273. 
     2            .AND. FACV(IL).NE.0. ) THEN
                AP1=-(1.+FACV(IL)*(RU(IL,L,ISO2)+AF1*XUA(IL,L,IO3)
     1              +AF2*(HPOT-SO2T)))/(2.*AF2*FACV(IL)*(
     2               1.+FACV(IL)*(RU(IL,L,ISO2)+AF1*XUA(IL,L,IO3))))
                AP2=SO2T/(AF2*FACV(IL)*(
     1               1.+FACV(IL)*(RU(IL,L,ISO2)+AF1*XUA(IL,L,IO3))))
                XU(IL,L,ISO2)=AP1+SQRT(AP1**2+AP2)
                XU(IL,L,IHPO)=HPOT/(1.+(RU(IL,L,IHPO)
     1                       +AF2*XU(IL,L,ISO2))*FACV(IL))
                SO4A=XU(IL,L,ISO2)*(AF1*XUA(IL,L,IO3)
     1              +AF2*XU(IL,L,IHPO))/EUO(IL,L)
                XU(IL,L,ISO4)=(X(IL,L,ISO4)+SO4A+(XU(IL,L+1,ISO4)
     1                        -X(IL,L,ISO4)-SO4A)*ATMP)
     2                        /(1.+RU(IL,L,ISO4)*FACV(IL))
                FACT=XU(IL,L,ISO2)*MBC(IL)*DZ(IL,L)*DSR(IL,L)/DP(IL,L)
                SC(IL,L,IHPO)=-AF2*XU(IL,L,IHPO)*FACT
                SC(IL,L,ISO2)=-(AF1*XUA(IL,L,IO3)+AF2*XU(IL,L,IHPO))
     1                         *FACT
                SC(IL,L,ISO4)=-SC(IL,L,ISO2)
             ELSE
                XU(IL,L,ISO2)=SO2T/(1.+RU(IL,L,ISO2)*FACV(IL))
                XU(IL,L,IHPO)=HPOT/(1.+RU(IL,L,IHPO)*FACV(IL))
                XU(IL,L,ISO4)=SO4T/(1.+RU(IL,L,ISO4)*FACV(IL))
             ENDIF
C
             XU(IL,L,ISO2)=MAX(XU(IL,L,ISO2),ZERO)
             XU(IL,L,ISO4)=MAX(XU(IL,L,ISO4),ZERO)
             XU(IL,L,IHPO)=MAX(XU(IL,L,IHPO),ZERO)
            ELSE                                               !IPAM
             FAC=MAX(T(IL,JB(IL))-.5*(T(IL,L+1)+T(IL,L)),ZERO)          
     1          /MAX(T(IL,JB(IL))-TFREEZ, ONE)                         
             C0=0.
             IF ( MBC(IL).GT.0. )  THEN
               C0=C0FAC*MIN(FAC**2, ONE)*MU(IL,L)/MBC(IL)
             ENDIF
             C0V(IL)=C0
            ENDIF                                              !IPAM
           ENDIF
  240    CONTINUE
C
C        *** UPDRAFT TRACER CONCENTRATIONS FOR TRACERS THAT
C        *** DO NOT UNDERGO CHEMICAL REACTIONS BUT ARE AFFECTED
C         *** BY SCAVENGING.
C
         DO 247 N=1,NTRAC
          IF ( IPAM.EQ.0 ) THEN
           KRIT1=ITRCCH(N).NE.0 .AND. 
     1      (N.NE.ISO2 .AND. N.NE.ISO4 .AND. N.NE.IHPO)
           KRIT2=ITRCCH(N).EQ.-1
          ELSE
           KRIT1=ITRCCH(N).NE.0
           KRIT2=ITRCCH(N).EQ.-1 .OR. N.EQ.IHPO
          ENDIF
          IF(KRIT1) THEN
           IF(KRIT2) THEN
            RUFAC=ASRPHOB
           ELSE
            RUFAC=ASRSO4
           ENDIF
C
           DO 245 IL=IL1G,IL2G
            IF(      L.GT.JT(IL) .AND. L.LT.JB(IL) .AND. EPS0(IL).GT.0.  
     1         .AND. MU(IL,L).GT.0. ) THEN
             RU(IL,L,N)=RUFAC*C0V(IL)
             XU(IL,L,N)=(X(IL,L,N)+(XU(IL,L+1,N)-X(IL,L,N))*ATMPV(IL))/
     1                  (1.+RU(IL,L,N)*FACV(IL))
             XU(IL,L,N)=MAX(XU(IL,L,N),ZERO)
            ENDIF
  245      CONTINUE
          ENDIF
  247    CONTINUE
C
C      *** DIAGNOSTIC CALCULATIONS FOR PH.
C
         IF(ISVCHEM.NE.0.AND.IPAM.EQ.0) THEN
          DO 250 IL=IL1G,IL2G
            IF ( AMHVAL(IL) > 1.E-20 ) THEN
              PHDROW(IDEEP(IL),L)=-LOG10(AMHVAL(IL))
            ELSE
              PHDROW(IDEEP(IL),L)=0.
            ENDIF
  250     CONTINUE
         ENDIF
  260 CONTINUE
C
C     *** CONVERT SCAVENGING RATES TO KG/KG/SEC.
C
      DO 285 N=1,NTRAC
        IF(ITRCCH(N).NE.0)                                    THEN
         DO 280 L=MSG+1,ILEV                                       
         DO 280 IL=IL1G,IL2G                                       
          RU(IL,L,N)=RU(IL,L,N)*MBC(IL)*DZ(IL,L)*XU(IL,L,N)*DSR(IL,L)
     1              /DP(IL,L)
  280    CONTINUE
        ENDIF
  285 CONTINUE
C
C     *** FLUX OF TRACERS INTO THE DOWNDRAFT DUE TO EVAPORATION
C     *** OF RAIN (UNITS KG/M**2/SEC).
C
      DO 320 N=1,NTRAC
        IF(ITRCCH(N).NE.0)                                    THEN
         DO 300 L=ILEV-1,MSG+2,-1                                       
         DO 300 IL=IL1G,IL2G
          IF(      L.GT.JT(IL) .AND. L.LT.JB(IL) .AND. EPS0(IL).GT.0.  
     1      .AND. MU(IL,L).GT.0.)            THEN
            SUMDX(IL,N)=SUMDX(IL,N)+WRK(IL)*RU(IL,L,N)*100.*DP(IL,L)
     1                             /(DSR(IL,L)*GRAV)
          ENDIF
 300     CONTINUE
        ENDIF
 320  CONTINUE
C
C     *** TRACER TENDENCY IN THE SUBCLOUD-LAYER DUE TO PRODUCTION
C     *** IN THE DOWNDRAFT (UNITS KG/KG/SEC). 
C
      DO 360 N=1,NTRAC
        IF(ITRCCH(N).NE.0)                                    THEN
         DO 350 IL=IL1G,IL2G
          L=MX(IL)
          SUMDX(IL,N)=SUMDX(IL,N)*GRAV*DSR(IL,L)/(100.*DP(IL,L))
 350     CONTINUE
        ENDIF
 360  CONTINUE
C
C     * TRACER DOWNDRAFT PROPERTIES.
C
      DO 910 N=1,NTRAC
      DO 910 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.)                                    THEN
           XD(IL,JD(IL),N)=X(IL,JD(IL)-1,N)
        ENDIF
  910 CONTINUE
C
      DO 922 N=1,NTRAC
       IF(ITRPHS(N).GT.0) THEN    
         DO 920 J=MSG+2,ILEV
         DO 920 IL=IL1G,IL2G
           IF((J.GT.JD(IL) .AND. J.LE.JB(IL)) .AND. EPS0(IL).GT.0.      
     1                                     .AND. JD(IL).LT.JB(IL))  THEN
             XD(IL,J,N)=X(IL,J-1,N)+(XD(IL,J-1,N)-X(IL,J-1,N))
     1                             *(MDO(IL,J-1)/MDO(IL,J))
           ENDIF
  920    CONTINUE
       ENDIF
  922 CONTINUE
      DO 930 N=IU,IV
      DO 930 IL=IL1G,IL2G
        IF(EPS0(IL).GT.0.)                                    THEN
           XDA(IL,JD(IL),N)=XGA(IL,JD(IL)-1,N)
        ENDIF
  930 CONTINUE
C
      DO 940 N=IU,IV
      DO 940 J=MSG+2,ILEV
      DO 940 IL=IL1G,IL2G
        IF((J.GT.JD(IL) .AND. J.LE.JB(IL)) .AND. EPS0(IL).GT.0.      
     1                                     .AND. JD(IL).LT.JB(IL))  THEN
          ATMP=MDO(IL,J)/MDO(IL,J-1)
C
C         * UNPERTURBED WIND IN DOWNDRAFT.
C
          XDA(IL,J,N)=XGA(IL,J-1,N)+(XDA(IL,J-1,N)-XGA(IL,J-1,N))/ATMP
        ENDIF
  940 CONTINUE
      DO 941 N=IU,IV
      DO 941 J=MSG+2,ILEV
      DO 941 IL=IL1G,IL2G
        IF((J.GT.JD(IL) .AND. J.LE.JB(IL)) .AND. EPS0(IL).GT.0.      
     1                                     .AND. JD(IL).LT.JB(IL))  THEN
          ATMP=MDO(IL,J)/MDO(IL,J-1)
C
C         * WIND IN DOWNDRAFT FROM PRESSURE GRADIENT FORCE PERTURBATION.
C
          XDA(IL,J,N)=CCD*XGA(IL,J-1,N)+(1.-CCD)*XDA(IL,J,N)
C
C         * DIAGNOSE PRESSURE GRADIENT FORCE TERM UNDER THE ASSUMPTION
C         * OF A CONSTANT SOURCE/SINK IN EACH GRID CELL.
C
          IF ( ABS(1.-ATMP).GT.1.E-02 ) THEN
             SCAT=EDO(IL,J-1)*(XDA(IL,J-1,N)-XGA(IL,J-1,N)
     1           -ATMP*(XDA(IL,J,N)-XGA(IL,J-1,N)))/(1.-ATMP)
          ELSE
             SCAT=MDO(IL,J)*(XDA(IL,J-1,N)-XDA(IL,J,N))/DZ(IL,J-1)
          ENDIF
          SCA(IL,J-1,N)=SCA(IL,J-1,N)
     1                 +SCAT*MBC(IL)*DZ(IL,J-1)*DSR(IL,J-1)/DP(IL,J-1)
        ENDIF
  941 CONTINUE
C
C     * UPDATE TRACER TENDENCIES.
C
      CALL XTEND5(DUDT,DVDT,XGA,XFA,XUA,XDA,SCA,DXDT,X,XF,XU,XD,
     1            SC,RU,MU,MD,DSR,WRKX,WRKXA,SUMDX,WRKS,DP,PFG,P,
     2            PRESSG,MBC,AMBCOR,JB,JT,MX,ITRPHS,DT,GRAV,ROG,
     3            ILG,ILEV,IL1G,IL2G,MSG,IU,IV,NTRAC,NTRACA )
C
C     * ENSURE POSITIVENESS BY SCALING OF PRODUCTION AND LOSS RATES.
C
      AFS=0.999999999
      DO 520 N=1,NTRAC
      DO 520 IL=IL1G,IL2G
        SUMDX(IL,N)=1.
  520 CONTINUE
      DO 527 N=1,NTRAC
       IF (ITRPHS(N).GT.0)                        THEN
         DO 526 L=MSG+1,ILEV
         DO 526 IL=IL1G,IL2G
           IF ( DXDT(IL,L,N).LT.0. ) THEN
              ACORT=-(AFS*(X(IL,L,N)+DXDT(IL,L,N)*DT))
     1            /(DXDT(IL,L,N)*DT)
              SUMDX(IL,N)=MAX(MIN(SUMDX(IL,N),ACORT),ZERO)
           END IF
  526    CONTINUE
       END IF
  527 CONTINUE
      DO 600 N=1,NTRAC
       IF (ITRPHS(N).GT.0)                        THEN
         DO 620 L=MSG+1,ILEV
         DO 620 IL=IL1G,IL2G
            DXDT(IL,L,N)=DXDT(IL,L,N)*SUMDX(IL,N)
            RU  (IL,L,N)=RU  (IL,L,N)*SUMDX(IL,N)
            SC  (IL,L,N)=SC  (IL,L,N)*SUMDX(IL,N)
  620    CONTINUE
       END IF
  600 CONTINUE
C
      IF (KCHEM) THEN
C
C      *** CHEMICAL PRODUCTION RATES (IN KG-S/M**2/SEC).
C
       XSRC(IL1G:IL2G,:)=0.
       DO 370 N=1,NTRAC
       DO 370 L=MSG+1,ILEV                                       
       DO 370 IL=IL1G,IL2G
          FACT=100.*DP(IL,L)/(DSR(IL,L)*GRAV)
          XSRC(IL,N)=XSRC(IL,N)+SC(IL,L,N)*FACT
  370  CONTINUE
C
C      *** WET DEPOSITION (IN KG-S/M**2/SEC).
C
       XSNK(IL1G:IL2G,:)=0.
       DO 310 N=1,NTRAC
       DO 310 L=ILEV-1,MSG+2,-1                                      
       DO 310 IL=IL1G,IL2G
         IF(      L.GT.JT(IL) .AND. L.LT.JB(IL) .AND. EPS0(IL).GT.0.  
     1      .AND. MU(IL,L).GT.0. ) THEN 
            FACT=(1.-WRK(IL))*100.*DP(IL,L)/(DSR(IL,L)*GRAV)
            XSNK(IL,N)=XSNK(IL,N)+RU(IL,L,N)*FACT
         ENDIF
 310   CONTINUE
      END IF
C
      RETURN
      END
