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 cam_logfile, only: iulog use gas_wetdep_opts, only: gas_wetdep_cnt, gas_wetdep_method, gas_wetdep_list 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 :: 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 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 abortutils, only : endrun integer :: k, m do_wetdep = gas_wetdep_cnt>0 .and. gas_wetdep_method=='MOZ' if ( .not. do_wetdep) return 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' ) if (masterproc) then write(iulog,*) 'sethet_inti: new ndx ',so2_ndx,soa_ndx,so4_ndx,cb2_ndx,oc2_ndx, & nh3_ndx,nh4no3_ndx,sa1_ndx,sa2_ndx,sa3_ndx,sa4_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, modal ) !----------------------------------------------------------------------- ! ... compute rainout loss rates (1/s) !----------------------------------------------------------------------- use shr_kind_mod, only : r8 => shr_kind_r8 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 abortutils, only : endrun 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) logical, intent(in) :: modal real(r8), intent(out) :: het_rates(ncol,pver,gas_pcnst) ! rainout loss rates !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- real(r8), parameter :: boltz = 1.38044e-16_r8 ! erg/K real(r8), parameter :: avo = 6.023e23_r8 ! 1/mol 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 * 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( 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 ( modal ) then if( so2_ndx > 0 ) then het_rates(i,k,so2_ndx) = het_rates(i,k,h2o2_ndx) end if else if( 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) end if 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