module MO_SETSOX use shr_kind_mod, only : r8 => shr_kind_r8 use cam_logfile, only : iulog private public :: sox_inti, setsox public :: has_sox save logical :: inv_o3 integer :: id_msa integer :: id_so2, id_nh3, id_hno3, id_h2o2, id_o3, id_ho2 integer :: id_so4, id_h2so4 logical :: has_sox = .true. logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ox, inv_nh4no3, inv_ho2 logical :: cloud_borne = .false. logical :: modal_aerosols = .false. contains !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine sox_inti !----------------------------------------------------------------------- ! ... initialize the hetero sox routine !----------------------------------------------------------------------- use mo_chem_utls, only : get_spc_ndx, get_inv_ndx use spmd_utils, only : masterproc use phys_control, only : phys_getopts use sox_cldaero_mod, only : sox_cldaero_init implicit none call phys_getopts( & prog_modal_aero_out=modal_aerosols ) cloud_borne = modal_aerosols !----------------------------------------------------------------- ! ... get species indicies !----------------------------------------------------------------- if (cloud_borne) then id_h2so4 = get_spc_ndx( 'H2SO4' ) else id_so4 = get_spc_ndx( 'SO4' ) endif id_msa = get_spc_ndx( 'MSA' ) inv_so2 = .false. id_so2 = get_inv_ndx( 'SO2' ) inv_so2 = id_so2 > 0 if ( .not. inv_so2 ) then id_so2 = get_spc_ndx( 'SO2' ) endif inv_NH3 = .false. id_NH3 = get_inv_ndx( 'NH3' ) inv_NH3 = id_NH3 > 0 if ( .not. inv_NH3 ) then id_NH3 = get_spc_ndx( 'NH3' ) endif inv_HNO3 = .false. id_HNO3 = get_inv_ndx( 'HNO3' ) inv_HNO3 = id_hno3 > 0 if ( .not. inv_HNO3 ) then id_HNO3 = get_spc_ndx( 'HNO3' ) endif inv_H2O2 = .false. id_H2O2 = get_inv_ndx( 'H2O2' ) inv_H2O2 = id_H2O2 > 0 if ( .not. inv_H2O2 ) then id_H2O2 = get_spc_ndx( 'H2O2' ) endif inv_HO2 = .false. id_HO2 = get_inv_ndx( 'HO2' ) inv_HO2 = id_HO2 > 0 if ( .not. inv_HO2 ) then id_HO2 = get_spc_ndx( 'HO2' ) endif inv_o3 = get_inv_ndx( 'O3' ) > 0 if (inv_o3) then id_o3 = get_inv_ndx( 'O3' ) else id_o3 = get_spc_ndx( 'O3' ) endif inv_ho2 = get_inv_ndx( 'HO2' ) > 0 if (inv_ho2) then id_ho2 = get_inv_ndx( 'HO2' ) else id_ho2 = get_spc_ndx( 'HO2' ) endif has_sox = (id_so2>0) .and. (id_h2o2>0) .and. (id_o3>0) .and. (id_ho2>0) if (cloud_borne) then has_sox = has_sox .and. (id_h2so4>0) else has_sox = has_sox .and. (id_so4>0) .and. (id_nh3>0) endif if (masterproc) then write(iulog,*) 'sox_inti: has_sox = ',has_sox endif if( has_sox ) then if (masterproc) then write(iulog,*) '-----------------------------------------' write(iulog,*) 'mozart will do sox aerosols' write(iulog,*) '-----------------------------------------' endif else return end if call sox_cldaero_init() end subroutine sox_inti !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine SETSOX( & ncol, & lchnk, & loffset,& dtime, & press, & pdel, & tfld, & mbar, & lwc, & cldfrc, & cldnum, & xhnm, & invariants, & qcw, & qin, & xphlwc, & aqso4, & aqh2so4,& aqso4_h2o2, & aqso4_o3, & yph_in & ) !----------------------------------------------------------------------- ! ... Compute heterogeneous reactions of SOX ! ! (0) using initial PH to calculate PH ! (a) HENRYs law constants ! (b) PARTIONING ! (c) PH values ! ! (1) using new PH to repeat ! (a) HENRYs law constants ! (b) PARTIONING ! (c) REACTION rates ! (d) PREDICTION !----------------------------------------------------------------------- ! use ppgrid, only : pcols, pver use chem_mods, only : gas_pcnst, nfs use chem_mods, only : adv_mass use physconst, only : mwdry, gravit use mo_constants, only : pi use sox_cldaero_mod, only : sox_cldaero_update, sox_cldaero_create_obj, sox_cldaero_destroy_obj use cldaero_mod, only : cldaero_conc_t ! implicit none ! !----------------------------------------------------------------------- ! ... Dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: ncol ! num of columns in chunk integer, intent(in) :: lchnk ! chunk id integer, intent(in) :: loffset ! offset of chem tracers in the advected tracers array real(r8), intent(in) :: dtime ! time step (sec) real(r8), intent(in) :: press(:,:) ! midpoint pressure ( Pa ) real(r8), intent(in) :: pdel(:,:) ! pressure thickness of levels (Pa) real(r8), intent(in) :: tfld(:,:) ! temperature real(r8), intent(in) :: mbar(:,:) ! mean wet atmospheric mass ( amu ) real(r8), target, intent(in) :: lwc(:,:) ! cloud liquid water content (kg/kg) real(r8), target, intent(in) :: cldfrc(:,:) ! cloud fraction real(r8), intent(in) :: cldnum(:,:) ! droplet number concentration (#/kg) real(r8), intent(in) :: xhnm(:,:) ! total atms density ( /cm**3) real(r8), intent(in) :: invariants(:,:,:) real(r8), target, intent(inout) :: qcw(:,:,:) ! cloud-borne aerosol (vmr) real(r8), intent(inout) :: qin(:,:,:) ! transported species ( vmr ) real(r8), intent(out) :: xphlwc(:,:) ! pH value multiplied by lwc real(r8), intent(out) :: aqso4(:,:) ! aqueous phase chemistry real(r8), intent(out) :: aqh2so4(:,:) ! aqueous phase chemistry real(r8), intent(out) :: aqso4_h2o2(:) ! SO4 aqueous phase chemistry due to H2O2 (kg/m2) real(r8), intent(out) :: aqso4_o3(:) ! SO4 aqueous phase chemistry due to O3 (kg/m2) real(r8), intent(in), optional :: yph_in ! ph value !----------------------------------------------------------------------- ! ... Local variables ! ! xhno3 ... in mixing ratio !----------------------------------------------------------------------- integer, parameter :: itermax = 20 real(r8), parameter :: ph0 = 5.0_r8 ! INITIAL PH VALUES real(r8), parameter :: const0 = 1.e3_r8/6.023e23_r8 real(r8), parameter :: xa0 = 11._r8 real(r8), parameter :: xb0 = -.1_r8 real(r8), parameter :: xa1 = 1.053_r8 real(r8), parameter :: xb1 = -4.368_r8 real(r8), parameter :: xa2 = 1.016_r8 real(r8), parameter :: xb2 = -2.54_r8 real(r8), parameter :: xa3 = .816e-32_r8 real(r8), parameter :: xb3 = .259_r8 real(r8), parameter :: kh0 = 9.e3_r8 ! HO2(g) -> Ho2(a) real(r8), parameter :: kh1 = 2.05e-5_r8 ! HO2(a) -> H+ + O2- real(r8), parameter :: kh2 = 8.6e5_r8 ! HO2(a) + ho2(a) -> h2o2(a) + o2 real(r8), parameter :: kh3 = 1.e8_r8 ! HO2(a) + o2- -> h2o2(a) + o2 real(r8), parameter :: Ra = 8314._r8/101325._r8 ! universal constant (atm)/(M-K) real(r8), parameter :: xkw = 1.e-14_r8 ! water acidity ! real(r8) :: xdelso4hp(ncol,pver) integer :: k, i, iter, file real(r8) :: wrk, delta real(r8) :: xph0, aden, xk, xe, x2 real(r8) :: tz, xl, px, qz, pz, es, qs, patm real(r8) :: Eso2, Eso4, Ehno3, Eco2, Eh2o, Enh3 real(r8) :: so2g, h2o2g, co2g, o3g real(r8) :: hno3a, nh3a, so2a, h2o2a, co2a, o3a real(r8) :: rah2o2, rao3, pso4, ccc real(r8) :: cnh3, chno3, com, com1, com2, xra real(r8) :: hno3g(ncol,pver), nh3g(ncol,pver) ! !----------------------------------------------------------------------- ! for Ho2(g) -> H2o2(a) formation ! schwartz JGR, 1984, 11589 !----------------------------------------------------------------------- real(r8) :: kh4 ! kh2+kh3 real(r8) :: xam ! air density /cm3 real(r8) :: ho2s ! ho2s = ho2(a)+o2- real(r8) :: r1h2o2 ! prod(h2o2) by ho2 in mole/L(w)/s real(r8) :: r2h2o2 ! prod(h2o2) by ho2 in mix/s real(r8), dimension(ncol,pver) :: & xhno3, xh2o2, xso2, xso4, xno3, & xnh3, xnh4, xo3, & cfact, & xph, xho2, & xh2so4, xmsa, xso4_init, & hehno3, & ! henry law const for hno3 heh2o2, & ! henry law const for h2o2 heso2, & ! henry law const for so2 henh3, & ! henry law const for nh3 heo3 !!, & ! henry law const for o3 real(r8) :: patm_x real(r8), dimension(ncol) :: work1 logical :: converged real(r8), pointer :: xso4c(:,:) real(r8), pointer :: xnh4c(:,:) real(r8), pointer :: xno3c(:,:) type(cldaero_conc_t), pointer :: cldconc real(r8) :: fact1_hno3, fact2_hno3, fact3_hno3 real(r8) :: fact1_so2, fact2_so2, fact3_so2, fact4_so2 real(r8) :: fact1_nh3, fact2_nh3, fact3_nh3 real(r8) :: tmp_hp, tmp_hso3, tmp_hco3, tmp_nh4, tmp_no3 real(r8) :: tmp_oh, tmp_so3, tmp_so4 real(r8) :: tmp_neg, tmp_pos real(r8) :: yph, yph_lo, yph_hi real(r8) :: ynetpos, ynetpos_lo, ynetpos_hi !----------------------------------------------------------------- ! ... NOTE: The press array is in pascals and must be ! mutiplied by 10 to yield dynes/cm**2. !----------------------------------------------------------------- !================================================================== ! ... First set the PH !================================================================== ! ... Initial values ! The values of so2, so4 are after (1) SLT, and CHEM !----------------------------------------------------------------- xph0 = 10._r8**(-ph0) ! initial PH value do k = 1,pver cfact(:,k) = xhnm(:,k) & ! /cm3(a) * 1.e6_r8 & ! /m3(a) * 1.38e-23_r8/287._r8 & ! Kg(a)/m3(a) * 1.e-3_r8 ! Kg(a)/L(a) end do cldconc => sox_cldaero_create_obj( cldfrc,qcw,lwc, cfact, ncol, loffset ) xso4c => cldconc%so4c xnh4c => cldconc%nh4c xno3c => cldconc%no3c xso4(:,:) = 0._r8 xno3(:,:) = 0._r8 xnh4(:,:) = 0._r8 do k = 1,pver xph(:,k) = xph0 ! initial PH value if ( inv_so2 ) then xso2 (:,k) = invariants(:,k,id_so2)/xhnm(:,k) ! mixing ratio else xso2 (:,k) = qin(:,k,id_so2) ! mixing ratio endif if (id_hno3 > 0) then xhno3(:,k) = qin(:,k,id_hno3) else xhno3(:,k) = 0.0_r8 endif if ( inv_h2o2 ) then xh2o2 (:,k) = invariants(:,k,id_h2o2)/xhnm(:,k) ! mixing ratio else xh2o2 (:,k) = qin(:,k,id_h2o2) ! mixing ratio endif if (id_nh3 > 0) then xnh3 (:,k) = qin(:,k,id_nh3) else xnh3 (:,k) = 0.0_r8 endif if ( inv_o3 ) then xo3 (:,k) = invariants(:,k,id_o3)/xhnm(:,k) ! mixing ratio else xo3 (:,k) = qin(:,k,id_o3) ! mixing ratio endif if ( inv_ho2 ) then xho2 (:,k) = invariants(:,k,id_ho2)/xhnm(:,k)! mixing ratio else xho2 (:,k) = qin(:,k,id_ho2) ! mixing ratio endif if (cloud_borne) then xh2so4(:,k) = qin(:,k,id_h2so4) else xso4 (:,k) = qin(:,k,id_so4) ! mixing ratio endif if (id_msa > 0) xmsa (:,k) = qin(:,k,id_msa) end do !----------------------------------------------------------------- ! ... Temperature dependent Henry constants !----------------------------------------------------------------- ver_loop0: do k = 1,pver !! pver loop for STEP 0 col_loop0: do i = 1,ncol if (cloud_borne .and. cldfrc(i,k)>0._r8) then xso4(i,k) = xso4c(i,k) / cldfrc(i,k) xnh4(i,k) = xnh4c(i,k) / cldfrc(i,k) xno3(i,k) = xno3c(i,k) / cldfrc(i,k) endif xl = cldconc%xlwc(i,k) if( xl >= 1.e-8_r8 ) then work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8 !----------------------------------------------------------------- ! 21-mar-2011 changes by rce ! ph calculation now uses bisection method to solve the electro-neutrality equation ! 3-mode aerosols (where so4 is assumed to be nh4hso4) ! old code set xnh4c = so4c ! new code sets xnh4c = 0, then uses a -1 charge (instead of -2) ! for so4 when solving the electro-neutrality equation !----------------------------------------------------------------- !----------------------------------------------------------------- ! calculations done before iterating !----------------------------------------------------------------- !----------------------------------------------------------------- pz = .01_r8*press(i,k) !! pressure in mb tz = tfld(i,k) patm = pz/1013._r8 xam = press(i,k)/(1.38e-23_r8*tz) !air density /M3 !----------------------------------------------------------------- ! ... hno3 !----------------------------------------------------------------- ! previous code ! hehno3(i,k) = xk*(1._r8 + xe/xph(i,k)) ! px = hehno3(i,k) * Ra * tz * xl ! hno3g = xhno3(i,k)/(1._r8 + px) ! Ehno3 = xk*xe*hno3g *patm ! equivalent new code ! hehno3 = xk + xk*xe/hplus ! hno3g = xhno3/(1 + px) ! = xhno3/(1 + hehno3*ra*tz*xl) ! = xhno3/(1 + xk*ra*tz*xl*(1 + xe/hplus) ! ehno3 = hno3g*xk*xe*patm ! = xk*xe*patm*xhno3/(1 + xk*ra*tz*xl*(1 + xe/hplus) ! = ( fact1_hno3 )/(1 + fact2_hno3 *(1 + fact3_hno3/hplus) ! [hno3-] = ehno3/hplus xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) ) xe = 15.4_r8 fact1_hno3 = xk*xe*patm*xhno3(i,k) fact2_hno3 = xk*ra*tz*xl fact3_hno3 = xe !----------------------------------------------------------------- ! ... so2 !----------------------------------------------------------------- ! previous code ! heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k))) ! px = heso2(i,k) * Ra * tz * xl ! so2g = xso2(i,k)/(1._r8+ px) ! Eso2 = xk*xe*so2g *patm ! equivalent new code ! heso2 = xk + xk*xe/hplus * xk*xe*x2/hplus**2 ! so2g = xso2/(1 + px) ! = xso2/(1 + heso2*ra*tz*xl) ! = xso2/(1 + xk*ra*tz*xl*(1 + (xe/hplus)*(1 + x2/hplus)) ! eso2 = so2g*xk*xe*patm ! = xk*xe*patm*xso2/(1 + xk*ra*tz*xl*(1 + (xe/hplus)*(1 + x2/hplus)) ! = ( fact1_so2 )/(1 + fact2_so2 *(1 + (fact3_so2/hplus)*(1 + fact4_so2/hplus) ! [hso3-] + 2*[so3--] = (eso2/hplus)*(1 + 2*x2/hplus) xk = 1.23_r8 *EXP( 3120._r8*work1(i) ) xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) ) x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) ) fact1_so2 = xk*xe*patm*xso2(i,k) fact2_so2 = xk*ra*tz*xl fact3_so2 = xe fact4_so2 = x2 !----------------------------------------------------------------- ! ... nh3 !----------------------------------------------------------------- ! previous code ! henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw) ! px = henh3(i,k) * Ra * tz * xl ! nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) ! Enh3 = xk*xe*nh3g/xkw *patm ! equivalent new code ! henh3 = xk + xk*xe*hplus/xkw ! nh3g = xnh34/(1 + px) ! = xnh34/(1 + henh3*ra*tz*xl) ! = xnh34/(1 + xk*ra*tz*xl*(1 + xe*hplus/xkw) ! enh3 = nh3g*xk*xe*patm/xkw ! = ((xk*xe*patm/xkw)*xnh34)/(1 + xk*ra*tz*xl*(1 + xe*hplus/xkw) ! = ( fact1_nh3 )/(1 + fact2_nh3 *(1 + fact3_nh3*hplus) ! [nh4+] = enh3*hplus xk = 58._r8 *EXP( 4085._r8*work1(i) ) xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) ) fact1_nh3 = (xk*xe*patm/xkw)*(xnh3(i,k)+xnh4(i,k)) fact2_nh3 = xk*ra*tz*xl fact3_nh3 = xe/xkw !----------------------------------------------------------------- ! ... h2o effects !----------------------------------------------------------------- Eh2o = xkw !----------------------------------------------------------------- ! ... co2 effects !----------------------------------------------------------------- co2g = 330.e-6_r8 !330 ppm = 330.e-6 atm xk = 3.1e-2_r8*EXP( 2423._r8*work1(i) ) xe = 4.3e-7_r8*EXP(-913._r8 *work1(i) ) Eco2 = xk*xe*co2g *patm !----------------------------------------------------------------- ! ... so4 effect !----------------------------------------------------------------- Eso4 = xso4(i,k)*xhnm(i,k) & ! /cm3(a) *const0/xl !----------------------------------------------------------------- ! now use bisection method to solve electro-neutrality equation ! ! during the iteration loop, ! yph_lo = lower ph value that brackets the root (i.e., correct ph) ! yph_hi = upper ph value that brackets the root (i.e., correct ph) ! yph = current ph value ! yposnet_lo and yposnet_hi = net positive ions for ! yph_lo and yph_hi !----------------------------------------------------------------- do iter = 1,itermax if (.not. present(yph_in)) then if (iter == 1) then ! 1st iteration ph = lower bound value yph_lo = 2.0_r8 yph_hi = yph_lo yph = yph_lo else if (iter == 2) then ! 2nd iteration ph = upper bound value yph_hi = 7.0_r8 yph = yph_hi else ! later iteration ph = mean of the two bracketing values yph = 0.5_r8*(yph_lo + yph_hi) end if else yph = yph_in end if ! calc current [H+] from ph xph(i,k) = 10.0_r8**(-yph) !----------------------------------------------------------------- ! ... hno3 !----------------------------------------------------------------- Ehno3 = fact1_hno3/(1.0_r8 + fact2_hno3*(1.0_r8 + fact3_hno3/xph(i,k))) !----------------------------------------------------------------- ! ... so2 !----------------------------------------------------------------- Eso2 = fact1_so2/(1.0_r8 + fact2_so2*(1.0_r8 + (fact3_so2/xph(i,k)) & *(1.0_r8 + fact4_so2/xph(i,k)))) !----------------------------------------------------------------- ! ... nh3 !----------------------------------------------------------------- Enh3 = fact1_nh3/(1.0_r8 + fact2_nh3*(1.0_r8 + fact3_nh3*xph(i,k))) tmp_nh4 = Enh3 * xph(i,k) tmp_hso3 = Eso2 / xph(i,k) tmp_so3 = tmp_hso3 * 2.0_r8*fact4_so2/xph(i,k) tmp_hco3 = Eco2 / xph(i,k) tmp_oh = Eh2o / xph(i,k) tmp_no3 = Ehno3 / xph(i,k) tmp_so4 = cldconc%so4_fact*Eso4 tmp_pos = xph(i,k) + tmp_nh4 tmp_neg = tmp_oh + tmp_hco3 + tmp_no3 + tmp_hso3 + tmp_so3 + tmp_so4 ynetpos = tmp_pos - tmp_neg ! yposnet = net positive ions/charge ! if the correct ph is bracketed by yph_lo and yph_hi (with yph_lo < yph_hi), ! then you will have yposnet_lo > 0 and yposnet_hi < 0 converged = .false. if (iter > 2) then if (ynetpos == 0.0_r8) then ! the exact solution was found (very unlikely) tmp_hp = xph(i,k) converged = .true. exit else if (ynetpos >= 0.0_r8) then ! net positive ions are >= 0 for both yph and yph_lo ! so replace yph_lo with yph yph_lo = yph ynetpos_lo = ynetpos else ! net positive ions are <= 0 for both yph and yph_hi ! so replace yph_hi with yph yph_hi = yph ynetpos_hi = ynetpos end if if (abs(yph_hi - yph_lo) .le. 0.005_r8) then ! |yph_hi - yph_lo| <= convergence criterion, so set ! final ph to their midpoint and exit ! (.005 absolute error in pH gives .01 relative error in H+) tmp_hp = xph(i,k) yph = 0.5_r8*(yph_hi + yph_lo) xph(i,k) = 10.0_r8**(-yph) converged = .true. exit else ! do another iteration converged = .false. end if else if (iter == 1) then if (ynetpos <= 0.0_r8) then ! the lower and upper bound ph values (2.0 and 7.0) do not bracket ! the correct ph, so use the lower bound tmp_hp = xph(i,k) converged = .true. exit end if ynetpos_lo = ynetpos else ! (iter == 2) if (ynetpos >= 0.0_r8) then ! the lower and upper bound ph values (2.0 and 7.0) do not bracket ! the correct ph, so use they upper bound tmp_hp = xph(i,k) converged = .true. exit end if ynetpos_hi = ynetpos end if end do ! iter if( .not. converged ) then write(iulog,*) 'SETSOX: pH failed to converge @ (',i,',',k,'), % change=', & 100._r8*delta end if else xph(i,k) = 1.e-7_r8 end if end do col_loop0 end do ver_loop0 ! end pver loop for STEP 0 !============================================================== ! ... Now use the actual PH !============================================================== ver_loop1: do k = 1,pver col_loop1: do i = 1,ncol work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8 tz = tfld(i,k) xl = cldconc%xlwc(i,k) patm = press(i,k)/101300._r8 ! press is in pascal xam = press(i,k)/(1.38e-23_r8*tz) ! air density /M3 !----------------------------------------------------------------------- ! ... hno3 !----------------------------------------------------------------------- xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) ) xe = 15.4_r8 hehno3(i,k) = xk*(1._r8 + xe/xph(i,k)) !----------------------------------------------------------------- ! ... h2o2 !----------------------------------------------------------------- xk = 7.4e4_r8 *EXP( 6621._r8*work1(i) ) xe = 2.2e-12_r8 *EXP(-3730._r8*work1(i) ) heh2o2(i,k) = xk*(1._r8 + xe/xph(i,k)) !----------------------------------------------------------------- ! ... so2 !----------------------------------------------------------------- xk = 1.23_r8 *EXP( 3120._r8*work1(i) ) xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) ) x2 = 6.0e-8_r8*EXP( 1120._r8*work1(i) ) wrk = xe/xph(i,k) heso2(i,k) = xk*(1._r8 + wrk*(1._r8 + x2/xph(i,k))) !----------------------------------------------------------------- ! ... nh3 !----------------------------------------------------------------- xk = 58._r8 *EXP( 4085._r8*work1(i) ) xe = 1.7e-5_r8*EXP(-4325._r8*work1(i) ) henh3(i,k) = xk*(1._r8 + xe*xph(i,k)/xkw) !----------------------------------------------------------------- ! ... o3 !----------------------------------------------------------------- xk = 1.15e-2_r8 *EXP( 2560._r8*work1(i) ) heo3(i,k) = xk !------------------------------------------------------------------------ ! ... for Ho2(g) -> H2o2(a) formation ! schwartz JGR, 1984, 11589 !------------------------------------------------------------------------ kh4 = (kh2 + kh3*kh1/xph(i,k)) / ((1._r8 + kh1/xph(i,k))**2) ho2s = kh0*xho2(i,k)*patm*(1._r8 + kh1/xph(i,k)) ! ho2s = ho2(a)+o2- r1h2o2 = kh4*ho2s*ho2s ! prod(h2o2) in mole/L(w)/s if ( cloud_borne ) then r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s / const0*1.e+6_r8 & ! correct a bug here ???? / xam else r2h2o2 = r1h2o2*xl & ! mole/L(w)/s * L(w)/fm3(a) = mole/fm3(a)/s * const0 & ! mole/fm3(a)/s * 1.e-3 = mole/cm3(a)/s / xam ! /cm3(a)/s / air-den = mix-ratio/s endif if ( .not. modal_aerosols ) then xh2o2(i,k) = xh2o2(i,k) + r2h2o2*dtime ! updated h2o2 by het production endif !----------------------------------------------- ! ... Partioning !----------------------------------------------- !----------------------------------------------------------------- ! ... hno3 !----------------------------------------------------------------- px = hehno3(i,k) * Ra * tz * xl hno3g(i,k) = (xhno3(i,k)+xno3(i,k))/(1._r8 + px) !------------------------------------------------------------------------ ! ... h2o2 !------------------------------------------------------------------------ px = heh2o2(i,k) * Ra * tz * xl h2o2g = xh2o2(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... so2 !------------------------------------------------------------------------ px = heso2(i,k) * Ra * tz * xl so2g = xso2(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... o3 !------------------------------------------------------------------------ px = heo3(i,k) * Ra * tz * xl o3g = xo3(i,k)/(1._r8+ px) !------------------------------------------------------------------------ ! ... nh3 !------------------------------------------------------------------------ px = henh3(i,k) * Ra * tz * xl if (id_nh3>0) then nh3g(i,k) = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px) else nh3g(i,k) = 0._r8 endif !----------------------------------------------- ! ... Aqueous phase reaction rates ! SO2 + H2O2 -> SO4 ! SO2 + O3 -> SO4 !----------------------------------------------- !------------------------------------------------------------------------ ! ... S(IV) (HSO3) + H2O2 !------------------------------------------------------------------------ rah2o2 = 8.e4_r8 * EXP( -3650._r8*work1(i) ) & / (.1_r8 + xph(i,k)) !------------------------------------------------------------------------ ! ... S(IV)+ O3 !------------------------------------------------------------------------ rao3 = 4.39e11_r8 * EXP(-4131._r8/tz) & + 2.56e3_r8 * EXP(-996._r8 /tz) /xph(i,k) !----------------------------------------------------------------- ! ... Prediction after aqueous phase ! so4 ! When Cloud is present ! ! S(IV) + H2O2 = S(VI) ! S(IV) + O3 = S(VI) ! ! reference: ! (1) Seinfeld ! (2) Benkovitz !----------------------------------------------------------------- !............................ ! S(IV) + H2O2 = S(VI) !............................ IF (XL .ge. 1.e-8_r8) THEN !! WHEN CLOUD IS PRESENTED if (cloud_borne) then patm_x = patm else patm_x = 1._r8 endif if (modal_aerosols) then pso4 = rah2o2 * 7.4e4_r8*EXP(6621._r8*work1(i)) * h2o2g * patm_x & * 1.23_r8 *EXP(3120._r8*work1(i)) * so2g * patm_x else pso4 = rah2o2 * heh2o2(i,k) * h2o2g * patm_x & * heso2(i,k) * so2g * patm_x ! [M/s] endif pso4 = pso4 & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / const0 & ! [/L(a)/s] / xhnm(i,k) ccc = pso4*dtime ccc = max(ccc, 1.e-30_r8) xso4_init(i,k)=xso4(i,k) IF (xh2o2(i,k) .gt. xso2(i,k)) THEN if (ccc .gt. xso2(i,k)) then xso4(i,k)=xso4(i,k)+xso2(i,k) if (cloud_borne) then xh2o2(i,k)=xh2o2(i,k)-xso2(i,k) xso2(i,k)=1.e-20_r8 else ! ???? bug ???? xso2(i,k)=1.e-20_r8 xh2o2(i,k)=xh2o2(i,k)-xso2(i,k) endif else xso4(i,k) = xso4(i,k) + ccc xh2o2(i,k) = xh2o2(i,k) - ccc xso2(i,k) = xso2(i,k) - ccc end if ELSE if (ccc .gt. xh2o2(i,k)) then xso4(i,k)=xso4(i,k)+xh2o2(i,k) xso2(i,k)=xso2(i,k)-xh2o2(i,k) xh2o2(i,k)=1.e-20_r8 else xso4(i,k) = xso4(i,k) + ccc xh2o2(i,k) = xh2o2(i,k) - ccc xso2(i,k) = xso2(i,k) - ccc end if END IF if (modal_aerosols) then xdelso4hp(i,k) = xso4(i,k) - xso4_init(i,k) endif !........................... ! S(IV) + O3 = S(VI) !........................... pso4 = rao3 * heo3(i,k)*o3g*patm_x * heso2(i,k)*so2g*patm_x ! [M/s] pso4 = pso4 & ! [M/s] = [mole/L(w)/s] * xl & ! [mole/L(a)/s] / const0 & ! [/L(a)/s] / xhnm(i,k) ! [mixing ratio/s] ccc = pso4*dtime ccc = max(ccc, 1.e-30_r8) xso4_init(i,k)=xso4(i,k) if (ccc .gt. xso2(i,k)) then xso4(i,k) = xso4(i,k) + xso2(i,k) xso2(i,k) = 1.e-20_r8 else xso4(i,k) = xso4(i,k) + ccc xso2(i,k) = xso2(i,k) - ccc end if END IF !! WHEN CLOUD IS PRESENTED end do col_loop1 end do ver_loop1 call sox_cldaero_update( & ncol, lchnk, loffset, dtime, mbar, pdel, press, tfld, cldnum, cldfrc, cfact, cldconc%xlwc, & xdelso4hp, xh2so4, xso4, xso4_init, nh3g, hno3g, xnh3, xhno3, xnh4c, xno3c, xmsa, xso2, xh2o2, qcw, qin, & aqso4, aqh2so4, aqso4_h2o2, aqso4_o3 ) xphlwc(:,:) = 0._r8 do k = 1, pver do i = 1, ncol if (cldfrc(i,k)>=1.e-5_r8 .and. lwc(i,k)>=1.e-8_r8) then xphlwc(i,k) = -1._r8*log10(xph(i,k)) * lwc(i,k) endif end do end do call sox_cldaero_destroy_obj(cldconc) end subroutine SETSOX end module MO_SETSOX