MODULE STRAT_CHEM_MOD
! automatic conversion to free f90 compatible form 
! free.pl strat_chem_model.f  
! linewidth: 72
! file names: strat_chem_model.f 
!
!23456789-123456789-123456789-123456789-123456789-123456789-123456789-23
!/
!/ THIS MODIFICATION SET FOR THE UNIFIED MODEL IS A 
!/ REAL STRATOSPHERIC CHEMICAL MODEL         
!/ WRITTEN BY J AUSTIN (10/02/94); MODIFIED 26/7/95 FOR 49 LEVEL MODEL  
!/ REVISED FOR CLIMATE LENGTH INTEGRATIONS 24/9/96
!/ REVISED FOR PARALLEL CODE 13/11/97
!/ REVISED CHEMICAL MODEL 14/1/98
!/ REVISED CHEMICAL MODEL WITH MAINZ ANALYTIC SCHEME 1/5/99
!/ REVISED INTEGRATION SCHEME AFTER STEILET AL. 1998 14/6/99
!/ UPDATED TO VN4.5 64 LEVELS J. AUSTIN 27/1/00
!/ UPDATED CH4 OXIDATION, TROPOSPHERIC CHEMISTRY 14/2/02
!/ UPDATED FOR USE IN THE FMS, 15/12/03
!/ REVISED, 13/10/04: SIMPLIFIED PSCs, IMPROVED LONG-LIVED TRACERS,
!/      AGE OF AIR PARAMETERISATION FOR HALOGEN RATES OF CHANGE
!/

implicit none
private

public CHEMISTRY, ZEN2, DCLY_DT, SEDIMENT

contains

      SUBROUTINE CHEMISTRY (ALON, ALAT, JLZ, KL, DY, H2O, DAGESQ, OZCOL, &
        PRESS, RHO, TEMP, VTEMP, COZEN, CDX, CHLB, OZB, BNAT, BICE, PHOTO, &
        SOLARDATA, OZON, COSP, COSPHC, ANOY, AEROSOL, DT, MERIDS, CH_TEND, &
        OZONE, O3_PROD, H2O_TEND, MYPE, ITIME)
! """" ----------------------------------------------------------------*
! VERSION NUMBER   : 1   DATE:  6/10/1989.
!   NOTES          : CODE REWRITTEN TO SOLVE EQUATIONS MORE DIRECTLY
!                  : IN FAMILY MODE TO COMPARE MODEL RESULTS
! VERSION NUMBER   : 2   DATE:  1/06/1990.
!   NOTES          : REACTION SET SIMPLIFIED TO REDUCE CPU TIME
! VERSION NUMBER   : 3   DATE: 02/04/1993.
!   NOTES          : REVISED HETEROGENEOUS REACTION SCHEME
! VERSION NUMBER   : 4   DATE: 09/02/1994ET SEQ.
!   NOTES          : ADAPTED FOR USE IN THE UNIFIED MODEL
! VERSION NUMBER   : 5   DATE: 15/12/2003
!   NOTES          : ADAPTED FOR USE IN THE FMS               
! VERSION NUMBER   : 6   DATE: 23/03/2006
!   NOTES          : CONVERTED TO f90 FOR USE IN THE FMS               
! """" ----------------------------------------------------------------*
! PURPOSE          : TO INTEGRATE A ZERO-DIMENSION(POINT) CHEMICAL      
!                    MODEL ALONG A PREDETERMINED TRAJECTORY FOR A       
!                    TIME PERIOD IDT
! METHOD           : USES THE EULER BACKWARD INTEGRATION SCHEME         
! INPUT/OUTPUT LIST:
!     JLX          = LAT.  NUMBER (NEEDED FORERROR DIAGNOSTICS)        
!     KL           = LEVEL NUMBER (NEEDED FOR PHOTOLYSIS RATES)         
!     DT           = INTEGRATION TIMESTEP
!     DY           = ARRAY OF FAMILY CONCENTRATIONS
!     OZCOL        = OZONE COLUMN ABOVE PARTICLE
!     FACT         = RATIO OF OZONE TO OY (FOR PHOTOLYSIS RATES)        
!     RHO          = AIR DENSITY AT THE PARTICLE
!     TEMP         = TEMPERATURE AT THE PARTICLE
!     CHEMINC      = INCREMENTS TO FAMILY CONCENTRATIONS      
!                    AFTER THIS CALL TO OAERO
!
!  SUBPROGRAM CALLS:
!     REACTION_RATES  - CALCULATES REACTION RATES (CALLED ONCE) 
!     PSC/GAMMA_AER/HETRATES - CALCULATE HETEROGENEOUS REACTION RATES
!     PHOTO   - CALCULATES PHOTODISSOCIATION RATES (CALLED ONCE)        
!----------------------------------------------------------------------*
      IMPLICIT NONE
      

      integer ::  nchem, jimpl, ihet, nhet, irx, irx4
      PARAMETER (NCHEM=21,JIMPL=18,IHET=72,NHET=9)
      PARAMETER (IRX=33,IRX4=IRX*4) 

      integer, intent(in) :: MERIDS, JLZ, KL, MYPE
      integer, intent(in) :: ITIME(6)
      real, intent(in)    :: CDX, DT
      real, intent(in), dimension(MERIDS) :: ALON, ALAT,   &
                           DAGESQ, OZCOL, PRESS, RHO, &
                           TEMP, VTEMP, COZEN, &
                           AEROSOL
      real, intent(in), dimension(MERIDS,NCHEM) :: DY
      real, intent(in), dimension(90,15) :: CHLB
      real, intent(in), dimension(144,90) :: OZB
      real, intent(in) :: PHOTO(132,14,11,48)
      real, intent(in) :: SOLARDATA(1801)
      real, intent(in) :: OZON(11,48)
      real, intent(inout) :: COSP(14)
      real, intent(inout) :: COSPHC(48)
      real, intent(in) :: ANOY(90,48)

      real, intent(inout), dimension(MERIDS)       :: H2O, BNAT, BICE, OZONE, O3_PROD, H2O_TEND
      real, intent(inout), dimension(MERIDS,NCHEM) :: CH_TEND
              
!----------------------------------------------------------------------*
! VARIABLES PASSED TO OR FROM HIGHER SUBROUTINES
!----------------------------------------------------------------------*


      INTEGER :: ISTEP_CHEM, IC, IL, ILX, JLX, IZ, NC, KLX, IR, &
                 ITX, ITNO, ITIME_LEFT, IT, IM, IDX, ICYC, IDAYS, LEAPD,&
                 JMOD, ISC
      REAL :: AGESQ
      INTEGER :: IMON(12)
!----------------------------------------------------------------------*
! VARIABLES FOR THIS SUBROUTINE
!----------------------------------------------------------------------*
#define _ALLOC

#ifdef _ALLOC
     REAL, ALLOCATABLE, DIMENSION(:,:) :: RATES, FOTO
     REAL, ALLOCATABLE, DIMENSION(:,:) :: DYY, DYZ, DQT, DPT, RAT, &
              DYZ1, DYY0, DYY1, COND, GAMMA
#else
     REAL :: RATES(MERIDS,IHET+NHET-1)
     REAL :: FOTO(MERIDS,IRX4)
     REAL :: DYY(MERIDS,NCHEM)
     REAL :: DYZ(MERIDS,JIMPL)
     REAL :: DQT(MERIDS,NCHEM)
     REAL :: DPT(MERIDS,NCHEM)
     REAL :: RAT(MERIDS,3)
     REAL :: DYZ1(MERIDS,JIMPL)
     REAL :: DYY0(MERIDS,NCHEM)
     REAL :: DYY1(MERIDS,NCHEM)
     REAL :: COND(MERIDS,3)
     REAL :: GAMMA(MERIDS,NHET)
#endif
     REAL, DIMENSION(MERIDS) ::  CLOX, BROX, ANOZ, H2O1, H2O0, DP1, DQ1, &
            H2SO4, TICE, WH2SO4, AM, AW, ALIQ, RMEAN, ASAT, RNAT, RICE, RTT

     REAL :: DVDT, CH4, AN2O, CH4TR, AN2OTR, ANAT, AICE, VTAU10, VTAU100, &
            VTAUYR, FACT, DELT, SUM, RATIO1, RATIO2, ANUM, ADEN, ANOX, DAX, &
            DA1, DA2, DD1, DB1, DB2, DD2, DDNUM, DA3, DC2, DC3, DD3, R1, R2, &
            DD13, VDDNUM, DAC2, DA12, DD12, CLX, TAU, RAT1, BX1, AX, BX, CX, &
            XYZ, DQT1, DQT2, DQT3, T0, T1, BRX, DIFF, VRHO, BX2, &
            RATIO, X1, X2, SOL, FRAC, SOLARTIME, SOLAR27, DX1, D11, D27, F11, &
            F27, PI, O2, H2, minusT0
      DATA  PI/3.141592653589793/
      DATA  O2/0.2095/,H2/0.5E-6/                                
      DATA  ISTEP_CHEM/900/ ! Chemical timestep (s)
      DATA  ANAT,AICE/10.0, 10.0/
      DATA  IMON/0,31,59,90,120,151,181,212,243,273,304,334/
!
!  ISC IS THE SRES SCENARIO FOR CH4 AND N2O
!  1 = A1B, 2 = A2, 3 = B1, 4 = timeslice
!
      DATA  ISC/1/
! ---------------------------------------------------------------------*
!                                                                       
!   THE ARRAYS DYY AND DYZ CONTAIN THE SPECIES AMOUNTS.  THE SPECIES    
!   DYY ARE CALCULATEDEXPLICITLY, WHILE THE SPECIES DYZ ARE FOUND      
!   IMPLICITLY USING PHOTOCHEMICALEQUILIBRIUM AND RELATED APPROACHES.  
!   THE VALUES IN DYZ ARE RECALCULATED AFRESH FOREVERY CALL TO OAERO.  
!   THE ORDER OF THE SPECIES CONTAINED IN THE ARRAYS DYY AND DYZ IS:    
!      DYY                                                              
!    1. HNO3     2. N2O5   3. H2O2    4. HCL    5. HOCL  6. CLONO2      
!    7. H2CO     8. OY     9. HOBR   10. HNO4  11. HBR  12. BRONO2      
!   13. CH3OOH  14. CO    15. NAT+NY 16. CLY   17. BRY  18. CH4  
!   19. Strat H2O 20. N2O 21. AGE  
!                                                                       
!      DYZ                                                              
!    1.   O3      2.  O3P    3.  O1D    4.  NO      5.  NO2    6. NO3   
!    7.   OH      8.  HO2    9.  H     10.  CL     11.  CLO   12. BR    
!   13.   BRO    14.  CL2O2 15.  BRCL  16.  CH3O2  17.  HONO  18. N
!----------------------------------------------------------------------*

#ifdef _ALLOC
     allocate(RATES(MERIDS,IHET+NHET-1))
     allocate(FOTO(MERIDS,IRX4))
     allocate(DYY(MERIDS,NCHEM))
     allocate(DYZ(MERIDS,JIMPL))
     allocate(DQT(MERIDS,NCHEM))
     allocate(DPT(MERIDS,NCHEM))
     allocate(RAT(MERIDS,3))
     allocate(DYZ1(MERIDS,JIMPL))
     allocate(DYY0(MERIDS,NCHEM))
     allocate(DYY1(MERIDS,NCHEM))
     allocate(COND(MERIDS,3))
     allocate(GAMMA(MERIDS,NHET))
#endif
!
! Set PSCs initially to zero
!
     DO IL = 1,MERIDS
       BNAT(IL) = 0.0
       BICE(IL) = 0.0
     ENDDO
!
! FIND CONCENTRATIONS AND TRENDS OF CH4 AND N2O, ACCORDING TO THE SRES SCENARIO
!
      CALL GHGS(ITIME,ISC,CH4,AN2O,CH4TR,AN2OTR)
!   
!  RELAX TOWARDS SPECIFIC LOWER BOUNDARY VALUES WITH A TIMESCALE OF 10 DAYS
!  (SPECIES 1-15) SPECIES 16-21 ARE TREATED DIFFERENTLY
!                   
      VTAU10 = 1.0/(10.0*86400.0)
      VTAU100 = 0.1*VTAU10
      VTAUYR = 1.0/(365.25*86400.0)
      IF(KL.EQ.48) THEN                                                  
        DO IC = 1,15                                            
          IF(IC.EQ.8) THEN
            DO IL = 1,MERIDS                                            
              ILX = 1 + ALON(IL)*72/PI
              JLX = 1 + (PI*0.5 + ALAT(IL))*89.0/PI
              IF(ILX.GT.144) ILX = ILX - 144 
              CH_TEND(IL,IC) = (OZB(ILX,JLX) - DY(IL,IC))*VTAU10      
            ENDDO
          ELSE             
            DO IL = 1,MERIDS                                            
              JLX = 1 + (PI*0.5 + ALAT(IL))*89.0/PI
              CH_TEND(IL,IC) = (CHLB(JLX,IC) - DY(IL,IC))*VTAU10     
            ENDDO
          ENDIF
        ENDDO                                             
        DO IL = 1,MERIDS
          CH_TEND(IL,16) = - DY(IL,16)*VTAU10     
          CH_TEND(IL,17) = - DY(IL,17)*VTAU10     
          CH_TEND(IL,18) = (CH4 - DY(IL,18))*VTAU10     
          CH_TEND(IL,19) = (3.0E-6 - DY(IL,19))*VTAU10     
          CH_TEND(IL,20) = (AN2O - DY(IL,20))*VTAU10     
        ENDDO
        RETURN
      ENDIF             
!----------------------------------------------------------------------*
!   NO CHANGE TO CHEMICAL CONCENTRATIONS IN THE TOP LEVELS,EXCEPT 
!   FOR AGE
!
!----------------------------------------------------------------------*
!      IF(PRESS(1).LT.8.0) THEN
!        DO IL = 1,MERIDS                                            
!        CH_TEND(IL,21) = VTAUYR     
!       ENDDO
!        RETURN
!     ENDIF
!----------------------------------------------------------------------*
!   RESTRICT THE ACTUAL CHEMICAL MODEL TO STRATOSPHERIC AND 
!   MESOSPHERIC AIR
!----------------------------------------------------------------------*
      DO IL = 1,MERIDS                                               
        DO IZ = 1,NCHEM                                                
          DYY(IL,IZ) = DY(IL,IZ)*RHO(IL)                                    
        ENDDO                  
        H2O(IL) = H2O(IL)*RHO(IL)
        DYZ(IL,1) = OZONE(IL)*RHO(IL)
        DYZ1(IL,1) = OZONE(IL)*RHO(IL)
!        IF(PRESS(IL).LT.10.0.AND.OZONE(IL).GT.1.0E-5) THEN
!          DYZ(IL,1) = 1.0E-5*RHO(IL)
!          DYZ1(IL,1) = DYZ(IL,1)
!        ENDIF                                            
        DO IZ = 2,JIMPL                                                
          DYZ(IL,IZ) = 0.                                                   
          DYZ1(IL,IZ) = 0.                                                  
        ENDDO
      ENDDO                                                   
!----------------------------------------------------------------------*
!   CALL SUBROUTINES FOR TEMPERATURE DEPENDENT RATES, HETEROGENEOUS RATES 
!   AND  FOR PHOTOLYSIS RATES (ONCE PER CALL TO CHEMISTRY)
!   PASS THE PRESSURE IN PA INTO THE OTHER CHEMISTRY S/RS
!----------------------------------------------------------------------*
      CALL REACTION_RATES (RHO,TEMP,VTEMP,PRESS,RATES,MERIDS,IHET,NHET)
      IF(KL.GT.10.AND.KL.LT.28) THEN
        DO IL = 1,MERIDS
!
!  Set H2SO4 as an analytical function of age, peaking at 0.5 ppbv
!
          X1 = 3.0*(ALOG10(PRESS(IL)) - 3.5)**2
          X2 = 3.0/(10.0*DYY(IL,21) + 1.0)*(DYY(IL,21) - 1.0)**2
          H2SO4(IL) = (0.01 + 0.49*EXP(-X1-X2))*1.0E-9*RHO(IL)
        ENDDO
        CALL PSC(TEMP,PRESS,RHO,DYY,H2O,H2SO4,ANAT,AICE,COND,TICE,      &
           WH2SO4,AM,AW,ALIQ,RMEAN,ASAT,RNAT,RICE,MERIDS,NCHEM,MYPE)
        DO IL = 1,MERIDS
          BNAT(IL) = COND(IL,2)/RHO(IL)
          BICE(IL) = COND(IL,3)/RHO(IL)
        ENDDO
        CALL GAMMA_AER(TEMP,PRESS,RHO,DYY,H2O,WH2SO4,AM,AW,RMEAN,GAMMA, &
           RTT,MERIDS,NCHEM,NHET,MYPE)
        CALL HETRATES(TEMP,PRESS,RHO,DYY,H2O,TICE,ANAT,AICE,AEROSOL,    &
           RMEAN,RNAT,RICE,GAMMA,RTT,RATES,IHET,NHET,MERIDS,NCHEM,MYPE)       
      ENDIF
!
      DO NC=1,NCHEM                                                  
        DO IL = 1,MERIDS                                               
          DYY0(IL,NC) = DYY(IL,NC) ! save amounts after natice subtract 
          H2O0(IL) = H2O(IL) 
        ENDDO
      ENDDO                                                
      KLX = KL
      CALL PHOTO_RATES(VTEMP,OZCOL,COZEN,KLX,PHOTO,OZON,COSP,COSPHC,         &
                        MERIDS,IRX,IRX4,FOTO,MYPE)
