ecearth3/sources/tm5mp/proj/cb05/chem_rates.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!------------------------------------------------------------------------------
!                    TM5                                                      !
!------------------------------------------------------------------------------
!BOP
!
! !MODULE:      CHEM_RATES
!
! !DESCRIPTION: routines that calculate rates. Also used if there is no
!                chemistry (without_chemistry) to calculate Henry coefficients
!                for wet deposition. 
!\
!\
! !INTERFACE: 
!
MODULE CHEM_RATES
  !
  ! !USES:
  !
  use GO, only : gol, goPr, goErr

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: rates       ! rate constants & Henry law
#ifndef without_chemistry
  public :: calrates
#endif
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter  ::  mname = 'chem_rates'
  !
  ! !REVISION HISTORY: 
  !   26 Mar 2010 - P. Le Sager - took off NO+XO2 three rates
  !   15 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------------

contains

  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RATES
  !
  ! !DESCRIPTION:
  !     
  !      calculation of look up tables rate constants/henry coefficients
  !
  !      method
  !      ------
  !      use known array of temperatures (integers) to calculate rate constants
  !      3 body reactions are explicitly calculated in chemistry      
  !\
  !\
  ! !INTERFACE:
  !
  SUBROUTINE RATES( status )
    !
    ! !USES:
    !
    use toolbox, only : zfarr
    use chem_param
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    :: status
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2010 - P. Le Sager - added protex
    !   16 Nov 2010 - P. Le Sager - bug fix : kh1 & kh2 are real now.
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    integer :: itemp, k,  i
    real    :: kh1,   kh2, br
    real    :: ztrec, zt3rec, temp

    ! --- const ------------------------------
    character(len=*), parameter ::  rname = mname//'/rates'

    ! start
    ! JEW : updated rates JPL(2006), incl. Evaluation Number 15 (March, 2006)
    !    
    do k=1,ntemp
       itemp=k+ntlow
       temp=float(itemp)
       ztrec=1./float(itemp)
       zt3rec=300./float(itemp)     
       !JEW: changed to IUPAC 2013           
       rates_lut(knoo3,k)=zfarr(2.07e-12,-1400.,ztrec)     
       rates_lut(kho2no,k)=zfarr(3.45e-12,270.,ztrec)
       !JEW: changed to IUPAC 2012   
       rates_lut(kno2oha,k)=3.2e-30*(temp/300.)**(-4.5)
       rates_lut(kno2ohb,k)=3.0e-11
       ! JEW : reference IUPAC
       rates_lut(kohhno3a,k)=zfarr(2.41e-14,460.,ztrec)
       rates_lut(kohhno3b,k)=zfarr(6.51e-34,1335.,ztrec)
       rates_lut(kohhno3c,k)=zfarr(2.69e-17,2199.,ztrec) 
       rates_lut(kno2o3,k)=zfarr(1.4e-13,-2470.,ztrec)
       rates_lut(knono3,k)=zfarr(1.8e-11,110.,ztrec)
       !JEW: changed to IUPAC 2012      
       rates_lut(kno2no3a,k)=3.6e-30*(temp/300.)**(-4.1)
       rates_lut(kno2no3b,k)=1.9e-12*(temp/300.)**0.2     
       rates_lut(kn2o5a,k)=1.3e-3*(temp/300.)**(-3.5)
       rates_lut(kn2o5b,k)=9.7e14*(temp/300.)**0.1 
       rates_lut(kohhono,k)=zfarr(2.5e-12,260.,ztrec)      
       rates_lut(khno4oh,k)=zfarr(3.2e-12,690.,ztrec)       
       rates_lut(kno2ho2a,k)=1.4e-31*(temp/300.)**(-3.1)
       rates_lut(kno2ho2b,k)=4.0e-12 ! independant of T
       rates_lut(khno4a,k)=zfarr(4.1e-5,-10650.,ztrec)
       rates_lut(khno4b,k)=zfarr(6.0e15,-11170.,ztrec)
       rates_lut(khonoa,k)=7.4e-31*(temp/300.)**(-2.4)
       rates_lut(khonob,k)=3.3e-11*(temp/300.)**(-0.3)
       !JEW: changed to JPL2006       
       rates_lut(kodm,k)=.2095*zfarr(3.3e-11,55.,ztrec)+ &
            .7808*zfarr(2.15e-11,110.,ztrec)
       !JEW: changed to JPL2006	    
       rates_lut(kh2ood,k)=zfarr(1.63e-10,60.,ztrec)
       rates_lut(ko3ho2,k)=zfarr(1.0e-14,-490.,ztrec)
       rates_lut(ko3oh,k)=zfarr(1.7e-12,-940.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(khpoh,k)=1.8e-12      
       rates_lut(kfrmoh,k)=zfarr(5.5e-12,125.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(kch4oh,k)=zfarr(2.45e-12,-1775.,ztrec)
       ! TvN: JPL2015 gives slightly different exponents:
       ! KCOOHA: 1.4 -> 1.0
       ! KCOOHC: -0.6 -> 0.0
       rates_lut(kcooha,k)=5.9e-33*zt3rec**1.4
       rates_lut(kcoohb,k)=1.1e-12*zt3rec**(-1.3)
       rates_lut(kcoohc,k)=1.5e-13*zt3rec**(-0.6)
       rates_lut(kcoohd,k)=2.1e9*zt3rec**(-6.1)         
       !JEW: changed to JPL2006       
       rates_lut(kohmper,k)=zfarr(3.8e-12,200.,ztrec)
       rates_lut(kohrooh,k)=zfarr(3.01e-12,190.,ztrec) ! CB05
       rates_lut(kho2oh,k)=zfarr(4.8e-11,250.,ztrec)
       rates_lut(kh2oh,k)=zfarr(2.8e-12,-1800.,ztrec)       
       !JEW: changed to IUPAC 2009 
       br=1./(1.+498.*exp(-1160./temp))    
       rates_lut(kmo2ho2a,k)=zfarr(3.8e-13,750.,ztrec)
       rates_lut(kmo2ho2a,k)=rates_lut(kmo2ho2a,k)*(1.0-br)
       rates_lut(kmo2ho2b,k)=zfarr(3.8e-13,750.,ztrec)
       rates_lut(kmo2ho2b,k)=rates_lut(kmo2ho2b,k)*br
       
       rates_lut(kmo2no,k)=zfarr(2.3e-12,360.,ztrec)      
       rates_lut(kmo2mo2,k)=zfarr(9.5e-14,390.,ztrec)     
       !JEW: changed to JPL2006       
       ! TvN: It looks like the first two reactions below
       !      have updated coefficients in JPL2015:
       rates_lut(kho2ho2a,k)=zfarr(3.5e-13,430.,ztrec)
       rates_lut(kho2ho2b,k)=zfarr(1.7e-33,1000.,ztrec)
       rates_lut(kho2ho2c,k)=zfarr(1.4e-21,2200.,ztrec)
       rates_lut(kc41,k)=5.8e-16
       rates_lut(kc46,k)=zfarr(8.1e-12,270.,ztrec)
       !
       ! Reformulated as defined in IUPAC recommendations
       !      
       rates_lut(kc47a,k)=3.28e-28*(temp/300.)**(-6.87)
       rates_lut(kc47b,k)=1.125e-11*(temp/300.)**(-1.105)
       rates_lut(kc48a,k)=zfarr(1.1e-5,-10100.,ztrec)
       rates_lut(kc48b,k)=zfarr(1.9e17,-14100.,ztrec)             
       !JEW: changed to JPL2006      
       rates_lut(kc49,k)=zfarr(2.9e-12,500.,ztrec)
       rates_lut(kc50,k)=zfarr(4.3e-13,1040.,ztrec)
       !------------------------------------------------------  
       rates_lut(kc52,k)=8.1e-13
       rates_lut(kc53,k)=zfarr(1.e15,-8000.,ztrec)
       rates_lut(kc54,k)=1.6e3
       !JEW: changed to JPL2006
       rates_lut(kc57,k)=zfarr(5.4e-12,-610.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(kc58,k)=zfarr(8.5e-16,1520.,ztrec)
       rates_lut(kc59,k)=zfarr(4.6e-14,400.,ztrec)
       !JEW: changed to JPL2006       
       ! TvN: JPL2015 recommends different coefficients
       !      for OH+C2H4 (KC61):
       ! (1.e-28,-4.5)   -> (1.1e-28,-3.5)
       ! (8.8e-12,-0.85) -> (8.4e-12,-1.75)
       rates_lut(kc61a,k)=1.e-28*(temp/300.)**(-4.5)
       rates_lut(kc61b,k)=8.8e-12*(temp/300.)**(-0.85)
       !JEW: changed to JPL2006       
       rates_lut(kc62,k)=zfarr(1.2e-14,-2630.,ztrec)
       !JEW: changed to IUPAC2004       
       rates_lut(kc73,k)=1.5e-11 ! IUPAC
       rates_lut(kc76,k)=zfarr(2.7e-11,390.,ztrec) ! IUPAC
       rates_lut(kc77,k)=zfarr(1.03e-14,-1995.,ztrec) ! IUPAC
       rates_lut(kc78,k)=zfarr(3.15e-12,-450.,ztrec) ! IUPAC
       !JEW: changed to JPL2006       
       rates_lut(kc79,k)=zfarr(2.6e-12,365.,ztrec)
       rates_lut(kc80,k)=zfarr(1.6e-12,-2200.,ztrec)
       rates_lut(kc81,k)=zfarr(2.6e-12,365.,ztrec) ! bug        
       rates_lut(kc82,k)=zfarr(7.5e-13,700.,ztrec) ! CB05
       rates_lut(kc83,k)=8.e-11
       rates_lut(kc84,k)=zfarr(5.9e-13,-360.,ztrec) ! CB05 temp dep
       rates_lut(kc85,k)=zfarr(7.5e-13,700.,ztrec) ! CB05
       rates_lut(KO3PO2,k)=6.0e-34*(temp/300.)**(-2.4)
       rates_lut(KO3PO3,k)=zfarr(8.0e-12,-2060.,ztrec)
       rates_lut(KXO2XO2N,k)=6.8e-14
       rates_lut(KXO2N,k)=6.8e-14
       ! sulfur and ammonia gas phase reactions
       ! the hynes et al. parameterisation is replaced by chin et al. 1996

!JEW: changed to JPL2006
       rates_lut(kdmsoha,k)= 1.11e-11*exp(-240./temp)
       rates_lut(kdmsohb,k)=1.0e-39*exp(5820./temp)
       rates_lut(kdmsohc,k)=5.0e-30*exp(6280./temp)
       rates_lut(kdmsno3,k)=zfarr(1.9e-13,520.,ztrec)!at 298 1.01e-12
!JEW: changed to JPL2006       
       rates_lut(kso2oha,k)=3.3e-31*(temp/300.)**(-4.3)
       rates_lut(kso2ohb,k)= 1.6e-12*(temp/300.)
       !
       ! fate of ammonia
       !
       rates_lut(knh3oh,k)=  zfarr(1.7e-12,-710.,ztrec)!1.56e-13 at 298k
       rates_lut(knh2oh,k)= 3.4e-11
       rates_lut(knh2no,k)=  zfarr(4.0e-12,450.,ztrec)!1.72e-11
       rates_lut(knh2no2,k)= zfarr(2.1e-12,650.,ztrec)!1.86e-11
       rates_lut(knh2ho2,k)= 3.4e-11
       ! TvN: JPL2015 only says value is smaller than 6.0e-21
       rates_lut(knh2o2,k)= 6.0e-21
       rates_lut(knh2o3,k)= zfarr(4.3e-12,-930.,ztrec)!1.89e-13 at 298k
       !
       ! duplicate rates for NH2O2
       !
       rates_lut(knh2o2o3,k)= rates_lut(knh2o3,k)
       rates_lut(knh2o2no,k)= rates_lut(knh2no,k)
       !
       ! for higher organics
       rates_lut(kohmcho,k) = zfarr(4.4e-12,365.,ztrec) ! IUPAC
       rates_lut(kohmch2cho,k) = zfarr(4.9e-12,405.,ztrec)
       
       rates_lut(kno3mcho,k) = zfarr(1.4e-12,-1860.,ztrec) 
       rates_lut(kno3mch2cho,k) = 6.4e-15
       
       rates_lut(kohole,k) = zfarr(8.2e-12,610.,ztrec)   ! IUPAC
       
       rates_lut(ko3ole,k) = 1.0e-17
      
       ! Corrected to expression in Williams et al. (ACP,2013) 
       !rates_lut(kno3ole,k) = zfarr(4.0e-14,-400.,ztrec)
       rates_lut(kno3ole,k) = zfarr(4.6e-14,400.,ztrec)
!
! the rates for additional BVOC's
!

       ! TvN: JPL2015 only gives two digits: 2.9e-12
       rates_lut(kohch3oh,k) = zfarr(2.85e-12,-345.,ztrec)
       rates_lut(kohhcooh,k) = 4.0e-13
       rates_lut(kohmcooh,k) = zfarr(4.2e-14,855.,ztrec)
       ! TvN: JPL2015 recommends different coefficients:
       ! (6.9e-12,-1000.) -> (7.66e-12,-1020.)
       ! (3.0e-12,20.)    -> (3.35e-12,0.)
       ! (7.6e-12,-585.)  -> (8.7e-12,-615.)
       rates_lut(kohc2h6,k) = zfarr(6.9e-12,-1000.,ztrec)
       rates_lut(kohethoh,k) = zfarr(3.0e-12,20.,ztrec)
       rates_lut(kohc3h8,k) = zfarr(7.6e-12,-585.,ztrec)
       ! TvN: JPL2015 recommends different coefficients
       ! (8.0e-27,-3.5)   -> (4.6e-27,-4.)
       ! (3.0e-11,-1.0)   -> (2.6e-11,-1.3)
       rates_lut(kohc3h6a,k) = 8.0e-27*(temp/300.)**(-3.5)
       rates_lut(kohc3h6b,k) = 3.0e-11*(temp/300.) **(-1.0)
       rates_lut(ko3c3h6,k) = zfarr(5.5e-15,-1880.,ztrec) ! IUPAC
       rates_lut(kno3c3h6,k) = zfarr(4.6e-13,-1155.,ztrec) ! IUPAC
!
! Terpenes
! 
       rates_lut(kohterp,k) =  zfarr(1.2e-11,440.,ztrec) ! IUPAC
       rates_lut(ko3terp,k) =  zfarr(6.3e-16,-580.,ztrec) ! IUPAC
       rates_lut(kno3terp,k) = zfarr(1.2e-12,490.,ztrec) ! IUPAC
!
! Acetone
! 
       rates_lut(kohaceta,k) = zfarr(8.8e-12,-1320.,ztrec) ! IUPAC
       rates_lut(kohacetb,k) = zfarr(1.7e-14,423.,ztrec) ! IUPAC
       
       rates_lut(kaco2ho2,k) = 1.0e-11
       rates_lut(kaco2mo2,k) = 3.8e-12 ! IUPAC
       rates_lut(kaco2no,k) =  8.0e-12
       
       rates_lut(kaco2xo2,k) = 6.8e-14
       rates_lut(kxo2xo2n,k) = 6.8e-14
       rates_lut(kxo2n,k) = 6.8e-14    
!
!  MACR, MVK
!                     
       rates_lut(kohmvk,k)=zfarr(2.6e-12,610.,ztrec)
       rates_lut(kohmacr,k)=zfarr(8.0e-12,380.,ztrec)	  
       rates_lut(ko3mvk,k)=zfarr(8.5e-16,-1520.,ztrec)
       rates_lut(ko3macr,k)=zfarr(1.4e-15,-2100.,ztrec)        
!
! Additional peroxy from C3 VOC
!	                                 
       rates_lut(kc3h7o2no,k)=zfarr(7.6e-12,-585.,ztrec)
       rates_lut(kc3h7o2ho2,k)=zfarr(7.5e-13,700.,ztrec)	  
       rates_lut(khypo2no,k)=zfarr(4.2e-12,180.,ztrec)
       rates_lut(khypo2ho2,k)=zfarr(7.5e-13,700.,ztrec)
!
!   methyl peroxy nitrate
!       
       rates_lut(kmo2no2a,k)=2.5e-30*(temp/300.)**(-5.5)       	    
       rates_lut(kmo2no2b,k)=1.8e-11

       rates_lut(kch3o2no2a,k)=9.0E-5*exp(-9690./temp)       	    
       rates_lut(kch3o2no2b,k)=1.1e16*exp(-10560./temp)             
       !
       ! **** solubility Henry equilibrium
       !    HNO3/so4/nh4 just a very high number to take H and 
       !      dissociation into account
       !
       henry(:,k)=0.
       henry(ih2o2,k)=1.0e5*exp(6300*ztrec)*6.656e-10
       henry(ihno3,k)=1e7
       henry(ich3o2h,k)=310.*exp(5200*ztrec)*2.664e-8
       henry(ich2o,k)=3000.*exp(7200*ztrec)*3.253e-11
       henry(irooh,k)=340.*exp(6000.*ztrec)*1.821e-9
       henry(ich2o,k)=henry(ich2o,k)*37
       henry(ich3o2no2,k)=2.0*exp(4700.*ztrec)*exp(-4700.*(1/298.15))
       henry(iorgntr,k)=zfarr(1.8e-6,6000.,ztrec)
       henry(iso4,k)=1.e7
       henry(inh4,k)=1.e7
       henry(imsa,k)=1.e7
       henry(iso2,k) =1.2*exp(3120.*ztrec)*2.85e-5
       henry(inh3,k) =75.0*exp(3400.*ztrec)*1.10e-5
       henry(io3,k)=1.1e-2*exp(2300.*ztrec)*4.45e-4

       ! JEW add two new scavenging rates for CH3COCHO and ALD2
       ! need KH(eff) due to hydration steps for both species
       henry(imgly,k) = 3.2e4*48.6
       kh1=17.*exp(5000.*ztrec)*exp(-5000.*(1/298.15))
       kh2=13.*exp(5700.*ztrec)*exp(-5700.*(1/298.15))
       henry(iald2,k) = ((kh1+kh2)/2.)*1.0246
       henry(ihcooh,k) = 8900.*exp(6100.*ztrec)*exp(-6100.*(1/298.15))
       henry(ich3oh,k) = 220.0*exp(5200.*ztrec)*2.66e-8
       henry(imcooh,k)= 4.1e3*exp(6300.*ztrec)*exp(-6300.*(1/298.15))
       henry(iethoh,k) = 190.0*exp(6600.*ztrec)*exp(-6600.*(1/298.15))
       henry(iacet,k) = 35.*exp(3800.*ztrec)*exp(-3800.*(1/298.15))
    end do !k temperature loop
    !
    ! marked tracers:
    !
    henry(io3s,:)      = henry(io3,:)

    !ok
    status = 0

  end subroutine rates
  !EOC

#ifndef without_chemistry
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  calrates
  !
  ! !DESCRIPTION:
  !          	calculate rate constants using lookup table rates_lut
  !		calculate third bodies 
  !		calculate heterogeneous removal on aerosols
  !
  ! External: CALCHET
  !
  !\
  !\
  ! !INTERFACE:
  !
  SUBROUTINE CALRATES( region, level, is, js, rjx, rr, reff_cloud, ye, dt_chem, sad_cld_out, sad_ice_out, sad_aer_out, status )
    !
    ! !USES:
    !
    use global_data,     only : region_dat, mass_dat
    use Dims,            only : CheckShape, idate, jm
    use chem_param
    use binas,           only : Avog, pi
    use tm5_distgrid,    only : dgrid, Get_DistGrid
    use reaction_data,   only : nreac
#ifdef with_m7
    use m7_data,         only : rw_mode, dens_mode
    use mo_aero_m7,      only : nmod, cmr2ras
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)    :: region, level, is, js
    real,    intent(in)    :: reff_cloud(is:,js:)
    real,    intent(in)    :: dt_chem
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out)   :: rr(is:,js:,:)  ! reaction rates
    integer, intent(out)   :: status
    real,    intent(out)   :: sad_cld_out(is:,js:), sad_ice_out(is:,js:), sad_aer_out(is:,js:)
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real, intent(inout)    :: rjx(is:,js:,:) ! photolysis rate
    real, intent(inout)    ::  ye(is:,js:,:) ! extra 2D fields
    !
    ! !REVISION HISTORY: 
    !   26 Mar 2010 - P. Le Sager - added protex
    !   16 Jun 2011 - P. Le Sager - bug fix from JEW : g_n2o5 values
    !   15 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !   12 Jun 2012 - P. Le Sager - fix O3 concentration
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/calc_rates' 
    
    real, dimension(:,:,:,:), pointer :: rm   ! tracer mass

    ! help variables
    integer                          :: itemp,i,j, i1, j1, i2, j2, imode, iaer,n,latindex
    real                             :: tr, temp, wv, airp, rx1, rx2, rx3
    real                             :: dum, h2ox, aird, air_vol, o2, o3
    real                             :: x1, x2, xice, xliq
    !
    ! cloud chemistry of n2o5
    real                             :: dg, kt_liq_n2o5, kt_ice_n2o5, kt_liq_ho2, kt_ice_ho2
    real                             :: g_n2o5_l_temp
    real                             :: r_ice, r_cloud       ! cm
    real                             :: sad_ice, sad_cloud, iwc, lwc, o3_p
    ! Aerosol heterogeneous chemistry
    real                             :: kt_aer_n2o5, kt_aer_no3, kt_aer_ho2
    real                             :: sad_aer, sad_aert, aer_conc, aer_dens, aer_rad
    real                             :: smr_aer
    real, parameter                  :: sad_factor = 4.*pi*1.e-3 ! 4*pi as prefactor for area, 1e-3 to convert air_vol to cm3
    ! Standard aerosol properties density and radius for NO3_A, MSA,NH4, and -if not with_m7- SO4
    real, parameter                  :: aer_dens_ref = 1.7      ! assumed particle density  [gr/cm3]
    real, parameter                  :: aer_rad_ref_so4_LRH  = 0.18e-4  ! assumed particle radius [cm]
    real, parameter                  :: aer_rad_ref_no3_LRH  = 0.2e-4   ! assumed particle radius [cm]
    real, parameter                  :: aer_rad_ref_so4_HRH  = 0.25e-4  ! assumed particle radius accounting for abs. H2O[cm]
    real, parameter                  :: aer_rad_ref_no3_HRH  = 0.27e-4  ! assumed particle radius accounting for abs. H2O[cm]    
    
    ! Parameters to describe subgrid scale mixing
    real                             :: zcc,nn
    real, parameter                  :: M_HO2  =  (1+2*16.)/Avog *1e-3 ! HO2 mass, kg/molec
    real, parameter                  :: M_N2O5 =  (2*14+5*16.)/Avog *1e-3 ! N2O5 mass, kg/molec
    real, parameter                  :: M_NO3  =  (14+3*16.)/Avog *1e-3 ! no3 mass, kg/molec
    real	                     :: C_Thermal_N2O5, C_Thermal_HO2,C_Thermal_NO3, gamma, ttemp, factor, g_n2o5_aer_so4
    real                             :: aer_rad_ref, tot, no3frac, so4frac, r_no3, r_so4
    real                             :: zsgs_mix_n2o5
    real,dimension(180)              :: h2_surface
    
    ! --- begin --------------------------------

    CALL GET_DISTGRID( DGRID(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    call CheckShape( (/i2-i1+1, j2-j1+1, nj     /), shape(rjx), status )
    IF_ERROR_RETURN(status = 1)
    call CheckShape( (/i2-i1+1, j2-j1+1, nreac  /), shape(rr ), status )
    IF_ERROR_RETURN(status = 1)
    call CheckShape( (/i2-i1+1, j2-j1+1, n_extra/), shape(ye ), status )
    IF_ERROR_RETURN(status = 1)

    rm => mass_dat(region)%rm
    h2_surface(1:180)=h2_lat(1:180,idate(2))
    !$OMP parallel &
    !$OMP  default (none) &
    !$OMP  shared  (level, jsr, jer, isr, ier, region, ye, rates_lut, rr, rjx, rm, fscale) &
    !$OMP  private (i, j, jss, jee, rx3, temp, itemp, airp, h2ox, aird, &
    !$OMP           tr, wv, o2, o3, o3_p, rx1, rx2)

    rx3=0.6
    rr(i1:i2,j1:j2,1:nreac)=0.0     
    do j=j1,j2
       do i=i1,i2
          temp=ye(i,j,i_temp)
          itemp=nint(temp-float(ntlow))
          itemp=min(max(itemp,1),ntemp) !limit temperatures in loop-up table range
          airp=ye(i,j,i_pres)
          !
          ! Calculation of relative humidity [%]
          ! Richardson's approximation for water vapor pressure
          ! should be calculated in subroutine rates
          !
          h2ox = ye(i,j,ih2on)
          aird = ye(i,j,iairn)
          tr=1.-373.15/temp
          wv=exp((((-.1299*tr-.6445)*tr-1.976)*tr+13.3185)*tr)
          ye(i,j,irh)=h2ox*temp/(1013.25*wv*7.24291e16)
          ye(i,j,irh) = max(min(ye(i,j,irh),100.0),0.0)   !limit rh between 0-100%

          o2=0.209476*aird          
          o3 = rm(i,j,level,io3) / ye(i,j,iairm) * aird * fscale(io3)   !kg ----> #/cm3
          !
          !**** calculate three body reaction rates
	  nn=0.75-(1.27*LOG10(0.41))          !	  
          rx1=rates_lut(KNO2OHA,itemp)*aird
          rx2=rates_lut(KNO2OHB,itemp)
          rr(i,j,kno2oh)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.41)/(1+(log10(RX1/RX2)/nn)**2))
          rx1=rates_lut(KOHHNO3C,itemp)
          rx2=rates_lut(KOHHNO3B,itemp)*aird
          rr(i,j,kohhno3)=rates_lut(KOHHNO3A,itemp)+rx1*rx2/(rx1+rx2)
	  !
	  ! reformulated following IUPAC 2013
	  !
	  nn=0.75-(1.27*LOG10(0.35))
          rx1=rates_lut(KNO2NO3A,itemp)*aird
          rx2=rates_lut(KNO2NO3B,itemp)
          rr(i,j,kno2no3)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.35)/(1+(log10(RX1/RX2)/nn)**2))
          rx1=rates_lut(KN2O5A,itemp)*aird*exp(-11000/temp)
          rx2=rates_lut(KN2O5B,itemp)*exp(-11080/temp)	
	  rr(i,j,kn2o5)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.35)/(1+(log10(RX1/RX2)/nn)**2))	  
	  nn=0.75-(1.27*LOG10(0.4))
	  rx1=rates_lut(KNO2HO2A,itemp)*aird
          rx2=rates_lut(KNO2HO2B,itemp)
          rr(i,j,kno2ho2)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.4)/(1+(log10(RX1/RX2)/nn)**2))
	  rx1=rates_lut(KHNO4A,itemp)*aird
          rx2=rates_lut(KHNO4B,itemp)
	  !
	  !=((4.1e-5*exp(-10650/T))*M*(6.0e15*exp(-11170/T)))/((4.1e-5*exp(-10650/T))*M+(6.0e15*exp(-11170/T)))*10^(log10(0.4)/(1+(log10((4.1e-5*exp(-10650/T))*M/(6.0e15*exp(-11170/T)))/(0.75-1.27*log10(0.4)))^2))	  
          !
	  rr(i,j,khno4m)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.4)/(1+(log10(RX1/RX2)/nn)**2))
	  rx1=rates_lut(khonoa,itemp)*aird
	  rx2=rates_lut(khonob,itemp)
	  rr(i,j,khono)=rx1/(1.+rx1/rx2)*0.81**(1./(1.+log10(rx1/rx2)**2))	  
	  nn=0.75-(1.27*LOG10(0.36))
	  rx1=rates_lut(KMO2NO2A,itemp)*aird
          rx2=rates_lut(KMO2NO2B,itemp)
	  rr(i,j,kmo2no2)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.36)/(1+(log10(RX1/RX2)/nn)**2))
	  rx1=rates_lut(KCH3O2NO2A,itemp)*aird
          rx2=rates_lut(KCH3O2NO2B,itemp)
	  rr(i,j,kch3o2no2m)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.36)/(1+(log10(RX1/RX2)/nn)**2))	  	  
	  rx1=rates_lut(kcooha,itemp)*aird
          rx2=rates_lut(kcoohb,itemp)
          rr(i,j,kcooh)=rx1/(1.+rx1/rx2)*rx3**(1./(1.+log10(rx1/rx2)**2))
          !JEW: now add the second term for CO + OH
          rx1=rates_lut(kcoohc,itemp)
          rx2=rates_lut(kcoohd,itemp)/aird
          rr(i,j,kcooh)=rr(i,j,kcooh)+(rx1/(1.+rx1/rx2)*rx3**(1./(1.+log10(rx1/rx2)**2)))
          rx1=rates_lut(KC61A,itemp)*aird
          rx2=rates_lut(KC61B,itemp)
          rr(i,j,kc61)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**2))
	  rx1=rates_lut(kohc3h6a,itemp)*aird
          rx2=rates_lut(kohc3h6b,itemp)
          rr(i,j,kohc3h6)=rx1/(1.+rx1/rx2)*0.5**(1./(1.+log10(rx1/rx2)**2))
          !
          ! photolysis rates and 2 body reactions    
          !		  
          rr(i,j,knoo3)=rates_lut(KNOO3,itemp)  
          rr(i,j,kho2no)=rates_lut(KHO2NO,itemp)
          rr(i,j,kmo2no)=rates_lut(KMO2NO,itemp)
          rr(i,j,kno2o3)=rates_lut(KNO2O3,itemp)
          rr(i,j,knono3)=rates_lut(KNONO3,itemp)
          rr(i,j,khno4oh)=rates_lut(KHNO4OH,itemp)
	  rr(i,j,kohhono)=rates_lut(KOHHONO,itemp) 
          rr(i,j,kodm)=rates_lut(KODM,itemp)
          rr(i,j,kh2ood)=rates_lut(KH2OOD,itemp)
          rr(i,j,ko3ho2)=rates_lut(KO3HO2,itemp)
          !old formulation          rr(i,j,kcooh)=1.5E-13+9E-14*airp/101325. 
          rr(i,j,ko3oh)=rates_lut(KO3OH,itemp)
          rr(i,j,khpoh)=rates_lut(KHPOH,itemp)
          rr(i,j,kfrmoh)=rates_lut(KFRMOH,itemp)
          rr(i,j,kch4oh)=rates_lut(KCH4OH,itemp)
          latindex=j*180/jm(region)
          rr(i,j,kh2oh)=rates_lut(kh2oh,itemp)*h2_surface(latindex)*1.0e-9*aird !h2 latitudinal constraint JEW update
          rr(i,j,kohmper)=rates_lut(KOHMPER,itemp)
          rr(i,j,kohrooh)=rates_lut(KOHROOH,itemp)
          rr(i,j,kmo2ho2a)=rates_lut(KMO2HO2a,itemp)
	  rr(i,j,kmo2ho2b)=rates_lut(KMO2HO2b,itemp)
          rr(i,j,kmo2mo2)=rates_lut(KMO2MO2,itemp)
          rr(i,j,kho2oh)=rates_lut(KHO2OH,itemp)
          rr(i,j,kho2ho2)=(rates_lut(KHO2HO2A,itemp) + &
               rates_lut(KHO2HO2B,itemp)*aird)*(1.+rates_lut(KHO2HO2C,itemp)*h2ox)
          rr(i,j,kc41)=rates_lut(kc41,itemp)      
          ! JEW for ALD take the average of the C2 and C3 rate co-efficients; measurements suggest
          ! that CH3CHO only comprises 50% of [higher aldehydes] - grossman et al, JGR, 2003.	       
          rr(i,j,kc43)=(rates_lut(kohmcho,itemp)+rates_lut(kohmch2cho,itemp))/2.
          rr(i,j,kc44)=(rates_lut(kno3mcho,itemp)+ rates_lut(kno3mch2cho,itemp))/2.
          rr(i,j,kc46)=rates_lut(KC46,itemp)
 	  !
	  !  PAN now reformulated (2014)
	  !
	  nn=0.75-(1.27*LOG10(0.3))
          rx1=rates_lut(KC47A,itemp)*aird
          rx2=rates_lut(KC47B,itemp)
          !  rr(i,j,kc47)=RX1/(1.+RX1/RX2)*0.3**(1./(1.+LOG10(RX1/RX2)**2))
	  !
	  rr(i,j,kc47)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.3)/(1+(log10(RX1/RX2)/nn)**2))
          rx1=rates_lut(KC48A,itemp)*aird
          rx2=rates_lut(KC48B,itemp)
          !  rr(i,j,kc48)=RX1/(1.+RX1/RX2)*0.3**(1./(1.+LOG10(RX1/RX2)**2)
	  !
	  rr(i,j,kc48)=(RX1*RX2)/(RX1+RX2)*10**(log10(0.3)/(1+(log10(RX1/RX2)/nn)**2))	 
          rr(i,j,kc49)=rates_lut(KC49,itemp)
          rr(i,j,kc50)=rates_lut(KC50,itemp)
          rr(i,j,kc52)=rates_lut(KC52,itemp)    
          rr(i,j,kc53)=rates_lut(KC53,itemp)
          rr(i,j,kc54)=rates_lut(KC54,itemp)
          ! JEW updated average for OLE oxidation reactions	  
          rr(i,j,kc57)=(rates_lut(kohmvk,itemp)+rates_lut(kohole,itemp)+rates_lut(kohmacr,itemp))/3.
          rr(i,j,kc58)=(rates_lut(ko3mvk,itemp)+rates_lut(ko3ole,itemp)+rates_lut(ko3macr,itemp))/3.
          rr(i,j,kc59)=(rates_lut(kno3ole,itemp)+6.0e-16+3.5e-15)/3.
          rr(i,j,kc62)=rates_lut(KC62,itemp)
          rr(i,j,kc73)=rates_lut(KC73,itemp)
          rr(i,j,kc76)=rates_lut(KC76,itemp)
          rr(i,j,kc77)=rates_lut(KC77,itemp)
          rr(i,j,kc78)=rates_lut(KC78,itemp)
          ! JEW use rates in cb05 for XO2 reactions	  
          rr(i,j,kc79)=rates_lut(KC79,itemp)
          rr(i,j,kc80)=rates_lut(KC80,itemp)
          rr(i,j,kc81)=rates_lut(KC81,itemp)
          rr(i,j,kc82)=rates_lut(KC82,itemp)
          rr(i,j,kc83)=rates_lut(KC83,itemp)
          rr(i,j,kc84)=rates_lut(kc84,itemp)
          rr(i,j,kc85)=rates_lut(kc85,itemp)
#ifndef without_photolysis

          rjx(i,j,jo3d)=rjx(i,j,jo3d)*rr(i,j,kh2ood)*h2ox / &
               ( rr(i,j,kodm)*aird + rr(i,j,kh2ood)*h2ox ) 
          ! this is fraction rjo3d=>OH 
#endif
          RX1=rates_lut(kso2oha,itemp)*aird
          RX2=rates_lut(kso2ohb,itemp)
          rr(i,j,kso2oh)=RX1/(1.+RX1/RX2)*0.6**(1./(1.+LOG10(RX1/RX2)**2))
          !
          ! dmsoha => so2
          ! dmsohb => 0.75 SO2 + 0.25 MSA
          !
          rr(i,j,kdmsoha)=rates_lut(kdmsoha,itemp)
          rr(i,j,kdmsohb)=rates_lut(kdmsohb,itemp)*o2/ &
               (1.+rates_lut(kdmsohc,itemp)*o2)
          rr(i,j,kdmsno3)=rates_lut(kdmsno3,itemp)	
          ! ammonia chemistry
          rr(i,j,knh3oh)=rates_lut(knh3oh,itemp)
          rr(i,j,knh2oh)=rates_lut(knh2oh,itemp)
          rr(i,j,knh2no)=rates_lut(knh2no,itemp)
          rr(i,j,knh2no2)=rates_lut(knh2no2,itemp)
          rr(i,j,knh2ho2)=rates_lut(knh2ho2,itemp)
          rr(i,j,knh2o2)=rates_lut(knh2o2,itemp)*o2
          rr(i,j,knh2o3)=rates_lut(knh2o3,itemp)
	  rr(i,j,knh2o2no)=rates_lut(knh2no,itemp)
	  rr(i,j,knh2o2o3)=rates_lut(knh2o3,itemp)
          rr(i,j,knh2o2ho2)=rates_lut(knh2ho2,itemp)
          rr(i,j,krn222)=2.11e-6               ! s-1 decay time Rn222 to Pb210
          rr(i,j,ko3po2)   = rates_lut(ko3po2,itemp)*aird
          rr(i,j,ko3po3)   = rates_lut(ko3po3,itemp)
	  rr(i,j,kxo2xo2n) = 6.8e-14
	  rr(i,j,kxo2n)    = 6.8e-14
          !
          ! O3P in molecules cm3
          !
          ! JEW: Reformulated June 2012
          !
          ! Runaway O3 above 50 hPa due to missing stratopsheric chemistry
          ! therefore use pressure as an index
          !
          o3_p=0.

          if(airp/100 > 50. .and. airp/100. < 350.) o3_p=(2*rjx(i,j,jo2)*o2)
          !
          !  [O2] not used in the ebi
          !
          rjx(i,j,jo2)=(o3_p/(rr(i,j,ko3po2)*o2+rr(i,j,ko3po3)*o3)) * o2 * rr(i,j,ko3po2) 
          rr(i,j,ko3po3)=rr(i,j,ko3po3)*o3_p
	  !
          ! additional biogenic reactions     
          rr(i,j,kohhcooh)=4.0e-13
          rr(i,j,kohch3oh)=rates_lut(kohch3oh,itemp)
          
	  ! additional no3 peroxy reactions
          rr(i,j,kno3ho2)=4.0e-12
          rr(i,j,kno3mo2)=1.2e-12
          rr(i,j,kno3c2o3)=4.0e-12
          rr(i,j,kno3xo2)=2.5e-12 ! Zaveri and Peters
          
	  ! additional C2 compounds          
          rr(i,j,kohmcooh)=rates_lut(kohmcooh,itemp)
          rr(i,j,kohc2h6)=rates_lut(kohc2h6,itemp)
          rr(i,j,kohethoh)=rates_lut(kohethoh,itemp)
          rr(i,j,kohc3h8)=rates_lut(kohc3h8,itemp)
          rr(i,j,ko3c3h6)=rates_lut(ko3c3h6,itemp)
          rr(i,j,kno3c3h6)=rates_lut(kno3c3h6,itemp)
          !
          !  TERPENE reactions
          !              
          rr(i,j,kohterp) = rates_lut(kohterp,itemp)
          rr(i,j,ko3terp) = rates_lut(ko3terp,itemp)
          rr(i,j,kno3terp) = rates_lut(kno3terp,itemp)    	  
          !
          !  ISPD reactions
          !
          rr(i,j,kohispd)=(rates_lut(kohmvk,itemp)+rates_lut(kohmacr,itemp))/2
          rr(i,j,ko3ispd)=(rates_lut(ko3mvk,itemp)+rates_lut(ko3macr,itemp))/2
          rr(i,j,kno3ispd)=(6.0e-16+3.4e-15)/2
	  !
          ! Radon decay
          rr(i,j,krn222)=2.11e-6               ! s-1 decay time Rn222 to Pb210
	  !	  
	  ! acetone
	  rr(i,j,kohacet)=rates_lut(kohaceta,itemp)+rates_lut(kohacetb,itemp)
	  rr(i,j,kaco2ho2)=1.0e-11
	  rr(i,j,kaco2mo2)=3.8e-12
	  rr(i,j,kaco2no)=8.0e-12
	  rr(i,j,kaco2xo2)=6.8e-14 ! CB05  peroxy-loss
	  rr(i,j,kxo2xo2n)=6.8e-14
	  rr(i,j,kxo2n)=6.8e-14
          !
	  !
	  rr(i,j,kc3h7o2no)=rates_lut(kc3h7o2no,itemp)
	  rr(i,j,kc3h7o2ho2)=rates_lut(kc3h7o2ho2,itemp)
	  rr(i,j,khypo2no)=rates_lut(khypo2no,itemp)
	  rr(i,j,khypo2ho2)=rates_lut(khypo2ho2,itemp)
	  !
	  !                 	  
       end do
    end do   !_lat/lon loop
    !$OMP END PARALLEL


    !
    !  heterogeneous reaction of N2O5 and H2O -> 2 HNO3 on cloud and aerosol
    !  included in gas phase chemistry.
    !  VH, August-December 2013: Methodology updated, and extended with reactions for NO3 and HO2. 
    !  rw_mode and dens_mode contain the average radius [m] and density [kg/m3] for each modes,
    !  as computed in subroutine m7_averageproperties (pm6rp and prhop)
    !
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (jsr, jer, isr, ier, region, ye, rr, reff_cloud,level, &
!$OMP           g_ho2,g_n2o5,g_no3, sad_cld_out, sad_ice_out,sad_aer_out, &
!$OMP           mode_nm, mode_tracers, rm,dt_chem, m_n2o5, m_ho2, m_no3) &
    !$OMP  private (i, j, jss, jee, airp, dg, aird, xliq, xice, &
    !$OMP           temp, sad_ice,sad_cloud,iwc,lwc,r_ice, zcc,&
    !$OMP           r_cloud, sad_aer,sad_aert, aer_conc, air_vol, &
    !$OMP           g_n2o5_l_temp,kt_aer_n2o5, kt_aer_ho2,kt_aer_no3, &
    !$OMP           kt_liq_n2o5, kt_ice_n2o5, kt_liq_ho2, kt_ice_ho2  ) &
#ifdef with_m7
    !$OMP shared  (rw_mode, dens_mode) &
    !$OMP private (aer_rad, aer_dens) 
#endif
    
    do j = j1, j2
       do i=i1, i2
          sad_cloud = 0.
          sad_ice = 0.
          aird = ye(i,j,iairn) ! #/cm3
          
          ! typically the optical reff is somewhat larger than the physical r by 1-2um
          ! therefore downsize reff by 2.0uM for droplets 9-13uM and 1.0 for those between
          ! 6-9um
          ! http://www-das.uwyo.edu/~geerts/cwx/notes/chap08/moist_cloud.html
          ! set cloud droplet radius [cm]
          if(reff_cloud(i,j)>=9.0) then
             r_cloud = (reff_cloud(i,j)-2.0)/1e4
          elseif(reff_cloud(i,j)>=6.0) then
             r_cloud = (reff_cloud(i,j)-1.0)/1e4
          elseif(reff_cloud(i,j)>= 4.0) then  
             r_cloud = (reff_cloud(i,j)-0.25)/1e4
          else
             r_cloud = 4.0e-4 ! ensure a minimum - sometime reff_cloud is zero, i.e. not defined...
          endif

          ! SAD for ice particles
          ! JEW: The surface area density for ice particles in now linked to
          ! the IWC by the parameterization of Heymsfield and McFarquar (1996)
          ! and the effective radii from Fu (1996)
          ! VH the computation of sad_ice and r_ice is optimized, taking account of units.
          ! VH compute SAD representative for cloudy part only (scale with cloud cover)
          r_ice=5e-3              ! adopt as the default [cm]
          if(ye(i,j,iciwc)>1.0e-15 .and. ye(i,j,icc) > 1e-5) then
             ! convert iwc from units [kg/kg] to  [gr/m3]
             iwc=(ye(i,j,iciwc)*aird*xmair*1e6/avog) / ye(i,j,icc) 
             sad_ice=1.0e-4*iwc**0.9       ! [cm2/cm3]
             ! calculate the r_eff using the relationship of Fu (1996)
             ! r_ice=(1.73205/(3*0.917))*((iwc/1e6)/sad_ice) filling constants:
	     !
	     ! Bug fix by VH: corrected 3/(4*rho_ice)
	     !
             r_ice=0.8179*((iwc/1e6)/sad_ice)  ! [cm]
             ! the value adopted in von Kuhlmann and Lawrence is too low
             sad_ice=10*sad_ice          ! [cm2/cm3]
          endif
  
          if(ye(i,j,iclwc)>1.0e-15 .and. ye(i,j,icc)>1e-5 ) then
             ! lwc convert from units [kg/kg] to [kg/m3]
             lwc=( ye(i,j,iclwc)*aird*xmair*1e3/avog ) / ye(i,j,icc)
             lwc=lwc*1e-3 ! convert [kg/m3] to [gr/cm3]
             !VH no_cloud=lwc/sphere_vol
             !VH sad_cloud=4*pi*r_cloud**2*no_cloud
             sad_cloud=3.*lwc/(r_cloud)
          endif
          
          ! Assume that loss on cloud particles is not dominating the full loss budget which can
          ! remove the full trace gas concentration within one time step, within grid cel that is 
          ! partly covered with cloud 
          zcc = min(0.99,ye(i,j,icc))
          !
          ! We use the original formulation in dentener and crutzen
          ! of : kt=(r/Dg + 4/vgamma)-1 * A(cm2/cm3)
          airp=ye(i,j,i_pres)
          dg=0.1*1e5/airp  !simple approximation for diffusion coeff. [cm2/s]

          !n2o5... (see IUPAC)
          temp=ye(i,j,i_temp)
	  C_Thermal_N2O5 = SQRT (8*1.38E-23 *temp / (PI*M_N2O5 ) )  * 1e2 ! m/s -> cm/sec
          C_Thermal_HO2  = SQRT (8*1.38E-23 *temp / (PI*M_HO2 ) )  * 1e2 ! m/s -> cm/sec
          C_Thermal_NO3  = SQRT (8*1.38E-23 *temp / (PI*M_NO3 ) )  * 1e2 ! m/s -> cm/sec	  
	  
	  !=====================================================================================
	  !
	  ! Full parameterization of gamma N2O5 variability of SO4= aerosols from Evans and Jacob, 2005
	  ! Implemented by JEW July 2014
	  !
	  !  RH dependence from Kane et al., Heterogeneous uptake of gaseous N2O5 by (NH4)2SO4
	  !  NH4HSO4 and H2SO4 aerosols, J. Phys. Chem A, 2001, 105, 6465-6470
	  !
	  !=====================================================================================
	  
	  gamma=2.79e-4 + ye(i,j,irh)*(1.3e-4 + ye(i,j,irh)*(-3.43e-6 + ye(i,j,irh)*7.52e-8))
	  ttemp=MAX(temp, 282.)
	  factor=10.**(-2.0-4.0e-2*(ttemp-294.0))/0.01
	  g_n2o5_aer_so4=gamma*factor
	  
	  if(zcc > 0.0) then	  
	  	  
           g_n2o5_l_temp = 2.7e-5*exp(1800./temp)
           kt_liq_n2o5=1./((r_cloud/dg)+(4./(C_Thermal_N2O5*g_n2o5_l_temp)))
           kt_ice_n2o5=1./((r_ice/dg)  +(4./(C_Thermal_N2O5*g_n2o5_i)))
	  	  !
	  ! Scaling factor accounting for subgrid-scale mixing within time step
          ! Assume a lagrangian time scale ZDT_LAG of 3h (i.e.time scale for mixing)
          ! ZSGS_MIX=0 if CC=0, and 1 if CC=1; ZSGS_MIX is lower than CC for 0<CC<1 ,reducing the effective reaction rate
	  	  
	    zsgs_mix_n2o5=zcc*(1.-exp(-0.01/(1.-zcc))) 
	   
	  !
           rr(i,j,kn2o5l)=(kt_liq_n2o5*sad_cloud+kt_ice_n2o5*sad_ice)*zsgs_mix_n2o5    !liquid cloud & ice n2o5

          ! ho2...
           kt_liq_ho2=1./((r_cloud/dg)+(4./(C_Thermal_HO2*g_ho2_l)))
           kt_ice_ho2=1./((r_ice/dg)  +(4./(C_Thermal_HO2*g_ho2_i)))
	   
	  endif 
 
          ! HO2 uptake on clouds (set to 0 until a better description is available)
          !rr(i,j,kho2_l) =(kt_liq_ho2*sad_cloud +kt_ice_ho2*sad_ice)*zcc     !liquid cloud and ice ho2
	  rr(i,j,kho2_l) = 0.
          
          ! output
          sad_cld_out(i,j)=sad_cloud
          sad_ice_out(i,j)=sad_ice
          
          ! Aerosol uptake
          rr(i,j,kn2o5_aer)=0.
          rr(i,j,kno3_aer) =0.
          rr(i,j,kho2_aer) =0.

	  air_vol=ye(i,j,iairm)*Avog / (aird * xmair)     !air mass [kg] / air density in units 1e-3[kg/m3]  -> volume (units) 1e-3[cm3]
	  sad_aert=0.
#ifdef with_m7
          do imode=1,nmod
	     aer_conc = 0.
             ! Total concentration of aerosol within this mode. Changed to units: [kg]/1e-3[cm3] -> [gr/cm3]
             do iaer=1,mode_nm(imode)
                aer_conc=aer_conc+rm(i,j,level,mode_tracers(iaer,imode)) / air_vol   
             enddo
	     if (aer_conc .gt.1e-15) then 
                !>>> TvN
                ! Avoid using dens_mode, since it does not account for the presence of
                ! ammonium nitrate and associated water, and MSA- in the soluble accumulation mode.
                ! Since these components are taken into account in rw_mode,
                ! the use of dens_mode would introduce inconsistencies.
                !aer_dens=max(1.0,dens_mode(region,imode)%d3(i,j,level)) *1.e-3    ! particle density [kg/m3] -> [gr/cm3]
                aer_rad =max(1e-10,rw_mode  (region,imode)%d3(i,j,level)) *1e2    ! particle median radius [m]  -> [cm]
                ! To calculate the surface area of each mode, 
                ! first convert the number median radius to the surface mean radius:
                smr_aer = aer_rad * cmr2ras(imode)
                ! The suface area can then be calculated as the mean surface area per particle
                ! times the number of particles in the mode:
                !sad_aer=max(3.0 * aer_conc / ( aer_dens * aer_rad ),1.e-15)    ! units [cm2/cm3]
                sad_aer = rm(i,j,level,mode_start(imode)) * sad_factor * smr_aer**2 / air_vol
                sad_aert=sad_aert+sad_aer    ! total aerosol SAD for evaluation purposes
                ! In principle one should average the full expression for
                ! the rate constant over the lognormal size distributions.
                ! In practice, an effective value for the radius is used
                ! in the expression for kt (see Jacob, Atm. Env., 2000),
                ! in our case the number median radius.
                ! <<< TvN
                ! Ensure that g_xxx > 0. !!!
                kt_aer_n2o5=1./((aer_rad/dg)+(4./(C_Thermal_N2O5 * g_n2o5(imode))))
                kt_aer_no3 =1./((aer_rad/dg)+(4./(C_Thermal_NO3 * g_no3 (imode))))
                kt_aer_ho2 =1./((aer_rad/dg)+(4./(C_Thermal_HO2 * g_ho2 (imode))))
                ! Fill reaction rates
                rr(i,j,kn2o5_aer)=rr(i,j,kn2o5_aer)+kt_aer_n2o5 * sad_aer
                rr(i,j,kno3_aer) =rr(i,j,kno3_aer) +kt_aer_no3  * sad_aer
                rr(i,j,kho2_aer) =rr(i,j,kho2_aer) +kt_aer_ho2  * sad_aer
                !                 aer_rad_max(imode)=max(aer_rad_max(imode),aer_rad)
	     endif
          enddo
          !>>> TvN
          ! Commented the part below,
          ! because rw_mode already accounts for the presence of NO3_a, NH4, MSA, and aerosol water
          ! in the soluble accumulation mode (see chemistry.F90).
          !
          ! simple treatment for NO3_a, NH4 and MSA
          !imode = nmod+1
          !aer_conc=(rm(i,j,level,inh4)+rm(i,j,level,ino3_a)+rm(i,j,level,imsa)) / air_vol
          !if (aer_conc.gt.1e-15) then
             
             ! FIXME - need to CHECK  AER_RAD_REF - first iteration
          !   if(ye(i,j,irh) > 69.0) then
          !      aer_rad_ref = ( rm(i,j,level,ino3_a)*aer_rad_ref_no3_HRH + &
          !                      rm(i,j,level,imsa  )*aer_rad_ref_so4_HRH ) /&
          !           (rm(i,j,level,ino3_a) + rm(i,j,level,imsa))
          !   else
          !      aer_rad_ref = ( rm(i,j,level,ino3_a)*aer_rad_ref_no3_LRH + &
          !                      rm(i,j,level,imsa  )*aer_rad_ref_so4_LRH ) /&
          !           (rm(i,j,level,ino3_a) + rm(i,j,level,imsa))
          !   endif

          !   sad_aer=max(3.0 * aer_conc / ( aer_dens_ref * aer_rad_ref ),1.e-15) ! units [cm2/cm3]
          !   sad_aert=sad_aert+sad_aer    ! total aerosol SAD for evaluation purposes
             ! Ensure that g_xxx > 0. !!!
          !   kt_aer_n2o5=1./((aer_rad_ref/dg)+(4./(C_Thermal_N2O5 * g_n2o5(imode))))
          !   kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(C_Thermal_NO3 * g_no3 (imode))))
          !   kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(C_Thermal_HO2 * g_ho2 (imode))))
             ! Fill reaction rates
          !   rr(i,j,kn2o5_aer)=rr(i,j,kn2o5_aer)+kt_aer_n2o5 * sad_aer
          !   rr(i,j,kno3_aer) =rr(i,j,kno3_aer) +kt_aer_no3  * sad_aer
          !   rr(i,j,kho2_aer) =rr(i,j,kho2_aer) +kt_aer_ho2  * sad_aer
          !endif
          !<<< TvN
