ecearth3/sources/tm5mp/proj/cbm4/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
    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 JPL2006           
       rates_lut(knoo3,k)=zfarr(3.e-12,-1500.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(kho2no,k)=zfarr(3.5e-12,250.,ztrec)
       !JEW: changed to JPL2006     
       rates_lut(kno2oha,k)=1.8e-30*zt3rec**3.0
       rates_lut(kno2ohb,k)=2.8e-11
       rates_lut(kohhno3a,k)=zfarr(2.41e-14,460.,ztrec)    !!!wp!!! new ravi
       rates_lut(kohhno3b,k)=zfarr(6.51e-34,1335.,ztrec)   !!!wp!!! new ravi
       rates_lut(kohhno3c,k)=zfarr(2.69e-17,2199.,ztrec)   !!!wp!!! new ravi
       rates_lut(kno2o3,k)=zfarr(1.2e-13,-2450.,ztrec)
       rates_lut(knono3,k)=zfarr(1.5e-11,170.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(kno2no3a,k)=2.0e-30*zt3rec**4.4
       rates_lut(kno2no3b,k)=1.4e-12*zt3rec**0.7
       !JEW: changed to JPL2006      
       rates_lut(kn2o5,k)=zfarr(2.7e-27,11000.,ztrec)
       rates_lut(khno4oh,k)=zfarr(1.3e-12,380.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(kno2ho2a,k)=2.0e-31*zt3rec**3.4
       rates_lut(kno2ho2b,k)=2.9e-12*zt3rec**1.1
       rates_lut(khno4m,k)=zfarr(2.1e-27,10900.,ztrec)
       !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
       !JEW: changed to JPL2006       
       rates_lut(kfrmoh,k)=zfarr(5.5e-12,125.,ztrec)
       !JEW: changed to JPL2006       
       rates_lut(kch4oh,k)=zfarr(2.45e-12,-1775.,ztrec)
       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 JPL2006       
       rates_lut(kmo2ho2,k)=zfarr(4.1e-13,750.,ztrec)
       rates_lut(kmo2no,k)=zfarr(2.8e-12,300.,ztrec)      
       rates_lut(kmo2mo2,k)=zfarr(9.5e-14,390.,ztrec)     
       !JEW: changed to JPL2006       
       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)
       ! from IUPAC (Atkinson et al, 2006)       
       rates_lut(kc47a,k)=2.7e-28*zt3rec**7.1
       rates_lut(kc47b,k)=1.2e-11*zt3rec**0.9
       rates_lut(kc48a,k)=zfarr(4.9e-3,-12100.,ztrec)
       rates_lut(kc48b,k)=zfarr(5.4e16,-13830.,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       
       rates_lut(kc61a,k)=1.e-28*zt3rec**4.5
       rates_lut(kc61b,k)=8.8e-12*zt3rec**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.04e-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*zt3rec**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(-253./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,500.,ztrec)!at 298 1.01e-12
       !JEW: changed to JPL2006       
       rates_lut(kso2oha,k)=3.3e-31*zt3rec**4.3
       rates_lut(kso2ohb,k)= 1.6e-12
       !
       ! fate of ammonia
       !
       rates_lut(knh3oh,k)=  zfarr(1.7e-12,-710.,ztrec)!1.56e-13 at 298k
       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
       rates_lut(knh2o2,k)= 6.0e-21
       rates_lut(knh2o3,k)= zfarr(4.3e-12,-930.,ztrec)!1.89e-13 at 298k
       !
       ! 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.5e-15

       rates_lut(kohmvk,k) = zfarr(1.86e-11,175.,ztrec)  ! IUPAC
       rates_lut(kohole,k) = zfarr(8.2e-12,610.,ztrec)   ! IUPAC
       rates_lut(kohmacr,k) = zfarr(2.6e-12,610.,ztrec)  ! IUPAC 

       rates_lut(ko3mvk,k) = zfarr(8.5e-16,-1520.,ztrec) ! Poesch et al 2000
       rates_lut(ko3ole,k) = 1.0e-17
       rates_lut(ko3macr,k) = zfarr(1.4e-15,-2100.,ztrec)! IUPAC

       rates_lut(kno3ole,k) = zfarr(4.0e-14,-400.,ztrec)

       !
       ! **** 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(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

    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
    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
    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  = 0.18e-4  ! assumed particle radius [cm]
    
    ! Parameters to describe subgrid scale mixing
    real                             :: zcc

    
    ! --- 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

    !$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
          
!------ Fix O3 conc (PLS, June 12, 2012)
!PREV          
!	  o3=y(i,j,io3)
!NOW
          o3 = rm(i,j,level,io3) / ye(i,j,iairm) * aird * fscale(io3)   !kg ----> #/cm3
!------          
          !
          !**** calculate three body reaction rates
          !
          rx1=rates_lut(KNO2OHA,itemp)*aird
          rx2=rates_lut(KNO2OHB,itemp)
          rr(i,j,kno2oh)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**2))
          rx1=rates_lut(KOHHNO3C,itemp)
          rx2=rates_lut(KOHHNO3B,itemp)*aird
          rr(i,j,kohhno3)=rates_lut(KOHHNO3A,itemp)+rx1*rx2/(rx1+rx2)
          rx1=rates_lut(KNO2NO3A,itemp)*aird
          rx2=rates_lut(KNO2NO3B,itemp)
          rr(i,j,kno2no3)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**2))
          rx1=rates_lut(KNO2HO2A,itemp)*aird
          rx2=rates_lut(KNO2HO2B,itemp)
          rr(i,j,kno2ho2)=RX1/(1.+RX1/RX2)*RX3**(1./(1.+LOG10(RX1/RX2)**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))
          !
          ! 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,kn2o5)=rr(i,j,kno2no3)/rates_lut(KN2O5,itemp)
          rr(i,j,khno4oh)=rates_lut(KHNO4OH,itemp) 
          rr(i,j,khno4m)=rr(i,j,kno2ho2)/rates_lut(KHNO4M,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)
          rr(i,j,kh2oh)=rates_lut(kh2oh,itemp)*550.e-9*aird !h2=550ppbv JEW update
          rr(i,j,kohmper)=rates_lut(KOHMPER,itemp)
          rr(i,j,kohrooh)=rates_lut(KOHROOH,itemp)
          rr(i,j,kmo2ho2)=rates_lut(KMO2HO2,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)
          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))
          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,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,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,krn222)=2.11e-6               ! s-1 decay time Rn222 to Pb210