! ---------------------------------------------------------------------*
!  INCLUDE VARIATIONS IN 11-YR CYCLE AND 27-DAY SOLAR CYCLE AND MULTIPLY 
!  FOTO BY SUN-EARTH DISTANCE TERM (CDX). 
!     IM....NO. OF MONTHS SINCE DEC 1949: USED TO DETERMINE SOLAR F10.7 FLUX.
!     SOLARTIME....TIME AFTER 1 JAN 1960 FOR (ARBITRARY) MIN OF 27 DAY CYCLE
!     ASSUME 12EQUAL MONTHS PER YEAR OF 365.25 DAYS
!  PHOTOLYSIS RATES ARE TAKEN TO BE LINEARLY RELATED TO THE F10.7 FLUX 
!  SOLAR MAX AND SOLAR MIN PHOTO RATES REFER TO APPROX F10.7 FLUX OF 
!  208.1 AND 72.0, CORRESPONDING TO THE YEARS 1991 AND 1996.
! ------------------------------------------------------------------
!
      IF(ISC.EQ.4) THEN
        FACT = 0.0
        SOL = 140.05
      ELSE
        IDAYS = (ITIME(1) - 1960)*365 + IMON(ITIME(2)) + ITIME(3)
        LEAPD = (ITIME(1) - 1960)/4
        JMOD = 4*LEAPD + 1960
        IF(ITIME(2).LT.3.AND.JMOD.EQ.ITIME(1)) LEAPD = LEAPD - 1
        FRAC = FLOAT(ITIME(4)*60 + ITIME(5))/1440.0
        SOLARTIME = FLOAT(IDAYS + LEAPD) + FRAC 
        ICYC = SOLARTIME/27
        SOLAR27 = (SOLARTIME - 27.0*ICYC)/27.0
        FACT = SIN(SOLAR27*2*PI)
        DX1 =  SOLARTIME*12.0/365.25 + 0.5       
        IDX = DX1
        DX1 = DX1 - IDX
        IM = 120 + IDX
        IF(IM.LT.1) IM = 1
        IF(IM.GT.1800) IM = 1800
        SOL = 0.1*(SOLARDATA(IM+1)*DX1 + (1.0 - DX1)*SOLARDATA(IM))
      ENDIF
      DO IR = 1,IRX
        DO IL = 1,MERIDS
          D11 = FOTO(IL,IR) - FOTO(IL,IR+IRX*2)
          D27 = FOTO(IL,IR+IRX) - FOTO(IL,IR+IRX*3)
          F11 = FOTO(IL,IR+IRX*2) + (SOL - 72.0)*D11/136.1
          F27 = FOTO(IL,IR+IRX*3) + (SOL - 72.0)*D27/136.1
          FOTO(IL,IR) = CDX*(F11 + F27*FACT) 
          IF(FOTO(IL,IR).LT.0.0) FOTO(IL,IR) = 0.0
        ENDDO
      ENDDO
! ---------------------------------------------------------------------*
! --- ITX    : TIMESTEP NUMBER                                          
! ---   IT     : ITERATION LOOP COUNTER                                 
! ---                                                                   
! ---   FIRST COMPUTE N + NO + NO2 + NO3 (ANOZ), CL + CLO (CLOX) AND        
! ---   BR + BRO + BRCL (BROX)                                          
! ---                                                                   
! ---------------------------------------------------------------------*
!
      ITX = 0
      ITNO = 6
      ITIME_LEFT = NINT(DT)
      DO 1100
!
      IF(ITIME_LEFT .LE. ISTEP_CHEM) THEN
        DELT = REAL(ITIME_LEFT)
        ITIME_LEFT = 0
      ELSE
        DELT = REAL(ISTEP_CHEM)
        ITIME_LEFT = ITIME_LEFT - ISTEP_CHEM
      ENDIF
      DVDT = 1.0/DELT
      ITX = ITX + 1
      IF (ITX .GE. 2) ITNO = 2
      DO 1000 IT = 1,ITNO
        DO 200 IL = 1,MERIDS
          SUM = DYY(IL,1) + 2.*DYY(IL,2) + DYY(IL,6) + DYY(IL,12) +          &
                DYY(IL,10)
          ANOZ(IL) = DYY(IL,15) - SUM
          IF(ANOZ(IL).LT.0.)  ANOZ(IL) = 0.
          DYZ(IL,5) = ANOZ(IL)
          IF(DYY(IL,1).GT.DYY(IL,15)) DYY(IL,1) = DYY(IL,15)
          SUM = DYY(IL,4) + DYY(IL,5) + DYY(IL,6)
          CLOX(IL) = DYY(IL,16) - SUM
          IF(CLOX(IL).LT.0.)  CLOX(IL) = 0.
          SUM = DYY(IL,12) + DYY(IL,11) + DYY(IL,9)
          BROX(IL) = DYY(IL,17) - SUM
          IF(BROX(IL).LT.0.)  BROX(IL) = 0.
  200   ENDDO
!----------------------------------------------------------------------*
!---
!---         COMPUTE PHOTOCHEMICAL EQUILIBRIUM FOR THE O1D/O3 AND       
!---         O3P/O3 RATIOS
!----------------------------------------------------------------------*
        DO 900 IL=1,MERIDS
          ANUM = (FOTO(IL,2) + FOTO(IL,28)*DYY(IL,20)/DYY(IL,8))
          ADEN =                                                            &
               (RATES(IL,3) + RATES(IL,9)*H2O(IL) +                         &
               RATES(IL,30)*DYY(IL,18) + RATES(IL,66)*DYY(IL,20) +          &
               RATES(IL,70)*DYY(IL,20))
          RATIO2 = ANUM/ADEN
          ANUM = RATES(IL,3)*RATIO2 + FOTO(IL,3) +                          &
                 (2.*FOTO(IL,1)*O2*RHO(IL) + DYZ(IL,5)*FOTO(IL,4) +         &
                  DYZ(IL,4)*FOTO(IL,31) +                                   &
                  DYZ(IL,6)*FOTO(IL,6) + RATES(IL,51)*DYZ(IL,7)**2)/        &
                  DYZ(IL,1)
          ADEN = RATES(IL, 1)           + RATES(IL, 2)*DYZ(IL, 1) +          &
                 RATES(IL, 5)*DYZ(IL,5) + RATES(IL,11)*DYZ(IL, 8) +          &
                 RATES(IL,12)*DYZ(IL,7) + RATES(IL,21)*DYZ(IL,11) +          &
                 RATES(IL,24)*DYY(IL,6) + RATES(IL,44)*DYY(IL, 9)           
          RATIO1 = ANUM/ADEN
          SUM = 1. + RATIO1 + RATIO2
          RAT(IL,1) = 1./SUM
          DYZ(IL,1) = DYY(IL,8)*RAT(IL,1)
          DYZ(IL,2) = RATIO1*DYZ(IL,1)
          DYZ(IL,3) = RATIO2*DYZ(IL,1)
          RAT(IL,2) = RATIO1*RAT(IL,1)
          RAT(IL,3) = RATIO2*RAT(IL,1)
          IF(DYZ(IL,1).LT.1.0E-15*RHO(IL)) THEN
            XYZ = 1.0E-15*RHO(IL) - DYZ(IL,1)
            DYZ(IL,1) = 1.0E-15*RHO(IL)
            DYZ(IL,2) = DYZ(IL,2) + XYZ
            RAT(IL,1) = DYZ(IL,1)/DYY(IL,8)
            RAT(IL,2) = DYZ(IL,2)/DYY(IL,8)
            RAT(IL,3) = 1.0 - RAT(IL,1) - RAT(IL,2)
            IF(RAT(IL,3).LT.0.0) RAT(IL,3) = 0.0
          ENDIF
  900   ENDDO
!----------------------------------------------------------------------*
!--
!---         APPLY PHOTOCHEMICAL EQUILIBRIUM FOR NO3 SUBJECT TO
!---        ENOUGH NOZ BEING PRESENT
!----------------------------------------------------------------------*
        DO 910 IL = 1,MERIDS
          ADEN = RATES(IL,7)*DYZ(IL,5) + FOTO(IL,6) + FOTO(IL,7) +          &
                 RATES(IL,53)*DYZ(IL,4) + RATES(IL,58)*DYZ(IL,5) +          &
                 2.*RATES(IL,59)*DYZ(IL,6) + RATES(IL,62)*DYY(IL,7)            
          IF(ADEN.LT.DVDT.AND.FOTO(IL,6).NE.0.) ADEN = DVDT                 
          IF(ADEN.NE.0.) THEN
            ANUM = (FOTO(IL,8) + RATES(IL,8))*DYY(IL,2) +                   &
              RATES(IL,19)*DYZ(IL, 7)*DYY(IL,1) +                           &
              RATES(IL, 6)*DYZ(IL, 5)*DYZ(IL,1) +                           &
             (RATES(IL,24)*DYZ(IL, 2) + FOTO(IL,12))*DYY(IL,6) +            &
               FOTO(IL,17)*DYY(IL,12)
            DYZ(IL,6) = ANUM/ADEN
            IF(DYZ(IL,6).GT.(0.5*ANOZ(IL))) DYZ(IL,6) = 0.5*ANOZ(IL)
          ENDIF           
!----------------------------------------------------------------------*
!---        APPLY PHOTOCHEMICAL EQUILIBRIUM FOR N SUBJECT TO         
!---       ENOUGH NOZ BEING PRESENT
!---
!----------------------------------------------------------------------*
          ADEN = RATES(IL,67)*DYZ(IL,4) + RATES(IL,68)*O2*RHO(IL) +         &
                 RATES(IL,71)*DYZ(IL,5)
          ANUM = FOTO(IL,31)*DYZ(IL,4)
          DYZ(IL,18) = ANUM/ADEN
          IF(DYZ(IL,18).GT.(0.5*ANOZ(IL))) DYZ(IL,18) = 0.5*ANOZ(IL)
!----------------------------------------------------------------------*
!---
!---    COMPUTE NO/NO2 RATIO ASSUMING PHOTOCHEMICAL EQUILIBRIUM         
!---
! ---   SOLVE SIMULTANEOUSEQUATIONS
! ---      DA1*NO - DA2*NO2 = DD1
! ---    - DB1*NO + DB2*NO2 = DD2
! ---
! ---       DA1*NO......THE RATE OF LOSS TERM OF NO
! ---       DA2*NO2.....THE RATE OF PRODUCTION OF NO FROM NO2
! ---       DB1*NO......THE RATE OF PRODUCTION OF NO2 FROM NO
! ---       DB2*NO2.....THE RATE OF LOSS TERM OF NO2
! ---       DD1.........THE TOTAL PRODUCTION TERM (P) FOR NO LESS
! ---                   CONTRIBUTIONS FROM NO2
! ---       DD2.........THE TOTAL PRODUCTION TERM (P) FOR NO2 LESS
! ---                   CONTRIBUTIONS FROM NO
! ---------------------------------------------------------------------*
          ANOX = ANOZ(IL) - DYZ(IL,6) - DYZ(IL,18)
          DAX = RATES(IL, 4)*DYZ(IL, 1) + RATES(IL,18)*DYZ(IL, 8) +          &
                RATES(IL,22)*DYZ(IL,11) + RATES(IL,43)*DYZ(IL,13) +         &
                RATES(IL,53)*DYZ(IL,6) + RATES(IL,55)*DYZ(IL,16) +          &
                RATES(IL,64)*DYZ(IL,2)          
          DA1 = DAX + RATES(IL,65)*DYZ(IL,7) + RATES(IL,67)*DYZ(IL,18) +    &
                FOTO(IL,31)
          DA2 = FOTO(IL,4) + RATES(IL,5)*DYZ(IL,2) +                        &
                RATES(IL,58)*DYZ(IL,6)
          DD1 = FOTO(IL,7)*DYZ(IL,6) + FOTO(IL,26)*DYZ(IL,17) +             &
                RATES(IL,68)*DYZ(IL,18)*O2*RHO(IL) +                        &
                2.*RATES(IL,66)*DYY(IL,20)*DYZ(IL,3)
          DB1 = DAX + RATES(IL,53)*DYZ(IL,6)
          DB2 = DA2 + RATES(IL, 6)*DYZ(IL, 1) + RATES(IL, 7)*DYZ(IL, 6) +    &
                      RATES(IL,17)*DYZ(IL, 7) + RATES(IL,23)*DYZ(IL,11) +    &
                      RATES(IL,40)*DYZ(IL,13) + RATES(IL,49)*DYZ(IL, 8) +   &
                      RATES(IL,71)*DYZ(IL,18)
          DD2 = RATES(IL, 8)*DYY(IL,2) + FOTO(IL, 5)*DYY(IL,1) +             &
                 FOTO(IL, 6)*DYZ(IL,6) + FOTO(IL, 8)*DYY(IL,2) +             &
               (RATES(IL,50)*DYZ(IL,7) + FOTO(IL,25) + RATES(IL,63))*       &
              DYY(IL,10) + 2.*RATES(IL,59)*DYZ(IL,6)**2 +                   &
              RATES(IL,60)*DYZ(IL,7)*DYZ(IL,17)
          DDNUM = DD2*DA1 + DB1*DD1
          IF(DDNUM.NE.0.) THEN
            RATIO = (DB2*DD1 + DA2*DD2)/DDNUM
          ELSE
            RATIO = 0.
          ENDIF
          DYZ(IL,5) = ANOX/(1. + RATIO)
          DYZ(IL,4) = RATIO*DYZ(IL,5)
  910   ENDDO
!----------------------------------------------------------------------*
!---
!---    COMPUTE THE BR/BRO AND BR/BRCL RATIO FROM THE EQUATIONS         
!---
! ---   SOLVE THREE WAY SIMULTANEOUSEQUATIONS                          
! ---      DA1*BR - DA2*BRO - DA3*BRCL = DD1                            
! ---    - DB1*BR + DB2*BRO            = DD2                            
! ---             - DC2*BRO + DC3*BRCL = DD3                            
! ---
! ---------------------------------------------------------------------*
        DO 930 IL = 1,MERIDS
          DA1 = RATES(IL,39)*DYZ(IL, 1) + RATES(IL,45)*DYZ(IL,8) +           &
                RATES(IL,46)*DYY(IL, 7)
          DA2 = RATES(IL,41)*DYZ(IL,11) + RATES(IL,43)*DYZ(IL,4) +           &
                 FOTO(IL,19)
          DA3 =  FOTO(IL,24)
          DD1 = RATES(IL,    47)*DYZ(IL, 7)*DYY(IL,11) +                     &
             2.*RATES(IL,IHET+7)*DYY(IL,11)*DYY(IL, 9) +                     &
                 FOTO(IL,    17)*DYY(IL,12)            +                     &
                 FOTO(IL,    18)*DYY(IL, 9)
          DB1 = RATES(IL,    39)*DYZ(IL, 1)
          DB2 = DA2 + RATES(IL,40)*DYZ(IL, 5) + RATES(IL,42)*DYZ(IL,8) +     &
                      RATES(IL,48)*DYZ(IL,11)
          DD2 = RATES(IL,44)*DYZ(IL, 2)*DYY(IL,9)
          DC2 = RATES(IL,48)*DYZ(IL,11)
          DC3 = DA3
          DD3 = RATES(IL,IHET+5)*DYY(IL, 9)*DYY(IL,4) +                      &
                RATES(IL,IHET+6)*DYY(IL,11)*DYY(IL,5)
!
          R1 = 0.0
          R2 = 0.0
          DD13 = DD1*DC3 + DD3*DA3
          DDNUM = DD13*DB1 + DD2*DA1*DC3
          IF(DDNUM.GT.0.0) THEN
            VDDNUM = 1.0/DDNUM
            DAC2 = DA2*DC3 + DC2*DA3
            R1 = (DD13*DB2 + DD2*DAC2)*VDDNUM
            DA12 = DA1*DB2 - DA2*DB1
            IF(DA12.LT.0.0) DA12 = 0.0
            DD12 = DA1*DD2 + DB1*DD1
            R2 = (DD12*DC2 + DD3*DA12)*VDDNUM
          ENDIF
          IF (FOTO(IL,24).GT.0.0) THEN    
            DYZ(IL,13) = BROX(IL)/(1. + R1 + R2)
            DYZ(IL,12) = R1*DYZ(IL,13)
            DYZ(IL,15) = R2*DYZ(IL,13)
          ELSE
            DYZ(IL,13) = 0.0
            DYZ(IL,12) = 0.0
            DYZ(IL,15) = BROX(IL)
          ENDIF
  930   ENDDO
!----------------------------------------------------------------------*
!---
!---    COMPUTE CL/CLO RATIO ASSUMING PHOTOCHEMICAL EQUILIBRIUM         
!---
! ---   SOLVE SIMULTANEOUS EQUATIONS
! ---      DA1*CL - DA2*CLO = DD1
! ---    - DB1*CL + DB2*CLO = DD2
! ---
! ---------------------------------------------------------------------*
        DO 940 IL = 1,MERIDS
          DA1 = RATES(IL,20)*DYZ(IL, 1) + RATES(IL,25)*DYY(IL,18) +         &
                RATES(IL,27)*DYZ(IL, 8) + RATES(IL,33)*DYY(IL, 7)      
          DA2 = RATES(IL,21)*DYZ(IL, 2) + RATES(IL,22)*DYZ(IL, 4) +         &
                RATES(IL,29)*DYZ(IL, 7) + RATES(IL,41)*DYZ(IL,13)      
          DD1 = RATES(IL,26)*DYZ(IL, 7)*DYY(IL,4) + FOTO(IL,15)*DYY(IL,5) +  &
              2.*FOTO(IL,16)*DYZ(IL,14) +  FOTO(IL,12)*DYY(IL,6) +           &
                 FOTO(IL,24)*DYZ(IL,15) +  FOTO(IL,33)*DYY(IL,4)
          DB1 = RATES(IL,20)*DYZ(IL, 1)
          DB2 = RATES(IL,21)*DYZ(IL, 2) + RATES(IL,22)*DYZ(IL, 4)    +       &
                RATES(IL,23)*DYZ(IL, 5) + RATES(IL,28)*DYZ(IL, 8)    +       &
               (RATES(IL,29) + RATES(IL,69))*DYZ(IL, 7) +                   &
                RATES(IL,37)*DYZ(IL,11)*2. +                                &
               (RATES(IL,41) + RATES(IL,48))*DYZ(IL,13)
          DD2 = RATES(IL,24)*DYZ(IL, 2)*DYY(IL,6) +                          &
             2.*RATES(IL,38)*DYZ(IL,14)
          DDNUM = DD2*DA1 + DB1*DD1
          IF(DDNUM.NE.0.) THEN
            RATIO = (DB2*DD1 + DA2*DD2)/DDNUM
          ELSE
            RATIO = 0.
          ENDIF
          CLX = CLOX(IL) - DYZ(IL,15)
          IF(CLX.LT.0.0) CLX = 0.0
          IF(FOTO(IL,15).GT.0.0)  THEN
            TAU = RATES(IL,37)/(FOTO(IL,16) + RATES(IL,38))                
            RAT1 = TAU*DYZ(IL,11)
            DYZ(IL,11) = CLX/(1. + RATIO + 2.*RAT1)
            RAT1 = TAU*DYZ(IL,11)
            DYZ(IL,11) = CLX/(1. + RATIO + 2.*RAT1)
            DYZ(IL,10) = DYZ(IL,11)*RATIO
            DYZ(IL,14) = DYZ(IL,11)*RAT1
          ELSE
            DYZ(IL,10) = 0.
            DYZ(IL,11) = 0.
            DYZ(IL,14) = CLX*0.5
          ENDIF
  940   ENDDO
