!Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
#include "MODAL_AERO_CPP_DEFINES.h"
#define WRF_PORT
#define MODAL_AERO
module MO_SETSOX
use shr_kind_mod, only : r8 => shr_kind_r8
#ifndef WRF_PORT
use cam_logfile, only : iulog
#else
use module_cam_support, only: iulog
#endif
#if (defined MODAL_AERO)
use modal_aero_data
#endif
private
public :: sox_inti, setsox
public :: has_sox
save
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
integer, target :: spc_ids(18)
integer, pointer :: id_so2, id_so4_1, id_so4_2, id_so4_3, id_so4_4, id_so4_5, id_so4_6, &
id_nh3, id_nh4_1, id_nh4_2, id_nh4_3, id_nh4_4, id_nh4_5, id_nh4_6, &
id_h2o2, id_ox, id_ho2, id_h2so4
#elif ( defined MODAL_AERO_3MODE )
integer, target :: spc_ids(8)
integer, pointer :: id_so2, id_so4_1, id_so4_2, id_so4_3, &
id_h2o2, id_ox, id_ho2, id_h2so4
integer :: id_nh3
#endif
logical :: inv_o3
integer :: id_hno3, id_msa
#else
integer, target :: spc_ids(8)
integer, pointer :: id_so2, id_so4, id_nh3, id_hno3, id_h2o2, id_ox, id_nh4no3, id_ho2
#endif
logical :: has_sox = .true.
logical :: inv_so2, inv_nh3, inv_hno3, inv_h2o2, inv_ox, inv_nh4no3, inv_ho2
contains
subroutine sox_inti
!-----------------------------------------------------------------------
! ... initialize the hetero sox routine
!-----------------------------------------------------------------------
use mo_chem_utls, only : get_spc_ndx, get_inv_ndx
#ifndef WRF_PORT
use spmd_utils, only : masterproc
#else
use module_cam_support, only: masterproc
#endif
implicit none
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
id_so2 => spc_ids(1)
id_so4_1 => spc_ids(2)
id_so4_2 => spc_ids(3)
id_so4_3 => spc_ids(4)
id_so4_4 => spc_ids(5)
id_so4_5 => spc_ids(6)
id_so4_6 => spc_ids(7)
id_nh3 => spc_ids(8)
id_nh4_1 => spc_ids(9)
id_nh4_2 => spc_ids(10)
id_nh4_3 => spc_ids(11)
id_nh4_4 => spc_ids(12)
id_nh4_5 => spc_ids(13)
id_nh4_6 => spc_ids(14)
id_h2o2 => spc_ids(15)
id_ox => spc_ids(16)
id_ho2 => spc_ids(17)
id_h2so4 => spc_ids(18)
#elif ( defined MODAL_AERO_3MODE )
id_so2 => spc_ids(1)
id_so4_1 => spc_ids(2)
id_so4_2 => spc_ids(3)
id_so4_3 => spc_ids(4)
id_h2o2 => spc_ids(5)
id_ox => spc_ids(6)
id_ho2 => spc_ids(7)
id_h2so4 => spc_ids(8)
#endif
#else
id_so2 => spc_ids(1)
id_so4 => spc_ids(2)
id_nh3 => spc_ids(3)
id_hno3 => spc_ids(4)
id_h2o2 => spc_ids(5)
id_ox => spc_ids(6)
id_ho2 => spc_ids(7)
#endif
!-----------------------------------------------------------------
! ... get species indicies
!-----------------------------------------------------------------
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
id_so4_1 = lptr_so4_cw_amode(1)
id_so4_2 = lptr_so4_cw_amode(2)
id_so4_3 = lptr_so4_cw_amode(4) ! no so4_c3
id_so4_4 = lptr_so4_cw_amode(5)
id_so4_5 = lptr_so4_cw_amode(6)
id_so4_6 = lptr_so4_cw_amode(7)
#elif ( defined MODAL_AERO_3MODE )
id_so4_1 = lptr_so4_cw_amode(1)
id_so4_2 = lptr_so4_cw_amode(2)
id_so4_3 = lptr_so4_cw_amode(3)
#endif
id_h2so4 = get_spc_ndx( 'H2SO4' )
id_msa = get_spc_ndx( 'MSA' )
#else
id_so4 = get_spc_ndx( 'SO4' )
#endif
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
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
id_nh4_1 = lptr_nh4_cw_amode(1)
id_nh4_2 = lptr_nh4_cw_amode(2)
id_nh4_3 = lptr_nh4_cw_amode(4) ! no nh4_c3
id_nh4_4 = lptr_nh4_cw_amode(5)
id_nh4_5 = lptr_nh4_cw_amode(6)
id_nh4_6 = lptr_nh4_cw_amode(7)
#endif
#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
#if (defined MODAL_AERO)
inv_o3 = get_inv_ndx( 'O3' ) > 0
if (inv_o3) then
id_ox = get_inv_ndx( 'O3' )
else
id_ox = 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
#else
inv_OX = .false.
id_OX = get_inv_ndx( 'OX' )
if( id_ox < 1 ) then
id_ox = get_inv_ndx( 'O3' )
end if
inv_OX = id_OX > 0
if ( .not. inv_OX ) then
id_OX = get_spc_ndx( 'OX' )
if ( id_OX < 0 ) then
id_OX = get_spc_ndx( 'O3' )
endif
endif
#endif
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
has_sox = all( spc_ids(1:18) > 0 )
#elif ( defined MODAL_AERO_3MODE )
has_sox = all( spc_ids(1:8) > 0 )
#endif
#else
has_sox = all( spc_ids(1:7) > 0 )
#endif
if (masterproc) then
write(iulog,*) 'sox_inti: has_sox = ',has_sox
#ifdef WRF_PORT
call wrf_message(iulog)
#endif
endif
if( has_sox ) then
if (masterproc) then
write(iulog,*) '-----------------------------------------'
#ifdef WRF_PORT
call wrf_message(iulog)
#endif
write(iulog,*) 'mozart will do sox aerosols'
#ifdef WRF_PORT
call wrf_message(iulog)
#endif
write(iulog,*) '-----------------------------------------'
#ifdef WRF_PORT
call wrf_message(iulog)
#endif
endif
end if
end subroutine sox_inti
subroutine SETSOX( ncol, &
press, &
dtime, &
tfld, &
qfld, &
lwc, &
#if (defined MODAL_AERO)
lchnk, &
pdel, &
mbar, &
clwlrat,&
cldfrc, &
cldnum, &
qcw, &
loffset,&
#endif
xhnm, &
qin, invariants )
!-----------------------------------------------------------------------
! ... 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
!-----------------------------------------------------------------------
!
#ifndef WRF_PORT
use ppgrid, only : pcols, pver
use chem_mods, only : gas_pcnst, nfs
#else
use module_cam_support, only: pcols, pver, gas_pcnst => gas_pcnst_modal_aero, &
nfs
#endif
#if (defined MODAL_AERO)
#ifndef WRF_PORT
use chem_mods, only : adv_mass
#else
use module_data_cam_mam_asect, only: adv_mass => mw_q_mo_array
#endif
use physconst, only : mwdry, gravit
#ifndef WRF_PORT
use mo_constants, only : pi
use cam_history, only : outfld
#else
use physconst, only : pi
use module_cam_support, only: outfld
#endif
#endif
!
implicit none
!
!-----------------------------------------------------------------------
! ... Dummy arguments
!-----------------------------------------------------------------------
integer, intent(in) :: ncol
real(r8), intent(in) :: dtime ! time step (sec)
real(r8), dimension(ncol ,pver,gas_pcnst), intent(inout) :: qin ! xported species ( vmr )
real(r8), dimension(ncol ,pver), intent(in) :: xhnm ! total atms density ( /cm**3)
real(r8), dimension(pcols,pver), intent(in) :: tfld, & ! temperature
qfld, & ! specific humidity( kg/kg )
press ! midpoint pressure ( Pa )
real(r8), dimension(ncol ,pver), intent(in) :: lwc !!, & ! cloud liquid water content (kg/kg)
real(r8), intent(in) :: invariants(ncol,pver,nfs)
#if (defined MODAL_AERO)
integer, intent(in) :: lchnk ! chunk id
real(r8), dimension(pcols,pver), intent(in) :: pdel ! pressure thickness of levels (Pa)
real(r8), dimension(ncol ,pver), intent(in) :: mbar ! mean wet atmospheric mass ( amu )
real(r8), dimension(ncol ,pver), intent(in) :: cldfrc, & ! cloud fraction
clwlrat,& ! first-order loss rate of l.s. condensate to precip (1/s)
cldnum ! droplet number concentration (#/kg)
real(r8), intent(inout) :: qcw(ncol,pver,gas_pcnst) ! cloud-borne aerosol (vmr)
integer, intent(in) :: loffset ! offset applied to modal aero "ptrs"
#endif
!-----------------------------------------------------------------------
! ... 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, &
xb0 = -.1_r8, &
xa1 = 1.053_r8, &
xb1 = -4.368_r8,&
xa2 = 1.016_r8, &
xb2 = -2.54_r8, &
xa3 = .816e-32_r8, &
xb3 = .259_r8
real(r8), parameter :: kh0 = 9.e3_r8, & ! HO2(g) -> Ho2(a)
kh1 = 2.05e-5_r8, & ! HO2(a) -> H+ + O2-
kh2 = 8.6e5_r8, & ! HO2(a) + ho2(a) -> h2o2(a) + o2
kh3 = 1.e8_r8, & ! HO2(a) + o2- -> h2o2(a) + o2
Ra = 8314._r8/101325._r8, & ! universal constant (atm)/(M-K)
xkw = 1.e-14_r8 ! water acidity
real(r8), parameter :: small_value = 1.e-20_r8
!
#if (defined MODAL_AERO)
integer :: l, n, m
integer :: ntot_msa_c
integer :: id_so4_1a, id_so4_2a, id_so4_3a, id_so4_4a, id_so4_5a, id_so4_6a, &
id_nh4_1a, id_nh4_2a, id_nh4_3a, id_nh4_4a, id_nh4_5a, id_nh4_6a
real(r8) :: delso4_hprxn, delso4_o3rxn, &
dso4dt_aqrxn, dso4dt_hprxn, &
dso4dt_gasuptk, dmsadt_gasuptk, &
dmsadt_gasuptk_tomsa, dmsadt_gasuptk_toso4, &
dqdt_aq, dqdt_wr, dqdt, &
rad_cd, radxnum_cd, num_cd, &
gasdiffus, gasspeed, knudsen, uptkrate, &
fuchs_sutugin, volx34pi_cd, &
fwetrem, sumf, &
frso2_g, frso2_c, frh2o2_g, frh2o2_c
real(r8) :: delnh3, delnh4
real(r8) :: faqgain_msa(ntot_amode), faqgain_so4(ntot_amode), &
qnum_c(ntot_amode)
real(r8) :: dqdt_aqso4(ncol,pver,gas_pcnst), &
dqdt_aqh2so4(ncol,pver,gas_pcnst), &
dqdt_aqhprxn(ncol,pver), dqdt_aqo3rxn(ncol,pver), &
xphlwc(ncol,pver), &
sflx(1:ncol)
#endif
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) :: hno3g, nh3g, 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
!
!-----------------------------------------------------------------------
! 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,&
xnh3, xnh4, xo3, &
xlwc, cfact, &
xph, xho2, &
#if (defined MODAL_AERO)
xh2so4, xmsa, xso4_init, xso4c, xnh4c, &
#endif
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), dimension(ncol) :: work1
logical :: converged
#if (defined MODAL_AERO)
id_so4_1a = id_so4_1 - loffset
id_so4_2a = id_so4_2 - loffset
id_so4_3a = id_so4_3 - loffset
#if ( defined MODAL_AERO_7MODE )
id_so4_4a = id_so4_4 - loffset
id_so4_5a = id_so4_5 - loffset
id_so4_6a = id_so4_6 - loffset
id_nh4_1a = id_nh4_1 - loffset
id_nh4_2a = id_nh4_2 - loffset
id_nh4_3a = id_nh4_3 - loffset
id_nh4_4a = id_nh4_4 - loffset
id_nh4_5a = id_nh4_5 - loffset
id_nh4_6a = id_nh4_6 - loffset
#endif
! make sure dqdt is zero initially, for budgets
dqdt_aqso4(:,:,:) = 0.0_r8
dqdt_aqh2so4(:,:,:) = 0.0_r8
dqdt_aqhprxn(:,:) = 0.0_r8
dqdt_aqo3rxn(:,:) = 0.0_r8
#endif
!-----------------------------------------------------------------
! ... 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
do k = 1,pver
xph(:,k) = xph0 ! initial PH value
xlwc(:,k) = lwc(:,k) *cfact(:,k) ! cloud water L(water)/L(air)
if ( inv_so2 ) then
xso2 (:,k) = invariants(:,k,id_so2) ! mixing ratio
else
xso2 (:,k) = qin(:,k,id_so2) ! mixing ratio
endif
#if (defined MODAL_AERO)
if (id_hno3 > 0) then
xhno3(:,k) = qin(:,k,id_hno3)
else
xhno3(:,k) = 0.0_r8
endif
#else
if ( inv_hno3 ) then
xhno3 (:,k) = invariants(:,k,id_hno3) ! mixing ratio
else
xhno3 (:,k) = qin(:,k,id_hno3) ! mixing ratio
endif
#endif
if ( inv_h2o2 ) then
xh2o2 (:,k) = invariants(:,k,id_h2o2) ! mixing ratio
else
xh2o2 (:,k) = qin(:,k,id_h2o2) ! mixing ratio
endif
#if (defined MODAL_AERO)
if (id_nh3 > 0) then
xnh3 (:,k) = qin(:,k,id_nh3)
else
xnh3 (:,k) = 0.0_r8
endif
#else
if ( inv_nh3 ) then
xnh3 (:,k) = invariants(:,k,id_nh3) ! mixing ratio
else
xnh3 (:,k) = qin(:,k,id_nh3) ! mixing ratio
endif
#endif
#if (defined MODAL_AERO)
if ( inv_o3 ) then
xo3 (:,k) = invariants(:,k,id_ox)/xhnm(:,k) ! mixing ratio
else
xo3 (:,k) = qin(:,k,id_ox) ! 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
#else
if ( inv_ox ) then
xo3 (:,k) = invariants(:,k,id_ox) ! mixing ratio
else
xo3 (:,k) = qin(:,k,id_ox) ! mixing ratio
endif
if ( inv_ho2 ) then
xho2 (:,k) = invariants(:,k,id_ho2) ! mixing ratio
else
xho2 (:,k) = qin(:,k,id_ho2) ! mixing ratio
endif
#endif
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
xso4c(:,k) = qcw(:,k,id_so4_1a) &
+ qcw(:,k,id_so4_2a) &
+ qcw(:,k,id_so4_3a) &
+ qcw(:,k,id_so4_4a) &
+ qcw(:,k,id_so4_5a) &
+ qcw(:,k,id_so4_6a)
#elif ( defined MODAL_AERO_3MODE )
xso4c(:,k) = qcw(:,k,id_so4_1a) &
+ qcw(:,k,id_so4_2a) &
+ qcw(:,k,id_so4_3a)
#endif
xh2so4(:,k)= qin(:,k,id_h2so4)
if (id_msa > 0) xmsa (:,k) = qin(:,k,id_msa)
#else
xso4 (:,k) = qin(:,k,id_so4) ! mixing ratio
#endif
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
xnh4c(:,k) = qcw(:,k,id_nh4_1a) &
+ qcw(:,k,id_nh4_2a) &
+ qcw(:,k,id_nh4_3a) &
+ qcw(:,k,id_nh4_4a) &
+ qcw(:,k,id_nh4_5a) &
+ qcw(:,k,id_nh4_6a)
#elif ( defined MODAL_AERO_3MODE )
xnh4c(:,k) = xso4c(:,k) * 1.0_r8 !assume NH4HSO4
#endif
#endif
end do
!-----------------------------------------------------------------
! ... Temperature dependent Henry constants
!-----------------------------------------------------------------
do k = 1,pver !! pver loop for STEP 0
do i = 1,ncol
#if (defined MODAL_AERO)
if (cldfrc(i,k) >= 1.0e-5_r8) then
xso4(i,k) = xso4c(i,k) / cldfrc(i,k)
xnh4(i,k) = xnh4c(i,k) / cldfrc(i,k)
xl = xlwc(i,k) / cldfrc(i,k) ! liquid water in the cloudy fraction of cell
#else
xl = xlwc(i,k)
#endif
if( xl >= 1.e-8_r8 ) then
work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8
!-----------------------------------------------------------------------
! ... hno3
!-----------------------------------------------------------------------
do iter = 1,itermax
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)
!-----------------------------------------------------------------
! ... Partioning and effect of pH
!-----------------------------------------------------------------
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
!-----------------------------------------------------------------
px = hehno3(i,k) * Ra * tz * xl
hno3g = xhno3(i,k)/(1._r8 + px)
xk = 2.1e5_r8 *EXP( 8700._r8*work1(i) )
xe = 15.4_r8
Ehno3 = xk*xe*hno3g *patm
!-----------------------------------------------------------------
! ... so2
!-----------------------------------------------------------------
px = heso2(i,k) * Ra * tz * xl
so2g = xso2(i,k)/(1._r8+ px)
xk = 1.23_r8 *EXP( 3120._r8*work1(i) )
xe = 1.7e-2_r8*EXP( 2090._r8*work1(i) )
Eso2 = xk*xe*so2g *patm
!-----------------------------------------------------------------
! ... nh3
!-----------------------------------------------------------------
px = henh3(i,k) * Ra * tz * xl
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px)
#elif ( defined MODAL_AERO_3MODE )
nh3g = xnh4(i,k)/px
#endif
#else
nh3g = xnh3(i,k)/(1._r8+ px)
#endif
xk = 58._r8 *EXP( 4085._r8*work1(i) )
xe = 1.7e-5_r8*EXP( -4325._r8*work1(i) )
Enh3 = xk*xe*nh3g/xkw *patm
!-----------------------------------------------------------------
! ... 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
!-----------------------------------------------------------------
! ... PH cal
!-----------------------------------------------------------------
com2 = (Eh2o + Ehno3 + Eso2 + Eco2) &
/ (1._r8 + Enh3 )
com2 = MAX( com2,1.e-20_r8 )
xph(i,k) = SQRT( com2 )
!-----------------------------------------------------------------
! ... Add so4 effect
!-----------------------------------------------------------------
Eso4 = xso4(i,k)*xhnm(i,k) & ! /cm3(a)
*const0/xl
xph(i,k) = MIN( 1.e-2_r8,MAX( 1.e-7_r8,xph(i,k) + 2._r8*Eso4 ) )
if( iter > 1 ) then
delta = ABS( (xph(i,k) - delta)/delta )
converged = delta < .01_r8
if( converged ) then
exit
else
delta = xph(i,k)
end if
else
delta = xph(i,k)
end if
end do
if( .not. converged ) then
write(iulog,*) 'SETSOX: pH failed to converge @ (',i,',',k,'), % change=', &
100._r8*delta
#ifdef WRF_PORT
call wrf_message(iulog)
#endif
end if
else
xph(i,k) = 1.e-7_r8
end if
#if (defined MODAL_AERO)
end if ! if (cldfrc(i,k) >= 1.0e-5_r8)
#endif
end do
end do ! end pver loop for STEP 0
!==============================================================
! ... Now use the actual PH
!==============================================================
do k = 1,pver
do i = 1,ncol
work1(i) = 1._r8 / tfld(i,k) - 1._r8 / 298._r8
tz = tfld(i,k)
#if (defined MODAL_AERO)
if (cldfrc(i,k) >= 1.0e-5_r8) then
xl = xlwc(i,k) / cldfrc(i,k)
#else
xl = xlwc(i,k)
#endif
patm = press(i,k)/101300._r8 ! press is in pascal
xam = press(i,k)/(1.38e-23_r8*tz) ! air density /M3
!-----------------------------------------------------------------
! ... 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)))
#if (defined MODAL_AERO)
!-----------------------------------------------------------------
! ... 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)
#endif
!-----------------------------------------------------------------
! ... 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 (defined MODAL_AERO)
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
r2h2o2 = 0.0_r8
#else
r2h2o2 = r1h2o2*xlwc(i,k) & ! 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
xh2o2(i,k) = xh2o2(i,k) + r2h2o2*dtime ! updated h2o2 by het production
!-----------------------------------------------
! ... Partioning
!-----------------------------------------------
!------------------------------------------------------------------------
! ... 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)
#if (defined MODAL_AERO)
!------------------------------------------------------------------------
! ... nh3
!------------------------------------------------------------------------
px = henh3(i,k) * Ra * tz * xl
#if ( defined MODAL_AERO_7MODE )
nh3g = (xnh3(i,k)+xnh4(i,k))/(1._r8+ px)
#elif ( defined MODAL_AERO_3MODE )
if(xl .ge. 1.e-8_r8) then
nh3g = xnh4(i,k)/px
endif
#endif
#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 (defined MODAL_AERO)
pso4 = rah2o2 * 7.4e4_r8*EXP(6621._r8*work1(i)) * h2o2g * patm &
* 1.23_r8 *EXP(3120._r8*work1(i)) * so2g * patm
pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
* xl & ! [mole/L(a)/s]
/ const0 & ! [/L(a)/s]
/ xhnm(i,k)
#else
pso4 = rah2o2 * heh2o2(i,k)*h2o2g &
* heso2(i,k) *so2g ! [M/s]
pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
* xlwc(i,k) & ! [mole/L(a)/s]
/ const0 & ! [/L(a)/s]
/ xhnm(i,k)
#endif
ccc = pso4*dtime
ccc = max(ccc, 1.e-30_r8)
#if (defined MODAL_AERO)
xso4_init(i,k)=xso4(i,k)
#endif
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 (defined MODAL_AERO)
xh2o2(i,k)=xh2o2(i,k)-xso2(i,k)
xso2(i,k)=1.e-20_r8
#else
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 (defined MODAL_AERO)
delso4_hprxn = xso4(i,k) - xso4_init(i,k)
#endif
!...........................
! S(IV) + O3 = S(VI)
!...........................
#if (defined MODAL_AERO)
pso4 = rao3 * heo3(i,k)*o3g*patm * heso2(i,k)*so2g*patm ! [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]
#else
pso4 = rao3 * heo3(i,k)*o3g * heso2(i,k)*so2g ! [M/s]
pso4 = pso4 & ! [M/s] = [mole/L(w)/s]
* xlwc(i,k) & ! [mole/L(a)/s]
/ const0 & ! [/L(a)/s]
/ xhnm(i,k) ! [mixing ratio/s]
#endif
ccc = pso4*dtime
ccc = max(ccc, 1.e-30_r8)
#if (defined MODAL_AERO)
xso4_init(i,k)=xso4(i,k)
#endif
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
#if (defined MODAL_AERO)
delso4_o3rxn = xso4(i,k) - xso4_init(i,k)
#endif
#if (defined MODAL_AERO)
#if ( defined MODAL_AERO_7MODE )
delnh3 = nh3g - xnh3(i,k)
delnh4 = - delnh3
#endif
#endif
#if (defined MODAL_AERO)
!-------------------------------------------------------------------------
! compute factors for partitioning aerosol mass gains among modes
! the factors are proportional to the activated particle MR for each
! mode, which is the MR of cloud drops "associated with" the mode
! thus we are assuming the cloud drop size is independent of the
! associated aerosol mode properties (i.e., drops associated with
! Aitken and coarse sea-salt particles are same size)
!
! qnum_c(n) = activated particle number MR for mode n (these are just
! used for partitioning among modes, so don't need to divide by cldfrc)
do n = 1, ntot_amode
qnum_c(n) = 0.0
l = numptrcw_amode(n) - loffset
if (l > 0) qnum_c(n) = max( 0.0_r8, qcw(i,k,l) )
end do
! force qnum_c(n) to be positive for n=modeptr_accum or n=1
n = modeptr_accum
if (n <= 0) n = 1
qnum_c(n) = max( 1.0e-10_r8, qnum_c(n) )
! faqgain_so4(n) = fraction of total so4_c gain going to mode n
! these are proportional to the activated particle MR for each mode
sumf = 0.0
do n = 1, ntot_amode
faqgain_so4(n) = 0.0
if (lptr_so4_cw_amode(n) > 0) then
faqgain_so4(n) = qnum_c(n)
sumf = sumf + faqgain_so4(n)
end if
end do
if (sumf > 0.0) then
do n = 1, ntot_amode
faqgain_so4(n) = faqgain_so4(n) / sumf
end do
end if
! at this point (sumf <= 0.0) only when all the faqgain_so4 are zero
! faqgain_msa(n) = fraction of total msa_c gain going to mode n
ntot_msa_c = 0
sumf = 0.0
do n = 1, ntot_amode
faqgain_msa(n) = 0.0
if (lptr_msa_cw_amode(n) > 0) then
faqgain_msa(n) = qnum_c(n)
ntot_msa_c = ntot_msa_c + 1
end if
sumf = sumf + faqgain_msa(n)
end do
if (sumf > 0.0) then
do n = 1, ntot_amode
faqgain_msa(n) = faqgain_msa(n) / sumf
end do
end if
! at this point (sumf <= 0.0) only when all the faqgain_msa are zero
!-----------------------------------------------------------------------
! compute uptake of h2so4 and msa to cloud water
!
! first-order uptake rate is
! 4*pi*(drop radius)*(drop number conc)
! *(gas diffusivity)*(fuchs sutugin correction)
! num_cd = (drop number conc in 1/cm^3)
num_cd = 1.0e-3_r8*cldnum(i,k)*cfact(i,k)/cldfrc(i,k)
num_cd = max( num_cd, 0.0_r8 )
! rad_cd = (drop radius in cm), computed from liquid water and drop number,
! then bounded by 0.5 and 50.0 micrometers
! radxnum_cd = (drop radius)*(drop number conc)
! volx34pi_cd = (3/4*pi) * (liquid water volume in cm^3/cm^3)
volx34pi_cd = xl*0.75_r8/pi
! following holds because volx34pi_cd = num_cd*(rad_cd**3)
radxnum_cd = (volx34pi_cd*num_cd*num_cd)**0.3333333_r8
! apply bounds to rad_cd to avoid the occasional unphysical value
if (radxnum_cd .le. volx34pi_cd*4.0e4_r8) then
radxnum_cd = volx34pi_cd*4.0e4_r8
rad_cd = 50.0e-4_r8
else if (radxnum_cd .ge. volx34pi_cd*4.0e8_r8) then
radxnum_cd = volx34pi_cd*4.0e8_r8
rad_cd = 0.5e-4_r8
else
rad_cd = radxnum_cd/num_cd
end if
! gasdiffus = h2so4 gas diffusivity from mosaic code (cm^2/s)
! (pmid must be Pa)
gasdiffus = 0.557_r8 * (tfld(i,k)**1.75_r8) / press(i,k)
! gasspeed = h2so4 gas mean molecular speed from mosaic code (cm/s)
gasspeed = 1.455e4_r8 * sqrt(tfld(i,k)/98.0_r8)
! knudsen number
knudsen = 3.0_r8*gasdiffus/(gasspeed*rad_cd)
! following assumes accomodation coefficient = 0.65
! (Adams & Seinfeld, 2002, JGR, and references therein)
! fuchs_sutugin = (0.75*accom*(1. + knudsen)) /
! (knudsen*(1.0 + knudsen + 0.283*accom) + 0.75*accom)
fuchs_sutugin = (0.4875_r8*(1. + knudsen)) / &
(knudsen*(1.184_r8 + knudsen) + 0.4875_r8)
! instantaneous uptake rate
uptkrate = 12.56637_r8*radxnum_cd*gasdiffus*fuchs_sutugin
! average uptake rate over dtime
uptkrate = (1.0_r8 - exp(-min(100._r8,dtime*uptkrate))) / dtime
! dso4dt_gasuptk = so4_c tendency from h2so4 gas uptake (mol/mol/s)
! dmsadt_gasuptk = msa_c tendency from msa gas uptake (mol/mol/s)
dso4dt_gasuptk = xh2so4(i,k) * uptkrate
if (id_msa > 0) then
dmsadt_gasuptk = xmsa(i,k) * uptkrate
else
dmsadt_gasuptk = 0.0_r8
end if
! if no modes have msa aerosol, then "rename" scavenged msa gas to so4
dmsadt_gasuptk_toso4 = 0.0_r8
dmsadt_gasuptk_tomsa = dmsadt_gasuptk
if (ntot_msa_c == 0) then
dmsadt_gasuptk_tomsa = 0.0_r8
dmsadt_gasuptk_toso4 = dmsadt_gasuptk
end if
!-----------------------------------------------------------------------
! now compute TMR tendencies
! this includes the above aqueous so2 chemistry AND
! the uptake of highly soluble aerosol precursor gases (h2so4, msa, ...)
! AND the wetremoval of dissolved, unreacted so2 and h2o2
dso4dt_aqrxn = (delso4_o3rxn + delso4_hprxn) / dtime
dso4dt_hprxn = delso4_hprxn / dtime
! fwetrem = fraction of in-cloud-water material that is wet removed
! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*clwlrat(i,k)))) )
fwetrem = 0.0 ! don't have so4 & msa wet removal here
! compute TMR tendencies for so4 and msa aerosol-in-cloud-water
do n = 1, ntot_amode
l = lptr_so4_cw_amode(n) - loffset
if (l > 0) then
dqdt_aqso4(i,k,l) = faqgain_so4(n)*dso4dt_aqrxn*cldfrc(i,k)
dqdt_aqh2so4(i,k,l) = faqgain_so4(n)* &
(dso4dt_gasuptk + dmsadt_gasuptk_toso4)*cldfrc(i,k)
dqdt_aq = dqdt_aqso4(i,k,l) + dqdt_aqh2so4(i,k,l)
dqdt_wr = -fwetrem*dqdt_aq
dqdt= dqdt_aq + dqdt_wr
qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
end if
l = lptr_msa_cw_amode(n) - loffset
if (l > 0) then
dqdt_aq = faqgain_msa(n)*dmsadt_gasuptk_tomsa*cldfrc(i,k)
dqdt_wr = -fwetrem*dqdt_aq
dqdt = dqdt_aq + dqdt_wr
qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
end if
l = lptr_nh4_cw_amode(n) - loffset
if (l > 0) then
if (delnh4 > 0.0_r8) then
dqdt_aq = faqgain_so4(n)*delnh4/dtime*cldfrc(i,k)
dqdt = dqdt_aq
qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
else
dqdt = (qcw(i,k,l)/max(xnh4c(i,k),1.0e-35_r8)) &
*delnh4/dtime*cldfrc(i,k)
qcw(i,k,l) = qcw(i,k,l) + dqdt*dtime
endif
end if
end do
! For gas species, tendency includes
! reactive uptake to cloud water that essentially transforms the gas to
! a different species. Wet removal associated with this is applied
! to the "new" species (e.g., so4_c) rather than to the gas.
! wet removal of the unreacted gas that is dissolved in cloud water.
! Need to multiply both these parts by cldfrc
! h2so4 (g) & msa (g)
qin(i,k,id_h2so4) = qin(i,k,id_h2so4) - dso4dt_gasuptk * dtime * cldfrc(i,k)
if (id_msa > 0) qin(i,k,id_msa) = qin(i,k,id_msa) - dmsadt_gasuptk * dtime * cldfrc(i,k)
! frso2_g is fraction of total so2 as gas, etc.
frso2_g = 1.0_r8 / ( 1.0_r8 + heso2(i,k) * Ra * tz * xl )
frh2o2_g = 1.0_r8 / ( 1.0_r8 + heh2o2(i,k) * Ra * tz * xl )
! frso2_c is fraction of total so2 in cloud water, etc.
frso2_c = max( 0.0_r8, (1.0_r8-frso2_g) )
frh2o2_c = max( 0.0_r8, (1.0_r8-frh2o2_g) )
! so2 -- the first order loss rate for so2 is frso2_c*clwlrat(i,k)
! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frso2_c*clwlrat(i,k)))) )
fwetrem = 0.0 ! don't include so2 wet removal here
dqdt_wr = -fwetrem*xso2(i,k)/dtime*cldfrc(i,k)
dqdt_aq = -dso4dt_aqrxn*cldfrc(i,k)
dqdt = dqdt_aq + dqdt_wr
qin(i,k,id_so2) = qin(i,k,id_so2) + dqdt * dtime
! h2o2 -- the first order loss rate for h2o2 is frh2o2_c*clwlrat(i,k)
! fwetrem = max( 0.0_r8, (1.0_r8-exp(-min(100._r8,dtime*frh2o2_c*clwlrat(i,k)))) )
fwetrem = 0.0 ! don't include h2o2 wet removal here
dqdt_wr = -fwetrem*xh2o2(i,k)/dtime*cldfrc(i,k)
dqdt_aq = -dso4dt_hprxn*cldfrc(i,k)
dqdt = dqdt_aq + dqdt_wr
qin(i,k,id_h2o2) = qin(i,k,id_h2o2) + dqdt * dtime
! NH3
#if ( defined MODAL_AERO_7MODE )
dqdt_aq = delnh3/dtime*cldfrc(i,k)
dqdt = dqdt_aq
qin(i,k,id_nh3) = qin(i,k,id_nh3) + dqdt * dtime
#endif
! for SO4 from H2O2/O3 budgets
dqdt_aqhprxn(i,k) = dso4dt_hprxn*cldfrc(i,k)
dqdt_aqo3rxn(i,k) = (dso4dt_aqrxn - dso4dt_hprxn)*cldfrc(i,k)
#endif
END IF !! WHEN CLOUD IS PRESENTED
#if (defined MODAL_AERO)
end if !! if (cldfrc(i,k) >= 1.0e-5_r8) then
#endif
end do
end do
!==============================================================
! ... Update the mixing ratios
!==============================================================
do k = 1,pver
#if (defined MODAL_AERO)
qin(:,k,id_so2) = MAX( qin(:,k,id_so2), small_value )
qcw(:,k,id_so4_1a) = MAX( qcw(:,k,id_so4_1a), small_value )
qcw(:,k,id_so4_2a) = MAX( qcw(:,k,id_so4_2a), small_value )
qcw(:,k,id_so4_3a) = MAX( qcw(:,k,id_so4_3a), small_value )
#if ( defined MODAL_AERO_7MODE )
qcw(:,k,id_so4_4a) = MAX( qcw(:,k,id_so4_4a), small_value )
qcw(:,k,id_so4_5a) = MAX( qcw(:,k,id_so4_5a), small_value )
qcw(:,k,id_so4_6a) = MAX( qcw(:,k,id_so4_6a), small_value )
#endif
#if ( defined MODAL_AERO_7MODE )
qin(:,k,id_nh3) = MAX( qin(:,k,id_nh3), small_value )
qcw(:,k,id_nh4_1a) = MAX( qcw(:,k,id_nh4_1a), small_value )
qcw(:,k,id_nh4_2a) = MAX( qcw(:,k,id_nh4_2a), small_value )
qcw(:,k,id_nh4_3a) = MAX( qcw(:,k,id_nh4_3a), small_value )
qcw(:,k,id_nh4_4a) = MAX( qcw(:,k,id_nh4_4a), small_value )
qcw(:,k,id_nh4_5a) = MAX( qcw(:,k,id_nh4_5a), small_value )
qcw(:,k,id_nh4_6a) = MAX( qcw(:,k,id_nh4_6a), small_value )
#endif
#else
if (.not. inv_so2) then
qin(:,k,id_so2) = MAX( xso2(:,k), small_value )
endif
if (.not. inv_h2o2) then
qin(:,k,id_h2o2)= MAX( xh2o2(:,k), small_value )
endif
qin(:,k,id_so4) = MAX( xso4(:,k), small_value )
#endif
end do
#if (defined MODAL_AERO)
do n = 1, ntot_amode
m = lptr_so4_cw_amode(n)
l = m - loffset
if (l > 0) then
sflx(:)=0._r8
do k=1,pver
do i=1,ncol
sflx(i)=sflx(i)+dqdt_aqso4(i,k,l)*adv_mass(l)/mbar(i,k) &
*pdel(i,k)/gravit ! kg/m2/s
enddo
enddo
call outfld( trim(cnst_name_cw(m))//'AQSO4', sflx(:ncol), ncol, lchnk)
sflx(:)=0._r8
do k=1,pver
do i=1,ncol
sflx(i)=sflx(i)+dqdt_aqh2so4(i,k,l)*adv_mass(l)/mbar(i,k) &
*pdel(i,k)/gravit ! kg/m2/s
enddo
enddo
call outfld( trim(cnst_name_cw(m))//'AQH2SO4', sflx(:ncol), ncol, lchnk)
endif
end do
sflx(:)=0._r8
do k=1,pver
do i=1,ncol
sflx(i)=sflx(i)+dqdt_aqhprxn(i,k)*specmw_so4_amode/mbar(i,k) &
*pdel(i,k)/gravit ! kg SO4 /m2/s
enddo
enddo
call outfld( 'AQSO4_H2O2', sflx(:ncol), ncol, lchnk)
sflx(:)=0._r8
do k=1,pver
do i=1,ncol
sflx(i)=sflx(i)+dqdt_aqo3rxn(i,k)*specmw_so4_amode/mbar(i,k) &
*pdel(i,k)/gravit ! kg SO4 /m2/s
enddo
enddo
call outfld( 'AQSO4_O3', sflx(:ncol), ncol, lchnk)
xphlwc(:,:) = 0._r8
do k = 1,pver
do i = 1,ncol
if (cldfrc(i,k) >= 1.0e-5_r8) then
if( lwc(i,k) >= 1.0e-8_r8 ) then
xphlwc(i,k) = -1._r8*log10(xph(i,k)) * lwc(i,k)
endif
endif
enddo
enddo
call outfld( 'XPH_LWC', xphlwc(:ncol,:), ncol ,lchnk )
#endif
end subroutine SETSOX
end module MO_SETSOX