LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_sethet.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 3 | 443 | 0.7 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 2 | 50.0 %

              Line data    Source code
           1             :
           2             : module mo_sethet
           3             :
           4             : !
           5             : ! LKE (10/11/2010): added HCN, CH3CN, HCOOH  to cesm1_0_beta07_offline version
           6             : !                   HCN, CH3CN have new Henry's Law coefficients, HCOOH is set to CH3COOH
           7             : ! LKE (10/18/2010): SO2 washout corrected based on recommendation of R.Easter, PNNL
           8             : !
           9             :   use shr_kind_mod,    only: r8 => shr_kind_r8
          10             :   use cam_logfile,     only: iulog
          11             :   use gas_wetdep_opts, only: gas_wetdep_cnt, gas_wetdep_method, gas_wetdep_list
          12             :   use phys_control,    only: phys_getopts
          13             :
          14             :   private
          15             :   public :: sethet_inti, sethet
          16             :
          17             :   save
          18             :
          19             :   integer :: h2o2_ndx, hno3_ndx, ch2o_ndx, ch3ooh_ndx, ch3coooh_ndx, &
          20             :        ho2no2_ndx, ch3cocho_ndx, xooh_ndx, onitr_ndx, glyald_ndx, &
          21             :        ch3cho_ndx, mvk_ndx, macr_ndx, pooh_ndx, c2h5ooh_ndx, &
          22             :        c3h7ooh_ndx, rooh_ndx, isopno3_ndx, onit_ndx, Pb_ndx, &
          23             :        macrooh_ndx, isopooh_ndx, ch3oh_ndx, c2h5oh_ndx, hyac_ndx, hydrald_ndx
          24             :   integer :: spc_h2o2_ndx, spc_hno3_ndx
          25             :   integer :: spc_so2_ndx
          26             :   integer :: spc_sogm_ndx, spc_sogi_ndx, spc_sogt_ndx, spc_sogb_ndx, spc_sogx_ndx
          27             :
          28             :   integer :: alkooh_ndx, mekooh_ndx, tolooh_ndx, terpooh_ndx, ch3cooh_ndx
          29             :   integer :: so2_ndx, soa_ndx, so4_ndx, cb2_ndx, oc2_ndx, nh3_ndx, nh4no3_ndx, &
          30             :              sa1_ndx, sa2_ndx, sa3_ndx, sa4_ndx, nh4_ndx, h2so4_ndx
          31             :   integer :: xisopno3_ndx,xho2no2_ndx,xonitr_ndx,xhno3_ndx,xonit_ndx
          32             :   integer :: clono2_ndx, brono2_ndx, hcl_ndx, n2o5_ndx, hocl_ndx, hobr_ndx, hbr_ndx
          33             :   integer :: ch3cn_ndx, hcn_ndx, hcooh_ndx
          34             :   integer, allocatable :: wetdep_map(:)
          35             :   integer :: sogm_ndx, sogi_ndx, sogt_ndx, sogb_ndx, sogx_ndx
          36             :   logical :: do_wetdep
          37             :
          38             :   ! prognostic modal aerosols
          39             :   logical :: prog_modal_aero
          40             :
          41             : contains
          42             :
          43        1536 :   subroutine sethet_inti
          44             :     !-----------------------------------------------------------------------
          45             :     !       ... intialize the wet removal rate constants routine
          46             :     !-----------------------------------------------------------------------
          47             :
          48             :     use mo_chem_utls,     only : get_het_ndx, get_spc_ndx
          49             :     use spmd_utils,       only : masterproc
          50             :     use cam_abortutils,   only : endrun
          51             :
          52             :     integer :: k, m
          53             :
          54        1536 :     do_wetdep = gas_wetdep_cnt>0 .and. gas_wetdep_method=='MOZ'
          55        1536 :     if ( .not. do_wetdep) return
          56             :
          57           0 :     call phys_getopts( prog_modal_aero_out = prog_modal_aero )
          58             :
          59           0 :     allocate( wetdep_map(gas_wetdep_cnt))
          60             :
          61           0 :     do k=1,gas_wetdep_cnt
          62           0 :        m = get_het_ndx( trim(gas_wetdep_list(k)))
          63           0 :        if (m>0) then
          64           0 :           wetdep_map(k) = m
          65             :        else
          66           0 :           call endrun('sethet_inti: cannot map '//trim(gas_wetdep_list(k)))
          67             :        endif
          68             :     enddo
          69             :
          70           0 :     xisopno3_ndx = get_het_ndx( 'XISOPNO3' )
          71           0 :     xho2no2_ndx  = get_het_ndx( 'XHO2NO2' )
          72           0 :     xonitr_ndx   = get_het_ndx( 'XONITR' )
          73           0 :     xhno3_ndx    = get_het_ndx( 'XHNO3' )
          74           0 :     xonit_ndx    = get_het_ndx( 'XONIT' )
          75             :
          76           0 :     spc_h2o2_ndx = get_spc_ndx( 'H2O2' )
          77           0 :     spc_hno3_ndx = get_spc_ndx( 'HNO3' )
          78           0 :     spc_so2_ndx  = get_spc_ndx( 'SO2' )
          79             :
          80           0 :     clono2_ndx = get_het_ndx( 'CLONO2' )
          81           0 :     brono2_ndx = get_het_ndx( 'BRONO2' )
          82           0 :     hcl_ndx    = get_het_ndx( 'HCL' )
          83           0 :     n2o5_ndx   = get_het_ndx( 'N2O5' )
          84           0 :     hocl_ndx   = get_het_ndx( 'HOCL' )
          85           0 :     hobr_ndx   = get_het_ndx( 'HOBR' )
          86           0 :     hbr_ndx    = get_het_ndx( 'HBR' )
          87             :
          88           0 :     h2o2_ndx   = get_het_ndx( 'H2O2' )
          89           0 :     hno3_ndx   = get_het_ndx( 'HNO3' )
          90           0 :     ch2o_ndx   = get_het_ndx( 'CH2O' )
          91           0 :     ch3ooh_ndx = get_het_ndx( 'CH3OOH' )
          92           0 :     ch3coooh_ndx = get_het_ndx( 'CH3COOOH' )
          93           0 :     ho2no2_ndx  = get_het_ndx( 'HO2NO2' )
          94           0 :     ch3cocho_ndx = get_het_ndx( 'CH3COCHO' )
          95           0 :     xooh_ndx    = get_het_ndx( 'XOOH' )
          96           0 :     onitr_ndx   = get_het_ndx( 'ONITR' )
          97           0 :     glyald_ndx  = get_het_ndx( 'GLYALD' )
          98           0 :     ch3cho_ndx  = get_het_ndx( 'CH3CHO' )
          99           0 :     mvk_ndx     = get_het_ndx( 'MVK' )
         100           0 :     macr_ndx    = get_het_ndx( 'MACR' )
         101           0 :     pooh_ndx    = get_het_ndx( 'POOH' )
         102           0 :     c2h5ooh_ndx = get_het_ndx( 'C2H5OOH' )
         103           0 :     c3h7ooh_ndx = get_het_ndx( 'C3H7OOH' )
         104           0 :     rooh_ndx    = get_het_ndx( 'ROOH' )
         105           0 :     isopno3_ndx = get_het_ndx( 'ISOPNO3' )
         106           0 :     onit_ndx    = get_het_ndx( 'ONIT' )
         107           0 :     Pb_ndx      = get_het_ndx( 'Pb' )
         108           0 :     macrooh_ndx = get_het_ndx( 'MACROOH' )
         109           0 :     isopooh_ndx = get_het_ndx( 'ISOPOOH' )
         110           0 :     ch3oh_ndx   = get_het_ndx( 'CH3OH' )
         111           0 :     c2h5oh_ndx  = get_het_ndx( 'C2H5OH' )
         112           0 :     hyac_ndx    = get_het_ndx( 'HYAC' )
         113           0 :     hydrald_ndx = get_het_ndx( 'HYDRALD' )
         114           0 :     alkooh_ndx  = get_het_ndx( 'ALKOOH' )
         115           0 :     mekooh_ndx  = get_het_ndx( 'MEKOOH' )
         116           0 :     tolooh_ndx  = get_het_ndx( 'TOLOOH' )
         117           0 :     terpooh_ndx = get_het_ndx( 'TERPOOH' )
         118           0 :     ch3cooh_ndx = get_het_ndx( 'CH3COOH' )
         119           0 :     so2_ndx     = get_het_ndx( 'SO2' )
         120           0 :     soa_ndx     = get_het_ndx( 'SOA' )
         121           0 :     sogb_ndx    = get_het_ndx( 'SOGB' )
         122           0 :     sogi_ndx    = get_het_ndx( 'SOGI' )
         123           0 :     sogm_ndx    = get_het_ndx( 'SOGM' )
         124           0 :     sogt_ndx    = get_het_ndx( 'SOGT' )
         125           0 :     sogx_ndx    = get_het_ndx( 'SOGX' )
         126           0 :     so4_ndx     = get_het_ndx( 'SO4' )
         127           0 :     cb2_ndx     = get_het_ndx( 'CB2' )
         128           0 :     oc2_ndx     = get_het_ndx( 'OC2' )
         129           0 :     nh3_ndx     = get_het_ndx( 'NH3' )
         130           0 :     nh4no3_ndx  = get_het_ndx( 'NH4NO3' )
         131           0 :     nh4_ndx     = get_het_ndx( 'NH4' )
         132           0 :     h2so4_ndx   = get_het_ndx( 'H2SO4' )
         133           0 :     sa1_ndx     = get_het_ndx( 'SA1' )
         134           0 :     sa2_ndx     = get_het_ndx( 'SA2' )
         135           0 :     sa3_ndx     = get_het_ndx( 'SA3' )
         136           0 :     sa4_ndx     = get_het_ndx( 'SA4' )
         137           0 :     ch3cn_ndx   = get_het_ndx( 'CH3CN' )
         138           0 :     hcn_ndx     = get_het_ndx( 'HCN' )
         139           0 :     hcooh_ndx   = get_het_ndx( 'HCOOH' )
         140             :
         141           0 :     if (masterproc) then
         142           0 :        write(iulog,*) 'sethet_inti: new ndx ',so2_ndx,soa_ndx,so4_ndx,cb2_ndx,oc2_ndx, &
         143           0 :             nh3_ndx,nh4no3_ndx,sa1_ndx,sa2_ndx,sa3_ndx,sa4_ndx
         144           0 :        write(iulog,*) ' '
         145           0 :        write(iulog,*) 'sethet_inti: diagnotics '
         146           0 :        write(iulog,'(10i5)') h2o2_ndx, hno3_ndx, ch2o_ndx, ch3ooh_ndx, ch3coooh_ndx, &
         147           0 :             ho2no2_ndx, ch3cocho_ndx, xooh_ndx, onitr_ndx, glyald_ndx, &
         148           0 :             ch3cho_ndx, mvk_ndx, macr_ndx, pooh_ndx, c2h5ooh_ndx, &
         149           0 :             c3h7ooh_ndx, rooh_ndx, isopno3_ndx, onit_ndx, Pb_ndx, &
         150           0 :             macrooh_ndx, isopooh_ndx, ch3oh_ndx, c2h5oh_ndx, hyac_ndx, hydrald_ndx
         151             :     endif
         152             :
         153             :   end subroutine sethet_inti
         154             :
         155           0 :   subroutine sethet( het_rates, press, zmid,  phis, tfld, &
         156           0 :                      cmfdqr, nrain, nevapr, delt, xhnm, &
         157           0 :                      qin, ncol, lchnk )
         158             :     !-----------------------------------------------------------------------
         159             :     !       ... compute rainout loss rates (1/s)
         160             :     !-----------------------------------------------------------------------
         161             :
         162             :     use physconst,        only : rga,pi
         163             :     use chem_mods,        only : gas_pcnst
         164             :     use ppgrid,           only : pver, pcols
         165             :     use phys_grid,        only : get_rlat_all_p
         166             :     use cam_abortutils,   only : endrun
         167             :     use mo_constants,     only : avo => avogadro, boltz_cgs
         168             :
         169             :     implicit none
         170             :     !-----------------------------------------------------------------------
         171             :     !       ... dummy arguments
         172             :     !-----------------------------------------------------------------------
         173             :     integer, intent(in)   ::    ncol                        ! columns in chunk
         174             :     integer, intent(in)   ::    lchnk                       ! chunk index
         175             :     real(r8), intent(in)  ::    delt                        ! time step ( s )
         176             :     real(r8), intent(in)  ::    press(pcols,pver)           ! pressure in pascals
         177             :     real(r8), intent(in)  ::    cmfdqr(ncol,pver)           ! dq/dt for convection
         178             :     real(r8), intent(in)  ::    nrain(ncol,pver)            ! stratoform precip
         179             :     real(r8), intent(in)  ::    nevapr(ncol,pver)           ! evaporation
         180             :     real(r8), intent(in)  ::    qin(ncol,pver,gas_pcnst)    ! xported species ( vmr )
         181             :     real(r8), intent(in)  ::    zmid(ncol,pver)             ! midpoint geopot (km)
         182             :     real(r8), intent(in)  ::    phis(pcols)                 ! surf geopot
         183             :     real(r8), intent(in)  ::    tfld(pcols,pver)            ! temperature (k)
         184             :     real(r8), intent(in)  ::    xhnm(ncol,pver)             ! total atms density ( /cm^3)
         185             :     real(r8), intent(out) ::    het_rates(ncol,pver,gas_pcnst) ! rainout loss rates
         186             :
         187             :     !-----------------------------------------------------------------------
         188             :     !       ... local variables
         189             :     !-----------------------------------------------------------------------
         190             :     real(r8), parameter ::  xrm   = .189_r8             ! mean diameter of rain drop (cm)
         191             :     real(r8), parameter ::  xum   = 748._r8             ! mean rain drop terminal velocity (cm/s)
         192             :     real(r8), parameter ::  xvv   = 6.18e-2_r8          ! kinetic viscosity (cm^2/s)
         193             :     real(r8), parameter ::  xdg   = .112_r8             ! mass transport coefficient (cm/s)
         194             :     real(r8), parameter ::  t0    = 298._r8             ! reference temperature (k)
         195             :     real(r8), parameter ::  xph0  = 1.e-5_r8            ! cloud [h+]
         196             :     real(r8), parameter ::  satf_hno3  = .016_r8        ! saturation factor for hno3 in clouds
         197             :     real(r8), parameter ::  satf_h2o2  = .016_r8        ! saturation factor for h2o2 in clouds
         198             :     real(r8), parameter ::  satf_so2   = .016_r8        ! saturation factor for so2 in clouds
         199             :     real(r8), parameter ::  satf_ch2o  = .1_r8          ! saturation factor for ch2o in clouds
         200             :     real(r8), parameter ::  satf_sog  =  .016_r8        ! saturation factor for sog in clouds
         201             :     real(r8), parameter ::  const0   = boltz_cgs * 1.e-6_r8 ! (atmospheres/deg k/cm^3)
         202             :     real(r8), parameter ::  hno3_diss = 15.4_r8         ! hno3 dissociation constant
         203             :     real(r8), parameter ::  geo_fac  = 6._r8            ! geometry factor (surf area/volume = geo_fac/diameter)
         204             :     real(r8), parameter ::  mass_air = 29._r8           ! mass of background atmosphere (amu)
         205             :     real(r8), parameter ::  mass_h2o = 18._r8           ! mass of water vapor (amu)
         206             :     real(r8), parameter ::  h2o_mol  = 1.e3_r8/mass_h2o ! (gm/mol water)
         207             :     real(r8), parameter ::  km2cm    = 1.e5_r8          ! convert km to cm
         208             :     real(r8), parameter ::  m2km     = 1.e-3_r8         ! convert m to km
         209             :     real(r8), parameter ::  cm3_2_m3 = 1.e-6_r8         ! convert cm^3 to m^3
         210             :     real(r8), parameter ::  m3_2_cm3 = 1.e6_r8          ! convert m^3 to cm^3
         211             :     real(r8), parameter ::  liter_per_gram = 1.e-3_r8
         212             :     real(r8), parameter ::  avo2  = avo * liter_per_gram * cm3_2_m3 ! (liter/gm/mol*(m/cm)^3)
         213             :
         214             :     integer  ::      i, m, k, kk                 ! indicies
         215             :     real(r8) ::      xkgm                        ! mass flux on rain drop
         216             :     real(r8) ::      all1, all2                  ! work variables
         217             :     real(r8) ::      stay                        ! fraction of layer traversed by falling drop in timestep delt
         218             :     real(r8) ::      xeqca1, xeqca2, xca1, xca2, xdtm
         219             :     real(r8) ::      xxx1, xxx2, yhno3, yh2o2
         220             :     real(r8) ::      all3, xeqca3, xca3, xxx3, yso2, so2_diss
         221             :     real(r8) ::      all4, xeqca4, xca4, xxx4
         222             :     real(r8) ::      all5, xeqca5, xca5, xxx5
         223             :     real(r8) ::      all6, xeqca6, xca6, xxx6
         224             :     real(r8) ::      all7, xeqca7, xca7, xxx7
         225             :     real(r8) ::      all8, xeqca8, xca8, xxx8
         226             :     real(r8) ::      ysogm,ysogi,ysogt,ysogb,ysogx
         227             :
         228             :     real(r8), dimension(ncol)  :: &
         229           0 :          xk0, work1, work2, work3, zsurf
         230             :     real(r8), dimension(pver)  :: &
         231             :          xgas1, xgas2
         232             :     real(r8), dimension(pver)  :: xgas3, xgas4, xgas5, xgas6, xgas7, xgas8
         233             :     real(r8), dimension(ncol)  :: &
         234           0 :          tmp0_rates, tmp1_rates
         235             :     real(r8), dimension(ncol,pver)  :: &
         236           0 :          delz, &              ! layer depth about interfaces (cm)
         237           0 :          xhno3, &             ! hno3 concentration (molecules/cm^3)
         238           0 :          xh2o2, &             ! h2o2 concentration (molecules/cm^3)
         239           0 :          xso2, &              ! so2 concentration (molecules/cm^3)
         240           0 :          xsogm, &             ! sogm concentration (molecules/cm^3)
         241           0 :          xsogi, &             ! sogi concentration (molecules/cm^3)
         242           0 :          xsogt, &             ! sogt concentration (molecules/cm^3)
         243           0 :          xsogb, &             ! sogb concentration (molecules/cm^3)
         244           0 :          xsogx, &             ! sogx concentration (molecules/cm^3)
         245           0 :          xliq, &              ! liquid rain water content in a grid cell (gm/m^3)
         246           0 :          rain                 ! conversion rate of water vapor into rain water (molecules/cm^3/s)
         247             :     real(r8), dimension(ncol,pver)  :: &
         248           0 :          xhen_hno3, xhen_h2o2, xhen_ch2o, xhen_ch3ooh, xhen_ch3co3h, &
         249           0 :          xhen_ch3cocho, xhen_xooh, xhen_onitr, xhen_ho2no2, xhen_glyald, &
         250           0 :          xhen_ch3cho, xhen_mvk, xhen_macr,xhen_sog
         251             :     real(r8), dimension(ncol,pver)  :: &
         252           0 :          xhen_nh3, xhen_ch3cooh
         253           0 :     real(r8), dimension(ncol,pver,8) :: tmp_hetrates
         254           0 :     real(r8), dimension(ncol,pver)  :: precip
         255           0 :     real(r8), dimension(ncol,pver)  :: xhen_hcn, xhen_ch3cn, xhen_so2
         256             :
         257             :     integer    ::      ktop_all
         258           0 :     integer    ::      ktop(ncol)                  ! 100 mb level
         259             :
         260             :     real(r8) :: rlat(pcols)                       ! latitude in radians for columns
         261             :     real(r8) :: p_limit
         262             :     real(r8), parameter :: d2r = pi/180._r8
         263             : !
         264             : ! jfl : new variables for rescaling sum of positive values to actual amount
         265             : !
         266             :     real(r8) :: total_rain,total_pos
         267             :     character(len=3) :: hetratestrg
         268             :     real(r8), parameter :: MISSING = -999999._r8
         269             :     integer ::  mm
         270             :
         271             : !
         272             :     !-----------------------------------------------------------------
         273             :     !        note: the press array is in pascals and must be
         274             :     !              mutiplied by 10 to yield dynes/cm**2.
         275             :     !-----------------------------------------------------------------
         276             :     !       ... set wet deposition for
         277             :     !           1. h2o2         2. hno3
         278             :     !           3. ch2o         4. ch3ooh
         279             :     !           5. pooh         6. ch3coooh
         280             :     !           7. ho2no2       8. onit
         281             :     !           9. mvk         10. macr
         282             :     !          11. c2h5ooh     12. c3h7ooh
         283             :     !          13. rooh        14. ch3cocho
         284             :     !          15. pb          16. macrooh
         285             :     !          17. xooh        18. onitr
         286             :     !          19. isopooh     20. ch3oh
         287             :     !          21. c2h5oh      22. glyald
         288             :     !          23. hyac        24. hydrald
         289             :     !          25. ch3cho      26. isopno3
         290             :     !-----------------------------------------------------------------
         291             :
         292           0 :     het_rates(:,:,:) = 0._r8
         293             :
         294           0 :     if ( .not. do_wetdep) return
         295             :
         296           0 :     call get_rlat_all_p(lchnk, ncol, rlat)
         297             :
         298           0 :     do mm = 1,gas_wetdep_cnt
         299           0 :        m = wetdep_map(mm)
         300           0 :        if ( m>0 ) then
         301           0 :           het_rates(:,:,m) = MISSING
         302             :        endif
         303             :     end do
         304             :
         305             :     !-----------------------------------------------------------------
         306             :     !   ... the 2 and .6 multipliers are from a formula by frossling (1938)
         307             :     !-----------------------------------------------------------------
         308             :     xkgm = xdg/xrm * 2._r8 + xdg/xrm * .6_r8 * sqrt( xrm*xum/xvv ) * (xvv/xdg)**(1._r8/3._r8)
         309             :
         310             :     !-----------------------------------------------------------------
         311             :     !   ... Find 100 mb level
         312             :     !-----------------------------------------------------------------
         313           0 :     do i = 1,ncol
         314           0 :        if ( abs(rlat(i)) > 60._r8*d2r ) then
         315             :           p_limit = 300.e2_r8
         316             :        else
         317           0 :           p_limit = 100.e2_r8
         318             :        endif
         319           0 :        k_loop: do k = pver,1,-1
         320           0 :           if( press(i,k) < p_limit ) then
         321           0 :              ktop(i) = k
         322           0 :              exit k_loop
         323             :           end if
         324             :        end do k_loop
         325             :     end do
         326           0 :     ktop_all = minval( ktop(:) )
         327             : !
         328             : ! jfl
         329             : !
         330             : ! this is added to rescale the variable precip (which can only be positive)
         331             : ! to the actual vertical integral of positive and negative values.  This
         332             : ! removes point storms
         333             : !
         334           0 :     do i = 1,ncol
         335             :        total_rain = 0._r8
         336             :        total_pos  = 0._r8
         337           0 :        do k = 1,pver
         338           0 :           precip(i,k) = cmfdqr(i,k) + nrain(i,k) - nevapr(i,k)
         339           0 :           total_rain = total_rain + precip(i,k)
         340           0 :           if ( precip(i,k) < 0._r8 ) precip(i,k) = 0._r8
         341           0 :           total_pos  = total_pos  + precip(i,k)
         342             :        end do
         343           0 :        if ( total_rain <= 0._r8 ) then
         344           0 :           precip(i,:) = 0._r8
         345             :        else
         346           0 :           do k = 1,pver
         347           0 :              precip(i,k) = precip(i,k) * total_rain/total_pos
         348             :           end do
         349             :        end if
         350             :     end do
         351             :
         352           0 :     do k = 1,pver
         353             :        !jfl       precip(:ncol,k) = cmfdqr(:ncol,k) + nrain(:ncol,k) - nevapr(:ncol,k)
         354           0 :        rain(:ncol,k)   = mass_air*precip(:ncol,k)*xhnm(:ncol,k) / mass_h2o
         355           0 :        xliq(:ncol,k)   = precip(:ncol,k) * delt * xhnm(:ncol,k) / avo*mass_air * m3_2_cm3
         356           0 :        if( spc_hno3_ndx > 0 ) then
         357           0 :           xhno3(:ncol,k)  = qin(:ncol,k,spc_hno3_ndx) * xhnm(:ncol,k)
         358             :        else
         359           0 :           xhno3(:ncol,k)  = 0._r8
         360             :        end if
         361           0 :        if( spc_h2o2_ndx > 0 ) then
         362           0 :           xh2o2(:ncol,k)  = qin(:ncol,k,spc_h2o2_ndx) * xhnm(:ncol,k)
         363             :        else
         364           0 :           xh2o2(:ncol,k)  = 0._r8
         365             :        end if
         366           0 :        if( spc_sogm_ndx > 0 ) then
         367           0 :           xsogm(:ncol,k)  = qin(:ncol,k,spc_sogm_ndx) * xhnm(:ncol,k)
         368             :        else
         369           0 :           xsogm(:ncol,k)  = 0._r8
         370             :        end if
         371           0 :        if( spc_sogi_ndx > 0 ) then
         372           0 :           xsogi(:ncol,k)  = qin(:ncol,k,spc_sogi_ndx) * xhnm(:ncol,k)
         373             :        else
         374           0 :           xsogi(:ncol,k)  = 0._r8
         375             :        end if
         376           0 :        if( spc_sogt_ndx > 0 ) then
         377           0 :           xsogt(:ncol,k)  = qin(:ncol,k,spc_sogt_ndx) * xhnm(:ncol,k)
         378             :        else
         379           0 :           xsogt(:ncol,k)  = 0._r8
         380             :        end if
         381           0 :        if( spc_sogb_ndx > 0 ) then
         382           0 :           xsogb(:ncol,k)  = qin(:ncol,k,spc_sogb_ndx) * xhnm(:ncol,k)
         383             :        else
         384           0 :           xsogb(:ncol,k)  = 0._r8
         385             :        end if
         386           0 :        if( spc_sogx_ndx > 0 ) then
         387           0 :           xsogx(:ncol,k)  = qin(:ncol,k,spc_sogx_ndx) * xhnm(:ncol,k)
         388             :        else
         389           0 :           xsogx(:ncol,k)  = 0._r8
         390             :        end if
         391           0 :        if( spc_so2_ndx > 0 ) then
         392           0 :           xso2(:ncol,k)  = qin(:ncol,k,spc_so2_ndx) * xhnm(:ncol,k)
         393             :        else
         394           0 :           xso2(:ncol,k)  = 0._r8
         395             :        end if
         396             :     end do
         397             :
         398           0 :     zsurf(:ncol) = m2km * phis(:ncol) * rga
         399           0 :     do k = ktop_all,pver-1
         400           0 :        delz(:ncol,k) = abs( (zmid(:ncol,k) - zmid(:ncol,k+1))*km2cm )
         401             :     end do
         402           0 :     delz(:ncol,pver) = abs( (zmid(:ncol,pver) - zsurf(:ncol) )*km2cm )
         403             :
         404             :     !-----------------------------------------------------------------
         405             :     !       ... part 0b,  for temperature dependent of henrys
         406             :     !                     xxhe1 = henry con for hno3
         407             :     !                     xxhe2 = henry con for h2o2
         408             :     !lwh 10/00 -- take henry''s law constants from brasseur et al. [1999],
         409             :     !             appendix j. for hno3, also consider dissociation to
         410             :     !             get effective henry''s law constant; equilibrium
         411             :     !             constant for dissociation from brasseur et al. [1999],
         412             :     !             appendix k. assume ph=5 (set as xph0 above).
         413             :     !             heff = h*k/[h+] for hno3 (complete dissociation)
         414             :     !             heff = h for h2o2 (no dissociation)
         415             :     !             heff = h * (1 + k/[h+]) (in general)
         416             :     !-----------------------------------------------------------------
         417           0 :     do k = ktop_all,pver
         418           0 :        work1(:ncol) = (t0 - tfld(:ncol,k))/(t0*tfld(:ncol,k))
         419             :        !-----------------------------------------------------------------
         420             :        !        ... effective henry''s law constants:
         421             :        !        hno3, h2o2, ch2o, ch3ooh, ch3coooh (brasseur et al., 1999)
         422             :        !       xooh, onitr, macrooh (j.-f. muller; brocheton, 1999)
         423             :        !       isopooh (equal to hno3, as for macrooh)
         424             :        !       ho2no2 (mozart-1)
         425             :        !       ch3cocho, hoch2cho (betterton and hoffman, environ. sci. technol., 1988)
         426             :        !       ch3cho (staudinger and roberts, crit. rev. sci. technol., 1996)
         427             :        !       mvk, macr (allen et al., environ. toxicol. chem., 1998)
         428             :        !       soag_bg(0-4), soag_ff_bb(0-4) (Hodzic et al., 2014)
         429             :        !-----------------------------------------------------------------
         430           0 :        xk0(:)             = 2.1e5_r8 *exp( 8700._r8*work1(:) )
         431           0 :        xhen_hno3(:,k)     = xk0(:) * ( 1._r8 + hno3_diss / xph0 )
         432           0 :        xhen_h2o2(:,k)     = 7.45e4_r8 * exp( 6620._r8 * work1(:) )
         433           0 :        xhen_ch2o(:,k)     = 6.3e3_r8 * exp( 6460._r8 * work1(:) )
         434           0 :        xhen_ch3ooh(:,k)   = 2.27e2_r8 * exp( 5610._r8 * work1(:) )
         435           0 :        xhen_ch3co3h(:,k)  = 4.73e2_r8 * exp( 6170._r8 * work1(:) )
         436           0 :        xhen_ch3cocho(:,k) = 3.70e3_r8 * exp( 7275._r8 * work1(:) )
         437           0 :        xhen_xooh(:,k)     = 90.5_r8 * exp( 5607._r8 * work1(:) )
         438           0 :        xhen_onitr(:,k)    = 7.51e3_r8 * exp( 6485._r8 * work1(:) )
         439           0 :        xhen_ho2no2(:,k)   = 2.e4_r8
         440           0 :        xhen_glyald(:,k)   = 4.1e4_r8 * exp( 4600._r8 * work1(:) )
         441           0 :        xhen_ch3cho(:,k)   = 1.4e1_r8 * exp( 5600._r8 * work1(:) )
         442           0 :        xhen_mvk(:,k)      = 21._r8 * exp( 7800._r8 * work1(:) )
         443           0 :        xhen_macr(:,k)     = 4.3_r8 * exp( 5300._r8 * work1(:) )
         444           0 :        xhen_ch3cooh(:,k)  = 4.1e3_r8 * exp( 6300._r8 * work1(:) )
         445           0 :        xhen_sog(:,k)      = 5.e5_r8 * exp (12._r8 * work1(:) )
         446             :        !
         447             :        ! calculation for NH3 using the parameters in drydep_tables.F90
         448             :        !
         449           0 :        xhen_nh3 (:,k)     = 1.e6_r8
         450           0 :        xhen_ch3cn(:,k)     = 50._r8 * exp( 4000._r8 * work1(:) )
         451           0 :        xhen_hcn(:,k)       = 12._r8 * exp( 5000._r8 * work1(:) )
         452           0 :        do i = 1, ncol
         453           0 :           so2_diss        = 1.23e-2_r8 * exp( 1960._r8 * work1(i) )
         454           0 :           xhen_so2(i,k)   = 1.23_r8 * exp( 3120._r8 * work1(i) ) * ( 1._r8 + so2_diss / xph0 )
         455             :        end do
         456             :        !
         457           0 :        tmp_hetrates(:,k,:) = 0._r8
         458             :     end do
         459             :
         460             :     !-----------------------------------------------------------------
         461             :     !       ... part 1, solve for high henry constant ( hno3, h2o2)
         462             :     !-----------------------------------------------------------------
         463           0 :     col_loop :  do i = 1,ncol
         464           0 :        xgas1(:) = xhno3(i,:)                     ! xgas will change during
         465           0 :        xgas2(:) = xh2o2(i,:)                     ! different levels wash
         466           0 :        xgas3(:) = xso2 (i,:)
         467           0 :        xgas4(:) = xsogm(i,:)
         468           0 :        xgas5(:) = xsogi(i,:)
         469           0 :        xgas6(:) = xsogt(i,:)
         470           0 :        xgas7(:) = xsogb(i,:)
         471           0 :        xgas8(:) = xsogx(i,:)
         472           0 :        level_loop1  : do kk = ktop(i),pver
         473           0 :           stay = 1._r8
         474           0 :           if( rain(i,kk) /= 0._r8 ) then            ! finding rain cloud
         475           0 :              all1 = 0._r8                           ! accumulation to justisfy saturation
         476           0 :              all2 = 0._r8
         477           0 :              all3 = 0._r8
         478           0 :              all4 = 0._r8
         479           0 :              all5 = 0._r8
         480           0 :              all6 = 0._r8
         481           0 :              all7 = 0._r8
         482           0 :              all8 = 0._r8
         483           0 :              stay = ((zmid(i,kk) - zsurf(i))*km2cm)/(xum*delt)
         484           0 :              stay = min( stay,1._r8 )
         485             :              !-----------------------------------------------------------------
         486             :              !       ... calculate the saturation concentration eqca
         487             :              !-----------------------------------------------------------------
         488           0 :              do k = kk,pver                      ! cal washout below cloud
         489           0 :                 xeqca1 =  xgas1(k) &
         490           0 :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_hno3(i,k)*const0*tfld(i,k))) &
         491           0 :                      *  xliq(i,kk)*avo2
         492             :                 xeqca2 =  xgas2(k) &
         493             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_h2o2(i,k)*const0*tfld(i,k))) &
         494           0 :                      *  xliq(i,kk)*avo2
         495             :                 xeqca3 =  xgas3(k) &
         496             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_so2( i,k)*const0*tfld(i,k))) &
         497           0 :                      *  xliq(i,kk)*avo2
         498             :                 xeqca4 =  xgas4(k) &
         499             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
         500           0 :                      *  xliq(i,kk)*avo2
         501             :                 xeqca5 =  xgas5(k) &
         502             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
         503           0 :                      *  xliq(i,kk)*avo2
         504             :                 xeqca6 =  xgas6(k) &
         505             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
         506           0 :                      *  xliq(i,kk)*avo2
         507             :                 xeqca7 =  xgas7(k) &
         508             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
         509           0 :                      *  xliq(i,kk)*avo2
         510             :                 xeqca8 =  xgas8(k) &
         511             :                      / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
         512           0 :                      *  xliq(i,kk)*avo2
         513             :
         514             :                 !-----------------------------------------------------------------
         515             :                 !       ... calculate ca; inside cloud concentration in #/cm3(air)
         516             :                 !-----------------------------------------------------------------
         517           0 :                 xca1 = geo_fac*xkgm*xgas1(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         518           0 :                 xca2 = geo_fac*xkgm*xgas2(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         519           0 :                 xca3 = geo_fac*xkgm*xgas3(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         520           0 :                 xca4 = geo_fac*xkgm*xgas4(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         521           0 :                 xca5 = geo_fac*xkgm*xgas5(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         522           0 :                 xca6 = geo_fac*xkgm*xgas6(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         523           0 :                 xca7 = geo_fac*xkgm*xgas7(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         524           0 :                 xca8 = geo_fac*xkgm*xgas8(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
         525             :
         526             :                 !-----------------------------------------------------------------
         527             :                 !       ... if is not saturated
         528             :                 !               hno3(gas)_new = hno3(gas)_old - hno3(h2o)
         529             :                 !           otherwise
         530             :                 !               hno3(gas)_new = hno3(gas)_old
         531             :                 !-----------------------------------------------------------------
         532           0 :                 all1 = all1 + xca1
         533           0 :                 all2 = all2 + xca2
         534           0 :                 if( all1 < xeqca1 ) then
         535           0 :                    xgas1(k) = max( xgas1(k) - xca1,0._r8 )
         536             :                 end if
         537           0 :                 if( all2 < xeqca2 ) then
         538           0 :                    xgas2(k) = max( xgas2(k) - xca2,0._r8 )
         539             :                 end if
         540           0 :                 all3 = all3 + xca3
         541           0 :                 if( all3 < xeqca3 ) then
         542           0 :                    xgas3(k) = max( xgas3(k) - xca3,0._r8 )
         543             :                 end if
         544           0 :                 all4 = all4 + xca4
         545           0 :                 all5 = all5 + xca5
         546           0 :                 all6 = all6 + xca6
         547           0 :                 all7 = all7 + xca7
         548           0 :                 all8 = all8 + xca8
         549           0 :                 if( all4 < xeqca4 ) then
         550           0 :                    xgas4(k) = max( xgas4(k) - xca4,0._r8 )
         551             :                 end if
         552           0 :                 if( all5 < xeqca5 ) then
         553           0 :                    xgas5(k) = max( xgas5(k) - xca5,0._r8 )
         554             :                 end if
         555           0 :                 if( all6 < xeqca6 ) then
         556           0 :                    xgas6(k) = max( xgas6(k) - xca6,0._r8 )
         557             :                 end if
         558           0 :                 if( all7 < xeqca7 ) then
         559           0 :                    xgas7(k) = max( xgas7(k) - xca7,0._r8 )
         560             :                 end if
         561           0 :                 if( all8 < xeqca8 ) then
         562           0 :                    xgas8(k) = max( xgas8(k) - xca8,0._r8 )
         563             :                 end if
         564             :              end do
         565             :           end if
         566             :           !-----------------------------------------------------------------
         567             :           !       ... calculate the lifetime of washout (second)
         568             :           !             after all layers washout
         569             :           !             the concentration of hno3 is reduced
         570             :           !             then the lifetime xtt is calculated by
         571             :           !
         572             :           !                  xtt = (xhno3(ini) - xgas1(new))/(dt*xhno3(ini))
         573             :           !                  where dt = passing time (s) in vertical
         574             :           !                             path below the cloud
         575             :           !                        dt = dz(cm)/um(cm/s)
         576             :           !-----------------------------------------------------------------
         577           0 :           xdtm = delz(i,kk) / xum                     ! the traveling time in each dz
         578           0 :           xxx1 = (xhno3(i,kk) - xgas1(kk))
         579           0 :           xxx2 = (xh2o2(i,kk) - xgas2(kk))
         580           0 :           if( xxx1 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         581           0 :              yhno3  = xhno3(i,kk)/xxx1 * xdtm
         582             :           else
         583             :              yhno3  = 1.e29_r8
         584             :           end if
         585           0 :           if( xxx2 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         586           0 :              yh2o2  = xh2o2(i,kk)/xxx2 * xdtm
         587             :           else
         588             :              yh2o2  = 1.e29_r8
         589             :           end if
         590           0 :           tmp_hetrates(i,kk,1) = max( 1._r8 / yh2o2,0._r8 ) * stay
         591           0 :           tmp_hetrates(i,kk,2) = max( 1._r8 / yhno3,0._r8 ) * stay
         592           0 :           xxx3 = (xso2( i,kk) - xgas3(kk))
         593           0 :           if( xxx3 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         594           0 :              yso2  = xso2( i,kk)/xxx3 * xdtm
         595             :           else
         596             :              yso2  = 1.e29_r8
         597             :           end if
         598           0 :           tmp_hetrates(i,kk,3) = max( 1._r8 / yso2, 0._r8 ) * stay
         599           0 :           xxx4 = (xsogm(i,kk) - xgas4(kk))
         600           0 :           xxx5 = (xsogi(i,kk) - xgas5(kk))
         601           0 :           xxx6 = (xsogt(i,kk) - xgas6(kk))
         602           0 :           xxx7 = (xsogb(i,kk) - xgas7(kk))
         603           0 :           xxx8 = (xsogx(i,kk) - xgas8(kk))
         604           0 :           if( xxx4 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         605           0 :              ysogm  = xsogm(i,kk)/xxx4 * xdtm
         606             :           else
         607             :              ysogm  = 1.e29_r8
         608             :           end if
         609           0 :           if( xxx5 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         610           0 :              ysogi  = xsogi(i,kk)/xxx5 * xdtm
         611             :           else
         612             :              ysogi  = 1.e29_r8
         613             :           end if
         614           0 :           if( xxx6 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         615           0 :              ysogt  = xsogt(i,kk)/xxx6 * xdtm
         616             :           else
         617             :              ysogt  = 1.e29_r8
         618             :           end if
         619           0 :           if( xxx7 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         620           0 :              ysogb  = xsogb(i,kk)/xxx7 * xdtm
         621             :           else
         622             :              ysogb  = 1.e29_r8
         623             :           end if
         624           0 :           if( xxx8 /= 0._r8 ) then                       ! if no washout lifetime = 1.e29
         625           0 :              ysogx  = xsogx(i,kk)/xxx8 * xdtm
         626             :           else
         627             :              ysogx  = 1.e29_r8
         628             :           end if
         629           0 :           tmp_hetrates(i,kk,4) = max( 1._r8 / ysogm,0._r8 ) * stay
         630           0 :           tmp_hetrates(i,kk,5) = max( 1._r8 / ysogi,0._r8 ) * stay
         631           0 :           tmp_hetrates(i,kk,6) = max( 1._r8 / ysogt,0._r8 ) * stay
         632           0 :           tmp_hetrates(i,kk,7) = max( 1._r8 / ysogb,0._r8 ) * stay
         633           0 :           tmp_hetrates(i,kk,8) = max( 1._r8 / ysogx,0._r8 ) * stay
         634             :        end do level_loop1
         635             :     end do col_loop
         636             :
         637             :     !-----------------------------------------------------------------
         638             :     !       ... part 2, in-cloud solve for low henry constant
         639             :     !                   hno3 and h2o2 have both in and under cloud
         640             :     !-----------------------------------------------------------------
         641           0 :     level_loop2 : do k = ktop_all,pver
         642           0 :        Column_loop2 : do i=1,ncol
         643           0 :           if ( rain(i,k) <= 0._r8 ) then
         644           0 :              het_rates(i,k,:) =  0._r8
         645             :              cycle
         646             :           endif
         647             :
         648           0 :           work1(i) = avo2 * xliq(i,k)
         649           0 :           work2(i) = const0 * tfld(i,k)
         650             :           work3(i) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch2o(i,k)*work2(i)))),0._r8 ) &
         651           0 :                * satf_ch2o
         652           0 :           if( ch2o_ndx > 0 ) then
         653           0 :              het_rates(i,k,ch2o_ndx)  = work3(i)
         654             :           end if
         655           0 :           if( isopno3_ndx > 0 ) then
         656           0 :              het_rates(i,k,isopno3_ndx) = work3(i)
         657             :           end if
         658           0 :           if( xisopno3_ndx > 0 ) then
         659           0 :              het_rates(i,k,xisopno3_ndx) = work3(i)
         660             :           end if
         661           0 :           if( hyac_ndx > 0 ) then
         662           0 :              het_rates(i,k,hyac_ndx) = work3(i)
         663             :           end if
         664           0 :           if( hydrald_ndx > 0 ) then
         665           0 :              het_rates(i,k,hydrald_ndx) = work3(i)
         666             :           end if
         667             :
         668           0 :           work3(i) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3ooh(i,k)*work2(i)))),0._r8 )
         669           0 :           if( ch3ooh_ndx > 0 ) then
         670           0 :              het_rates(i,k,ch3ooh_ndx)  = work3(i)
         671             :           end if
         672           0 :           if( pooh_ndx > 0 ) then
         673           0 :              het_rates(i,k,pooh_ndx)  = work3(i)
         674             :           end if
         675           0 :           if( c2h5ooh_ndx > 0 ) then
         676           0 :              het_rates(i,k,c2h5ooh_ndx) = work3(i)
         677             :           end if
         678           0 :           if( c3h7ooh_ndx > 0 ) then
         679           0 :              het_rates(i,k,c3h7ooh_ndx) = work3(i)
         680             :           end if
         681           0 :           if( rooh_ndx > 0 ) then
         682           0 :              het_rates(i,k,rooh_ndx) = work3(i)
         683             :           end if
         684           0 :           if( ch3oh_ndx > 0 ) then
         685           0 :              het_rates(i,k,ch3oh_ndx) = work3(i)
         686             :           end if
         687           0 :           if( c2h5oh_ndx > 0 ) then
         688           0 :              het_rates(i,k,c2h5oh_ndx) = work3(i)
         689             :           end if
         690           0 :           if( alkooh_ndx  > 0 ) then
         691           0 :              het_rates(i,k,alkooh_ndx) = work3(i)
         692             :           end if
         693           0 :           if( mekooh_ndx  > 0 ) then
         694           0 :              het_rates(i,k,mekooh_ndx) = work3(i)
         695             :           end if
         696           0 :           if( tolooh_ndx  > 0 ) then
         697           0 :              het_rates(i,k,tolooh_ndx) = work3(i)
         698             :           end if
         699           0 :           if( terpooh_ndx > 0 ) then
         700           0 :              het_rates(i,k,terpooh_ndx) = work3(i)
         701             :           end if
         702             :
         703           0 :           if( ch3coooh_ndx > 0 ) then
         704           0 :              het_rates(i,k,ch3coooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3co3h(i,k)*work2(i)))),0._r8 )
         705             :           end if
         706           0 :           if( ho2no2_ndx > 0 ) then
         707           0 :              het_rates(i,k,ho2no2_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ho2no2(i,k)*work2(i)))),0._r8 )
         708             :           end if
         709           0 :           if( xho2no2_ndx > 0 ) then
         710           0 :              het_rates(i,k,xho2no2_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ho2no2(i,k)*work2(i)))),0._r8 )
         711             :           end if
         712           0 :           if( ch3cocho_ndx > 0 ) then
         713           0 :              het_rates(i,k,ch3cocho_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cocho(i,k)*work2(i)))),0._r8 )
         714             :           end if
         715           0 :           if( xooh_ndx > 0 ) then
         716           0 :              het_rates(i,k,xooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_xooh(i,k)*work2(i)))),0._r8 )
         717             :           end if
         718           0 :           if( onitr_ndx > 0 ) then
         719           0 :              het_rates(i,k,onitr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_onitr(i,k)*work2(i)))),0._r8 )
         720             :           end if
         721           0 :           if( xonitr_ndx > 0 ) then
         722           0 :              het_rates(i,k,xonitr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_onitr(i,k)*work2(i)))),0._r8 )
         723             :           end if
         724           0 :           if( glyald_ndx > 0 ) then
         725           0 :              het_rates(i,k,glyald_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_glyald(i,k)*work2(i)))),0._r8 )
         726             :           end if
         727           0 :           if( ch3cho_ndx > 0 ) then
         728           0 :              het_rates(i,k,ch3cho_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cho(i,k)*work2(i)))),0._r8 )
         729             :           end if
         730           0 :           if( mvk_ndx > 0 ) then
         731           0 :              het_rates(i,k,mvk_ndx)  = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_mvk(i,k)*work2(i)))),0._r8 )
         732             :           end if
         733           0 :           if( macr_ndx > 0 ) then
         734           0 :              het_rates(i,k,macr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_macr(i,k)*work2(i)))),0._r8 )
         735             :           end if
         736           0 :           if( h2o2_ndx > 0 ) then
         737           0 :              work3(i) = satf_h2o2 * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_h2o2(i,k)*work2(i)))),0._r8 )
         738           0 :              het_rates(i,k,h2o2_ndx) =  work3(i) + tmp_hetrates(i,k,1)
         739             :           end if
         740           0 :           if ( prog_modal_aero .and. so2_ndx>0 .and. h2o2_ndx>0 ) then
         741           0 :              het_rates(i,k,so2_ndx) = het_rates(i,k,h2o2_ndx)
         742           0 :           elseif( so2_ndx > 0 ) then
         743           0 :              work3(i) = satf_so2 * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_so2( i,k)*work2(i)))),0._r8 )
         744           0 :              het_rates(i,k,so2_ndx ) =  work3(i) + tmp_hetrates(i,k,3)
         745             :           endif
         746             : !
         747           0 :           work3(i) = satf_sog * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_sog(i,k)*work2(i)))),0._r8 )
         748           0 :           if( sogm_ndx > 0 ) then
         749           0 :              het_rates(i,k,sogm_ndx) =  work3(i) + tmp_hetrates(i,k,4)
         750             :           end if
         751           0 :           if( sogi_ndx > 0 ) then
         752           0 :              het_rates(i,k,sogi_ndx) =  work3(i) + tmp_hetrates(i,k,5)
         753             :           end if
         754           0 :           if( sogt_ndx > 0 ) then
         755           0 :              het_rates(i,k,sogt_ndx) =  work3(i) + tmp_hetrates(i,k,6)
         756             :           end if
         757           0 :           if( sogb_ndx > 0 ) then
         758           0 :              het_rates(i,k,sogb_ndx) =  work3(i) + tmp_hetrates(i,k,7)
         759             :           end if
         760           0 :           if( sogx_ndx > 0 ) then
         761           0 :              het_rates(i,k,sogx_ndx) =  work3(i) + tmp_hetrates(i,k,8)
         762             :           end if
         763             : !
         764             :           work3(i) = tmp_hetrates(i,k,2) + satf_hno3 * &
         765           0 :                max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_hno3(i,k)*work2(i)))),0._r8 )
         766           0 :           tmp0_rates(i)   = work3(i)
         767           0 :           tmp1_rates(i)   = .2_r8*work3(i)
         768           0 :           if( hno3_ndx > 0 ) then
         769           0 :              het_rates(i,k,hno3_ndx) = work3(i)
         770             :           end if
         771           0 :           if( xhno3_ndx > 0 ) then
         772           0 :              het_rates(i,k,xhno3_ndx) = work3(i)
         773             :           end if
         774           0 :           if( onit_ndx > 0 ) then
         775           0 :              het_rates(i,k,onit_ndx) = work3(i)
         776             :           end if
         777           0 :           if( xonit_ndx > 0 ) then
         778           0 :              het_rates(i,k,xonit_ndx) = work3(i)
         779             :           end if
         780           0 :           if( Pb_ndx > 0 ) then
         781           0 :              het_rates(i,k,Pb_ndx) = work3(i)
         782             :           end if
         783           0 :           if( macrooh_ndx > 0 ) then
         784           0 :              het_rates(i,k,macrooh_ndx) = work3(i)
         785             :           end if
         786           0 :           if( isopooh_ndx > 0 ) then
         787           0 :              het_rates(i,k,isopooh_ndx) = work3(i)
         788             :           end if
         789             :
         790           0 :           if( clono2_ndx > 0 ) then
         791           0 :              het_rates(i,k, clono2_ndx) = work3(i)
         792             :           end if
         793           0 :           if( brono2_ndx > 0 ) then
         794           0 :              het_rates(i,k, brono2_ndx) = work3(i)
         795             :           end if
         796           0 :           if( hcl_ndx > 0 ) then
         797           0 :              het_rates(i,k, hcl_ndx) = work3(i)
         798             :           end if
         799           0 :           if( n2o5_ndx > 0 ) then
         800           0 :              het_rates(i,k, n2o5_ndx) = work3(i)
         801             :           end if
         802           0 :           if( hocl_ndx > 0 ) then
         803           0 :              het_rates(i,k, hocl_ndx) = work3(i)
         804             :           end if
         805           0 :           if( hobr_ndx > 0 ) then
         806           0 :              het_rates(i,k, hobr_ndx) = work3(i)
         807             :           end if
         808           0 :           if( hbr_ndx > 0 ) then
         809           0 :              het_rates(i,k, hbr_ndx) = work3(i)
         810             :           end if
         811             :
         812           0 :           if( soa_ndx > 0 ) then
         813           0 :              het_rates(i,k,soa_ndx) = tmp1_rates(i)
         814             :           end if
         815           0 :           if( oc2_ndx > 0 ) then
         816           0 :              het_rates(i,k,oc2_ndx) = tmp1_rates(i)
         817             :           end if
         818           0 :           if( cb2_ndx > 0 ) then
         819           0 :              het_rates(i,k,cb2_ndx) = tmp1_rates(i)
         820             :           end if
         821           0 :           if( so4_ndx > 0 ) then
         822           0 :              het_rates(i,k,so4_ndx) = tmp1_rates(i)
         823             :           end if
         824           0 :           if( sa1_ndx > 0 ) then
         825           0 :              het_rates(i,k,sa1_ndx) = tmp1_rates(i)
         826             :           end if
         827           0 :           if( sa2_ndx > 0 ) then
         828           0 :              het_rates(i,k,sa2_ndx) = tmp1_rates(i)
         829             :           end if
         830           0 :           if( sa3_ndx > 0 ) then
         831           0 :              het_rates(i,k,sa3_ndx) = tmp1_rates(i)
         832             :           end if
         833           0 :           if( sa4_ndx > 0 ) then
         834           0 :              het_rates(i,k,sa4_ndx) = tmp1_rates(i)
         835             :           end if
         836             :
         837           0 :           if( h2so4_ndx > 0 ) then
         838           0 :              het_rates(i,k,h2so4_ndx) = tmp0_rates(i)
         839             :           end if
         840           0 :           if( nh4_ndx > 0 ) then
         841           0 :              het_rates(i,k,nh4_ndx) = tmp0_rates(i)
         842             :           end if
         843           0 :           if( nh4no3_ndx > 0 ) then
         844           0 :              het_rates(i,k,nh4no3_ndx ) = tmp0_rates(i)
         845             :           end if
         846           0 :           if( nh3_ndx > 0 ) then
         847           0 :              het_rates(i,k,nh3_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_nh3(i,k)*work2(i)))),0._r8 )
         848             :           end if
         849             :
         850           0 :           if( ch3cooh_ndx > 0 ) then
         851           0 :              het_rates(i,k,ch3cooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cooh(i,k)*work2(i)))),0._r8 )
         852             :           end if
         853           0 :           if( hcooh_ndx > 0 ) then
         854           0 :              het_rates(i,k,hcooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cooh(i,k)*work2(i)))),0._r8 )
         855             :           endif
         856           0 :           if ( hcn_ndx > 0 ) then
         857           0 :              het_rates(i,k,hcn_ndx     ) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_hcn(i,k)*work2(i)))),0._r8 )
         858             :           endif
         859           0 :           if ( ch3cn_ndx > 0 ) then
         860           0 :              het_rates(i,k,ch3cn_ndx   ) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cn(i,k)*work2(i)))),0._r8 )
         861             :           endif
         862             :
         863             :        end do Column_loop2
         864             :     end do level_loop2
         865             :
         866             :     !-----------------------------------------------------------------
         867             :     !   ... Set rates above tropopause = 0.
         868             :     !-----------------------------------------------------------------
         869           0 :     do mm = 1,gas_wetdep_cnt
         870           0 :        m = wetdep_map(mm)
         871           0 :        do i = 1,ncol
         872           0 :           do k = 1,ktop(i)
         873           0 :              het_rates(i,k,m) = 0._r8
         874             :           end do
         875             :        end do
         876           0 :        if ( any( het_rates(:ncol,:,m) == MISSING) ) then
         877           0 :           write(hetratestrg,'(I3)') m
         878           0 :           call endrun('sethet: het_rates (wet dep) not set for het reaction number : '//hetratestrg)
         879             :        endif
         880             :     end do
         881             :
         882           0 :   end subroutine sethet
         883             :
         884             : end module mo_sethet

Generated by: LCOV version 1.14