!----------------------------------------------------------------------*
!---
!---    COMPUTE HO2 AND OH ASSUMING PHOTOCHEMICAL EQUILIBRIUM           
!---
! ---   SOLVE SIMULTANEOUSEQUATIONS
! ---      DA1*OH - DA2*HO2 = DD1
! ---    - DB1*OH + DB2*HO2 = DD2
! ---
! ---------------------------------------------------------------------*
        DO 950 IL = 1,MERIDS
          BX1 = RATES(IL,17)*DYZ(IL, 5) +                                   &
                RATES(IL,19)*DYY(IL, 1) + RATES(IL,26)*DYY(IL, 4) +         &
                RATES(IL,31)*DYY(IL,18) + RATES(IL,50)*DYY(IL,10) +         &
                RATES(IL,57)*DYY(IL,13) + RATES(IL,60)*DYZ(IL,17) +         &
                RATES(IL,65)*DYZ(IL,4) + RATES(IL,69)*DYZ(IL,11)
          BX2 = RATES(IL,27)*DYZ(IL,10) + RATES(IL,28)*DYZ(IL,11) +         &
                RATES(IL,42)*DYZ(IL,13) + RATES(IL,45)*DYZ(IL,12) +         &
                RATES(IL,49)*DYZ(IL, 5) + RATES(IL,56)*DYZ(IL,16)
          DB1 = RATES(IL,10)*DYZ(IL, 1) + RATES(IL,12)*DYZ(IL, 2) +         &
                RATES(IL,29)*DYZ(IL,11) + RATES(IL,32)*DYY(IL, 7) +         &
                RATES(IL,36)*DYY(IL,3) +  RATES(IL,47)*DYY(IL,11) +         &
                RATES(IL,52)*DYY(IL,14) + RATES(IL,54)*H2*RHO(IL)
          DA1 = DB1 + BX1 + RATES(IL,16)*DYZ(IL, 8) +                       &
                 (RATES(IL,35) + RATES(IL,51))*DYZ(IL,7)*2.
          DA2 = RATES(IL,11)*DYZ(IL, 2) + RATES(IL,15)*DYZ(IL,1) +          &
                RATES(IL,18)*DYZ(IL, 4)
          DB2 = RATES(IL,16)*DYZ(IL, 7) + RATES(IL,34)*DYZ(IL,8)*2. +       &
                DA2 + BX2
          DD1 = DYZ(IL,3)*(2.*RATES(IL, 9)*H2O(IL) +                        &
                              RATES(IL,30)*DYY(IL,18)) +                    &
                RATES(IL,44)*DYZ(IL,2)*DYY(IL,9) +                          &
                2.*RATES(IL,61)*DYZ(IL,16)**2 +                             &
                 FOTO(IL, 5)*DYY(IL, 1) +  FOTO(IL,15)*DYY(IL, 5) +         &
              2.*FOTO(IL,11)*DYY(IL, 3) +  FOTO(IL,18)*DYY(IL, 9) +         &
                FOTO(IL,26)*DYZ(IL,17) +                                    &
! FOTO(IL,27) IS TAKEN AS A PROXY FOR CH3O2 + HV --> H2CO + OH      
                FOTO(IL,27)*(DYZ(IL,16)+DYY(IL,13)) +                       &
                FOTO(IL,32)*H2O(IL)
          DD2 = (RATES(IL,33)*DYZ(IL,10) +                                  &
                RATES(IL,46)*DYZ(IL,12) + FOTO(IL,13)*2.0)*DYY(IL,7) +      &
                RATES(IL,55)*DYZ(IL,4)*DYZ(IL,16) +                         &
                RATES(IL,62)*DYZ(IL,6)*DYY(IL,7) +                          &
                RATES(IL,63)*DYY(IL,10) + FOTO(IL,25)*DYY(IL,10) +          &
                FOTO(IL,27)*DYY(IL,13) +                                    &
                FOTO(IL,32)*H2O(IL)  +  FOTO(IL,33)*DYY(IL,4) +             &
                FOTO(IL,32)*DYY(IL,18)*1.53/1.83
          DDNUM = DB2*DD1 + DA2*DD2
          IF(DD1.NE.0.0.AND.DD2.NE.0.0) THEN
            RATIO2 = (DD2*DA1 + DB1*DD1)/DDNUM
            AX = 4.*(RATES(IL,16)*RATIO2 + RATES(IL,34)*RATIO2**2 +        &
                  RATES(IL,35) + RATES(IL,51))
            BX = BX1 + RATIO2*BX2
            CX = DD1 + DD2
            IF(CX.LT.0.) CX = 0.0
            IF(BX.LT.0.) BX = 0.0
            IF(BX.EQ.0.) THEN
              DYZ(IL,7) = SQRT(2.*CX/AX)
            ELSE
              DYZ(IL,7) = (-BX + SQRT(BX**2 + 2.*AX*CX))/AX
            ENDIF
            DYZ(IL,8) = RATIO2*DYZ(IL,7)
          ELSE
            DYZ(IL,7) = 0.0
            DYZ(IL,8) = 0.0
          ENDIF
          RATIO1 = RATES(IL,12)*DYZ(IL,2)/(RATES(IL,14) +                   &
                   RATES(IL,13)*DYZ(IL,1))
          DYZ(IL,9) = RATIO1*DYZ(IL,7)
  950   ENDDO
!----------------------------------------------------------------------*
!---
!---    COMPUTE CH3O2 AND HONO ASSUMING PHOTOCHEMICAL EQUILIBRIUM           
!---
!----------------------------------------------------------------------*
        DO 960 IL = 1,MERIDS
          AX = 4.*RATES(IL,61)
!
! FOTO(IL,27) IS TAKEN AS A PROXY FOR CH3O2 + HV --> H2CO + OH
!
          BX = RATES(IL,55)*DYZ(IL,4) + RATES(IL,56)*DYZ(IL,8) +            &
               FOTO(IL,27)
          CX = (RATES(IL,25)*DYZ(IL,10) + RATES(IL,30)*DYZ(IL,3) +          &
                RATES(IL,31)*DYZ(IL,7) + FOTO(IL,32)*1.53/1.83)*DYY(IL,18) + &
                RATES(IL,57)*DYZ(IL,7)*DYY(IL,13)
          IF(CX.LT.0.) CX = 0.0
          IF(BX.LT.0.) BX = 0.0
          IF(BX.EQ.0.) THEN
            DYZ(IL,16) = SQRT(2.*CX/AX)
          ELSE
            DYZ(IL,16) = (-BX + SQRT(BX**2 + 2.*AX*CX))/AX
          ENDIF
  960   ENDDO!CONTINUE
        DO 970 IL = 1,MERIDS
          DDNUM = RATES(IL,60)*DYZ(IL,7) + FOTO(IL,26)
          IF(DDNUM.NE.0.0)  THEN
            DYZ(IL,17) = RATES(IL,65)*DYZ(IL,7)*DYZ(IL,4)/DDNUM
          ELSE
            DYZ(IL,17) = RATES(IL,65)*DYZ(IL,4)/RATES(IL,60)
          ENDIF
  970   ENDDO!CONTINUE
        IF(IT.EQ.1.AND.ITX.EQ.1) THEN
          DO IC = 1,JIMPL
            DO IL = 1,MERIDS
              DYZ1(IL,IC) = DYZ(IL,IC)
            ENDDO
          ENDDO
        ELSE
          DO IC = 1,JIMPL
            DO IL = 1,MERIDS
              XYZ = (DYZ(IL,IC) + DYZ1(IL,IC))*0.5
              DYZ(IL,IC) = XYZ
              DYZ1(IL,IC) = XYZ
            ENDDO
          ENDDO
        ENDIF
!----------------------------------------------------------------------*
!   INTEGRATE MODEL FOR SINGLE TIMESTEP OF LENGTH DT
! --- THE CHEMISTRY SCHEME CALCULATES DESTRUCTION AND
! --- PRODUCTION TERMS READY FOR AN INTEGRATION OF THE FORM
! ---
! ---                 Y = P/Q + (Y - P/Q)*EXP(-Q.DT)
! ---
! ---EXCEPT WHERE Q.DT IS SMALL, WHERE 
! ---   
! ---                 Y = (Y/DT + P)/(Q + 1/DT)
! ---
! --- (FIRST ORDER EULER BACKWARD SCHEME) IS USED.
! ---
! ---   DQT : TOTAL DESTRUCTION RATE (FOR EACH EXPLICIT SPECIES)        
! ---   DPT : TOTAL PRODUCTION RATE (FOR EACH EXPLICIT SPECIES)         
! ---
! ---------------------------------------------------------------------*
        DO 990 IL = 1,MERIDS
          JLX = 1 + (PI*0.5 + ALAT(IL))*89.0/PI
!
!   HNO3
!
          DPT(IL,1) = RATES(IL,17)*DYZ(IL,5)*DYZ(IL,7) +                    &
                      RATES(IL,62)*DYZ(IL,6)*DYY(IL,7) +                    &
                      RATES(IL,IHET+1)*DYY(IL, 2)*DYY(IL, 4)    +           &
                      RATES(IL,IHET+2)*H2O(IL)*DYY(IL, 2)*2. +              &
                      RATES(IL,IHET+3)*H2O(IL)*DYY(IL, 6)    +              &
                      RATES(IL,IHET+4)*DYY(IL, 6)*DYY(IL, 4)
!
!  The HNO3 production term from heterogeneous reaction BrONO2 + H2O
!  is sometimes rather fast leading to too much HNO3 extrapolated to the 
!  full timestep. Hence this term is neglected since it should be small in 
!  any case.      
!     *            RATES(IL,IHET+8)*DYY(IL,12)*H2O(IL)  
!
          DQT(IL,1) = RATES(IL,19)*DYZ(IL,7) + FOTO(IL,5)            
!       
!   N2O5
!
          DPT(IL,2) = RATES(IL,7)*DYZ(IL,5)*DYZ(IL,6)
          DQT(IL,2) = RATES(IL,8) + FOTO(IL,8) +                            &
                      RATES(IL,IHET+1)*DYY(IL, 4) +                         &
                      RATES(IL,IHET+2)*H2O(IL)
!
!  H2O2
!
          DPT(IL,3) = RATES(IL,34)*DYZ(IL,8)**2 +                           &
                      RATES(IL,35)*DYZ(IL,7)**2
          DQT(IL,3) = FOTO(IL,11) + RATES(IL,36)*DYZ(IL,7)
!
!  HCl
!
          DPT(IL,4) = DYZ(IL,10)*(RATES(IL,27)*DYZ(IL, 8)  +                &
                                  RATES(IL,33)*DYY(IL, 7)  +                &
                                  RATES(IL,25)*DYY(IL,18)) +                &
                      RATES(IL,69)*DYZ(IL,11)*DYZ(IL,7)
          DQT(IL,4) = RATES(IL,    26)*DYZ(IL, 7)    +                      &
                      RATES(IL,  IHET)*DYY(IL, 5)    +                      &
                      RATES(IL,IHET+1)*DYY(IL, 2)    +                      &
                      RATES(IL,IHET+4)*DYY(IL, 6)    +                      &
                      RATES(IL,IHET+5)*DYY(IL, 9)    +  FOTO(IL,33)
!
!  HOCl
!
          DPT(IL,5) = RATES(IL,    28)*DYZ(IL,11)*DYZ(IL, 8) +              &
                      RATES(IL,IHET+3)*DYY(IL, 6)*H2O(IL)
          DQT(IL,5) = FOTO(IL,15)                 +                         &
                      RATES(IL,  IHET)*DYY(IL, 4) +                         &
                      RATES(IL,IHET+6)*DYY(IL,11)
!
!  ClONO2
!
          DPT(IL,6) = RATES(IL,23)*DYZ(IL,11)*DYZ(IL,5)
          DQT(IL,6) = RATES(IL,24)*DYZ(IL, 2) + FOTO(IL,12) +               &
                      RATES(IL,IHET+3)*H2O(IL) +                            &
                      RATES(IL,IHET+4)*DYY(IL, 4)
!
!  H2CO
!
          DPT(IL,7) = DYZ(IL,16)*(RATES(IL,55)*DYZ(IL,4) +                  &
                    2.*RATES(IL,61)*DYZ(IL,16)) + FOTO(IL,27)*DYY(IL,13)
          DQT(IL,7) = FOTO(IL,13) + FOTO(IL,14) +                           &
                      RATES(IL,32)*DYZ(IL, 7) + RATES(IL,33)*DYZ(IL,10) +   &
                      RATES(IL,46)*DYZ(IL,12) + RATES(IL,62)*DYZ(IL,6)          
!
! Oy
!
          DPT(IL,8) = RATES(IL,51)*DYZ(IL,7)**2 +                           &
                   2.*FOTO(IL,1)*O2*RHO(IL) + FOTO(IL, 4)*DYZ(IL, 5) +      &
                      FOTO(IL,19)*DYZ(IL,13) + FOTO(IL,28)*DYY(IL,20)    
          DQT1 = RAT(IL,1)*                                                 &
                (RATES(IL, 4)*DYZ(IL, 4) +                                  &
                RATES(IL,10)*DYZ(IL, 7) + RATES(IL,13)*DYZ(IL, 9) +         &
                RATES(IL,15)*DYZ(IL, 8) + RATES(IL,20)*DYZ(IL,10) +         &
                RATES(IL,39)*DYZ(IL,12))
!     &      RATES(IL,15)*DYZ(IL, 8) + RATES(IL,20)*DYZ(IL,10))
          DQT2 = RAT(IL,2)*                                                 &
            (2.*RATES(IL, 2)*DYZ(IL, 1) + RATES(IL, 5)*DYZ(IL, 5) +         &
                RATES(IL,11)*DYZ(IL, 8) + RATES(IL,12)*DYZ(IL, 7) +         &
                RATES(IL,21)*DYZ(IL,11) + RATES(IL,24)*DYY(IL, 6) +         &
                RATES(IL,44)*DYY(IL,9)  + RATES(IL,64)*DYZ(IL, 4))
          DQT3 = RAT(IL,3)*                                                 &
               (RATES(IL, 9)*H2O(IL) +                                      &
                RATES(IL,30)*DYY(IL,18))
!      DQT4 = RATES(IL,39)*OZONE(IL)*ISTEP_CHEM
!      IF(DQT4.GT.1.0) DQT4 = 1.0 
!      DQT(IL,8) = DQT1 + DQT2 + DQT3 + DQT4*DYZ(IL,12)/(ISTEP_CHEM*OZONE(IL))
          DQT(IL,8) = DQT1 + DQT2 + DQT3 
!
!  HOBr
!
          DPT(IL,9) = RATES(IL,    42)*DYZ(IL,13)*DYZ(IL, 8) +              &
                      RATES(IL,IHET+8)*H2O(IL)*DYY(IL,12)
          DQT(IL,9) = FOTO(IL,18) + RATES(IL,44)*DYZ(IL,2) +                &
                      RATES(IL,IHET+5)*DYY(IL, 4)          +                &
                      RATES(IL,IHET+7)*DYY(IL,11)
!  
!  HNO4
!
          DPT(IL,10) = RATES(IL,49)*DYZ(IL,8)*DYZ(IL,5)
          DQT(IL,10) = FOTO(IL,25) + RATES(IL,50)*DYZ(IL,7) +               &
                       RATES(IL,63)
!
!  HBr
!
          DPT(IL,11) = RATES(IL,45)*DYZ(IL,12)*DYZ(IL,8) +                  &
                       RATES(IL,46)*DYZ(IL,12)*DYY(IL,7)
          DQT(IL,11) = RATES(IL,    47)*DYZ(IL, 7) +                        &
                       RATES(IL,IHET+6)*DYY(IL, 5) +                        &
                       RATES(IL,IHET+7)*DYY(IL, 9)
!
!  BrONO2
!
          DPT(IL,12) = RATES(IL,40)*DYZ(IL,13)*DYZ(IL,5)
          DQT(IL,12) = FOTO(IL,17) + RATES(IL,IHET+8)*H2O(IL)
!
!  CH3OOH
!
          DPT(IL,13) = RATES(IL,56)*DYZ(IL,16)*DYZ(IL,8)
          DQT(IL,13) = FOTO(IL,27) + RATES(IL,57)*DYZ(IL,7)
!
! CO
!
          DPT(IL,14) = DYY(IL,7)*                                           &
                  (RATES(IL,32)*DYZ(IL,7) + RATES(IL,33)*DYZ(IL,10) +       &
                   RATES(IL,46)*DYZ(IL,12) + RATES(IL,62)*DYZ(IL,6) +       &
                   FOTO(IL,13) + FOTO(IL,14))
          DQT(IL,14) = RATES(IL,52)*DYZ(IL,7)