! - BEFORE -   25 Jun 2012 - P. Le Sager 
!	  rr(i,j,ko3po2)=rates_lut(ko3po2,itemp)*o2
!	  rr(i,j,ko3po3)=rates_lut(ko3po3,itemp)
!	  rr(i,j,kxo2xo2n)=6.8e-14
!	  rr(i,j,kxo2n)=6.8e-14
!	  rjx(i,j,jo2)=(2*rjx(i,j,jo2)*o2)*rr(i,j,ko3po2)/ &
!	               (rr(i,j,ko3po2) + rr(i,j,ko3po3)*o3)                 
!
! - NOW -
          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
! - - 
          
       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) &
    !$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)
          g_n2o5_l_temp = 2.7e-5*exp(1800./temp)
          kt_liq_n2o5=1./((r_cloud/dg)+(4./(2.4e4*g_n2o5_l_temp)))
          kt_ice_n2o5=1./((r_ice/dg)  +(4./(2.4e4*g_n2o5_i)))
          rr(i,j,kn2o5l)=(kt_liq_n2o5*sad_cloud+kt_ice_n2o5*sad_ice)*zcc    !liquid cloud & ice n2o5

          ! ho2...
          kt_liq_ho2=1./((r_cloud/dg)+(4./(2.4e4*g_ho2_l)))
          kt_ice_ho2=1./((r_ice/dg)  +(4./(2.4e4*g_ho2_i)))

          ! 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./(2.4e4 * g_n2o5(imode))))
                kt_aer_no3 =1./((aer_rad/dg)+(4./(2.4e4 * g_no3 (imode))))
                kt_aer_ho2 =1./((aer_rad/dg)+(4./(2.4e4 * 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
          !   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./(2.4e4 * g_n2o5(imode))))
          !   kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(2.4e4 * g_no3 (imode))))
          !   kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(2.4e4 * 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

             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./((aer_rad_ref/dg)+(4./(2.4e4 * g_n2o5_aer)))
             kt_aer_no3 =1./((aer_rad_ref/dg)+(4./(2.4e4 * g_no3_aer)))
             kt_aer_ho2 =1./((aer_rad_ref/dg)+(4./(2.4e4 * 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)

    ! ok
    status = 0
    
  end subroutine calrates
  !EOC
#endif /*CHEM*/
  
END MODULE CHEM_RATES