#else
          ! Uptake on on available aerosol: SO4, NO3, NH4, and MSA.
          ! Total concentration of sulfate aerosol  (include others) [gr/cm3]
          aer_conc=(rm(i,j,level,iso4)+rm(i,j,level,inh4)+rm(i,j,level,ino3_a)+rm(i,j,level,imsa)) / air_vol
          
          if (aer_conc.gt.1e-20) then
             
             no3frac=rm(i,j,level,ino3_a)/air_vol    !IGNORE NH4?
             so4frac=(rm(i,j,level,iso4)+rm(i,j,level,imsa))/air_vol
             !
             ! JEW (2014) : Introduce a crude fix for swelling of particles at high RH
             ! NH4+ exists in either SO4-- or NO3- state, so don't double count
             ! for aerosol volumes.
             !  
             if(ye(i,j,irh) < 70.0) then
                r_no3=(no3frac/aer_conc)*aer_rad_ref_no3_LRH
                r_so4=(so4frac/aer_conc)*aer_rad_ref_so4_LRH
                aer_rad_ref=r_no3+r_so4
             endif

             if(ye(i,j,irh) > 69.0) then
                r_no3=(no3frac/aer_conc)*aer_rad_ref_no3_HRH
                r_so4=(so4frac/aer_conc)*aer_rad_ref_so4_HRH
                aer_rad_ref=r_no3+r_so4
             endif

             sad_aer  = 3.0 * aer_conc / ( aer_dens_ref * aer_rad_ref )        ! units [cm2/cm3]
             sad_aert=sad_aer

             ! Ensure that g_xxx > 0. !!!
             kt_aer_n2o5=1./((r_no3/dg)+(4./(C_Thermal_N2O5 * g_n2o5_aer)))
             kt_aer_n2o5=kt_aer_n2o5+(1./((r_so4/dg)+(4./(C_Thermal_N2O5 * g_n2o5_aer_so4))))
             kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(C_Thermal_NO3 * g_no3_aer)))
             kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(C_Thermal_HO2 * g_ho2_aer)))

             ! Fill reaction rates
             rr(i,j,kn2o5_aer)=kt_aer_n2o5 * sad_aer
             rr(i,j,kno3_aer) =kt_aer_no3  * sad_aer
             rr(i,j,kho2_aer) =kt_aer_ho2  * sad_aer
             
          endif
#endif /* M7 */

          ! output
          sad_aer_out(i,j)=sad_aert
          
       end do
    end do
    !$OMP END PARALLEL

    nullify(rm)
    status = 0
    
  end subroutine calrates
  !EOC
#endif /*CHEM*/
  
END MODULE CHEM_RATES