!
!  NOY + NAT: Normal Stratospheric terms, relax towards ANOY in the 
!  troposphere
!
          AGESQ = 1.0E-2
          IF(DAGESQ(IL).LT.AGESQ) THEN
            DPT(IL,15) = ANOY(JLX,KL)*RHO(IL)*VTAU10
            DQT(IL,15) = VTAU10
          ELSE
            DPT(IL,15) = 2.*RATES(IL,66)*DYZ(IL,3)*DYY(IL,20)
            DQT(IL,15) = 2.*DYZ(IL,18)*                                     &
                   (RATES(IL,67)*DYZ(IL,4)+RATES(IL,71)*DYZ(IL,5))/DYY(IL,15)
          ENDIF
!
!  Cly: solved separately after call to chemistry; set to zero here
!
          DPT(IL,16) = 0.0
          DQT(IL,16) = 0.0
!
!  Bry: solved separately after call to chemistry; set to zero here
!
          DPT(IL,17) = 0.0
          DQT(IL,17) = 0.0
!
!  CH4: Production rate assumes relaxation to a uniform tropospheric value 
!  with a 10 day time scale
!  DAGESQ < 0.01 is used to denote the troposphere. 
!
          IF(DAGESQ(IL).LT.AGESQ) THEN
            DPT(IL,18) = CH4*VTAU10*RHO(IL)
            DQT(IL,18) = VTAU10 
          ELSE
            DPT(IL,18) = 0.0
!  
! Mesospheric loss of CH4 by lyman alpha photolysis is taken to be the H2O 
! photolysis rate weighted by the ratio of the cross sections
! 
            DQT(IL,18) = RATES(IL,25)*DYZ(IL,10) + RATES(IL,30)*DYZ(IL,3) +   &
                         RATES(IL,31)*DYZ(IL,7) + FOTO(IL,32)*1.53/1.83
          ENDIF
!
!  h2ostrat: not used
!
          DPT(IL,19) = 0.0
          DQT(IL,19) = 0.0
!
!  N2O: Production rate assumes relaxation to a uniform tropospheric value 
!  with a 10 day time scale
!  DAGESQ < 0.01 is used to denote the troposphere. 
!
          IF(DAGESQ(IL).LT.AGESQ) THEN
            DPT(IL,20) = AN2O*VTAU10*RHO(IL)
            DQT(IL,20) = VTAU10
          ELSE
            DPT(IL,20) = RATES(IL,71)*DYZ(IL,5)*DYZ(IL,18)
            DQT(IL,20) = (RATES(IL,66)+RATES(IL,70))*DYZ(IL,3) +           &
                          FOTO(IL,28)
          ENDIF
!
!  AGE OF AIR: Increase at the rate of 1 per second outside the troposphere.
!  Results expressed in years. Relax towards zero with a 10 
!  day timescale in the troposphere, denoted by DAGESQ less than 0.01 
!
          IF(DAGESQ(IL).LT.AGESQ) THEN
            DPT(IL,21) = 0.0
            DQT(IL,21) = VTAU10
          ELSE
            DPT(IL,21) = VTAUYR*RHO(IL)
            DQT(IL,21) = 0.0
          ENDIF
!
!  H2O + ICE
!
          DP1(IL) = DYZ(IL,7)*( RATES(IL,16)*DYZ(IL,8) +                    &
             RATES(IL,19)*DYY(IL,1) + RATES(IL,26)*DYY(IL,4) +              &
             RATES(IL,31)*DYY(IL,18) + RATES(IL,32)*DYY(IL,7) +             &
             RATES(IL,36)*DYY(IL,3) + RATES(IL,50)*DYY(IL,10) +             &
             RATES(IL,47)*DYY(IL,11) +                                      &
             RATES(IL,51)*DYZ(IL,7) + RATES(IL,54)*H2*RHO(IL) +             &
             RATES(IL,57)*DYY(IL,13) + RATES(IL,60)*DYZ(IL,17) )
          DQ1(IL) = FOTO(IL,32) + RATES(IL,9)*DYZ(IL,3)
  990   ENDDO
! ---------------------------------------------------------------------*
! ---
! ---  INTEGRATE NCHEM SPECIES EXPLICITLY
! ---
! ---------------------------------------------------------------------*
        DO 980 NC = 1,NCHEM
          DO  IL = 1,MERIDS
            T0 = DQT(IL,NC)/DVDT
            IF(T0 .LT. 1.0E-6) THEN 
              DYY1(IL,NC) = (DYY(IL,NC)*DVDT + DPT(IL,NC))/(DQT(IL,NC)+DVDT)
! The following line changes answers but should be equivalent to the above..
!              DYY1(IL,NC) = (DYY(IL,NC) + DPT(IL,NC)/DVDT)/(T0+1.0)
            ELSE
              minusT0 = -1.0*T0
              T1 = DPT(IL,NC)/DQT(IL,NC)
!              DYY1(IL,NC) = T1 + (DYY(IL,NC) - T1)*EXP(-1.0*T0)
              DYY1(IL,NC) = T1 + (DYY(IL,NC) - T1)*EXP(minusT0)
            ENDIF
          ENDDO
  980   ENDDO!CONTINUE 
        DO  IL = 1,MERIDS
          H2O1(IL) = (H2O(IL)*DVDT + DP1(IL))/(DQ1(IL)+DVDT)
!
! SPECIAL CODE TO ENSURE THAT THE BROMINE SPECIES ARE UNDER CONTROL DURING
! THE RAPIDLY VARYING PERIOD AT DAWN AND DUSK
!
          BRX = DYY1(IL,9) + DYY1(IL,11) + DYY1(IL,12) + DYZ(IL,12) +       &
                DYZ(IL,13) + DYZ(IL,15)
          IF(DYY1(IL,17).EQ.0.0) THEN
            RATIO = 0.0
          ELSE
            RATIO = BRX/DYY1(IL,17)
          ENDIF
          IF(RATIO.GT.1.0) THEN
            DYY1(IL,9) = DYY1(IL,9)/RATIO
            DYY1(IL,11) = DYY1(IL,11)/RATIO
            DYY1(IL,12) = DYY1(IL,12)/RATIO
            DYZ(IL,12) = DYZ(IL,12)/RATIO
            DYZ(IL,13) = DYZ(IL,13)/RATIO
            DYZ(IL,15) = DYZ(IL,15)/RATIO
          ENDIF
!
!  SPECIAL CODE TO ENSURE THAT THE HETEROGENEOUS SCHEME DOES NOT PRODUCE
! EXCESSIVE HNO3. IF HNO3 EXCEEDS NOy - 2.*N2O5 SET HNO3 TO  NOy - 2.*N2O5
!
          DIFF = DYY1(IL,15) - 2.0*DYY1(IL,2)
          IF (DYY1(IL,1).GT.DIFF) DYY1(IL,1) = DIFF
        ENDDO
 1000 ENDDO
!
!  COMPUTE INCREMENT TO O3_CHEM DIAGNOSTIC TRACER
!
      DO IL = 1,MERIDS
        O3_PROD(IL) = O3_PROD(IL) +                                     &
             (DPT(IL,8) - DQT(IL,8)*DYY(IL,8))*RAT(IL,1)*DELT/RHO(IL)
      ENDDO
      IF(ITIME_LEFT .GT. 0) THEN
        DO NC = 1,NCHEM
          DO IL = 1,MERIDS
            DYY(IL,NC) = DYY1(IL,NC)
            H2O(IL) = H2O1(IL)
          ENDDO
        ENDDO
        CYCLE
      ELSE
        EXIT
      ENDIF
 1100 ENDDO
! ---------------------------------------------------------------------*
! ---
! --- COMPUTE CHEMICAL INCREMENTS FOR THIS SUPER-TIMESTEP   
! ---
! ---------------------------------------------------------------------*
!
!  GET CHEM TENDENCY AND CONVERT TO MIXING RATIOS
!
      DO IL = 1,MERIDS
        VRHO = 1.0/(RHO(IL)*DT)
        DO NC = 1,NCHEM
          CH_TEND(IL,NC) = CH_TEND(IL,NC) +                             &
                           (DYY1(IL,NC) - DYY0(IL,NC))*VRHO
       ENDDO
          H2O_TEND(IL) = H2O_TEND(IL) + (H2O1(IL) - H2O0(IL))*VRHO
      ENDDO
!
! SET OZONE CONCENTRATION
!
      DO IL = 1,MERIDS
        OZONE(IL) = DYZ(IL,1)/RHO(IL)
      ENDDO
!
#ifdef _ALLOC
      deallocate(RATES)
      deallocate(FOTO)
      deallocate(DYY)
      deallocate(DYZ)
      deallocate(DQT)
      deallocate(DPT)
      deallocate(RAT)
      deallocate(DYZ1)
      deallocate(DYY0)
      deallocate(DYY1)
      deallocate(COND)
      deallocate(GAMMA)
#endif
      RETURN
END SUBROUTINE CHEMISTRY
!         
SUBROUTINE DCLY_DT(age,dfdage,tropc,tracer1,cly,bry,dclydt,dbrydt, &
           agefact2,merids,lats,levels,itime)       
!----------------------------------------------------------------------*
!      age:  Age of air since stratospheric entry in years.                
!
!      The order of dfdage array elements is:
!
!      1. CFC11 as a fraction of its tropospheric concentration.        
!      2.-7.  CFC12, CFC113, CCL4, CH3CL, CH3CCL3, HCFC22 
!      8  Bry as a fraction of its tropospheric concentration
!
!      TROPC contains the tropospheric concentrations for 1950.0-2100.0 
!      at annual intervals. The species are ordered as dfdage elements 1-7 
!      and Total organic Cl and Br.
!----------------------------------------------------------------------*
      IMPLICIT NONE
      
      integer, intent(in) :: MERIDS,LATS,LEVELS,itime(6) 
      REAL, intent(in) :: age(merids,lats,levels),dfdage(90,48,8),TROPC(151,9), &
            cly(merids,lats,levels),bry(merids,lats,levels),         &
            tracer1(merids,lats,levels)
      REAL, intent(inout) :: dclydt(merids,lats,levels), &
            dbrydt(merids,lats,levels)
      REAL, intent(in) :: AGEFACT2
            
      
      INTEGER MDAY(12)
      INTEGER IANN,IMON,IDAY,IMON1,IT1,IT2,IC,KL,JL, &
              IL,IY,IY1,IM,MX
      REAL  CLWEIGHT(7),                    &
            dfdtau(merids,lats,8),          &
            cfc(merids,lats,9),cfct(9,2)

      REAL  YSTART,TFACT,AGEFACT1,DX1,DX2,TIME,DT1,DT2,SUM1,   &
            SUM2,FACTOR,YY,Y1,Y2,clytot,brytot
      DATA MDAY/31,28,31,30,31,30,31,31,30,31,30,31/
      DATA YSTART/1950.0/
      DATA CLWEIGHT/ 3.0, 2.0, 3.0, 4.0, 1.0, 3.0, 1.0 / 
!
!      Interpolate between months for TAB array
!
      iann = itime(1)
      imon = itime(2)
      iday = itime(3)
      tfact = 1.0/(365.25*86400.0)
!
! agefact1 and agefact2 are fudge factors, to correct for low model 
! atmospheric age in the computation of the halogen rates of change.
!
      agefact1 = 1.25
      DX1 = IDAY/MDAY(IMON)
      IF(IMON.EQ.2.AND.4*(IANN/4).EQ.IANN) DX1 = IDAY/29.0
      DX2 = 1.0 - DX1
      IMON1 = IMON + 1
      IF(IMON1.EQ.13) IMON1 = 1
!
!  Compute multiplying factor for missing CFCs, and include factor for 
!  conversion of rates to a per second rate. 
!
      TIME = REAL(IANN) + (REAL(imon)-1.0)/12.0
      IT1 = INT(TIME - YSTART) + 1
      IF(IT1 .LT. 1) IT1 = 1
      IF(IT1 .GT. 150) IT1 = 150
      IT2 = IT1 + 1
      DT1 = TIME - REAL(IT1) - (YSTART - 1.0)
      DT2 = 1.0 - DT1
      sum1 = 0.0
      sum2 = 0.0
      do ic = 1,7
        sum1 = sum1 + clweight(ic)*tropc(it1,ic)
        sum2 = sum2 + clweight(ic)*tropc(it2,ic)
      enddo
      factor = (dt2*tropc(it1,8) + dt1*tropc(it2,8))*tfact/             &
                (sum1*dt2 + sum2*dt1)
      do 500 kl = 1,levels
!----------------------------------------------------------------------*
!  determine age and dfdage at air parcel positions using tracer
!----------------------------------------------------------------------*
        do jl = 1,lats
          do il = 1,merids
            yy = tracer1(il,jl,kl)
            IF(YY.LT.1.)  YY =  1.01
            IF(YY.GT.89.99) YY = 89.99
            IY = YY
            Y1 = YY - IY
            Y2 = 1. - Y1
            IY1 = IY + 1
            do ic = 1,8
            dfdtau(il,jl,ic) = y2*dfdage(iy,kl,ic) + y1*dfdage(iy1,kl,ic)
            enddo
          enddo
        enddo
!----------------------------------------------------------------------*
!  compute CFCs at time t - age
!----------------------------------------------------------------------*
!      DO 100 JL = 1,lats
!      DO 100 IL = 1,merids
        DO JL = 1,lats
          DO IL = 1,merids
            DO IM = 1,2
              MX = IM + IMON - 1
              IF(MX.EQ.13) MX = 1
              TIME = REAL(IANN) + (REAL(MX)-1.0)/12.0 - age(IL,JL,kl)*agefact1  
!
!      Interpolate the tropospheric gas concs in time
!
              IT1 = INT(TIME - YSTART) + 1
              IF(IT1 .LT. 1) IT1 = 1
              IF(IT1 .GT. 150) IT1 = 150
              IT2 = IT1 + 1
              DT1 = TIME - REAL(IT1) - (YSTART - 1.0)
              DT2 = 1.0 - DT1
              DO IC=1,9
                CFCT(IC,IM) = TROPC(IT1,IC)*DT2 + TROPC(IT2,IC)*DT1                
              ENDDO! IC
            ENDDO! IM
            DO IC = 1,9
              CFC(IL,JL,IC) =  CFCT(IC,1)*DX2 +  CFCT(IC,2)*DX1  
            ENDDO !IC
          ENDDO
        ENDDO
! 100  CONTINUE
!----------------------------------------------------------------------*
!  Finally compute dclydt and dbrydt
!----------------------------------------------------------------------*
        do jl = 1,lats
          do il = 1,merids
            dclydt(il,jl,kl) = 0.0
            do ic = 1,7
              dclydt(il,jl,kl) = dclydt(il,jl,kl) +             &
                  1.0e-12*factor*dfdtau(il,jl,ic)*clweight(ic)* &
                  cfc(il,jl,ic)*agefact2
            enddo
            clytot = cfc(il,jl,8)*1.0e-12
            if(cly(il,jl,kl).ge.clytot) dclydt(il,jl,kl) = 0.0
            dbrydt(il,jl,kl) = 1.0e-12*tfact*dfdtau(il,jl,8)*cfc(il,jl,9)
            brytot = 1.0e-12*cfc(il,jl,9)
            if(bry(il,jl,kl).ge.brytot) dbrydt(il,jl,kl) = 0.0
          enddo
        enddo
 500  ENDDO
      RETURN
END SUBROUTINE DCLY_DT
!
SUBROUTINE REACTION_RATES (RHO,TEMP,VTEMP,PRESS,RATES,            &
                           MERIDS,IHET,NHET)
!-------------------------------------------------------------------------
!
!   THIS SUBROUTINE DETERMINES THE HOMOGENEOUS GAS PHASE REACTION RATES
!
!-------------------------------------------------------------------------
IMPLICIT NONE
INTEGER  , intent(in) :: MERIDS,IHET,NHET
REAL, intent(in), dimension(MERIDS) ::  RHO,TEMP,VTEMP,PRESS
REAL, intent(inout) :: RATES(MERIDS,IHET+NHET-1)
!
INTEGER  IL,IC
REAL TEMP3(MERIDS)
!
REAL  O2,R0,R1,RX,RX1,RX2
!
  DATA O2/0.2095/      
!
  DO IL = 1,MERIDS
    TEMP3(IL) = TEMP(IL)/300.0
  ENDDO   
!
!  RATES TAKEN FROM JPL 2006 
!
  DO IL = 1,MERIDS
!   REACTION 1 O+O2+M --> O3+M
      RATES(IL,1) = 6.0E-34*TEMP3(IL)**(-2.40)*O2*RHO(IL)*RHO(IL)       
!   REACTION 2 O+O3 --> O2+O2
      RATES(IL,2) = 8.0E-12*EXP(-2060*VTEMP(IL))
!   REACTION 3 OSD+M --> O+M N2 & O2 3RD BODY RATES ALLOWED FOR       
      RATES(IL,3) = RHO(IL)*(6.93E-12*EXP(70*VTEMP(IL)) +              &
                               1.677E-11*EXP(110*VTEMP(IL)))
!   REACTION 4 NO+O3 --> NO2+O2
      RATES(IL,4) = 3.0E-12*EXP(-1500*VTEMP(IL))
!   REACTION 5 NO2+O --> NO+O2
      RATES(IL,5) = 5.1E-12*EXP(210*VTEMP(IL))
!   REACTION 6 NO2+O3 --> NO3+O2
      RATES(IL,6) = 1.2E-13*EXP(-2450*VTEMP(IL))
!   REACTION 7 NO2+NO3 + M --> N2O5 + M
      R0 = 2.0E-30*RHO(IL)*TEMP3(IL)**(-4.4)
      R1 = 1.4E-12*TEMP3(IL)**(-0.7)
      RX = ALOG10(R0/R1)
      RATES(IL,7) =(R0/(1.0 + R0/R1))*                                &
                          (0.6**(1.0/(1.0 + (RX*RX))))
