
module mo_sethet

!
! LKE (10/11/2010): added HCN, CH3CN, HCOOH  to cesm1_0_beta07_offline version
!                   HCN, CH3CN have new Henry's Law coefficients, HCOOH is set to CH3COOH
! LKE (10/18/2010): SO2 washout corrected based on recommendation of R.Easter, PNNL
!
  use shr_kind_mod,    only: r8 => shr_kind_r8
  use cam_logfile,     only: iulog
  use gas_wetdep_opts, only: gas_wetdep_cnt, gas_wetdep_method, gas_wetdep_list
  use phys_control,    only: phys_getopts

  private
  public :: sethet_inti, sethet

  save

  integer :: h2o2_ndx, hno3_ndx, ch2o_ndx, ch3ooh_ndx, ch3coooh_ndx, &
       ho2no2_ndx, ch3cocho_ndx, xooh_ndx, onitr_ndx, glyald_ndx, &
       ch3cho_ndx, mvk_ndx, macr_ndx, pooh_ndx, c2h5ooh_ndx, &
       c3h7ooh_ndx, rooh_ndx, isopno3_ndx, onit_ndx, Pb_ndx, &
       macrooh_ndx, isopooh_ndx, ch3oh_ndx, c2h5oh_ndx, hyac_ndx, hydrald_ndx
  integer :: spc_h2o2_ndx, spc_hno3_ndx
  integer :: spc_so2_ndx
  integer :: spc_sogm_ndx, spc_sogi_ndx, spc_sogt_ndx, spc_sogb_ndx, spc_sogx_ndx

  integer :: alkooh_ndx, mekooh_ndx, tolooh_ndx, terpooh_ndx, ch3cooh_ndx
  integer :: so2_ndx, soa_ndx, so4_ndx, cb2_ndx, oc2_ndx, nh3_ndx, nh4no3_ndx, &
             sa1_ndx, sa2_ndx, sa3_ndx, sa4_ndx, nh4_ndx, h2so4_ndx
  integer :: soag0_ndx, soag15_ndx, soag24_ndx, soag31_ndx
  integer :: soag32_ndx, soag33_ndx, soag34_ndx, soag35_ndx
  integer :: bsoag1_ndx,bsoag2_ndx,bsoag3_ndx,bsoag4_ndx,bsoag5_ndx
  integer :: bsoag6_ndx,bsoag7_ndx,bsoag8_ndx,bsoag9_ndx,bsoag10_ndx
  integer :: bsoag11_ndx,bsoag12_ndx,bsoag13_ndx,bsoag14_ndx,bsoag15_ndx
  integer :: bsoag16_ndx,bsoag17_ndx
  integer :: xisopno3_ndx,xho2no2_ndx,xonitr_ndx,xhno3_ndx,xonit_ndx
  integer :: clono2_ndx, brono2_ndx, hcl_ndx, n2o5_ndx, hocl_ndx, hobr_ndx, hbr_ndx 
  integer :: ch3cn_ndx, hcn_ndx, hcooh_ndx
  integer, allocatable :: wetdep_map(:)
  integer :: sogm_ndx, sogi_ndx, sogt_ndx, sogb_ndx, sogx_ndx
  logical :: do_wetdep

  ! prognostic modal aerosols
  logical :: prog_modal_aero