!   REACTION 8 N2O5+M --> NO3+NO2+M
      RATES(IL,8) = RATES(IL,7)/(2.7E-27*EXP(11000.0*VTEMP(IL)))
!   REACTION 9 OSD+H2O --> OH+OH
      RATES(IL,9) = 1.63E-10*EXP(60.0*VTEMP(IL))
!   REACTION 10 OH+O3 --> HO2+O2
      RATES(IL,10) = 1.7E-12*EXP(-940.0*VTEMP(IL))
!   REACTION 11 O+HO2 --> OH+O2
      RATES(IL,11) = 3.0E-11*EXP(200.0*VTEMP(IL))
!   REACTION 12 O+OH --> H+O2
      RATES(IL,12) = 2.2E-11*EXP(120.0*VTEMP(IL))
!   REACTION 13 H+O3 --> OH+O2
      RATES(IL,13) = 1.4E-10*EXP(-470.0*VTEMP(IL))
!   REACTION 14 H+O2+M --> HO2+M : MULTIPLIED BY O2 NUMBER DENSITY    
      R0 = 4.4E-32*RHO(IL)*TEMP3(IL)**(-1.3)
      R1 = 4.7E-11*TEMP3(IL)**(-0.2)
      RX = ALOG10(R0/R1)
      RATES(IL,14) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
      RATES(IL,14) = RATES(IL,14)*O2*RHO(IL)
!   REACTION 15 HO2+O3 --> OH+O2+O2
      RATES(IL,15) = 1.0E-14*EXP(-490.0*VTEMP(IL))
!   REACTION 16 OH+HO2 --> H2O+O2
      RATES(IL,16) = 4.8E-11*EXP(250.0*VTEMP(IL))
!   REACTION 17 OH+NO2+M --> HONO2+M
!   INCLUDE ALSO OH+NO2+M --> HOONO + M
      R0 = 1.8E-30*RHO(IL)*TEMP3(IL)**(-3.0)
      R1 = 2.8E-11
      RX = ALOG10(R0/R1)
      RATES(IL,17) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
      R0 = 9.1E-32*RHO(IL)*TEMP3(IL)**(-3.9)
      R1 = 4.2E-11*TEMP3(IL)**(-0.5)
      RX = ALOG10(R0/R1)
      RATES(IL,17) = RATES(IL,17) + (R0/(1.0 + R0/R1))*           &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 18 NO+HO2 --> NO2+OH
      RATES(IL,18) = 3.5E-12*EXP(250.0*VTEMP(IL))
!   REACTION 19 OH+HNO3 --> H2O+NO3 : SPECIAL TREATMENT               
      RX1 = 6.5E-34*RHO(IL)*EXP(1335.0*VTEMP(IL))            
      RX2 = 2.7E-17*EXP(2199.0*VTEMP(IL))            
      RATES(IL,19) =  2.4E-14*EXP(460.0*VTEMP(IL)) +                        &
                      RX1/(1.0 + RX1/RX2)           
!   REACTION 20 CL+O3 --> CLO+O2
      RATES(IL,20) = 2.3E-11*EXP(-200.0*VTEMP(IL))
!   REACTION 21 O+CLO --> CL+O2
      RATES(IL,21) = 2.8E-11*EXP(85.0*VTEMP(IL))
!   REACTION 22 CLO+NO --> NO2+CL
      RATES(IL,22) = 6.4E-12*EXP(290.0*VTEMP(IL))
!   REACTION 23 CLO+NO2+M --> CLONO2+M
      R0 = 1.8E-31*RHO(IL)*TEMP3(IL)**(-3.4)
      R1 = 1.5E-11*TEMP3(IL)**(-1.9)
      RX = ALOG10(R0/R1)
      RATES(IL,23) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 24 O+CLONO2--> CLO+NO3
      RATES(IL,24) = 2.9E-12*EXP(-800.0*VTEMP(IL))
!   REACTION 25 CL+CH4 (+O2) --> HCL+CH3O2
      RATES(IL,25) = 7.3E-12*EXP(-1280.0*VTEMP(IL))
!   REACTION 26 OH+HCL --> H2O+CL
      RATES(IL,26) = 2.6E-12*EXP(-350.0*VTEMP(IL))
!   REACTION 27 CL+HO2 --> HCL+O2
      RATES(IL,27) = 1.8E-11*EXP(170.0*VTEMP(IL))
!   REACTION 28 CLO+HO2 --> HOCL+O2
      RATES(IL,28) = 2.7E-12*EXP(220.0*VTEMP(IL))
!   REACTION 29 CLO+OH --> HO2+CL
      RATES(IL,29) = 7.4E-12*EXP(270.0*VTEMP(IL))
!   REACTION 30 OSD+CH4 --> OH + CH3O2
      RATES(IL,30) = 1.5E-10*EXP(0.0*VTEMP(IL))
!   REACTION 31 OH+CH4 --> H2O + CH3O2
      RATES(IL,31) = 2.45E-12*EXP(-1775.0*VTEMP(IL))
!   REACTION 32 H2CO+OH (+O2) --> H2O + HO2 + CO
      RATES(IL,32) = 5.5E-12*EXP(125.0*VTEMP(IL))
!   REACTION 33 H2CO+CL (+O2) --> HCL+ HO2 + CO        
      RATES(IL,33) = 8.1E-11*EXP(-30.0*VTEMP(IL))
!   REACTION 34 HO2+HO2 --> H2O2+O2
      RATES(IL,34) = 3.5E-13*EXP(430.0*VTEMP(IL)) +                        &
                     1.7E-33*EXP(1000.0*VTEMP(IL))*RHO(IL)                
!   REACTION 35 OH+OH+M --> H2O2+M
      R0 = 6.9E-31*RHO(IL)*TEMP3(IL)**(-1.0)
      R1 = 2.6E-11*TEMP3(IL)**0.0
      RX = ALOG10(R0/R1)
      RATES(IL,35) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 36 OH+H2O2 --> HO2+H2O
      RATES(IL,36) = 2.9E-12*EXP(-160.0*VTEMP(IL))
!   REACTION 37 CLO+CLO+M --> CL2O2 + M
      R0 = 1.6E-32*RHO(IL)*TEMP3(IL)**(-4.5)
      R1 = 2.0E-12*TEMP3(IL)**(-2.4)
      RX = ALOG10(R0/R1)
      RATES(IL,37) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 38 CL2O2+M --> 2CLO + M
      RATES(IL,38) = RATES(IL,37)/(9.3E-28*EXP(8835.0*VTEMP(IL)))  
!   REACTION 39 BR + O3 --> BRO + O2
      RATES(IL,39) = 1.7E-11*EXP(-800.0*VTEMP(IL))
!   REACTION 40 BRO + NO2 + M --> BRONO2 + M
      R0 = 5.2E-31*RHO(IL)*TEMP3(IL)**(-3.2)
      R1 = 6.9E-12*TEMP3(IL)**(-2.9)
      RX = ALOG10(R0/R1)
      RATES(IL,40) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 41 BRO + CLO --> BR + CL + 02
      RATES(IL,41) =  2.3E-12*EXP(260.0*VTEMP(IL))   
!   REACTION 42 BRO + HO2 --> HOBR + O2
      RATES(IL,42) = 4.5E-12*EXP(460.0*VTEMP(IL))
!   REACTION 43 BRO + NO --> BR + NO2
      RATES(IL,43) = 8.8E-12*EXP(260.0*VTEMP(IL))
!   REACTION 44 HOBR + O --> BRO + OH
      RATES(IL,44) = 1.2E-10*EXP(-430.0*VTEMP(IL))
!   REACTION 45 BR + HO2 --> HBR + O2
      RATES(IL,45) = 4.8E-12*EXP(-310.0*VTEMP(IL))
!   REACTION 46 BR + H2CO (+O2) --> HBR + HO2 + CO
      RATES(IL,46) = 1.7E-11*EXP(-800.0*VTEMP(IL))
!   REACTION 47 HBR + OH --> BR + H2O
      RATES(IL,47) = 5.5E-12*EXP(200.0*VTEMP(IL))
!   REACTION 48 BRO + CLO --> BRCL + O2
      RATES(IL,48) = 4.1E-13*EXP(290.0*VTEMP(IL))
!   REACTION 49 HO2 + NO2 + M --> HNO4 + M
      R0 = 2.0E-31*RHO(IL)*TEMP3(IL)**(-3.4)
      R1 = 2.9E-12*TEMP3(IL)**(-1.1)
      RX = ALOG10(R0/R1)
      RATES(IL,49) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 50 HNO4 + OH --> H2O + NO2 + O2
      RATES(IL,50) = 1.3E-12*EXP(380.0*VTEMP(IL))
!
! REACTIONS ADDED FOR TROPOSPHERIC CHEMISTRY
!
!   REACTION 51 OH + OH --> H2O + O
      RATES(IL,51) = 1.8E-12
!   REACTION 52 CO + OH (+O2) --> HO2 +CO2
      R0 = 5.9E-33*RHO(IL)*TEMP3(IL)**(-1.4)
      R1 = 1.1E-12*TEMP3(IL)**(1.3)
      RX = ALOG10(R0/R1)
      RATES(IL,52) =    (R0/(1.0 + R0/R1))*                           &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
      R0 = 1.5E-13*TEMP3(IL)**0.6
      R1 = 2.1E9*TEMP3(IL)**(6.1)/RHO(IL)       
      RX = ALOG10(R0/R1)       
      RATES(IL,52) = RATES(IL,52)+ (R0/(1.0 + R0/R1))*                 &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 53 NO3 + NO --> 2NO2
      RATES(IL,53) = 1.5E-11*EXP(170.0*VTEMP(IL))
!   REACTION 54 H2 + OH (+O2) --> HO2 + H2O
      RATES(IL,54) = 5.5E-12*EXP(-2000.0*VTEMP(IL))
!   REACTION 55 CH3O2 + NO (+O2) --> H2CO + NO2 + HO2
      RATES(IL,55) = 2.8E-12*EXP(300.0*VTEMP(IL))
!   REACTION 56 CH3O2 + HO2 --> CH3OOH + O2
      RATES(IL,56) = 4.1E-13*EXP(750.0*VTEMP(IL))
!   REACTION 57 CH3OOH + OH --> CH3O2 + H2O 
      RATES(IL,57) = 3.8E-12*EXP(200.0*VTEMP(IL))
!   REACTION 58 NO2 + NO3 --> NO + NO2 + O2
      RATES(IL,58) = 4.5E-14*EXP(-1260.0*VTEMP(IL))
!   REACTION 59 NO3 + NO3 --> 2NO2 + O2 
      RATES(IL,59) = 8.5E-13*EXP(-2450.0*VTEMP(IL))
!   REACTION 60 OH + HONO --> H2O + NO2
      RATES(IL,60) = 1.8E-11*EXP(-390.0*VTEMP(IL))
!   REACTION 61 CH3O2 + CH3O2 --> 2H2CO + 2OH
      RATES(IL,61) = 9.5E-14*EXP(390.0*VTEMP(IL))
!   REACTION 62 NO3 + HCHO (+O2) --> HNO3 + CO + HO2
      RATES(IL,62) = 5.8E-16
!   REACTION 63 HO2NO2 + M --> HO2 + NO2
      RATES(IL,63) = RATES(IL,49)*4.762E26*EXP(-10900*VTEMP(IL))
!   REACTION 64 NO + O + M --> NO2 + M
      R0 = 9.0E-31*RHO(IL)*TEMP3(IL)**(-1.5)
      R1 = 3.0E-11
      RX = ALOG10(R0/R1)
      RATES(IL,64) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 65 OH + NO + M --> HONO + M
      R0 = 7.0E-31*RHO(IL)*TEMP3(IL)**(-2.6)
      R1 = 3.6E-11*TEMP3(IL)**(-0.1)
      RX = ALOG10(R0/R1)
      RATES(IL,65) =(R0/(1.0 + R0/R1))*                               &
                          (0.6**(1.0/(1.0 + (RX*RX)))) 
!   REACTION 66 N2O + OSD --> 2NO
      RATES(IL,66) = 6.7E-11*EXP(20.0*VTEMP(IL))
!   REACTION 67 N + NO --> N2 + O
      RATES(IL,67) = 2.1E-11*EXP(100.0*VTEMP(IL))
!   REACTION 68 N + O2 --> NO + O
      RATES(IL,68) = 1.5E-11*EXP(-3600.0*VTEMP(IL))
!   REACTION 69 CLO+OH --> --> HCL+O2
      RATES(IL,69) = 6.0E-13*EXP(-230.0*VTEMP(IL))
!   REACTION 70 N2O + OSD --> N2 + O2
      RATES(IL,70) = 4.7E-11*EXP(20.0*VTEMP(IL))
!   REACTION 71 N + NO2 --> N2O + O
      RATES(IL,71) = 5.8E-12*EXP(220.0*VTEMP(IL))
  ENDDO
  DO IC=IHET,IHET+NHET-1
    DO IL=1,MERIDS
      RATES(IL,IC) = 0.0                                            
    ENDDO
  ENDDO
  RETURN
END SUBROUTINE REACTION_RATES

SUBROUTINE PSC(TEMP,PRESS,RHO,DYY,H2O,H2SO4,ANAT,AICE,COND,TICE,  &
           WH2SO4,AM,AW,ALIQ,RMEAN,ASAT,RNAT,RICE,MERIDS,NCHEM,MYPE)
!---------------------------------------------------------------------------
!  CALCULATES IMPORTANT TERMS FOR BINARY AEROSOL/SAT/NAT/ICE SCHEME.
! THE MAIN OUTPUT IS THE SURFACE AREA OR PARTICLE RADIUS
!----------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, intent(in)                           :: MERIDS,NCHEM,MYPE
REAL, intent(in)                              :: ANAT, AICE
REAL, intent(in),    dimension(MERIDS)        :: PRESS, RHO, TEMP, H2SO4

REAL, intent(inout), dimension(MERIDS)        :: H2O, TICE, WH2SO4, AM, AW, &
                                                 ALIQ, RMEAN, ASAT, RNAT, RICE
REAL, intent(inout), dimension(MERIDS,3)      :: COND
REAL, intent(inout), dimension(MERIDS, NCHEM) :: DYY

INTEGER IL
REAL DENS(MERIDS),VHET(MERIDS,4),WF(MERIDS,2),PH2O(MERIDS)
REAL AVGDR,RR,PI,ADROP,SIGMAL,TSAT,ABT,AMT,BMT,PX,C2,   &
     A1,A2,A3,A4,C3,C4,CLIMIT,P0H2O,Y1,Y2,T1,RMODE,RMODESAT,SANAT
!----------------------------------------------------------------------------
!
! AVGDR IS THE RECIPROCAL OF THE AVOGADRO CONSTANT
! RR IS THE GAS CONSTANT
!
!----------------------------------------------------------------------------
  DATA AVGDR,RR/1.66056547E-24, 8.3144/,PI/3.1415926536/
!----------------------------------------------------------------------------
!
!   ADROP, ANAT AND AICE ARE THE (FIXED) NUMBER OF DROPS OR PARTICLES PER CC
!   IN THE POLAR STRATOSPHERC CLOUDS. A LOG NORMAL SIZE DISTRIBUTION IS ASSUMED
!   WITH SIGMA = 1.8 GIVING LOG(SIGMA)**2 = 0.3455  
!
!----------------------------------------------------------------------------
  DATA ADROP/10.0/,SIGMAL/0.34549316/      
!          
    DO 40 IL=1,MERIDS
!----------------------------------------------------------------------------
!
! COMPUTE SOLID CONDENSED MATTER: COND(IL,1:3) -- SAT, NAT, ICE
!               
!----------------------------------------------------------------------------
      PH2O(IL) = H2O(IL)*PRESS(IL)/(RHO(IL)*101325.0)
      COND(IL,1) = 0.0
      COND(IL,2) = 0.0
      COND(IL,3) = 0.0
      TSAT = 3236.0/(11.502 - ALOG10(PH2O(IL)*760.0))           
      IF(TEMP(IL).LT.TSAT) COND(IL,1) = H2SO4(IL)        
      ABT = 39.1104 - 11397.0/TEMP(IL) + 9.179E-3*TEMP(IL)                
      AMT = -2.7836 - 8.8E-4*TEMP(IL)
      BMT = -2.1249 + ALOG10(PH2O(IL)*101325.0)
      PX = AMT*BMT + ABT
      C2 = DYY(IL,1) - (RHO(IL)*100.0/PRESS(IL))*10.0**(PX)
      IF(C2.GT.0.0) COND(IL,2) = C2
      TICE(IL) = 2668.70/(10.4310 - ALOG10(760.0*PH2O(IL)))               
      A1 = 7.5502 - 2668.7/TEMP(IL)
      C3 = H2O(IL) - (RHO(IL)*101325.0/PRESS(IL))*10.0**(A1)          
      C4 = (RHO(IL)*101325.0/PRESS(IL))*10.0**(A1)          
      IF(C3.GT.0.0) COND(IL,3) = C3
      DYY(IL,1) = DYY(IL,1) - COND(IL,2)
      DYY(IL,15) = DYY(IL,15) - COND(IL,2) 
      CLIMIT = RHO(IL)*1.E-15
      IF(DYY(IL,15).LT.CLIMIT) DYY(IL,15) = CLIMIT
      H2O(IL) = H2O(IL) - COND(IL,3) 
 40 ENDDO!CONTINUE
!-------------------------------------------------------------------------
!
!  COMPUTE WEIGHT % H2SO4. FROM TABAZADEHET AL., GRL, 24, 1931-1934, 1997.
!  TABLE A1 OF SHIAET AL. JGR, 106, 24,529-24,274, 2001
!
!-------------------------------------------------------------------------
    DO 50 IL = 1,MERIDS
      P0H2O =EXP(18.452406985 - 3505.1578807/TEMP(IL) -                &
             330918.55082/(TEMP(IL)**2) + 12725068.262/(TEMP(IL)**3))    
      AW(IL) = H2O(IL)*PRESS(IL)*0.01/(RHO(IL)*P0H2O)
      IF(AW(IL).LE.0.05) THEN
        Y1 = 12.37208932*AW(IL)**(-0.16125516114) - 30.490657554*AW(IL)  &
            -  2.1133114241
        Y2 = 13.455394705*AW(IL)**(-0.1921312255) - 34.285174607*AW(IL)        &
            - 1.7620073078
      ELSEIF(AW(IL).LT.0.85.AND.AW(IL).GT.0.05) THEN
        Y1 = 11.820654354*AW(IL)**(-0.20786404244) - 4.807306373*AW(IL)   &
            -  5.1727540348
        Y2 = 12.891938068*AW(IL)**(-0.23233847708) - 6.4261237757*AW(IL)   &
            - 4.9005471319
      ELSE
        Y1 = -180.06541028*AW(IL)**(-0.38601102592)                     &
            - 93.317846778*AW(IL) + 273.88132245
        Y2 = -176.95814097*AW(IL)**(-0.36257048154)                     &
             - 90.469744201*AW(IL) + 267.45509988
      ENDIF              
      AM(IL) = Y1 + (TEMP(IL) - 190.0)*(Y2 - Y1)/70.0
      WH2SO4(IL) = 9800.0*AM(IL)/(98.0*AM(IL) + 1000.0)
      WF(IL,1) = 0.01*WH2SO4(IL) 
      WF(IL,2) = 0.0
 50 ENDDO!CONTINUE
!---------------------------------------------------------------------------
!
!  COMPUTE DENSITY OF BINARY AEROSOL        
!
!---------------------------------------------------------------------------
    CALL DENSITY(WF,TEMP,DENS,MERIDS)
!---------------------------------------------------------------------------
!
!  COMPUTE VOLUME OF BINARY AEROSOL/SAT/NAT/ICE 
!  1.6, 1.35 and 0.928  are the densities of SAT, NAT and ICE     
!
!---------------------------------------------------------------------------
    DO 100 IL=1,MERIDS
      T1 = H2SO4(IL)*PRESS(IL)/(RHO(IL)*TEMP(IL)*RR)
      VHET(IL,1) = T1*98.076E-6/(WF(IL,1)*DENS(IL))             
      VHET(IL,2) = COND(IL,1)*170.1*AVGDR/1.6
      VHET(IL,3) = COND(IL,2)*117.1*AVGDR/1.35
      VHET(IL,4) = COND(IL,3)*18.02*AVGDR/0.928
100 ENDDO
!---------------------------------------------------------------------------
!
!  COMPUTE PARTICLE PARAMETERS FROM WHICH THE HETEROGENEOUS REACTION RATES
!  ARE DETERMINED; ASSUME SURFACE AREA FROM SAT IS LIMITED BY NAT AMOUNT   
!
!---------------------------------------------------------------------------
    A1 =EXP(-4.5*SIGMAL)
    A2 =EXP(0.5*SIGMAL)
    A3 =EXP(2.0*SIGMAL)
    A4 = 1.33333333*PI*ADROP
    DO 150 IL=1,MERIDS
      RMODE = (VHET(IL,1)*A1/A4)**0.33333333
      RMEAN(IL) = RMODE*A2
      ALIQ(IL) = 3.0*A4*(RMODE**2)*A3
      IF(RMEAN(IL).LT.1.0E-12) RMEAN(IL) = 1.0E-12
      RMODESAT = (VHET(IL,2)*A1/A4)**0.33333333
      ASAT(IL) = 3.0*A4*(RMODESAT**2)*A3
      RNAT(IL) = (VHET(IL,3)/(1.33333333*PI*ANAT))**0.33333333
      RICE(IL) = (VHET(IL,4)/(1.33333333*PI*AICE))**0.33333333
      SANAT = 4.0*PI*ANAT*RNAT(IL)**2
      ASAT(IL) = ASAT(IL) - SANAT
      ALIQ(IL) = ALIQ(IL) - SANAT
      IF(ASAT(IL).LT.0.0) ASAT(IL) = 0.0
      IF(ALIQ(IL).LT.0.0) ALIQ(IL) = 0.0
150 ENDDO
    RETURN
END SUBROUTINE PSC
!
SUBROUTINE GAMMA_AER(TEMP,PRESS,RHO,DYY,H2O,WH2SO4,AM,AW,RMEAN,   &
                     GAMMA,RTT,MERIDS,NCHEM,NHET,MYPE)
!-------------------------------------------------------------------------
!
! SUBROUTINE TO CALCULATE REACTION PROBABILITIES ON SULPHATE AEROSOL
! BASED ON JPL'03 RECOMMENDATION
!  
!-------------------------------------------------------------------------
  IMPLICIT NONE
  INTEGER, intent(in) ::  MERIDS,NCHEM,NHET,MYPE
  REAL, intent(in) :: TEMP(MERIDS),PRESS(MERIDS),RHO(MERIDS),DYY(MERIDS,NCHEM), &
        H2O(MERIDS),WH2SO4(MERIDS),AW(MERIDS),AM(MERIDS),           &
        RMEAN(MERIDS)
  REAL, intent(inout) :: RTT(MERIDS), GAMMA(MERIDS,NHET)
!
  INTEGER  IL
  REAL  AMH2SO4(MERIDS),XH2SO4(MERIDS),VISC(MERIDS),AACID(MERIDS),  &
        TEMP2(MERIDS)
!-------------------------------------------------------------------------
!  
! CALCULATE H2SO4 MOLARITY (AMH2SO4), MOLE FRACTION (XH2SO4), 
!  VISCOSITY (VISC) AND ACID ACTIVITY (AACID)
!  TABLE A2 OF SHIAET AL. JGR, 106, 24,529-24,274, 2001.
! 
!-------------------------------------------------------------------------
  REAL  T2,Z1,Z2,Z3,RHOX,AA,X,T1,T3,AKH,AKH2O,AKHYDR,DIFF,SCLONO2,  &
         CCLONO2,GAMMAB1,HHCL,AKHCL,Q1,RQ,A1,FCLONO2,GAMMARXN,      &
         GAMMABHCL,GAMMAS,FHCL,GAMMASP,GAMMABHCLP,GAMMAB,GCLONO2,   &
         SHOCL,HHOCL,FHOCL,WT,AK0,AK1,AK2,T0,HCLONO2,   &
         AMHCL,AKHOCL,CHOCL
!
!  The parameterisations used here break down below about 185K, so the
!  temperature is here limited to 185K and above (TEMP2).
!
    DO IL = 1,MERIDS
      TEMP2(IL) = TEMP(IL)
      IF(TEMP2(IL).LT.185.0) TEMP2(IL) = 185.0
      RTT(IL) = SQRT(TEMP(IL))
    ENDDO
    DO 20 IL = 1,MERIDS
      T2 = TEMP2(IL)**2
      Z1 = 0.12364 - 5.6E-7*T2
      Z2 = -0.02954 + 1.814E-7*T2
      Z3 = 2.343E-3 - 1.487E-6*TEMP2(IL) - 1.324E-8*T2
      RHOX = 1.0 + Z1*AM(IL) + Z2*AM(IL)**1.5 + Z3*AM(IL)**2
      AMH2SO4(IL) = RHOX*WH2SO4(IL)/9.8
      XH2SO4(IL) = WH2SO4(IL)/                                          &
         (WH2SO4(IL) + (100.0 - WH2SO4(IL))*98.0/18.0)
      AA = 169.5 + 5.18*WH2SO4(IL) - 0.0825*WH2SO4(IL)**2 +             &
           3.27E-3*WH2SO4(IL)**3 
      T0 = 144.11 + 0.166*WH2SO4(IL) - 0.015*WH2SO4(IL)**2 +            &
           2.18E-4*WH2SO4(IL)**3
      X = TEMP2(IL)**(-1.43) 
      VISC(IL) = AA*X*EXP(448.0/(TEMP2(IL) - T0))
      T1 = 60.51 - 0.095*WH2SO4(IL) + 0.0077*WH2SO4(IL)**2              &
           - 1.61E-5*WH2SO4(IL)**3
      T2 =  (-805.89 + 253.05*WH2SO4(IL)**0.076)/RTT(IL)  
      T3 =   (1.76 + 2.52E-4*WH2SO4(IL)**2)*RTT(IL)
      AACID(IL) =EXP(T1 + T2 - T3)
 20 ENDDO!CONTINUE
!-------------------------------------------------------------------------
!
!  CALCULATE REACTION PROBABILITES FOR CLONO2 + H2O AND CLONO2 + HCL AND 
!  HENRY'S LAW COEFFICIENTS.
!  TABLE A3 OF SHIAET AL. JGR, 106, 24,529-24,274, 2001.
!
!  The following formulation for the water activity is from Shi et al.,
!  but the differences between their parameterisation and that calculated
!  using the actua model H2O is not large.
!
!      awx = exp((-69.775*xh2so4(il) - 18253.7*xh2so4(il)**2 + 
!     2     31072.2*xh2so4(il)**3 - 25668.8*xh2so4(il)**4)*
!     3     (1.0/temp(il) - 26.9033/(temp(il)**2)))      
!      AKHYDR = AWx*(AKH2O + AKH*AACID(IL))
!
!-------------------------------------------------------------------------
    DO 30 IL = 1,MERIDS
      AKH = 1.22E12*EXP(-6200.0/TEMP2(IL))
      AKH2O = 1.95E10*EXP(-2800.0/TEMP2(IL))     
      AKHYDR = AW(IL)*(AKH2O + AKH*AACID(IL))
      DIFF = 5.0E-8*TEMP2(IL)/VISC(IL)
      SCLONO2 = 0.306 + 24.0/TEMP2(IL)
      HCLONO2 = 1.6E-6*EXP(4710.0/TEMP2(IL) - SCLONO2*AMH2SO4(IL))
      CCLONO2 = 1474.*TEMP2(IL)**0.5
      GAMMAB1 = (4.0*HCLONO2*0.082*TEMP2(IL)/CCLONO2)*                  &
                (DIFF*AKHYDR)**0.5
!
      HHCL = (0.094 - 0.61*XH2SO4(IL) + 1.2*XH2SO4(IL)**2)*             &
        EXP(-8.68 + (8515. - 10718.*XH2SO4(IL)**0.7)/TEMP2(IL))
      AMHCL = HHCL*DYY(IL,4)*PRESS(IL)/(RHO(IL)*101325.0)
      AKHCL = 7.9E11*AACID(IL)*DIFF*AMHCL
      Q1 = (DIFF/(AKHYDR + AKHCL))**0.5
      RQ = RMEAN(IL)/Q1
      A1 = RQ + 0.312*RQ**2
      FCLONO2 = A1/(3.0 + A1)
      GAMMARXN = FCLONO2*GAMMAB1*(1.0 + AKHCL/AKHYDR)**0.5
      GAMMABHCL = GAMMARXN*AKHCL/(AKHCL + AKHYDR)
      GAMMAS = 66.12*EXP(-1374./TEMP2(IL))*HCLONO2*AMHCL
      IF(DYY(IL,4).NE.0.0) THEN
        FHCL = 1.0/(1.0 + 0.612*(GAMMAS + GAMMABHCL)*DYY(IL,6)/DYY(IL,4))
      ELSE
        FHCL = 0.0
      ENDIF
      GAMMASP = FHCL*GAMMAS
      GAMMABHCLP = FHCL*GAMMABHCL     
      GAMMAB = GAMMABHCLP + GAMMARXN*AKHYDR/(AKHCL + AKHYDR)
      GCLONO2 = 1.0/(1.0 + 1.0/(GAMMASP + GAMMAB))
      GAMMA(IL,5) = GCLONO2*(GAMMASP + GAMMABHCLP)/                     &
            (GAMMASP + GAMMAB)
      GAMMA(IL,4) = GCLONO2 - GAMMA(IL,5)
!-------------------------------------------------------------------------
!
!  CALCULATE REACTION PROBABILITES FOR HOCL + HCL AND HENRY'S LAW COEFFICIENTS.
!  TABLE A4 OF SHIAET AL. JGR, 106, 24,529-24,274, 2001.
!
!-------------------------------------------------------------------------
      SHOCL = 0.0776 + 59.18/TEMP2(IL)
      HHOCL = 1.91E-6*EXP(5862.4/TEMP2(IL) - SHOCL*AMH2SO4(IL))
      DIFF = 6.4E-8 *TEMP2(IL)/VISC(IL) 
      AKHOCL = 1.25E9*AACID(IL)*DIFF*AMHCL
      CHOCL = 2009.*RTT(IL)
      Q1 = (DIFF/AKHOCL)**0.5
      RQ = RMEAN(IL)/Q1
      A1 = RQ + 0.312*RQ**2
      FHOCL = A1/(3.0 + A1)
      GAMMARXN = (FHOCL*4.0*HHOCL*0.082*TEMP2(IL)/CHOCL)*               &
                 (DIFF*AKHOCL)**0.5
      GAMMA(IL,1) = 1.0/(1.0 + 1.0/(GAMMARXN*FHCL))
 30 ENDDO!CONTINUE
!-------------------------------------------------------------------------
!
!  CALCULATE REACTION PROBABILITES FOR N2O5 + H2O 
!  ROBINSONET AL. JGR, 102, 3583-3601, 1997.
!
!-------------------------------------------------------------------------
    DO 40 IL = 1,MERIDS
      WT = WH2SO4(IL)
      IF(WH2SO4(IL).GT.80.0) WT = 80.0
      AK0 = -25.5265 - 0.133188*WT + 0.00930846*WT**2 -                 &
              9.0194E-5*WT**3  
      AK1 = 9283.76 + 115.345*WT - 5.19258*WT**2 +                      &
           0.0483464*WT**3  
      AK2 = -851801. - 22191.2*WT + 766.916*WT**2 -                     &
           6.85427*WT**3
      GAMMA(IL,3) =EXP(AK0 + AK1/TEMP2(IL) + AK2/(TEMP2(IL)**2))  
 40 ENDDO!CONTINUE
!-------------------------------------------------------------------------
!
!  REACTION PROBABILITES FOR 
!       N2O5 + HCL
!       HOBR + HCL
!       HOCL + HBR 
!  NO RECOMMENDATION IN JPL '03, ASSUMED ZERO
!
!-------------------------------------------------------------------------
    DO IL = 1,MERIDS
      GAMMA(IL,2) = 0.0
      GAMMA(IL,6) = 0.0
      GAMMA(IL,7) = 0.0
    ENDDO
!-------------------------------------------------------------------------
!
!  REACTION PROBABILITES FOR HOBR + HBR 
!  ABBATT, JGR, 100, 14009-14017, 1995. 
!
!-------------------------------------------------------------------------
    DO IL = 1,MERIDS
      GAMMA(IL,8) = 0.25
    ENDDO
!-------------------------------------------------------------------------
!
!  REACTION PROBABILITES FOR BRONO2 + H2O 
!  USE JPL '03 RECOMMENDATION (HANSON PERS. COMM.)
!
!-------------------------------------------------------------------------
    DO IL = 1,MERIDS
      GAMMA(IL,9) = 1.0/(1.2422 + 1.0/(0.114 +EXP(29.24 -              &
                    0.396*WH2SO4(IL))))
    ENDDO
    RETURN
END SUBROUTINE GAMMA_AER
!
SUBROUTINE HETRATES(TEMP,PRESS,RHO,DYY,H2O,TICE,ANAT,AICE,ALIQ,        &
        RMEAN,RNAT,RICE,GAMMA,RTT,RATES,IHET,NHET,MERIDS,NCHEM,MYPE)

IMPLICIT NONE
!------------------------------------------------------------------------
!
!  This subroutine computes the equivalent 2nd order reaction rates for 
!  the heterogeneous reactions on aerosol, nat and ice.
!
!------------------------------------------------------------------------
      
  INTEGER, intent(in) :: IHET, NHET, MERIDS, NCHEM, MYPE
  REAL, intent(in) ::  TEMP(MERIDS),PRESS(MERIDS),RHO(MERIDS), &
       DYY(MERIDS,NCHEM), H2O(MERIDS),TICE(MERIDS),            &
       ALIQ(MERIDS),RMEAN(MERIDS),RNAT(MERIDS),RICE(MERIDS),   &
       GAMMA(MERIDS,NHET), RTT(MERIDS)
  REAL, intent(inout) :: RATES(MERIDS,IHET+NHET-1)
!
  INTEGER INN(9),IC,IL
  REAL AMW(9),GNAT(9),GICE(9),                          &
       CHEMC(MERIDS,NHET),CONST(NHET),                              &
       G2NAT(MERIDS,NHET),DELT(MERIDS),SICE(MERIDS)
!
  REAL  ANAT,AICE,PI,ANUM,ADEN
  DATA AMW/52.45, 108.00, 108.00, 97.45, 97.45, 96.91,              &
           52.45, 96.91, 141.91/               
!------------------------------------------------------------------------
!  AMW = MOLECULAR WEIGHT OF GAS PHASE SPECIES
!------------------------------------------------------------------------
  DATA PI/3.141592653589793/    
  DATA GNAT/0.1, 3.0E-3, 4.0E-4, 4.0E-3, 0.2,                    &
            0.0, 0.0, 0.0, 0.0/
  DATA GICE/0.2, 0.03, 0.02, 0.3, 0.3, 0.3,                         &
            0.03, 0.1, 0.3/
  DATA INN/4,4,0,0,4,4,11,11,0/