contains

  subroutine sethet_inti
    !-----------------------------------------------------------------------      
    !       ... intialize the wet removal rate constants routine
    !-----------------------------------------------------------------------      

    use mo_chem_utls, only : get_het_ndx, get_spc_ndx
    use spmd_utils,   only : masterproc
    use cam_abortutils,   only : endrun

    integer :: k, m
    
    do_wetdep = gas_wetdep_cnt>0 .and. gas_wetdep_method=='MOZ'
    if ( .not. do_wetdep) return

    call phys_getopts( prog_modal_aero_out = prog_modal_aero )

    allocate( wetdep_map(gas_wetdep_cnt))

    do k=1,gas_wetdep_cnt
       m = get_het_ndx( trim(gas_wetdep_list(k))) 
       if (m>0) then
          wetdep_map(k) = m
       else
          call endrun('sethet_inti: cannot map '//trim(gas_wetdep_list(k)))
       endif
    enddo

    xisopno3_ndx = get_het_ndx( 'XISOPNO3' )
    xho2no2_ndx  = get_het_ndx( 'XHO2NO2' )
    xonitr_ndx   = get_het_ndx( 'XONITR' )
    xhno3_ndx    = get_het_ndx( 'XHNO3' )
    xonit_ndx    = get_het_ndx( 'XONIT' )

    spc_h2o2_ndx = get_spc_ndx( 'H2O2' )
    spc_hno3_ndx = get_spc_ndx( 'HNO3' )
    spc_so2_ndx  = get_spc_ndx( 'SO2' )

    clono2_ndx = get_het_ndx( 'CLONO2' )
    brono2_ndx = get_het_ndx( 'BRONO2' )
    hcl_ndx    = get_het_ndx( 'HCL' )
    n2o5_ndx   = get_het_ndx( 'N2O5' )
    hocl_ndx   = get_het_ndx( 'HOCL' )
    hobr_ndx   = get_het_ndx( 'HOBR' )
    hbr_ndx    = get_het_ndx( 'HBR' )

    h2o2_ndx   = get_het_ndx( 'H2O2' )
    hno3_ndx   = get_het_ndx( 'HNO3' )
    ch2o_ndx   = get_het_ndx( 'CH2O' )
    ch3ooh_ndx = get_het_ndx( 'CH3OOH' )
    ch3coooh_ndx = get_het_ndx( 'CH3COOOH' )
    ho2no2_ndx  = get_het_ndx( 'HO2NO2' )
    ch3cocho_ndx = get_het_ndx( 'CH3COCHO' )
    xooh_ndx    = get_het_ndx( 'XOOH' )
    onitr_ndx   = get_het_ndx( 'ONITR' )
    glyald_ndx  = get_het_ndx( 'GLYALD' )
    ch3cho_ndx  = get_het_ndx( 'CH3CHO' )
    mvk_ndx     = get_het_ndx( 'MVK' )
    macr_ndx    = get_het_ndx( 'MACR' )
    pooh_ndx    = get_het_ndx( 'POOH' )
    c2h5ooh_ndx = get_het_ndx( 'C2H5OOH' )
    c3h7ooh_ndx = get_het_ndx( 'C3H7OOH' )
    rooh_ndx    = get_het_ndx( 'ROOH' )
    isopno3_ndx = get_het_ndx( 'ISOPNO3' )
    onit_ndx    = get_het_ndx( 'ONIT' )
    Pb_ndx      = get_het_ndx( 'Pb' )
    macrooh_ndx = get_het_ndx( 'MACROOH' )
    isopooh_ndx = get_het_ndx( 'ISOPOOH' )
    ch3oh_ndx   = get_het_ndx( 'CH3OH' )
    c2h5oh_ndx  = get_het_ndx( 'C2H5OH' )
    hyac_ndx    = get_het_ndx( 'HYAC' )
    hydrald_ndx = get_het_ndx( 'HYDRALD' )
    alkooh_ndx  = get_het_ndx( 'ALKOOH' )
    mekooh_ndx  = get_het_ndx( 'MEKOOH' )
    tolooh_ndx  = get_het_ndx( 'TOLOOH' )
    terpooh_ndx = get_het_ndx( 'TERPOOH' )
    ch3cooh_ndx = get_het_ndx( 'CH3COOH' )
    so2_ndx     = get_het_ndx( 'SO2' )
    soa_ndx     = get_het_ndx( 'SOA' )
    sogb_ndx    = get_het_ndx( 'SOGB' )
    sogi_ndx    = get_het_ndx( 'SOGI' )
    sogm_ndx    = get_het_ndx( 'SOGM' )
    sogt_ndx    = get_het_ndx( 'SOGT' )
    sogx_ndx    = get_het_ndx( 'SOGX' )
    so4_ndx     = get_het_ndx( 'SO4' )
    cb2_ndx     = get_het_ndx( 'CB2' )
    oc2_ndx     = get_het_ndx( 'OC2' )
    nh3_ndx     = get_het_ndx( 'NH3' )
    nh4no3_ndx  = get_het_ndx( 'NH4NO3' )
    nh4_ndx     = get_het_ndx( 'NH4' )
    h2so4_ndx   = get_het_ndx( 'H2SO4' )
    sa1_ndx     = get_het_ndx( 'SA1' )
    sa2_ndx     = get_het_ndx( 'SA2' )
    sa3_ndx     = get_het_ndx( 'SA3' )
    sa4_ndx     = get_het_ndx( 'SA4' )
    ch3cn_ndx   = get_het_ndx( 'CH3CN' )
    hcn_ndx     = get_het_ndx( 'HCN' )
    hcooh_ndx   = get_het_ndx( 'HCOOH' )
    soag0_ndx   = get_het_ndx( 'SOAG0' )
    soag15_ndx  = get_het_ndx( 'SOAG15' )
    soag24_ndx  = get_het_ndx( 'SOAG24' )
    soag31_ndx  = get_het_ndx( 'SOAG31' )
    soag32_ndx  = get_het_ndx( 'SOAG32' )
    soag33_ndx  = get_het_ndx( 'SOAG33' )
    soag34_ndx  = get_het_ndx( 'SOAG34' )
    soag35_ndx  = get_het_ndx( 'SOAG35' )
    bsoag1_ndx  = get_het_ndx( 'BSOAG1' )
    bsoag2_ndx  = get_het_ndx( 'BSOAG2' )
    bsoag3_ndx  = get_het_ndx( 'BSOAG3' )
    bsoag4_ndx  = get_het_ndx( 'BSOAG4' )
    bsoag5_ndx  = get_het_ndx( 'BSOAG5' )
    bsoag6_ndx  = get_het_ndx( 'BSOAG6' )
    bsoag7_ndx  = get_het_ndx( 'BSOAG7' )
    bsoag8_ndx  = get_het_ndx( 'BSOAG8' )
    bsoag9_ndx  = get_het_ndx( 'BSOAG9' )
    bsoag10_ndx = get_het_ndx( 'BSOAG10' )
    bsoag11_ndx = get_het_ndx( 'BSOAG11' )
    bsoag12_ndx = get_het_ndx( 'BSOAG12' )
    bsoag13_ndx = get_het_ndx( 'BSOAG13' )
    bsoag14_ndx = get_het_ndx( 'BSOAG14' )
    bsoag15_ndx = get_het_ndx( 'BSOAG15' )
    bsoag16_ndx = get_het_ndx( 'BSOAG16' )
    bsoag17_ndx = get_het_ndx( 'BSOAG17' )


    if (masterproc) then
       write(iulog,*) 'sethet_inti: new ndx ',so2_ndx,soa_ndx,so4_ndx,cb2_ndx,oc2_ndx, &
            soag0_ndx,soag15_ndx,soag24_ndx,soag31_ndx,soag32_ndx,soag33_ndx,soag34_ndx,soag35_ndx, &
            nh3_ndx,nh4no3_ndx,sa1_ndx,sa2_ndx,sa3_ndx,sa4_ndx, &
            bsoag1_ndx,bsoag2_ndx,bsoag3_ndx,bsoag4_ndx,bsoag5_ndx, &
            bsoag6_ndx,bsoag7_ndx,bsoag8_ndx,bsoag9_ndx,bsoag10_ndx, &
            bsoag11_ndx,bsoag12_ndx,bsoag13_ndx,bsoag14_ndx,bsoag15_ndx,bsoag16_ndx,bsoag17_ndx
       write(iulog,*) ' '
       write(iulog,*) 'sethet_inti: diagnotics '
       write(iulog,'(10i5)') h2o2_ndx, hno3_ndx, ch2o_ndx, ch3ooh_ndx, ch3coooh_ndx, &
            ho2no2_ndx, ch3cocho_ndx, xooh_ndx, onitr_ndx, glyald_ndx, &
            ch3cho_ndx, mvk_ndx, macr_ndx, pooh_ndx, c2h5ooh_ndx, &
            c3h7ooh_ndx, rooh_ndx, isopno3_ndx, onit_ndx, Pb_ndx, &
            macrooh_ndx, isopooh_ndx, ch3oh_ndx, c2h5oh_ndx, hyac_ndx, hydrald_ndx
    endif

  end subroutine sethet_inti

  subroutine sethet( het_rates, press, zmid,  phis, tfld, &
                     cmfdqr, nrain, nevapr, delt, xhnm, &
                     qin, ncol, lchnk )
    !-----------------------------------------------------------------------      
    !       ... compute rainout loss rates (1/s)
    !-----------------------------------------------------------------------      

    use physconst,    only : rga,pi
    use chem_mods,    only : gas_pcnst
    use ppgrid,       only : pver, pcols
    use phys_grid,    only : get_rlat_all_p
    use cam_abortutils,   only : endrun
    use mo_constants, only : avo => avogadro, boltz_cgs

    implicit none
    !-----------------------------------------------------------------------      
    !       ... dummy arguments
    !-----------------------------------------------------------------------      
    integer, intent(in)   ::    ncol                        ! columns in chunk
    integer, intent(in)   ::    lchnk                       ! chunk index
    real(r8), intent(in)  ::    delt                        ! time step ( s )
    real(r8), intent(in)  ::    press(pcols,pver)           ! pressure in pascals
    real(r8), intent(in)  ::    cmfdqr(ncol,pver)           ! dq/dt for convection
    real(r8), intent(in)  ::    nrain(ncol,pver)            ! stratoform precip
    real(r8), intent(in)  ::    nevapr(ncol,pver)           ! evaporation
    real(r8), intent(in)  ::    qin(ncol,pver,gas_pcnst)    ! xported species ( vmr )
    real(r8), intent(in)  ::    zmid(ncol,pver)             ! midpoint geopot (km)
    real(r8), intent(in)  ::    phis(pcols)                 ! surf geopot
    real(r8), intent(in)  ::    tfld(pcols,pver)            ! temperature (k)
    real(r8), intent(in)  ::    xhnm(ncol,pver)             ! total atms density ( /cm^3)
    real(r8), intent(out) ::    het_rates(ncol,pver,gas_pcnst) ! rainout loss rates

    !-----------------------------------------------------------------------      
    !       ... local variables
    !-----------------------------------------------------------------------      
    real(r8), parameter ::  xrm   = .189_r8             ! mean diameter of rain drop (cm)
    real(r8), parameter ::  xum   = 748._r8             ! mean rain drop terminal velocity (cm/s)
    real(r8), parameter ::  xvv   = 6.18e-2_r8          ! kinetic viscosity (cm^2/s)
    real(r8), parameter ::  xdg   = .112_r8             ! mass transport coefficient (cm/s)
    real(r8), parameter ::  t0    = 298._r8             ! reference temperature (k)
    real(r8), parameter ::  xph0  = 1.e-5_r8            ! cloud [h+]
    real(r8), parameter ::  satf_hno3  = .016_r8        ! saturation factor for hno3 in clouds 
    real(r8), parameter ::  satf_h2o2  = .016_r8        ! saturation factor for h2o2 in clouds 
    real(r8), parameter ::  satf_so2   = .016_r8        ! saturation factor for so2 in clouds 
    real(r8), parameter ::  satf_ch2o  = .1_r8          ! saturation factor for ch2o in clouds 
    real(r8), parameter ::  satf_sog  =  .016_r8        ! saturation factor for sog in clouds
    real(r8), parameter ::  const0   = boltz_cgs * 1.e-6_r8 ! (atmospheres/deg k/cm^3)
    real(r8), parameter ::  hno3_diss = 15.4_r8         ! hno3 dissociation constant
    real(r8), parameter ::  geo_fac  = 6._r8            ! geometry factor (surf area/volume = geo_fac/diameter)
    real(r8), parameter ::  mass_air = 29._r8           ! mass of background atmosphere (amu)
    real(r8), parameter ::  mass_h2o = 18._r8           ! mass of water vapor (amu)
    real(r8), parameter ::  h2o_mol  = 1.e3_r8/mass_h2o ! (gm/mol water)
    real(r8), parameter ::  km2cm    = 1.e5_r8          ! convert km to cm
    real(r8), parameter ::  m2km     = 1.e-3_r8         ! convert m to km
    real(r8), parameter ::  cm3_2_m3 = 1.e-6_r8         ! convert cm^3 to m^3
    real(r8), parameter ::  m3_2_cm3 = 1.e6_r8          ! convert m^3 to cm^3
    real(r8), parameter ::  liter_per_gram = 1.e-3_r8
    real(r8), parameter ::  avo2  = avo * liter_per_gram * cm3_2_m3 ! (liter/gm/mol*(m/cm)^3)

    integer  ::      i, m, k, kk                 ! indicies
    real(r8) ::      xkgm                        ! mass flux on rain drop
    real(r8) ::      all1, all2                  ! work variables
    real(r8) ::      stay                        ! fraction of layer traversed by falling drop in timestep delt
    real(r8) ::      xeqca1, xeqca2, xca1, xca2, xdtm
    real(r8) ::      xxx1, xxx2, yhno3, yh2o2
    real(r8) ::      all3, xeqca3, xca3, xxx3, yso2, so2_diss
    real(r8) ::      all4, xeqca4, xca4, xxx4
    real(r8) ::      all5, xeqca5, xca5, xxx5
    real(r8) ::      all6, xeqca6, xca6, xxx6
    real(r8) ::      all7, xeqca7, xca7, xxx7
    real(r8) ::      all8, xeqca8, xca8, xxx8
    real(r8) ::      ysogm,ysogi,ysogt,ysogb,ysogx

    real(r8), dimension(ncol)  :: &
         xk0, work1, work2, work3, zsurf
    real(r8), dimension(pver)  :: &
         xgas1, xgas2
    real(r8), dimension(pver)  :: xgas3, xgas4, xgas5, xgas6, xgas7, xgas8
    real(r8), dimension(ncol)  :: &
         tmp0_rates, tmp1_rates
    real(r8), dimension(ncol,pver)  :: &
         delz, &              ! layer depth about interfaces (cm)
         xhno3, &             ! hno3 concentration (molecules/cm^3)
         xh2o2, &             ! h2o2 concentration (molecules/cm^3)
         xso2, &              ! so2 concentration (molecules/cm^3)
         xsogm, &             ! sogm concentration (molecules/cm^3)
         xsogi, &             ! sogi concentration (molecules/cm^3)
         xsogt, &             ! sogt concentration (molecules/cm^3)
         xsogb, &             ! sogb concentration (molecules/cm^3)
         xsogx, &             ! sogx concentration (molecules/cm^3)
         xliq, &              ! liquid rain water content in a grid cell (gm/m^3)
         rain                 ! conversion rate of water vapor into rain water (molecules/cm^3/s)
    real(r8), dimension(ncol,pver)  :: &
         xhen_hno3, xhen_h2o2, xhen_ch2o, xhen_ch3ooh, xhen_ch3co3h, &
         xhen_ch3cocho, xhen_xooh, xhen_onitr, xhen_ho2no2, xhen_glyald, &
         xhen_ch3cho, xhen_mvk, xhen_macr,xhen_sog
    real(r8), dimension(ncol,pver)  :: &
         xhen_nh3, xhen_ch3cooh
    real(r8), dimension(ncol,pver,8) :: tmp_hetrates
    real(r8), dimension(ncol,pver)  :: precip
    real(r8), dimension(ncol,pver)  :: xhen_hcn, xhen_ch3cn, xhen_so2

    integer    ::      ktop_all       
    integer    ::      ktop(ncol)                  ! 100 mb level

    real(r8) :: rlat(pcols)                       ! latitude in radians for columns
    real(r8) :: p_limit
    real(r8), parameter :: d2r = pi/180._r8
!
! jfl : new variables for rescaling sum of positive values to actual amount
!
    real(r8) :: total_rain,total_pos
    character(len=3) :: hetratestrg
    real(r8), parameter :: MISSING = -999999._r8
    integer ::  mm

!
    !-----------------------------------------------------------------
    !        note: the press array is in pascals and must be
    !              mutiplied by 10 to yield dynes/cm**2.
    !-----------------------------------------------------------------
    !       ... set wet deposition for
    !           1. h2o2         2. hno3
    !           3. ch2o         4. ch3ooh
    !           5. pooh         6. ch3coooh
    !           7. ho2no2       8. onit
    !           9. mvk         10. macr
    !          11. c2h5ooh     12. c3h7ooh
    !          13. rooh        14. ch3cocho
    !          15. pb          16. macrooh
    !          17. xooh        18. onitr
    !          19. isopooh     20. ch3oh
    !          21. c2h5oh      22. glyald
    !          23. hyac        24. hydrald
    !          25. ch3cho      26. isopno3
    !-----------------------------------------------------------------

    het_rates(:,:,:) = 0._r8

    if ( .not. do_wetdep) return

    call get_rlat_all_p(lchnk, ncol, rlat)

    do mm = 1,gas_wetdep_cnt
       m = wetdep_map(mm)
       if ( m>0 ) then
          het_rates(:,:,m) = MISSING
       endif
    end do

    !-----------------------------------------------------------------
    !	... the 2 and .6 multipliers are from a formula by frossling (1938)
    !-----------------------------------------------------------------
    xkgm = xdg/xrm * 2._r8 + xdg/xrm * .6_r8 * sqrt( xrm*xum/xvv ) * (xvv/xdg)**(1._r8/3._r8) 

    !-----------------------------------------------------------------
    !	... Find 100 mb level
    !-----------------------------------------------------------------
    do i = 1,ncol
       if ( abs(rlat(i)) > 60._r8*d2r ) then
          p_limit = 300.e2_r8
       else
          p_limit = 100.e2_r8 
       endif
       k_loop: do k = pver,1,-1
          if( press(i,k) < p_limit ) then
             ktop(i) = k
             exit k_loop
          end if
       end do k_loop
    end do
    ktop_all = minval( ktop(:) )
!
! jfl
!
! this is added to rescale the variable precip (which can only be positive)
! to the actual vertical integral of positive and negative values.  This
! removes point storms
!
    do i = 1,ncol
       total_rain = 0._r8
       total_pos  = 0._r8
       do k = 1,pver
          precip(i,k) = cmfdqr(i,k) + nrain(i,k) - nevapr(i,k)
          total_rain = total_rain + precip(i,k)
          if ( precip(i,k) < 0._r8 ) precip(i,k) = 0._r8
          total_pos  = total_pos  + precip(i,k)
       end do
       if ( total_rain <= 0._r8 ) then
          precip(i,:) = 0._r8
       else
          do k = 1,pver
             precip(i,k) = precip(i,k) * total_rain/total_pos
          end do
       end if
    end do

    do k = 1,pver
       !jfl       precip(:ncol,k) = cmfdqr(:ncol,k) + nrain(:ncol,k) - nevapr(:ncol,k)
       rain(:ncol,k)   = mass_air*precip(:ncol,k)*xhnm(:ncol,k) / mass_h2o
       xliq(:ncol,k)   = precip(:ncol,k) * delt * xhnm(:ncol,k) / avo*mass_air * m3_2_cm3
       if( spc_hno3_ndx > 0 ) then
          xhno3(:ncol,k)  = qin(:ncol,k,spc_hno3_ndx) * xhnm(:ncol,k)
       else
          xhno3(:ncol,k)  = 0._r8
       end if
       if( spc_h2o2_ndx > 0 ) then
          xh2o2(:ncol,k)  = qin(:ncol,k,spc_h2o2_ndx) * xhnm(:ncol,k)
       else
          xh2o2(:ncol,k)  = 0._r8
       end if
       if( spc_sogm_ndx > 0 ) then
          xsogm(:ncol,k)  = qin(:ncol,k,spc_sogm_ndx) * xhnm(:ncol,k)
       else
          xsogm(:ncol,k)  = 0._r8
       end if
       if( spc_sogi_ndx > 0 ) then
          xsogi(:ncol,k)  = qin(:ncol,k,spc_sogi_ndx) * xhnm(:ncol,k)
       else
          xsogi(:ncol,k)  = 0._r8
       end if
       if( spc_sogt_ndx > 0 ) then
          xsogt(:ncol,k)  = qin(:ncol,k,spc_sogt_ndx) * xhnm(:ncol,k)
       else
          xsogt(:ncol,k)  = 0._r8
       end if
       if( spc_sogb_ndx > 0 ) then
          xsogb(:ncol,k)  = qin(:ncol,k,spc_sogb_ndx) * xhnm(:ncol,k)
       else
          xsogb(:ncol,k)  = 0._r8
       end if
       if( spc_sogx_ndx > 0 ) then
          xsogx(:ncol,k)  = qin(:ncol,k,spc_sogx_ndx) * xhnm(:ncol,k)
       else
          xsogx(:ncol,k)  = 0._r8
       end if
       if( spc_so2_ndx > 0 ) then
          xso2(:ncol,k)  = qin(:ncol,k,spc_so2_ndx) * xhnm(:ncol,k)
       else
          xso2(:ncol,k)  = 0._r8
       end if
    end do

    zsurf(:ncol) = m2km * phis(:ncol) * rga
    do k = ktop_all,pver-1
       delz(:ncol,k) = abs( (zmid(:ncol,k) - zmid(:ncol,k+1))*km2cm ) 
    end do
    delz(:ncol,pver) = abs( (zmid(:ncol,pver) - zsurf(:ncol) )*km2cm ) 

    !-----------------------------------------------------------------
    !       ... part 0b,  for temperature dependent of henrys
    !                     xxhe1 = henry con for hno3
    !                     xxhe2 = henry con for h2o2
    !lwh 10/00 -- take henry''s law constants from brasseur et al. [1999],
    !             appendix j. for hno3, also consider dissociation to
    !             get effective henry''s law constant; equilibrium
    !             constant for dissociation from brasseur et al. [1999],
    !             appendix k. assume ph=5 (set as xph0 above).
    !             heff = h*k/[h+] for hno3 (complete dissociation)
    !             heff = h for h2o2 (no dissociation)
    !             heff = h * (1 + k/[h+]) (in general)
    !-----------------------------------------------------------------
    do k = ktop_all,pver
       work1(:ncol) = (t0 - tfld(:ncol,k))/(t0*tfld(:ncol,k))
       !-----------------------------------------------------------------
       ! 	... effective henry''s law constants:
       !	hno3, h2o2, ch2o, ch3ooh, ch3coooh (brasseur et al., 1999)
       !       xooh, onitr, macrooh (j.-f. muller; brocheton, 1999)
       !       isopooh (equal to hno3, as for macrooh)
       !       ho2no2 (mozart-1)
       !       ch3cocho, hoch2cho (betterton and hoffman, environ. sci. technol., 1988)
       !       ch3cho (staudinger and roberts, crit. rev. sci. technol., 1996)
       !       mvk, macr (allen et al., environ. toxicol. chem., 1998)
       !-----------------------------------------------------------------
       xk0(:)             = 2.1e5_r8 *exp( 8700._r8*work1(:) )
       xhen_hno3(:,k)     = xk0(:) * ( 1._r8 + hno3_diss / xph0 )
       xhen_h2o2(:,k)     = 7.45e4_r8 * exp( 6620._r8 * work1(:) )
       xhen_ch2o(:,k)     = 6.3e3_r8 * exp( 6460._r8 * work1(:) )
       xhen_ch3ooh(:,k)   = 2.27e2_r8 * exp( 5610._r8 * work1(:) )
       xhen_ch3co3h(:,k)  = 4.73e2_r8 * exp( 6170._r8 * work1(:) )
       xhen_ch3cocho(:,k) = 3.70e3_r8 * exp( 7275._r8 * work1(:) )
       xhen_xooh(:,k)     = 90.5_r8 * exp( 5607._r8 * work1(:) )
       xhen_onitr(:,k)    = 7.51e3_r8 * exp( 6485._r8 * work1(:) )
       xhen_ho2no2(:,k)   = 2.e4_r8
       xhen_glyald(:,k)   = 4.1e4_r8 * exp( 4600._r8 * work1(:) )
       xhen_ch3cho(:,k)   = 1.4e1_r8 * exp( 5600._r8 * work1(:) )
       xhen_mvk(:,k)      = 21._r8 * exp( 7800._r8 * work1(:) )
       xhen_macr(:,k)     = 4.3_r8 * exp( 5300._r8 * work1(:) )
       xhen_ch3cooh(:,k)  = 4.1e3_r8 * exp( 6300._r8 * work1(:) )
       xhen_sog(:,k)      = 5.e5_r8 * exp (12._r8 * work1(:) )
       !
       ! calculation for NH3 using the parameters in drydep_tables.F90
       !
       xhen_nh3 (:,k)     = 1.e6_r8
       xhen_ch3cn(:,k)     = 50._r8 * exp( 4000._r8 * work1(:) )
       xhen_hcn(:,k)       = 12._r8 * exp( 5000._r8 * work1(:) )
       do i = 1, ncol
          so2_diss        = 1.23e-2_r8 * exp( 1960._r8 * work1(i) )
          xhen_so2(i,k)   = 1.23_r8 * exp( 3120._r8 * work1(i) ) * ( 1._r8 + so2_diss / xph0 )
       end do
       !
       tmp_hetrates(:,k,:) = 0._r8
    end do

    !-----------------------------------------------------------------
    !       ... part 1, solve for high henry constant ( hno3, h2o2)
    !-----------------------------------------------------------------
    col_loop :  do i = 1,ncol
       xgas1(:) = xhno3(i,:)                     ! xgas will change during 
       xgas2(:) = xh2o2(i,:)                     ! different levels wash 
       xgas3(:) = xso2 (i,:)
       xgas4(:) = xsogm(i,:)
       xgas5(:) = xsogi(i,:)
       xgas6(:) = xsogt(i,:)
       xgas7(:) = xsogb(i,:)
       xgas8(:) = xsogx(i,:)
       level_loop1  : do kk = ktop(i),pver
          stay = 1._r8
          if( rain(i,kk) /= 0._r8 ) then            ! finding rain cloud           
             all1 = 0._r8                           ! accumulation to justisfy saturation
             all2 = 0._r8 
             all3 = 0._r8 
             all4 = 0._r8 
             all5 = 0._r8 
             all6 = 0._r8 
             all7 = 0._r8 
             all8 = 0._r8 
             stay = ((zmid(i,kk) - zsurf(i))*km2cm)/(xum*delt)
             stay = min( stay,1._r8 )
             !-----------------------------------------------------------------
             !       ... calculate the saturation concentration eqca
             !-----------------------------------------------------------------
             do k = kk,pver                      ! cal washout below cloud
                xeqca1 =  xgas1(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_hno3(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca2 =  xgas2(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_h2o2(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca3 =  xgas3(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_so2( i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca4 =  xgas4(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca5 =  xgas5(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca6 =  xgas6(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca7 =  xgas7(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2
                xeqca8 =  xgas8(k) &
                     / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
                     *  xliq(i,kk)*avo2

                !-----------------------------------------------------------------
                !       ... calculate ca; inside cloud concentration in #/cm3(air)
                !-----------------------------------------------------------------
                xca1 = geo_fac*xkgm*xgas1(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca2 = geo_fac*xkgm*xgas2(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca3 = geo_fac*xkgm*xgas3(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca4 = geo_fac*xkgm*xgas4(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca5 = geo_fac*xkgm*xgas5(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca6 = geo_fac*xkgm*xgas6(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca7 = geo_fac*xkgm*xgas7(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
                xca8 = geo_fac*xkgm*xgas8(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3

                !-----------------------------------------------------------------
                !       ... if is not saturated
                !               hno3(gas)_new = hno3(gas)_old - hno3(h2o)
                !           otherwise
                !               hno3(gas)_new = hno3(gas)_old
                !-----------------------------------------------------------------
                all1 = all1 + xca1
                all2 = all2 + xca2
                if( all1 < xeqca1 ) then
                   xgas1(k) = max( xgas1(k) - xca1,0._r8 )
                end if
                if( all2 < xeqca2 ) then
                   xgas2(k) = max( xgas2(k) - xca2,0._r8 )
                end if
                all3 = all3 + xca3
                if( all3 < xeqca3 ) then
                   xgas3(k) = max( xgas3(k) - xca3,0._r8 )
                end if
                all4 = all4 + xca4
                all5 = all5 + xca5
                all6 = all6 + xca6
                all7 = all7 + xca7
                all8 = all8 + xca8
                if( all4 < xeqca4 ) then
                   xgas4(k) = max( xgas4(k) - xca4,0._r8 )
                end if
                if( all5 < xeqca5 ) then
                   xgas5(k) = max( xgas5(k) - xca5,0._r8 )
                end if
                if( all6 < xeqca6 ) then
                   xgas6(k) = max( xgas6(k) - xca6,0._r8 )
                end if
                if( all7 < xeqca7 ) then
                   xgas7(k) = max( xgas7(k) - xca7,0._r8 )
                end if
                if( all8 < xeqca8 ) then
                   xgas8(k) = max( xgas8(k) - xca8,0._r8 )
                end if
             end do
          end if
          !-----------------------------------------------------------------
          !       ... calculate the lifetime of washout (second)
          !             after all layers washout 
          !             the concentration of hno3 is reduced 
          !             then the lifetime xtt is calculated by
          !
          !                  xtt = (xhno3(ini) - xgas1(new))/(dt*xhno3(ini))
          !                  where dt = passing time (s) in vertical
          !                             path below the cloud
          !                        dt = dz(cm)/um(cm/s)
          !-----------------------------------------------------------------
          xdtm = delz(i,kk) / xum                     ! the traveling time in each dz
          xxx1 = (xhno3(i,kk) - xgas1(kk))
          xxx2 = (xh2o2(i,kk) - xgas2(kk))
          if( xxx1 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             yhno3  = xhno3(i,kk)/xxx1 * xdtm    
          else
             yhno3  = 1.e29_r8
          end if
          if( xxx2 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             yh2o2  = xh2o2(i,kk)/xxx2 * xdtm     
          else
             yh2o2  = 1.e29_r8
          end if
          tmp_hetrates(i,kk,1) = max( 1._r8 / yh2o2,0._r8 ) * stay
          tmp_hetrates(i,kk,2) = max( 1._r8 / yhno3,0._r8 ) * stay
          xxx3 = (xso2( i,kk) - xgas3(kk))
          if( xxx3 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             yso2  = xso2( i,kk)/xxx3 * xdtm     
          else
             yso2  = 1.e29_r8
          end if
          tmp_hetrates(i,kk,3) = max( 1._r8 / yso2, 0._r8 ) * stay
          xxx4 = (xsogm(i,kk) - xgas4(kk))
          xxx5 = (xsogi(i,kk) - xgas5(kk))
          xxx6 = (xsogt(i,kk) - xgas6(kk))
          xxx7 = (xsogb(i,kk) - xgas7(kk))
          xxx8 = (xsogx(i,kk) - xgas8(kk))
          if( xxx4 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             ysogm  = xsogm(i,kk)/xxx4 * xdtm
          else
             ysogm  = 1.e29_r8
          end if
          if( xxx5 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             ysogi  = xsogi(i,kk)/xxx5 * xdtm
          else
             ysogi  = 1.e29_r8
          end if
          if( xxx6 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             ysogt  = xsogt(i,kk)/xxx6 * xdtm
          else
             ysogt  = 1.e29_r8
          end if
          if( xxx7 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             ysogb  = xsogb(i,kk)/xxx7 * xdtm
          else
             ysogb  = 1.e29_r8
          end if
          if( xxx8 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
             ysogx  = xsogx(i,kk)/xxx8 * xdtm
          else
             ysogx  = 1.e29_r8
          end if
          tmp_hetrates(i,kk,4) = max( 1._r8 / ysogm,0._r8 ) * stay
          tmp_hetrates(i,kk,5) = max( 1._r8 / ysogi,0._r8 ) * stay
          tmp_hetrates(i,kk,6) = max( 1._r8 / ysogt,0._r8 ) * stay
          tmp_hetrates(i,kk,7) = max( 1._r8 / ysogb,0._r8 ) * stay
          tmp_hetrates(i,kk,8) = max( 1._r8 / ysogx,0._r8 ) * stay
       end do level_loop1
    end do col_loop

    !-----------------------------------------------------------------
    !       ... part 2, in-cloud solve for low henry constant
    !                   hno3 and h2o2 have both in and under cloud
    !-----------------------------------------------------------------
    level_loop2 : do k = ktop_all,pver
       Column_loop2 : do i=1,ncol
          if ( rain(i,k) <= 0._r8 ) then
             het_rates(i,k,:) =  0._r8 
             cycle
          endif

          work1(i) = avo2 * xliq(i,k)
          work2(i) = const0 * tfld(i,k)
          work3(i) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch2o(i,k)*work2(i)))),0._r8 ) &
               * satf_ch2o
          if( ch2o_ndx > 0 ) then
             het_rates(i,k,ch2o_ndx)  = work3(i)
          end if
          if( isopno3_ndx > 0 ) then
             het_rates(i,k,isopno3_ndx) = work3(i)
          end if
          if( xisopno3_ndx > 0 ) then
             het_rates(i,k,xisopno3_ndx) = work3(i)
          end if
          if( hyac_ndx > 0 ) then
             het_rates(i,k,hyac_ndx) = work3(i)
          end if
          if( hydrald_ndx > 0 ) then
             het_rates(i,k,hydrald_ndx) = work3(i)
          end if

          work3(i) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3ooh(i,k)*work2(i)))),0._r8 )
          if( ch3ooh_ndx > 0 ) then
             het_rates(i,k,ch3ooh_ndx)  = work3(i)
          end if
          if( soag0_ndx > 0 ) then
             het_rates(i,k,soag0_ndx)  = work3(i)
          endif
          if( soag15_ndx > 0 ) then
             het_rates(i,k,soag15_ndx)  = work3(i)
          endif
          if( soag24_ndx > 0 ) then
             het_rates(i,k,soag24_ndx)  = work3(i)
          endif
          if( soag31_ndx > 0 ) then
             het_rates(i,k,soag31_ndx)  = work3(i)
          endif
          if( soag32_ndx > 0 ) then
             het_rates(i,k,soag32_ndx)  = work3(i)
          endif
          if( soag33_ndx > 0 ) then
             het_rates(i,k,soag33_ndx)  = work3(i)
          endif
          if( soag34_ndx > 0 ) then
             het_rates(i,k,soag34_ndx)  = work3(i)
          endif
          if( soag35_ndx > 0 ) then
             het_rates(i,k,soag35_ndx)  = work3(i)
          endif
          if( bsoag1_ndx > 0 ) then
             het_rates(i,k,bsoag1_ndx)  = work3(i)
          endif
          if( bsoag2_ndx > 0 ) then
             het_rates(i,k,bsoag2_ndx)  = work3(i)
          endif
          if( bsoag3_ndx > 0 ) then
             het_rates(i,k,bsoag3_ndx)  = work3(i)
          endif
          if( bsoag4_ndx > 0 ) then
             het_rates(i,k,bsoag4_ndx)  = work3(i)
          endif
          if( bsoag5_ndx > 0 ) then
             het_rates(i,k,bsoag5_ndx)  = work3(i)
          endif
          if( bsoag6_ndx > 0 ) then
             het_rates(i,k,bsoag6_ndx)  = work3(i)
          endif
          if( bsoag7_ndx > 0 ) then
             het_rates(i,k,bsoag7_ndx)  = work3(i)
          endif
          if( bsoag8_ndx > 0 ) then
             het_rates(i,k,bsoag8_ndx)  = work3(i)
          endif
          if( bsoag9_ndx > 0 ) then
             het_rates(i,k,bsoag9_ndx)  = work3(i)
          endif
          if( bsoag10_ndx > 0 ) then
             het_rates(i,k,bsoag10_ndx)  = work3(i)
          endif
          if( bsoag11_ndx > 0 ) then
             het_rates(i,k,bsoag11_ndx)  = work3(i)
          endif
          if( bsoag12_ndx > 0 ) then
             het_rates(i,k,bsoag12_ndx)  = work3(i)
          endif
          if( bsoag13_ndx > 0 ) then
             het_rates(i,k,bsoag13_ndx)  = work3(i)
          endif
          if( bsoag14_ndx > 0 ) then
             het_rates(i,k,bsoag14_ndx)  = work3(i)
          endif
          if( bsoag15_ndx > 0 ) then
             het_rates(i,k,bsoag15_ndx)  = work3(i)
          endif
          if( bsoag16_ndx > 0 ) then
             het_rates(i,k,bsoag16_ndx)  = work3(i)
          endif
          if( bsoag17_ndx > 0 ) then
             het_rates(i,k,bsoag17_ndx)  = work3(i)
          endif
          if( pooh_ndx > 0 ) then
             het_rates(i,k,pooh_ndx)  = work3(i)
          end if
          if( c2h5ooh_ndx > 0 ) then
             het_rates(i,k,c2h5ooh_ndx) = work3(i)
          end if
          if( c3h7ooh_ndx > 0 ) then
             het_rates(i,k,c3h7ooh_ndx) = work3(i)
          end if
          if( rooh_ndx > 0 ) then
             het_rates(i,k,rooh_ndx) = work3(i)
          end if
          if( ch3oh_ndx > 0 ) then
             het_rates(i,k,ch3oh_ndx) = work3(i)
          end if
          if( c2h5oh_ndx > 0 ) then
             het_rates(i,k,c2h5oh_ndx) = work3(i)
          end if
          if( alkooh_ndx  > 0 ) then
             het_rates(i,k,alkooh_ndx) = work3(i)
          end if
          if( mekooh_ndx  > 0 ) then
             het_rates(i,k,mekooh_ndx) = work3(i)
          end if
          if( tolooh_ndx  > 0 ) then
             het_rates(i,k,tolooh_ndx) = work3(i)
          end if
          if( terpooh_ndx > 0 ) then
             het_rates(i,k,terpooh_ndx) = work3(i)
          end if

          if( ch3coooh_ndx > 0 ) then
             het_rates(i,k,ch3coooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3co3h(i,k)*work2(i)))),0._r8 )
          end if
          if( ho2no2_ndx > 0 ) then
             het_rates(i,k,ho2no2_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ho2no2(i,k)*work2(i)))),0._r8 )
          end if
          if( xho2no2_ndx > 0 ) then
             het_rates(i,k,xho2no2_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ho2no2(i,k)*work2(i)))),0._r8 )
          end if
          if( ch3cocho_ndx > 0 ) then
             het_rates(i,k,ch3cocho_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cocho(i,k)*work2(i)))),0._r8 )
          end if
          if( xooh_ndx > 0 ) then
             het_rates(i,k,xooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_xooh(i,k)*work2(i)))),0._r8 )
          end if
          if( onitr_ndx > 0 ) then
             het_rates(i,k,onitr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_onitr(i,k)*work2(i)))),0._r8 )
          end if
          if( xonitr_ndx > 0 ) then
             het_rates(i,k,xonitr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_onitr(i,k)*work2(i)))),0._r8 )
          end if
          if( glyald_ndx > 0 ) then
             het_rates(i,k,glyald_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_glyald(i,k)*work2(i)))),0._r8 )
          end if
          if( ch3cho_ndx > 0 ) then
             het_rates(i,k,ch3cho_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cho(i,k)*work2(i)))),0._r8 )
          end if
          if( mvk_ndx > 0 ) then
             het_rates(i,k,mvk_ndx)  = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_mvk(i,k)*work2(i)))),0._r8 )
          end if
          if( macr_ndx > 0 ) then
             het_rates(i,k,macr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_macr(i,k)*work2(i)))),0._r8 )
          end if
          if( h2o2_ndx > 0 ) then
             work3(i) = satf_h2o2 * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_h2o2(i,k)*work2(i)))),0._r8 )    
             het_rates(i,k,h2o2_ndx) =  work3(i) + tmp_hetrates(i,k,1)
          end if
          if ( prog_modal_aero .and. so2_ndx>0 .and. h2o2_ndx>0 ) then
             het_rates(i,k,so2_ndx) = het_rates(i,k,h2o2_ndx)
          elseif( so2_ndx > 0 ) then
             work3(i) = satf_so2 * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_so2( i,k)*work2(i)))),0._r8 )    
             het_rates(i,k,so2_ndx ) =  work3(i) + tmp_hetrates(i,k,3)
          endif