!------------------------------------------------------------------------
!  INN = INDEX NUMBER OF LIQUID/SOLID PHASE SPECIES (0= H2O)
!-------------------------------------------------------------------------
!     REACTION 70 HOCL + HCL --> H2O + CL2 (HETEROGENEOUS)              
!     REACTION 71 N2O5 + HCL --> HNO3 + CLNO2 (HETEROGENEOUS)           
!     REACTION 72 N2O5+H2O --> 2HNO3 (HETEROGENEOUS)
!     REACTION 73 CLONO2+H2O --> HOCL+HNO3 (HETEROGENEOUS)              
!     REACTION 74 CLONO2+HCL --> CL2+HNO3 (HETEROGENEOUS)               
!     REACTION 75 HOBR + HCL --> BRCL + H2O (HETEROGENEOUS)             
!     REACTION 76 HOCL + HBR --> BRCL + H2O (HETEROGENEOUS)             
!     REACTION 77 HOBR + HBR --> 2BR + H2O (HETEROGENEOUS)              
!     REACTION 78 BRONO2 + H2O --> HOBR + HNO3 (HETEROGENEOUS)
!     (THE FIRST MOLECULE ON THE LEFT HAND SIDE IS IN THE GAS PHASE,
!     THE SECOND MOLECULE IS IN THE LIQUID/SOLID PHASE)         
!
!-------------------------------------------------------------------------
!
!    aliq is the liquid surface area. const is sqrt(8R/(pi*mw)) for the
!    gaseous phase species, with the mean molecular
!    speed equal to const*sqrt(temp)
!
!-------------------------------------------------------------------------
    DO 60 IC = 1,NHET
      CONST(IC) = SQRT(8.0*8.3144E7/(PI*AMW(IC)))
      DO 50 IL = 1,MERIDS
        IF(INN(IC).EQ.0) THEN
          CHEMC(IL,IC) = H2O(IL)
        ELSE
          CHEMC(IL,IC) = DYY(IL,INN(IC))
        ENDIF
        RATES(IL,IHET-1+IC) =  0.0 
 50   ENDDO
 60 ENDDO
!-------------------------------------------------------------------------
!
!    Reactions IHET to IHET+NHET-1 on NAT
!
!-------------------------------------------------------------------------
    DO IC = 1,NHET
      DO IL = 1,MERIDS
        G2NAT(IL,IC) = GNAT(IC)
      ENDDO
    ENDDO
    DO 110 IL=1,MERIDS
      DELT(IL) = TEMP(IL) - TICE(IL)
      SICE(IL) =  10**(2668.70*(1.0/TEMP(IL) - 1.0/TICE(IL)))      
      IF (SICE(IL) .GT. 3.0) SICE(IL) = 3.0
      G2NAT(IL,4) =EXP(-9.03 + 2.81*SICE(IL))
      G2NAT(IL,5) = 1.0/(4.3478 + 1.4241*EXP(0.518*DELT(IL)))             
110 ENDDO
!
    DO IC=1,NHET
      DO IL=1,MERIDS
        ANUM = PI*CONST(IC)*RTT(IL)*G2NAT(IL,IC)*ANAT*RNAT(IL)**2
        ADEN = CHEMC(IL,IC)
        IF(ADEN.GT.0.0.AND.ANUM.GT.0.0)                               &
            RATES(IL,IHET-1+IC) =  ANUM/ADEN
      ENDDO
    ENDDO
!-------------------------------------------------------------------------
!
!    Reactions IHET to IHET+NHET-1 on ICE
!
!------------------------------------------------------------------------ 
    DO IC=1,NHET
      DO IL=1,MERIDS
        ANUM = PI*CONST(IC)*RTT(IL)*GICE(IC)*AICE*RICE(IL)**2  
        ADEN = CHEMC(IL,IC)                
        IF(ADEN.GT.0.0.AND.ANUM.GT.0.0)                               &
           RATES(IL,IHET-1+IC) = ANUM/ADEN 
      ENDDO
    ENDDO
!-------------------------------------------------------------------------
!
!    Reactions IHET to IHET+NHET-1 on AEROSOL
!    aliq is the liquid surface area. const is sqrt(8R/(pi*mw)) for the
!    gaseous phase species, with the mean molecular speed equal to 
!        const*sqrt(temp)
!
!-------------------------------------------------------------------------
    DO IC = 1,NHET
      DO IL = 1,MERIDS
        IF(CHEMC(IL,IC).GT.0.0)                                       &
           RATES(IL,IHET+IC-1) =  RATES(IL,IHET+IC-1) +               &
           0.25*GAMMA(IL,IC)*CONST(IC)*RTT(IL)*ALIQ(IL)/CHEMC(IL,IC)
      ENDDO
    ENDDO
!    
    RETURN
END SUBROUTINE HETRATES
!
SUBROUTINE DENSITY(WF,T,DENS,MERIDS)
IMPLICIT NONE
!
!    Density of ternary solution in g cm-3
!
  INTEGER, intent(in)    :: MERIDS
  REAL,    intent(in )   :: WF(MERIDS,2), T(MERIDS)
  REAL,    intent(inout) :: DENS(MERIDS)

  INTEGER IL
  REAL X(22),AMR(3)
  REAL W,WH,T2,V1,A1,A2,VS,VN,VMCAL
  DATA X/2.393284E-02,-4.359335E-05,7.961181E-08,0.0,-0.198716351,   &
         1.39564574E-03,-2.020633E-06,0.51684706,-3.0539E-03,        &
         4.505475E-06,-0.30119511,1.840408E-03,-2.7221253742E-06,    &
        -0.11331674116,8.47763E-04,-1.22336185E-06,0.3455282,        &
        -2.2111E-03,3.503768245E-06,-0.2315332,1.60074E-03,          &
        -2.5827835E-06/
  DATA AMR/0.05550622,0.01019576,0.01586899/
    DO 100 IL=1,MERIDS
      W = WF(IL,1) + WF(IL,2)
      WH = 1.0 - W
      T2 = T(IL)**2
      V1 = X(1) + X(2)*T(IL) + X(3)*T2 + X(4)*T2*T(IL)
      A1 = X(8) + X(9)*T(IL) + X(10)*T2
      A2 = X(11) + X(12)*T(IL) + X(13)*T2
      VS = X(5) + X(6)*T(IL) + X(7)*T2 + A1*W + A2*W**2
      A1 = X(17) + X(18)*T(IL) + X(19)*T2
      A2 = X(20) + X(21)*T(IL) + X(22)*T2
      VN = X(14) + X(15)*T(IL) + X(16)*T2 + A1*W + A2*W**2
      VMCAL = WH*V1*AMR(1) + VS*WF(IL,1)*AMR(2) + VN*WF(IL,2)*AMR(3)
      DENS(IL) = 1.0E-3/VMCAL
100 ENDDO
!
  RETURN
END SUBROUTINE DENSITY
     
SUBROUTINE PHOTO_RATES (VTEMP,OZONX,COSPHI,KL,PHOTO,OZON,COSP,    &
                        COSPHC,MERIDS,IRX,IRX4,FOTO,MYPE)