!
          work3(i) = satf_sog * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_sog(i,k)*work2(i)))),0._r8 )
          if( sogm_ndx > 0 ) then
             het_rates(i,k,sogm_ndx) =  work3(i) + tmp_hetrates(i,k,4)
          end if
          if( sogi_ndx > 0 ) then
             het_rates(i,k,sogi_ndx) =  work3(i) + tmp_hetrates(i,k,5)
          end if
          if( sogt_ndx > 0 ) then
             het_rates(i,k,sogt_ndx) =  work3(i) + tmp_hetrates(i,k,6)
          end if
          if( sogb_ndx > 0 ) then
             het_rates(i,k,sogb_ndx) =  work3(i) + tmp_hetrates(i,k,7)
          end if
          if( sogx_ndx > 0 ) then
             het_rates(i,k,sogx_ndx) =  work3(i) + tmp_hetrates(i,k,8)
          end if
!
          work3(i) = tmp_hetrates(i,k,2) + satf_hno3 * &
               max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_hno3(i,k)*work2(i)))),0._r8 )    
          tmp0_rates(i)   = work3(i)
          tmp1_rates(i)   = .2_r8*work3(i)
          if( hno3_ndx > 0 ) then
             het_rates(i,k,hno3_ndx) = work3(i)
          end if
          if( xhno3_ndx > 0 ) then
             het_rates(i,k,xhno3_ndx) = work3(i)
          end if
          if( onit_ndx > 0 ) then
             het_rates(i,k,onit_ndx) = work3(i)
          end if
          if( xonit_ndx > 0 ) then
             het_rates(i,k,xonit_ndx) = work3(i)
          end if
          if( Pb_ndx > 0 ) then
             het_rates(i,k,Pb_ndx) = work3(i)
          end if
          if( macrooh_ndx > 0 ) then
             het_rates(i,k,macrooh_ndx) = work3(i)
          end if
          if( isopooh_ndx > 0 ) then
             het_rates(i,k,isopooh_ndx) = work3(i)
          end if

          if( clono2_ndx > 0 ) then
             het_rates(i,k, clono2_ndx) = work3(i)
          end if
          if( brono2_ndx > 0 ) then
             het_rates(i,k, brono2_ndx) = work3(i)
          end if
          if( hcl_ndx > 0 ) then
             het_rates(i,k, hcl_ndx) = work3(i)
          end if
          if( n2o5_ndx > 0 ) then
             het_rates(i,k, n2o5_ndx) = work3(i)
          end if
          if( hocl_ndx > 0 ) then
             het_rates(i,k, hocl_ndx) = work3(i)
          end if
          if( hobr_ndx > 0 ) then
             het_rates(i,k, hobr_ndx) = work3(i)
          end if
          if( hbr_ndx > 0 ) then
             het_rates(i,k, hbr_ndx) = work3(i)
          end if

          if( soa_ndx > 0 ) then
             het_rates(i,k,soa_ndx) = tmp1_rates(i)
          end if
          if( oc2_ndx > 0 ) then
             het_rates(i,k,oc2_ndx) = tmp1_rates(i)
          end if
          if( cb2_ndx > 0 ) then
             het_rates(i,k,cb2_ndx) = tmp1_rates(i)
          end if
          if( so4_ndx > 0 ) then
             het_rates(i,k,so4_ndx) = tmp1_rates(i)
          end if
          if( sa1_ndx > 0 ) then
             het_rates(i,k,sa1_ndx) = tmp1_rates(i)
          end if
          if( sa2_ndx > 0 ) then
             het_rates(i,k,sa2_ndx) = tmp1_rates(i)
          end if
          if( sa3_ndx > 0 ) then
             het_rates(i,k,sa3_ndx) = tmp1_rates(i)
          end if
          if( sa4_ndx > 0 ) then
             het_rates(i,k,sa4_ndx) = tmp1_rates(i)
          end if

          if( h2so4_ndx > 0 ) then
             het_rates(i,k,h2so4_ndx) = tmp0_rates(i)
          end if
          if( nh4_ndx > 0 ) then
             het_rates(i,k,nh4_ndx) = tmp0_rates(i)
          end if
          if( nh4no3_ndx > 0 ) then
             het_rates(i,k,nh4no3_ndx ) = tmp0_rates(i)
          end if
          if( nh3_ndx > 0 ) then
             het_rates(i,k,nh3_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_nh3(i,k)*work2(i)))),0._r8 )
          end if

          if( ch3cooh_ndx > 0 ) then
             het_rates(i,k,ch3cooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cooh(i,k)*work2(i)))),0._r8 )
          end if
          if( hcooh_ndx > 0 ) then
             het_rates(i,k,hcooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cooh(i,k)*work2(i)))),0._r8 )
          endif
          if ( hcn_ndx > 0 ) then
             het_rates(i,k,hcn_ndx     ) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_hcn(i,k)*work2(i)))),0._r8 )
          endif
          if ( ch3cn_ndx > 0 ) then
             het_rates(i,k,ch3cn_ndx   ) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cn(i,k)*work2(i)))),0._r8 )
          endif

       end do Column_loop2
    end do level_loop2

    !-----------------------------------------------------------------
    !	... Set rates above tropopause = 0.
    !-----------------------------------------------------------------
    do mm = 1,gas_wetdep_cnt
       m = wetdep_map(mm)
       do i = 1,ncol
          do k = 1,ktop(i)
             het_rates(i,k,m) = 0._r8
          end do
       end do
       if ( any( het_rates(:ncol,:,m) == MISSING) ) then
          write(hetratestrg,'(I3)') m
          call endrun('sethet: het_rates (wet dep) not set for het reaction number : '//hetratestrg)
       endif
    end do

  end subroutine sethet

end module mo_sethet