!---------------------------------------------------------------------- 
!   THIS SUBROUTINE CALCULATES THE PHOTOLYSIS RATES USING A LOOK-UP     
!   TABLE. THE RATES ARE PUT INTO FOTO AND ARE AS FOLLOWS               
!   1.  O2 + HV --> 2O
!   2.  O3 + HV --> OSD + O2
!   3.  O3 + HV --> O3P + O2
!   4.  NO2 + HV --> NO + O
!   5.  HNO3 + HV --> OH + NO2 (300K) (LATER OVERWRITTEN WITH VALUE     
!   6.  NO3 + HV --> NO2 + O
!   7.  NO3 + HV --> NO + O2
!   8.  N2O5 + HV --> NO2 + NO3  (300K)  (LATER OVERWRITTEN WITH VALUE  
!   9.  N2O5 + HV --> NO2 + NO3  (250K)               AT TEMP)          
!  10.  N2O5 + HV --> NO2 + NO3  (200K)
!  11.  H2O2 + HV --> 2OH
!  12.  CLONO2 + HV --> CL + NO3 (300K) (LATER OVERWRITTEN WITH VALUE   
!  13.  H2CO + HV --> H + HCO --> 2H + CO
!  14.  H2CO + HV --> H2 + CO
!  15.  HOCL + HV --> OH + CL
!  16.  CL2O2 + HV --> CLO2 + CL
!  17.  BRONO2 + HV --> BR + NO3
!  18.  HOBR + HV --> OH + BR
!  19.  BRO + HV --> BR + O
!  20.  HNO3 + HV --> OH + NO2 (250K)
!  21.  HNO3 + HV --> OH + NO2 (200K)
!  22.  CLONO2 + HV --> CL + NO3 (250K)
!  23.  CLONO2 + HV --> CL + NO3 (200K)
!  24.  BRCL + HV --> BR + CL
!  25.  HNO4 + HV --> HO2 + NO2
!  26.  HONO + HV --> OH + NO
!  27.  CH3OOH + HV --> H2CO + HO2 + OH
!  28.  N2O + HV --> N2 + OSD  (300K)  (LATER OVERWRITTEN WITH VALUE  
!  29.  N2O + HV --> N2 + OSD  (250K)               AT TEMP)          
!  30.  N2O + HV --> N2 + OSD  (200K)
!  31.  NO + HV --> N + O
!  32.  H2O + HV --> H + OH
!  33.  HCL + HV --> H + CL
!---------------------------------------------------------------------- 
IMPLICIT NONE
  INTEGER, intent(in) :: MERIDS,IRX,IRX4,MYPE, KL
  REAL, intent(in) :: PHOTO(IRX4,14,11,48), OZON(11,48),      &
       VTEMP(MERIDS),OZONX(MERIDS),COSPHI(MERIDS)
  REAL, intent(inout) ::   COSP(14), FOTO(MERIDS,IRX4), COSPHC(48)


  REAL  X1(MERIDS),X2(MERIDS),Y1(MERIDS),Y2(MERIDS),Z1(MERIDS),       &
        Z2(MERIDS),OZON1(MERIDS),OZON2(MERIDS),           &
        COSP1(MERIDS),COSP2(MERIDS)  
  INTEGER IL,ICOX,IOX,IR,IM
  INTEGER IOZ(MERIDS),ICOZ(MERIDS)    
  LOGICAL LCOZ(MERIDS)
!
!      REAL  ALG,CX1,AX1,BX1,CX2,AX1,BX2,AX3,BX3,XX,ZZ,P1,P2,XX2,AX2,CX3  
  REAL  ALG,CX1,AX1,BX1,CX2,AX2,BX2,CX3,AX3,BX3,XX,ZZ,P1,P2,XX2  
!    
!
     IF(KL.EQ.48) COSPHC(KL) = COSPHC(2)
     COSP(1) = -COSPHC(KL)*0.5
     COSP(2) = -COSPHC(KL)*0.25
     DO 110 IL = 1,MERIDS
       LCOZ(IL) = COSPHI(IL).GT.(-COSPHC(KL))
 110 ENDDO
     DO 120 IL = 1,MERIDS
       ICOZ(IL) = 1
       COSP1(IL) = COSP(1)
       COSP2(IL) = COSP(2)
 120 ENDDO
     DO 140 ICOX =  2,13
       DO 130 IL = 1,MERIDS
         IF(COSPHI(IL).GT.COSP(ICOX))  THEN
           ICOZ(IL) = ICOX
           COSP1(IL) = COSP(ICOX)
           COSP2(IL) = COSP(ICOX+1)
         ENDIF
 130   ENDDO
 140 ENDDO
!
!  INTERPOLATE LINEARLY IN COSPHI AND LINEARLY IN
!  OZONE AMOUNT FOR THE LOG OF THE PHOTODISSOCIATION RATE.              
!
     DO 150 IL = 1,MERIDS
       IOZ(IL) = 1
       OZON1(IL) = OZON(1,KL)
       OZON2(IL) = OZON(2,KL)
 150 ENDDO
     DO 170 IOX = 2,10
       DO 160 IL = 1,MERIDS
         IF(OZONX(IL).GT.OZON(IOX,KL))  THEN
           IOZ(IL) = IOX
           OZON1(IL) = OZON(IOX,KL)
           OZON2(IL) = OZON(IOX+1,KL)
         ENDIF
 160   ENDDO
 170 ENDDO
     DO 190 IL = 1,MERIDS
       X1(IL) = OZONX(IL) - OZON1(IL)
       X2(IL) = OZON2(IL) - OZONX(IL)
       XX = X1(IL) + X2(IL)
       IF(XX.EQ.0.0) THEN
         WRITE(6,*) 'TEST 1 ', MYPE
         WRITE(6,*) IL,KL,X1(IL),X2(IL),OZONX(IL),OZON1(IL),             &
                    OZON2(IL),IOZ(IL)
       ENDIF
       Y1(IL) = 1./(X1(IL) + X2(IL))
       Z1(IL) = COSPHI(IL) - COSP1(IL)
       Z2(IL) = COSP2(IL) - COSPHI(IL)
       ZZ = Z1(IL) + Z2(IL)
       IF(ZZ.EQ.0.0) THEN
         WRITE(6,*) 'TEST 2 ', MYPE
         WRITE(6,*) IL,KL,Z1(IL),Z2(IL),COSPHI(IL),COSP1(IL),            &
                    COSP2(IL),ICOZ(IL)
       ENDIF
       Y2(IL) = 1./ZZ
 190 ENDDO
     DO 200 IL = 1,MERIDS
       X1(IL) = X1(IL)*Y1(IL)
       X2(IL) = X2(IL)*Y1(IL)
       Z1(IL) = Z1(IL)*Y2(IL)
       Z2(IL) = Z2(IL)*Y2(IL)
 200 ENDDO
     DO 220 IR = 1,IRX4
       DO 210 IL = 1,MERIDS
         FOTO(IL,IR) = 0.
 210   ENDDO
 220 ENDDO
     DO 240 IL = 1,MERIDS
       IF(LCOZ(IL))  THEN
         DO 230 IR = 1,IRX4
           P1 = PHOTO(IR,ICOZ(IL)+1,IOZ(IL)+1,KL)*X1(IL) +                 &
                X2(IL)*PHOTO(IR,ICOZ(IL)+1,IOZ(IL),KL)
           P2 = PHOTO(IR,ICOZ(IL),IOZ(IL)+1,KL)*X1(IL) +                   &
             X2(IL)*PHOTO(IR,ICOZ(IL),IOZ(IL),KL)
           ALG =  P1*Z1(IL) + P2*Z2(IL)
           FOTO(IL,IR) = ALG
 230     ENDDO
       ENDIF
 240 ENDDO
!
!  PHOTODISSOCIATION RATES FOR TEMPERATURE DEPENDENT N2O5 CROSS-SECTION 
!  FIT A QUADRATIC TO THE VALUES AT 300K, 250K AND 200K                 
!
     DO 250 IL = 1,MERIDS
       IF(LCOZ(IL))  THEN
         XX = 5.0 - VTEMP(IL)*1000.0
         XX2 = XX*XX
         DO 245 IM = 0,3
           CX1 = FOTO(IL,10+IRX*IM)
           AX1 = (9.*FOTO(IL,8+IRX*IM) - 15.*FOTO(IL,9+IRX*IM) +             &
                  6.*CX1)*0.1
           BX1 = FOTO(IL,9+IRX*IM) - AX1 - CX1
           FOTO(IL,8+IRX*IM) = AX1*XX2 + BX1*XX + CX1
           CX2 = FOTO(IL,21+IRX*IM)      
           AX2 = (9.*FOTO(IL,5+IRX*IM) - 15.*FOTO(IL,20+IRX*IM) +            &
                  6.*CX2)*0.1
           BX2 = FOTO(IL,20+IRX*IM) - AX2 - CX2
           FOTO(IL,5+IRX*IM) = AX2*XX2 + BX2*XX + CX2
           CX3 = FOTO(IL,23+IRX*IM)      
           AX3 = (9.*FOTO(IL,12+IRX*IM) - 15.*FOTO(IL,22+IRX*IM) +           &
                  6.*CX3)*0.1
           BX3 = FOTO(IL,22+IRX*IM) - AX3 - CX3
           FOTO(IL,12+IRX*IM) = AX3*XX2 + BX3*XX + CX3
           CX1 = FOTO(IL,30+IRX*IM)
           AX1 = (9.*FOTO(IL,28+IRX*IM) - 15.*FOTO(IL,29+IRX*IM) +           &
                  6.*CX1)*0.1
           BX1 = FOTO(IL,29+IRX*IM) - AX1 - CX1
           FOTO(IL,28+IRX*IM) = AX1*XX2 + BX1*XX + CX1
 245     ENDDO
       ENDIF
 250 ENDDO
     DO 300 IL = 1,MERIDS
       IF(LCOZ(IL))  THEN  
         DO IM = 0,3
           DO IR = 1,IRX
             IF(FOTO(IL,IR+IM*IRX).GT.10.0) FOTO(IL,IR+IM*IRX) = 10.0
             IF(FOTO(IL,IR+IM*IRX).LT.-100.0) FOTO(IL,IR+IM*IRX) = -100.0
             FOTO(IL,IR+IM*IRX) =EXP(FOTO(IL,IR+IM*IRX))  
           ENDDO
         ENDDO
       ENDIF 
 300 ENDDO
     RETURN
END SUBROUTINE PHOTO_RATES
!/ ----------------------------------------------------------------     
!      
SUBROUTINE ZEN2(ITIME,DT,DRAD,DARG,DELTA,CDX)         
!
!  THIS SUBROUTINE CALCULATES THE COSINE OF THE SOLAR ZENITH ANGLE,     
!  COZEN, ALLOWING FOR LEAP YEARS AND THEELLIPTICITY OF THEEARTH'S    
!  ORBIT.
! ELONG IS THE LONGITUDE IN DEG.E. ,XLAT THE LATITUDE NORTH. GMT IS   
!  THE G.M.T. AND DAY IS THE DAY IN THE MONTH. TIMING INFO. IS GIVEN BY 
!  ITIME.
!
IMPLICIT NONE
INTEGER, intent(in) :: ITIME(6)
REAL, intent(in)    :: DT,DRAD
REAL, intent(inout)   :: DARG,DELTA,CDX

REAL DGMT,DP,D1,D2,DF1,DTY,DR,DET,    &
      DAL1,DAL,CD1,DTD
!
INTEGER  IYEAR,IMON,IDAY
!
    IYEAR = ITIME(1)
    IMON = ITIME(2)
    IDAY = ITIME(3)
    DGMT = ITIME(4) + ITIME(5)/60.0 + (ITIME(6) + 0.5*DT)/3600.0
    IF(DGMT.LT.24.)  GOTO 5
    DGMT = DGMT - 24.
    IDAY = IDAY + 1.
  5 DP = 365.*IYEAR + IDAY + 31.*(IMON-1.)
    if ((IMON - 2) <= 0) then
      D1 = (IYEAR - 1.)*0.25
      D2 = (IYEAR - 1.)*1.E-2 + 1.
      DF1 = DP + INT(D1) - INT(0.75*INT(D2))
    else
      D1 = 4.E-1*IMON + 2.3
      D2 = IYEAR*1.E-2
      DF1 = DP - INT(D1) + INT(25.*D2) - INT(0.75*(INT(D2) + 1.))       
    endif
    DTY = (DF1 - 693960.)/36525.
    DTD = 100.*DTY
    DR = (2.7969668E02 + 3.6000768925E04*DTY + 3.03E-4*DTY*DTY)*    &
      DRAD
    DET = (-9.32701E01  - 1.42E-1*DTD)*SIN(DR) + (-4.3194497E02 +    &
      3.3E-2*DTD)*COS(DR) + 5.965E02*SIN(2.*DR) - 2.*COS(2.*DR)       &
    + 4.2E00*SIN(3.*DR) + 1.93E01*COS(3.*DR) - 1.28E01*SIN(4.*DR)   
    DAL1 = DET*DRAD/240.0
    DAL = DR - DAL1
!
!  DELTA IS THE SOLAR DECLINATION AND COZEN IS THE COSINE OF THE SOLAR  
!  ZENITH ANGLE
!
    DELTA = ATAN(4.336E-1*SIN(DAL))
    DARG = (DGMT - 12.0)*15.*DRAD + DAL1
    CD1 = DTD*DRAD*360.0
    CDX = 1.00011 + 0.034221*COS(CD1) + 0.001280*SIN(CD1) +            &
          0.000719*COS(2.*CD1) + 0.000077*SIN(2.*CD1)                 
    RETURN
END SUBROUTINE ZEN2
!
SUBROUTINE SEDIMENT(ANAT,AICE,ANOY,AH2O,AHNO3,P_FIELD,LEVELS,     &
                    PRESS,DT,MYPE)      
!
!  CALCULATES SEDIMENTATION RATES OF TYPE I AND TYPE 2 PARTICLES               
!  AND VERTICALLY ADVECTS MODEL NAT AND ICE
!
  INTEGER, intent(in) ::  P_FIELD, LEVELS, MYPE
  REAL, intent(inout) :: ANAT(P_FIELD,LEVELS),AICE(P_FIELD,LEVELS), &
       ANOY(P_FIELD,LEVELS),AH2O(P_FIELD,LEVELS),                   &
       AHNO3(P_FIELD,LEVELS)
  REAL, intent(in) :: PRESS(P_FIELD,LEVELS)
  REAL, intent(in) :: DT   


  REAL  SNATS(P_FIELD,LEVELS),SICES(P_FIELD,LEVELS),                 &
        F1(P_FIELD,LEVELS),F2(P_FIELD,LEVELS), &
        PNAT(P_FIELD),PICE(P_FIELD),PNAT2(P_FIELD),PICE2(P_FIELD),        &
        ANAT2(P_FIELD,LEVELS),AICE2(P_FIELD,LEVELS),                 &
        ANATMAX(P_FIELD),AICEMAX(P_FIELD)
!
!  DATA: V1 IS SEDIMENTATION VELOCITY (M/S)
!  OF ICE PARTICLES V2 IS SEDIMENTATION VELOCITY OF         
!  NAT PARTICLES, ASSUMED RADII R1, R2.
!  AM1, AM2 ARE THE MOLECULAR WEIGHTS AND RHO1, RHO2 ARE THE DENSITIES 
!  OF THE PSCs IN G/CM3
!
!      DATA  V1, V2/1.27E-2, 1.39E-4/
!      DATA  R1, R2/7.0E-6, 0.5E-6/
!      DATA AM1,AM2/18.0,117.0/
!      DATA RHO1,RHO2/0.928, 1.35/
  REAL ::   V1=1.27E-2,  V2=1.39E-4
  REAL ::   R1=7.0E-6,   R2=0.5E-6
  REAL ::  AM1=18.0,    AM2=117.0
  REAL :: RHO1=0.928,  RHO2=1.35
  REAL :: RATIO, TEMP, PFRAC, DZ, CONST, D1, D2, FIXNAT, FIXICE
  INTEGER :: KL, JL

!
!  CALCULATE FRACTION OF NAT PARTICLES USED AS TYPE 2 CORES (F1)               
!  AND FRACTION OF NAT PARTICLES THAT REMAIN AS TYPE 1 CORES (F2)              
!  DETERMINE MAXIMUM NAT AND ICE TO APPLY LIMITERS TO ADVECTED AMOUNTS
!
    ANATMAX(:) = 0.0
    AICEMAX(:) = 0.0
    RATIO = AM1*RHO2/(AM2*RHO1)*(R2/R1)**3
    DO 20 KL = 1,LEVELS 
      DO 10 JL = 1,P_FIELD 
        IF(AICE(JL,KL).LT.1.0E-18) THEN 
          AICE(JL,KL) = 0.       
        ENDIF
        IF(ANAT(JL,KL).LT.1.0E-18)  THEN
          ANAT(JL,KL) = 0.
          F1(JL,KL) = 0.
          F2(JL,KL) = 0.
        ELSE
          F1(JL,KL) = AICE(JL,KL)*RATIO/ANAT(JL,KL)
          IF(F1(JL,KL).GT.1.0)  F1(JL,KL) = 1.0
          F2(JL,KL) = 1.0 - F1(JL,KL)
        ENDIF
        IF(ANAT(JL,KL).GT.ANATMAX(JL)) ANATMAX(JL) = ANAT(JL,KL)
        IF(AICE(JL,KL).GT.AICEMAX(JL)) AICEMAX(JL) = AICE(JL,KL)
 10   ENDDO
 20 ENDDO
!
! VERTICALLY ADVECT NAT AND ICE. NOTE THAT PART OF NAT IS ADVECTED             
! AT TYPE 2 RATE AND THE REMAINDER AT TYPE 1 RATE. CALCULATE DESCENT IN
!  1 TIMESTEP; USE APPROXIMATE VERTICAL DISPLACEMENT BETWEEN LAYERS
!
    TEMP = 195.0
    PNAT(:) = 0.0
    PICE(:) = 0.0
    DO 50 KL = 2,LEVELS
      DO 40 JL = 1,P_FIELD
        PFRAC = PRESS(JL,KL)/PRESS(JL,KL-1)
        DZ = 29.26*TEMP*ALOG(PFRAC)
        CONST = DT/DZ
        D1 = ANAT(JL,KL) - ANAT(JL,KL-1)
        D2 = AICE(JL,KL) - AICE(JL,KL-1)         
        SNATS(JL,KL) = -CONST*D1*(V1*F1(JL,KL) + V2*F2(JL,KL))
        SICES(JL,KL) = -CONST*D2*V1
        PNAT(JL) = PNAT(JL) + PRESS(JL,KL)*ANAT(JL,KL)
        PICE(JL) = PICE(JL) + PRESS(JL,KL)*AICE(JL,KL)
 40   ENDDO!CONTINUE
 50 ENDDO!CONTINUE
!
!  set sedimented nat and ice to zero at top and bottom
!
    SNATS(:,1) = 0.0
    SICES(:,1) = 0.0 
    SNATS(:,LEVELS) = 0.0
    SICES(:,LEVELS) = 0.0
    DO 70 KL = 1,LEVELS
      DO 60 JL = 1,P_FIELD
        ANAT2(JL,KL) = ANAT(JL,KL) + SNATS(JL,KL)
        AICE2(JL,KL) = AICE(JL,KL) + SICES(JL,KL)
!
!  APPLY LIMITERS TO NEW NAT AND ICE
!
        IF(ANAT2(JL,KL).LT.0.0) ANAT2(JL,KL) = 0.0
        IF(AICE2(JL,KL).LT.0.0) AICE2(JL,KL) = 0.0
        IF(ANAT2(JL,KL).GT.ANATMAX(JL)) ANAT2(JL,KL) = ANATMAX(JL)
        IF(AICE2(JL,KL).GT.AICEMAX(JL)) AICE2(JL,KL) = AICEMAX(JL)
 60   ENDDO!CONTINUE
 70 ENDDO!CONTINUE
!
! APPLY MASS FIXER
!
    PNAT2(:) = 0.0
    PICE2(:) = 0.0
    DO 90 KL = 1,LEVELS
      DO 80 JL = 1,P_FIELD
        PNAT2(JL) = PNAT2(JL) + PRESS(JL,KL)*ANAT2(JL,KL)
        PICE2(JL) = PICE2(JL) + PRESS(JL,KL)*AICE2(JL,KL)
 80   ENDDO!CONTINUE
 90 ENDDO!CONTINUE
    DO 110 JL = 1,P_FIELD
      IF(PNAT2(JL).EQ.0.0) THEN 
        FIXNAT = 1.0
      ELSE 
        FIXNAT = PNAT(JL)/PNAT2(JL)
      ENDIF
      IF(PICE2(JL).EQ.0.0) THEN
        FIXICE = 1.0
      ELSE
        FIXICE = PICE(JL)/PICE2(JL)
      ENDIF
      DO 100 KL = 1,LEVELS
        ANAT2(JL,KL) = ANAT2(JL,KL)*FIXNAT 
        AICE2(JL,KL) = AICE2(JL,KL)*FIXICE 
100   ENDDO
110 ENDDO
!
!  ADJUST NOY AND H2O TENDENCY FIELDS
!
    DO KL = 1,LEVELS
      DO JL = 1,P_FIELD
        ANOY(JL,KL) = ANOY(JL,kl) + (ANAT2(JL,KL) - ANAT(JL,KL))/DT 
        AHNO3(JL,KL) = AHNO3(JL,kl) + (ANAT2(JL,KL) - ANAT(JL,KL))/DT 
        AH2O(JL,KL) = AH2O(JL,KL) + (AICE2(JL,KL) - AICE(JL,KL))/DT 
     ENDDO
   ENDDO
   RETURN
END SUBROUTINE SEDIMENT
!
SUBROUTINE GHGS(ITIME,ISC,CH4,AN2O,CH4TR,AN2OTR)
!
!  FINDS THE CONCENTRATIONS OF THE GHGS CH4 AND N2O ACCORDING TO THE TIME 
!  ITIME AND SCENARIO ISC, CHOSEN FROM SRES
!
  INTEGER, intent(in) :: ITIME(6),ISC
  REAL, intent(inout) :: CH4,AN2O,CH4TR,AN2OTR

  REAL TIMEN2O(22),AN2O_A1B(22),AN2O_A2(22),AN2O_B1(22)
  REAL TIMECH4(36),CH4_A1B(36),CH4_A2(36),CH4_B1(36)
  REAL TIME, SECS, X1,X2, Z1, C1, C2
  INTEGER IMON(12), IT
  DATA  IMON/0,31,59,90,120,151,181,212,243,273,304,334/
  DATA TIMEN2O/                                                     &
     1950.5, 1976.5, 1980.5, 1984.5, 1988.5, 1990.5,                &
     1992.5, 1994.5, 1996.5, 1998.5, 2000.5, 2010.5,                &
     2020.5, 2030.5, 2040.5, 2050.5, 2060.5, 2070.5,                &
     2080.5, 2090.5, 2100.5, 2500.5/
  DATA AN2O_A1B/                                                    &
     278.0,  299.0,  301.0,  303.0,  305.0,  309.6,                 &
     311.0,  312.0,  312.4,  314.0,  316.0,  324.0,                 &
     331.0,  338.0,  344.0,  350.0,  356.0,  360.0,                 &
     365.0,  368.0,  372.0,  372.0/
  DATA AN2O_A2/                                                     &
     278.0,  299.0,  301.0,  303.0,  305.0,  309.6,                 &
     311.0,  312.0,  312.4,  314.0,  316.0,  325.0,                 &
     335.0,  347.0,  360.0,  373.0,  387.0,  401.0,                 &
     416.0,  432.0,  447.0,  447.0/
  DATA AN2O_B1/                                                     &
     278.0,  299.0,  301.0,  303.0,  305.0,  309.6,                 &
     311.0,  312.0,  312.4,  314.0,  316.0,  324.0,                 &
     333.0,  341.0,  349.0,  357.0,  363.0,  368.0,                 &
     371.0,  374.0,  375.0,  375.0/
  DATA TIMECH4/                                                     &
     1950.5, 1952.5, 1954.5, 1956.5, 1958.5, 1960.5,                &
     1962.5, 1964.5, 1966.5, 1968.5, 1970.5, 1972.5,                &
     1974.5, 1976.5, 1978.5, 1980.5, 1982.5, 1984.5,                &
     1986.5, 1988.5, 1990.5, 1992.5, 1994.5, 1996.5,                &
     1998.5, 2000.5, 2010.5, 2020.5, 2030.5, 2040.5,                &
     2050.5, 2060.5, 2070.5, 2080.5, 2090.5, 2100.5/
  DATA CH4_A1B/                                                     &
     1147.5, 1163.8, 1182.1, 1202.4, 1224.2, 1247.5,                &
     1272.2, 1298.4, 1326.0, 1355.3, 1386.0, 1417.5,                &
     1449.3, 1481.5, 1514.0, 1547.1, 1580.9, 1614.2,                &
     1644.9, 1671.6, 1694.3, 1714.0, 1721.0, 1728.0,                &
     1745.0, 1760.0, 1871.0, 2026.0, 2202.0, 2337.0,                &
     2400.0, 2386.0, 2301.0, 2191.0, 2078.0, 1974.0/
  DATA CH4_A2/                                                      &
     1147.5, 1163.8, 1182.1, 1202.4, 1224.2, 1247.5,                &
     1272.2, 1298.4, 1326.0, 1355.3, 1386.0, 1417.5,                &
     1449.3, 1481.5, 1514.0, 1547.1, 1580.9, 1614.2,                &
     1644.9, 1671.6, 1694.3, 1714.0, 1721.0, 1728.0,                &
     1745.0, 1760.0, 1861.0, 1997.0, 2163.0, 2357.0,                &
     2562.0, 2779.0, 3011.0, 3252.0, 3493.0, 3731.0/
  DATA CH4_B1/                                                      &
     1147.5, 1163.8, 1182.1, 1202.4, 1224.2, 1247.5,                &
     1272.2, 1298.4, 1326.0, 1355.3, 1386.0, 1417.5,                &
     1449.3, 1481.5, 1514.0, 1547.1, 1580.9, 1614.2,                &
     1644.9, 1671.6, 1694.3, 1714.0, 1721.0, 1728.0,                &
     1745.0, 1760.0, 1827.0, 1891.0, 1927.0, 1919.0,                &
     1881.0, 1836.0, 1797.0, 1741.0, 1663.0, 1574.0/
!
    TIME = ITIME(1) + (IMON(ITIME(2)) + ITIME(3))/365.0
    IF(ISC.EQ.4) TIME = ITIME(1) + 0.5
    IF(TIME.LT.1950.5) TIME = 1950.5
    IF(TIME.GT.2100.5) TIME = 2100.5
    SECS = 365.25*86400.0
    DO IT = 2,22
      IF(TIMEN2O(IT).GT.TIME) THEN
        X1 = TIMEN2O(IT) - TIME
        X2 = TIME - TIMEN2O(IT-1)
        Z1 = X1 + X2
        IF(ISC.EQ.1) THEN
          C1 = AN2O_A1B(IT-1)
          C2 = AN2O_A1B(IT)
        ELSEIF(ISC.EQ.2) THEN
          C1 = AN2O_A2(IT-1)
          C2 = AN2O_A2(IT)
        ELSEIF(ISC.EQ.3) THEN
          C1 = AN2O_B1(IT-1)
          C2 = AN2O_B1(IT)
        ELSEIF(ISC.EQ.4) THEN
          C1 = AN2O_A1B(IT-1)
          C2 = AN2O_A1B(IT)
        ENDIF
        AN2O = (X1*C1 + X2*C2)/Z1
        AN2OTR = (C2 - C1)/(AN2O*SECS*Z1)
        GOTO 10
      ENDIF
    ENDDO
 10 DO IT = 2,36
      IF(TIMECH4(IT).GT.TIME) THEN
        X1 = TIMECH4(IT) - TIME
        X2 = TIME - TIMECH4(IT-1)
        Z1 = X1 + X2
        IF(ISC.EQ.1) THEN
          C1 = CH4_A1B(IT-1)
          C2 = CH4_A1B(IT)
        ELSEIF(ISC.EQ.2) THEN
          C1 = CH4_A2(IT-1)
          C2 = CH4_A2(IT)
        ELSEIF(ISC.EQ.3) THEN
          C1 = CH4_B1(IT-1)
          C2 = CH4_B1(IT)
        ELSEIF(ISC.EQ.4) THEN
          C1 = CH4_A1B(IT-1)
          C2 = CH4_A1B(IT)
        ENDIF
        CH4 = (X1*C1 + X2*C2)/Z1
        CH4TR = (C2 - C1)/(CH4*SECS*Z1)
        GOTO 20
      ENDIF
    ENDDO
 20 AN2O = AN2O*1.0E-9
    CH4 = CH4*1.0E-9
    RETURN
END SUBROUTINE GHGS
!/ ----------------------------------------------------------------            
END MODULE STRAT_CHEM_MOD
