module sulchem 2,4
!-----------------------------------------------------------------------
! Purpose:
! Contains the sulfur cycle chemistry codes.
!
! Author:
! Chemistry code is from Mary Barth.
! Module coded by B. Eaton.
!-----------------------------------------------------------------------
use shr_kind_mod
, only: r8 => shr_kind_r8
use ppgrid, only: pcols, pver
use pmgrid
use abortutils, only: endrun
implicit none
save
private
public :: &
inisulchem, &! Initialize sulfur cycle chemistry.
chemwdepdr ! Chemistry and wet deposition driver.
! private module data
integer, parameter :: DBLKIND = selected_real_kind(12)
integer, parameter ::&
nz=18, &
nza=8, &
nalb=2, &
nz1=nz-1, &
nza1=nza-1, &
nz2=51, &
nsiz =nz *nza*nalb, &
nsiz2=nz2*nza*nalb
real(r8) ::&
jper2(nsiz2), &! j values interpolated to new ht grid
delz(nz1), &! Difference between tabl. z's
delang(nza1), &! Difference between tabl. sec(zen.ang)
delalb, &! Difference in tabl. albedo = .5-.2
pie, &! 3.14159...
dayspy ! Number of days per 1 year
real(r8), parameter :: mwa2so2 = 28.966/64. ! ratio molecular wt of air to so2
real(r8), parameter :: mwa2h2o2 = 28.966/34. ! ratio molecular wt of air to h2o2
real(r8), parameter, dimension(nsiz) :: jper = &! Tabulated j values
(/ -11.63010025,-11.53059959,-11.47570038,-11.44130039,-11.40859985,-11.40690041,-11.42399979,-11.46020031, &
-11.48509979,-11.49100018,-11.46759987,-11.41129971,-11.32050037,-11.20450020,-10.96660042,-10.65079975, &
-10.31389999,-10.15680027,-11.93789959,-11.81239986,-11.73900032,-11.68900013,-11.62989998,-11.60620022, &
-11.60490036,-11.61820030,-11.62650013,-11.62040043,-11.58790016,-11.52509975,-11.43029976,-11.31110001, &
-11.06599998,-10.75179958,-10.39190006,-10.20240021,-12.20409966,-12.06019974,-11.97229958,-11.90939999, &
-11.82750034,-11.78409958,-11.76560020,-11.75730038,-11.74950027,-11.73139954,-11.68959999,-11.62010002, &
-11.52070045,-11.39859962,-11.14780045,-10.83170033,-10.46100044,-10.24470043,-12.51459980,-12.35309982, &
-12.24969959,-12.17240047,-12.06379986,-11.99619961,-11.95629978,-11.92070007,-11.89220047,-11.85890007, &
-11.80539989,-11.72710037,-11.62139988,-11.49520016,-11.23849964,-10.91720009,-10.54150009,-10.29699993, &
-13.14570045,-12.95919991,-12.83010006,-12.72649956,-12.56420040,-12.44439983,-12.35589981,-12.25739956, &
-12.18089962,-12.11240005,-12.03289986,-11.93480015,-11.81389999,-11.67720032,-11.40919971,-11.07330036, &
-10.69830036,-10.41100025,-14.34010029,-14.15550041,-14.01509953,-13.88949966,-13.65480042,-13.43579960, &
-13.23670006,-12.97910023,-12.77519989,-12.61489964,-12.47119999,-12.32730007,-12.17119980,-12.00699997, &
-11.70969963,-11.35060024,-10.96660042,-10.65600014,-15.41930008,-15.26990032,-15.16590023,-15.07689953, &
-14.89920044,-14.68239975,-14.41380024,-13.96300030,-13.55609989,-13.23820019,-12.98649979,-12.77019978, &
-12.56400013,-12.36279964,-12.01889992,-11.63490009,-11.22630024,-10.92059994,-17.23080063,-17.09090042, &
-16.99939919,-16.92959976,-16.82430077,-16.74169922,-16.66230011,-16.49410057,-16.09469986,-15.37609959, &
-14.65410042,-14.07880020,-13.62609959,-13.26119995,-12.75020027,-12.25749969,-11.79109955,-11.45209980, &
-11.17109966,-11.14850044,-11.14080048,-11.14089966,-11.15730000,-11.18900013,-11.23019981,-11.29240036, &
-11.33609962,-11.35649967,-11.34700012,-11.30399990,-11.22620010,-11.12189960,-10.90159988,-10.60239983, &
-10.27840042,-10.12619972,-11.52420044,-11.47599983,-11.44919968,-11.43290043,-11.42099953,-11.42879963, &
-11.44970036,-11.48649979,-11.51080036,-11.51659966,-11.49520016,-11.44289970,-11.35809994,-11.24800014, &
-11.01669979,-10.71510029,-10.36550045,-10.18050003,-11.81890011,-11.75230026,-11.71059990,-11.68099976, &
-11.64519978,-11.63220024,-11.63479996,-11.64830017,-11.65499973,-11.64729977,-11.61489964,-11.55410004, &
-11.46290016,-11.34809971,-11.10849953,-10.80249977,-10.44029999,-10.22789955,-12.15250015,-12.06820011, &
-12.01109982,-11.96700001,-11.90400028,-11.86610031,-11.84650040,-11.83150005,-11.81620026,-11.79199982, &
-11.74639988,-11.67520046,-11.57610035,-11.45580006,-11.20800018,-10.89459991,-10.52560043,-10.28439999, &
-12.80440044,-12.69659996,-12.61540031,-12.54609966,-12.43089962,-12.34109974,-12.27270031,-12.19359970, &
-12.12870026,-12.06789970,-11.99429989,-11.90139961,-11.78509998,-11.65240002,-11.39009953,-11.05939960, &
-10.68850040,-10.40359974,-13.98359966,-13.87979984,-13.79139996,-13.70489979,-13.52700043,-13.34549999, &
-13.17129993,-12.93599987,-12.74400043,-12.59039974,-12.45110035,-12.31060028,-12.15719986,-11.99520016, &
-11.70090008,-11.34430027,-10.96230030,-10.65270042,-15.04539967,-14.97089958,-14.91469955,-14.86209965, &
-14.74139977,-14.57019997,-14.33740044,-13.92109966,-13.53120041,-13.22140026,-12.97410011,-12.76060009, &
-12.55630016,-12.35649967,-12.01439953,-11.63179970,-11.22420025,-10.91909981,-16.85479927,-16.78750038, &
-16.74150085,-16.70509911,-16.64769936,-16.59880066,-16.54520035,-16.40889931,-16.04310036,-15.35249996, &
-14.64319992,-14.07279968,-13.62240028,-13.25860023,-12.74860001,-12.25650024,-11.79049969,-11.45160007 /)
real(r8), parameter, dimension(nz) :: zz = &! Tabulated heights
(/ 0., 1., 2., 3., 5., 7., 9., 12., 15., 18., 21., 24., &
27., 30., 35., 40., 45., 50. /)
real(r8), parameter, dimension(nza) :: zasec = &! Tabulated sec(zen.ang)
(/ 1., 1.3, 1.6, 2., 3., 6., 12., 50. /)
real(r8), parameter, dimension(nalb) :: alb = &! Tabulated albedo = 0.2 and 0.5
(/ 0.2, 0.5 /)
real(r8), parameter, dimension(121) :: yield = &! yield of DMS+OH rxn going to SO2
(/ 0.76905, 0.76953, 0.77002, 0.77050, 0.77098, 0.77146, &
0.77195, 0.77243, 0.77291, 0.77340, 0.77388, 0.77436, &
0.77485, 0.77533, 0.77581, 0.77629, 0.77678, 0.77726, &
0.77774, 0.77822, 0.77870, 0.77918, 0.77967, 0.78015, &
0.78063, 0.78111, 0.78159, 0.78207, 0.78255, 0.78303, &
0.78351, 0.78399, 0.78447, 0.78495, 0.78543, 0.78591, &
0.78640, 0.78688, 0.78737, 0.78786, 0.78835, 0.78884, &
0.78934, 0.78984, 0.79035, 0.79086, 0.79137, 0.79190, &
0.79243, 0.79297, 0.79353, 0.79409, 0.79467, 0.79526, &
0.79587, 0.79650, 0.79715, 0.79783, 0.79853, 0.79927, &
0.80004, 0.80084, 0.80169, 0.80258, 0.80352, 0.80452, &
0.80558, 0.80671, 0.80790, 0.80918, 0.81055, 0.81200, &
0.81356, 0.81522, 0.81700, 0.81890, 0.82093, 0.82309, &
0.82540, 0.82786, 0.83047, 0.83325, 0.83618, 0.83928, &
0.84255, 0.84598, 0.84957, 0.85332, 0.85722, 0.86126, &
0.86543, 0.86973, 0.87412, 0.87861, 0.88316, 0.88777, &
0.89240, 0.89705, 0.90169, 0.90631, 0.91087, 0.91537, &
0.91978, 0.92409, 0.92829, 0.93236, 0.93630, 0.94009, &
0.94373, 0.94721, 0.95054, 0.95370, 0.95670, 0.95954, &
0.96223, 0.96476, 0.96714, 0.96938, 0.97148, 0.97344, &
0.97528 /)
real(r8), parameter, dimension(121) :: dmsrate = &! rate coeff of DMS+OH
(/ 0.21261E-10,0.21112E-10,0.20966E-10,0.20823E-10,0.20683E-10, &
0.20546E-10,0.20411E-10,0.20280E-10,0.20151E-10,0.20024E-10, &
0.19900E-10,0.19779E-10,0.19660E-10,0.19543E-10,0.19428E-10, &
0.19316E-10,0.19206E-10,0.19098E-10,0.18991E-10,0.18887E-10, &
0.18785E-10,0.18684E-10,0.18585E-10,0.18488E-10,0.18393E-10, &
0.18299E-10,0.18207E-10,0.18116E-10,0.18027E-10,0.17939E-10, &
0.17852E-10,0.17766E-10,0.17682E-10,0.17598E-10,0.17516E-10, &
0.17434E-10,0.17354E-10,0.17274E-10,0.17194E-10,0.17115E-10, &
0.17036E-10,0.16958E-10,0.16880E-10,0.16801E-10,0.16722E-10, &
0.16643E-10,0.16563E-10,0.16482E-10,0.16401E-10,0.16317E-10, &
0.16233E-10,0.16146E-10,0.16057E-10,0.15966E-10,0.15871E-10, &
0.15774E-10,0.15673E-10,0.15567E-10,0.15458E-10,0.15343E-10, &
0.15223E-10,0.15097E-10,0.14965E-10,0.14826E-10,0.14680E-10, &
0.14526E-10,0.14365E-10,0.14195E-10,0.14016E-10,0.13828E-10, &
0.13632E-10,0.13426E-10,0.13211E-10,0.12987E-10,0.12754E-10, &
0.12514E-10,0.12265E-10,0.12009E-10,0.11747E-10,0.11479E-10, &
0.11207E-10,0.10932E-10,0.10655E-10,0.10376E-10,0.10098E-10, &
0.98217E-11,0.95482E-11,0.92787E-11,0.90145E-11,0.87565E-11, &
0.85057E-11,0.82630E-11,0.80289E-11,0.78042E-11,0.75893E-11, &
0.73845E-11,0.71899E-11,0.70058E-11,0.68321E-11,0.66687E-11, &
0.65155E-11,0.63722E-11,0.62384E-11,0.61140E-11,0.59985E-11, &
0.58915E-11,0.57926E-11,0.57013E-11,0.56173E-11,0.55402E-11, &
0.54694E-11,0.54047E-11,0.53456E-11,0.52917E-11,0.52426E-11, &
0.51981E-11,0.51578E-11,0.51214E-11,0.50886E-11,0.50591E-11, &
0.50327E-11 /)
!##############################################################################
contains
!##############################################################################
subroutine inisulchem() 1
!-----------------------------------------------------------------------
! Purpose:
! Initialize sulfur cycle chemistry module.
!
! Author: B. Eaton
!-----------------------------------------------------------------------
implicit none
! Local variables
integer ::&
i, l, ik, k
!-----------------------------------------------------------------------
! initialize variables for jh2o2 interpolation
! interpolate jval data onto 1 km height grid
do l = 1, nza*nalb
do ik = 1, nz2-1
k = 2
11 continue
if( zz(k) .gt. ik-1 ) then
jper2((l-1)*nz2+ik) = jper((l-1)*nz+k-1) &
+ ((ik-1)-zz(k-1))/(zz(k)-zz(k-1)) &
* (jper((l-1)*nz+k)-jper((l-1)*nz+k-1))
else
k = k + 1
go to 11
endif
end do
jper2(l*nz2) = jper(l*nz)
end do
do i = 1,nza1
delang(i) = 1. / (zasec(i+1) - zasec(i))
end do
delalb = 1. / (alb(2) - alb(1))
! Constants from radiation routines needed for dicor()
dayspy = 365.
pie = 4.*atan(1.)
end subroutine inisulchem
!##############################################################################
subroutine chemwdepdr( lchnk, nstep, lat, clat, lon, calday, dtime, & 1,44
pmid, pdel, zm, t, q, &
cldwat, cldf, concld, cldv, cmfdqr, &
prain, evapr, cme, rain, icwmr1, &
chtr, fracis, ncol )
!-----------------------------------------------------------------------
! Purpose:
! Sulfur chemistry and wet deposition driver. The prognostic species
! are DMS, SO2, SO4, and H2O2.
!
! Author: B. Eaton
!-----------------------------------------------------------------------
use acbnd, only: acbndget
use history, only: outfld
use wetdep, only: wetdepa, wetdepg
#ifdef SCYC_MASSBGT
use massbgt, only: gamwet
#endif
implicit none
integer, intent(in) :: lchnk
integer, intent(in) :: nstep
integer, intent(in) :: lat(pcols) ! latitude indices
integer, intent(in) :: lon(pcols) ! longitude indices
integer, intent(in) :: ncol
real(r8), intent(in) ::&
clat(pcols), &
calday, &
dtime, &! time step
pmid(pcols,pver), &! pressure at layer midpoints
pdel(pcols,pver), &! delta-p between layer interfaces
zm(pcols,pver), &! height of layer midpoints
t(pcols,pver), &! air temperature (K)
q(pcols,pver), &! specific humidity (kg/kg)
cldwat(pcols,pver), &! cloud water mixing ratio
cldf(pcols,pver), &! total cloud fraction
concld(pcols,pver), &! convective cloud fraction
cldv(pcols,pver), &! cloudy volume which is occupied by rain or cloud water
cmfdqr(pcols,pver), &! dq/dt due to moist convective rainout
prain(pcols,pver), &! rate of conversion of condensate to precip
evapr(pcols,pver), &! rate of evaporation of falling precip
cme(pcols,pver), &! rate of cond-evap within the cloud
rain(pcols,pver), &! total precip mixing ratio
icwmr1(pcols,pver) ! in cloud water mixing ratio for zhang+hack scheme
! real(r8) chtr(pcols,pver,4) ! chemical tracer array
real(r8), intent(inout) :: chtr(pcols,pver,4) ! chemical tracer array
! obuf(*) ! history output buffer
real(r8), intent(out) ::&
! real(r8) &
fracis(pcols,pver,4) ! fraction of transported species that are isouble
! Local variables:
integer :: i, k, m
real(r8), dimension(pcols,pver) ::&
so2, &! so2 array
so4, &! so4 array
dms, &! dms
h2o2, &! h2o2
so2new, &! so2 array
so2tmp, &! temporary so2 array
so4new, &! so4 array
so4tmp, &! temporary so4 array
dmsnew, &! dms
dmssnk, &! total sink of dms
so2srcg, &! total source of so2
so2srcg2, &! dms + oh source of so2
so4src, &! total source of so4
so4srcg, &! gas source of so4
so4srca, &! aq. source of so4
h2o2src, &! total source of h2o2
h2o2snkg, &! gas sink of h2o2
h2o2snka ! aqueous sink of h2o2
real(r8), dimension(pcols,pver) ::&
h2o2new, &! h2o2
h2o2tmp, &! temporary h2o2
h2o23d, &! h2o2 -- from IMAGES
oh, &! oh -- from IMAGES
no3, &! no3 -- from IMAGES
ho2im, &! ho2 -- from IMAGES
o3, &! o3 -- from IMAGES
ho2, &! ho2 = ho2im*dicorfac
ph, &! ph of drops
hion ! Hydrogen ion concentration in drops
real(r8) ::&
dicorfac(pcols) ! factors to increase HO2 diurnal averages
integer ::&
indcp(pcols,pver), &! Indices of cloudy grid points
ncldypts(pver) ! Number of cloudy grid points
real(r8), dimension(pcols,pver) ::&
jh2o2, &! photolysis rate of H2O2
iscavt, &
scavt
real(DBLKIND) :: &
ekh2o2(pcols,pver), &! H2O2 Henry's Law coefficient
ekso2(pcols,pver) ! SO2 effective Henry's Law coefficient
real(r8) o3int(pcols), const
! The chemistry and wet dep schemes currently only use the in cloud water mixing
! ratio as a sum of icwmr1+icwmr2. The value supplied by MATCH is that sum,
! and icwmr2 is set to zero.
real(r8), dimension(pcols,pver) ::&
icwmr2 ! in cloud water mixing ratio for hack scheme
real(r8) :: molwt
real(r8) :: sol_fact
!----------------------------------------------------------------------c
do k = 1,pver
do i = 1,ncol
icwmr2(i,k) = 0.0
end do
end do
do k = 1,pver
do i = 1,ncol
so2(i,k) = chtr(i,k,1)
so4(i,k) = chtr(i,k,2)
dms(i,k) = chtr(i,k,3)
h2o2(i,k) = chtr(i,k,4)
end do
end do
do m = 1, 4
do k = 1, pver
do i = 1, ncol
fracis(i,k,m) = 1.
end do
end do
end do
call cldychmini( t, pmid, cldwat, rain, cldf, &
cldv, icwmr1, icwmr2, so2, so4, &
ph, hion, ekso2, ekh2o2, indcp, &
ncldypts, ncol )
call outfld( 'PH', ph, pcols, lchnk )
! Nominal concentrations
call acbndget( 'H2O2', ncol, lchnk, h2o23d )
call acbndget( 'O3', ncol, lchnk, o3 )
call acbndget( 'HO2', ncol, lchnk, ho2im ) ! 1.52e8
call acbndget( 'NO3', ncol, lchnk, no3 ) ! 1.64e7
call acbndget( 'OH', ncol, lchnk, oh ) ! 1.e6
! alter HO2 from diurnal average to instantaneous value
! calculate integrate o3 (dobson units)
call dicor( lat, calday, clat, ncol, dicorfac)
const = 6.02e23 / (10.* 28.97 * 2.7e16*9.8)
o3int = 0.
do k = 1, pver
do i = 1, ncol
ho2(i,k) = ho2im(i,k)*dicorfac(i)
o3int(i) = o3int(i) + o3(i,k)*const*pdel(i,k)
end do
end do
call outfld( 'O3INT', o3int, pcols, lchnk )
call outfld( 'O3', o3, pcols, lchnk ) !"mol/mol in oxid_bnd.nc"
call getj( clat, lat, lon, calday, zm, ncol, jh2o2 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Toggle the chemistry and wet deposition routines !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if( mod( nstep/2, 2 ) .eq. 0 ) then
call aqchem( so2, so4, h2o2, t, pmid, &
cldwat, rain, so2new, so4new, h2o2new, &
o3, cldf, so4srca, hion, h2o2snka, &
ekso2, ekh2o2, dtime, indcp, ncldypts, &
lat, lon, cldv, icwmr1, icwmr2, ncol )
so2tmp = so2new
so4tmp = so4new
h2o2tmp = h2o2new
call gaschem( so2tmp, so4tmp, dms, h2o2tmp, t, &
pmid, so2new, so4new, dmsnew, h2o2new, &
h2o23d, oh, no3, ho2, q, &
jh2o2, so4srcg, dmssnk, so2srcg, so2srcg2, &
h2o2src, h2o2snkg, dtime, lat, ncol )
#ifdef SCYC_MASSBGT
call endrun ('CHEMWDEPDR: gamwet expects scalar for actual argument lat')
call gamwet( 'dmssnk', lat, pdel, dmssnk )
call gamwet( 'so2srcg', lat, pdel, so2srcg )
call gamwet( 'so4srca', lat, pdel, so4srca )
call gamwet( 'so4srcg', lat, pdel, so4srcg )
#endif
! Wet deposition of SO2.
molwt = 64._r8
call wetdepg( lat, lon, t, pmid, q, pdel, &
cldf, concld, cmfdqr, prain, evapr, &
rain, cldwat, so2new, dtime, molwt, &
ekso2, scavt, iscavt, cldv, icwmr1, &
icwmr2, fracis(1,1,1), ncol )
do k = 1,pver
do i = 1,ncol
so2new(i,k) = so2new(i,k) + scavt(i,k)*dtime
end do
end do
call outfld( 'SO2WET', scavt, pcols, lchnk )
#ifdef SCYC_MASSBGT
call endrun ('CHEMWDEPDR: gamwet expects scalar for actual argument lat')
call gamwet( 'so2wet', lat, pdel, scavt )
#endif
! Wet deposition of SO4.
sol_fact = 0.3_r8
call wetdepa( lat, t, pmid, q, pdel, &
cldf, concld, cmfdqr, icwmr1, prain, cme, &
evapr, cldwat, so4new, dtime, &
scavt, iscavt, cldv, fracis(1,1,2), sol_fact, ncol )
do k = 1,pver
do i = 1,ncol
so4new(i,k) = so4new(i,k) + scavt(i,k)*dtime
end do
end do
call outfld( 'SO4WET', scavt, pcols, lchnk )
#ifdef SCYC_MASSBGT
call endrun ('CHEMWDEPDR: gamwet expects scalar for actual argument lat')
call gamwet( 'so4wet', lat, pdel, scavt )
#endif
! Wet deposition of H2O2.
molwt = 34._r8
call wetdepg( lat, lon, t, pmid, q, pdel, &
cldf, concld, cmfdqr, prain, evapr, &
rain, cldwat, h2o2new, dtime, molwt, &
ekh2o2, scavt, iscavt, cldv, icwmr1, &
icwmr2, fracis(1,1,4), ncol )
do k = 1,pver
do i = 1,ncol
h2o2new(i,k) = h2o2new(i,k) + scavt(i,k)*dtime
end do
end do
call outfld( 'H2O2WET', scavt, pcols, lchnk )
!###############################
else !# toggle wetdep and chemistry #
!###############################
! Wet deposition of SO2.
molwt = 64._r8
call wetdepg( lat, lon, t, pmid, q, pdel, &
cldf, concld, cmfdqr, prain, evapr, &
rain, cldwat, so2, dtime, molwt, &
ekso2, scavt, iscavt, cldv, icwmr1, &
icwmr2, fracis(1,1,1), ncol )
do k = 1,pver
do i = 1,ncol
so2(i,k) = so2(i,k) + scavt(i,k)*dtime
end do
end do
call outfld( 'SO2WET', scavt, pcols, lchnk )
#ifdef SCYC_MASSBGT
call endrun ('CHEMWDEPDR: gamwet expects scalar for actual argument lat')
call gamwet( 'so2wet', lat, pdel, scavt )
#endif
! Wet deposition of SO4.
sol_fact = 0.3_r8
call wetdepa( lat, t, pmid, q, pdel, &
cldf, concld, cmfdqr, icwmr1, prain, cme, &
evapr, cldwat, so4, dtime, &
scavt, iscavt, cldv, fracis(1,1,2), sol_fact, ncol )
do k = 1,pver
do i = 1,ncol
so4(i,k) = so4(i,k) + scavt(i,k)*dtime
end do
end do
call outfld( 'SO4WET', scavt, pcols, lchnk )
#ifdef SCYC_MASSBGT
call endrun ('CHEMWDEPDR: gamwet expects scalar for actual argument lat')
call gamwet( 'so4wet', lat, pdel, scavt )
#endif
! Wet deposition of H2O2.
molwt = 34._r8
call wetdepg( lat, lon, t, pmid, q, pdel, &
cldf, concld, cmfdqr, prain, evapr, &
rain, cldwat, h2o2, dtime, molwt, &
ekh2o2, scavt, iscavt, cldv, icwmr1, &
icwmr2, fracis(1,1,4), ncol )
do k = 1,pver
do i = 1,ncol
h2o2(i,k) = h2o2(i,k) + scavt(i,k)*dtime
end do
end do
call outfld( 'H2O2WET', scavt, pcols, lchnk )
call aqchem( so2, so4, h2o2, t, pmid, &
cldwat, rain, so2new, so4new, h2o2new, &
o3, cldf, so4srca, hion, h2o2snka, &
ekso2, ekh2o2, dtime, indcp, ncldypts, &
lat, lon, cldv, icwmr1, icwmr2, ncol )
so2tmp = so2new
so4tmp = so4new
h2o2tmp = h2o2new
call gaschem( so2tmp, so4tmp, dms, h2o2tmp, t, &
pmid, so2new, so4new, dmsnew, h2o2new, &
h2o23d, oh, no3, ho2, q, &
jh2o2, so4srcg, dmssnk, so2srcg, so2srcg2, &
h2o2src, h2o2snkg, dtime, lat, ncol )
#ifdef SCYC_MASSBGT
call endrun ('CHEMWDEPDR: gamwet expects scalar for actual argument lat')
call gamwet( 'dmssnk', lat, pdel, dmssnk )
call gamwet( 'so2srcg', lat, pdel, so2srcg )
call gamwet( 'so4srca', lat, pdel, so4srca )
call gamwet( 'so4srcg', lat, pdel, so4srcg )
#endif
endif ! End of toggle
do k = 1, pver
do i = 1, ncol
chtr(i,k,1) = so2new(i,k)
chtr(i,k,2) = so4new(i,k)
chtr(i,k,3) = dmsnew(i,k)
chtr(i,k,4) = h2o2new(i,k)
so4src(i,k) = so4srca(i,k) + so4srcg(i,k)
end do
end do
call outfld( 'DMSSNK ', dmssnk, pcols, lchnk )
call outfld( 'SO2SRCG ', so2srcg, pcols, lchnk )
call outfld( 'SO2SRCG2', so2srcg2, pcols, lchnk )
call outfld( 'SO4SRC ', so4src, pcols, lchnk )
call outfld( 'SO4SRCG ', so4srcg, pcols, lchnk )
call outfld( 'H2O2SRC ', h2o2src, pcols, lchnk )
call outfld( 'H2O2SNKA', h2o2snka, pcols, lchnk )
call outfld( 'H2O2SNKG', h2o2snkg, pcols, lchnk )
end subroutine chemwdepdr
!##############################################################################
subroutine aqchem( so2, so4, h2o2, tair, pres, & 2
qc, qr, so2new, so4new, h2o2new, &
o3, cldf, so4srca, hion, h2o2snka, &
ekso2, ekh2o2, ztodt, ind, ncpts, &
lat, lon, cldv, icwmr1, icwmr2, ncol )
!-----------------------------------------------------------------------
! Purpose:
! Calculate the sources and sinks from the aqueous chemistry:
! S(IV) + H2O2 --> SO4
! S(IV) + O3 --> SO4
!
! where concentrations are concentrations in aqueous phase.
!
! Author: M. Barth
!-----------------------------------------------------------------------
implicit none
integer, intent(in) ::&
lat(pcols), &! latitude indices
lon(pcols), &! longititude indices
ind(pcols,pver), &! indices for cloudy grid points
ncpts(pver), &! number of cloudy grid points
ncol ! number of atmospheric columns used in chunk
real(r8), intent(in) ::&
so2(pcols,pver), &! so2 array
so4(pcols,pver), &! so4 from aqueous chemistry
h2o2(pcols,pver), &! h2o2 in kg/kg
o3(pcols,pver), &! o3
tair(pcols,pver), &! air temperature (K)
pres(pcols,pver), &! midpoint pressures
qc(pcols,pver), &! cloud water mixing ratio
qr(pcols,pver), &! rain mixing ratio
ztodt, &! time step
cldf(pcols,pver) ! fraction cloud cover
real(r8), intent(in) ::&
cldv(pcols,pver), &! estimate of local volume occupied by clouds
icwmr1(pcols,pver), &! in cloud water mixing ration for zhang scheme
icwmr2(pcols,pver) ! in cloud water mixing ration for hack scheme
real(r8), intent(inout) ::&
hion(pcols,pver) ! hydrogen ion concentration (mol/L)
real(DBLKIND) &
theff(pcols), &! ekso2 at cloudy grid points
sheff, &! temp array for theff
ekh2o2(pcols,pver), &! H2O2 Henry's Law coefficient
ekso2(pcols,pver) ! SO2 effective Henry's Law coefficient
real(r8), intent(out) ::&
so2new(pcols,pver), &! so2 array
so4new(pcols,pver), &! so4 from aqueous chemistry
h2o2new(pcols,pver), &! h2o2 from aqueous chemistry
so4srca(pcols,pver), &! so4 production rate from achem
h2o2snka(pcols,pver) ! aqueous sink of H2O2
! Local variables:
integer :: i, k, n
integer :: i2 ! indice for cloudy grid points
integer :: ilw ! number of cloudy grid points
integer, parameter :: ncyc = 20 ! number of times to loop through aqchem
! perform aqchem on smaller dt
real(r8), dimension(pcols,pver) ::&
cwat, &! cloud water
rwat, &! rain
twat ! total liquid water = cwat + rwat
real(r8) ::&
h2o2rate, &! rate of S(IV) + H2O2 reaction
o3rate, &! first order rate coefficient for S(IV) + O3
molarso4, &! total so4 in units mol SO4/L H2O
ahso3 ! (molar)**2 concentration of HSO3-
real(r8), dimension(pcols) ::&
tso2, &! so2 at cloudy grid points
tso4, &! so4 at cloudy grid points
th2o2, &! h2o2 at cloudy grid points
to3, &! o3 at cloudy grid points
ttwat, &! twat at cloudy grid points
temp, &! tair at cloudy grid points
tprs, &! pres at cloudy grid points
tso2n, &! so2new at cloudy grid points
tso4n, &! so4new at cloudy grid points
th2o2n, &! h2o2new at cloudy grid points
tso4src, &! so4srca at cloudy grid points
th2o2snk, &! h2o2snk at cloudy grid points
thion, &! hion at cloudy grid points
tkh2o2 ! ekh2o2 at cloudy grid points
real(r8) ::&
molso2, &! total so2 in units mol SO2/mol air
molh2o2, &! total h2o2 in units mol H2O2/mol air
molarh2o2, &! total h2o2 in units mol H2O2/L H2O
molo3, &! total o3 in units mol O3/mol air
eko3, &! O3 Henry's Law coefficient
so4src, &! source of SO4 from SO2 + H2O2
so2snk, &! sink of SO2 from SO2 + H2O2
h2o2snk, &! sink of H2O2 from SO2 + H2O2
siv, &! concentration of S(IV) (mol/L)
a1, a2, &! fraction of S(IV) that is HSO3, SO3
totsb, &! sum of sulfur at beg of aqchem for one point
totse, &! sum of sulfur at end of aqchem for one point
weight, &
tc ! temperature in deg C
real(r8) ::&
dtchem, &! chemistry timestep
patm, &! pressure (atm)
fa1, &! = Khs*K1s/[H+]
fa2, &! = Khs*K1s*K2s/[H+]**2
khs, &! Henry's Law constant for SO2
khsk1s, &! Khs * K1s
kh12s, &! Khs * K1s * K2s
kh2o2, &! rate constant for H2O2 + HSO3-
ko3hs, &! rate constant for O3 + HSO3-
ko3so3, &! rate constant for O3 + SO3=
qso2, &! temp array for so2
qso4, &! temp array for so4
qh2o2, &! temp array for h2o2
ahion, &! temp array for ahion
molrate, &! rate of rxn in mol/mol-s
tmprate, &! temp rate of rxn
sumrate, &! h2o2rate + o3rate
omsm
! Function:
!-drb real(r8) e298, dhr, tt, ek
real(r8) e298, dhr, tt
real (DBLKIND) ek
ek(e298, dhr, tt) = e298*exp(dhr*(1./tt - 1./298.))
!----------------------------------------------------------------------
! call t_startf ('aqchem')
dtchem = ztodt/ncyc
omsm = 1. - 1.e-10
do k=1,pver
do i=1,ncol
! Determine cloud water and rain water mixing ratios
tc = tair(i,k) - 273.16
weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
cwat(i,k) = (qc(i,k)/max(cldf(i,k), 0.01_r8) + &
! add suspended water from convection and do aqchem in only the liquid phase
icwmr1(i,k) + icwmr2(i,k)) * (1.-weight)
if(tair(i,k) .gt. 273.16) then
weight = 0. ! fraction of condensate that is ice
else
weight = 1.
endif
rwat(i,k) = qr(i,k)/max(cldv(i,k), 0.01_r8) &
*(1.-weight) ! aqchem in only the liquid phase
! Sum cwat and rwat
twat(i,k) = cwat(i,k) + rwat(i,k)
! if(twat(i,k) .gt. 0.5e-3) then
! tc = tair(i,k) - 273.16
! weight = max(0._r8,min(-tc*0.05,1.0_r8))
! write(6,'(a,3i4,5e12.5)') 'TWAT',i,lat,k, twat(i,k), cwat(i,k),
! _ rwat(i,k), icwmr1(i,k)*(1.-weight), icwmr2(i,k)*(1.-weight)
! write (6,*) ' rwat stuff, qr, cldv, cldf ', qr(i,k), cldv(i,k),
! $ cldf(i,k)
! endif
end do
end do
! Initialize arrays
do k=1,pver
do i=1,ncol
so2new(i,k) = so2(i,k) !keep concentration
so4new(i,k) = so4(i,k)
h2o2new(i,k) = h2o2(i,k)
so4srca(i,k) = 0.
h2o2snka(i,k) = 0.
end do
end do
! Initialize cloudy grid point arrays
! write(6,*) ' Beginning of aqchem'
do k=1,pver
ilw = ncpts(k)
do i2=1,ilw
tso2(i2) = so2(ind(i2,k),k)
tso4(i2) = so4(ind(i2,k),k)
th2o2(i2) = h2o2(ind(i2,k),k)
to3(i2) = o3(ind(i2,k),k)
ttwat(i2) = twat(ind(i2,k),k)
temp(i2) = tair(ind(i2,k),k)
tprs(i2) = pres(ind(i2,k),k)
thion(i2) = hion(ind(i2,k),k)
theff(i2) = ekso2(ind(i2,k),k)
tkh2o2(i2) = ekh2o2(ind(i2,k),k)
tso4src(i2) = 0.
th2o2snk(i2) = 0.
! write(6,'(3i5,4e12.5,f8.3,a)') ind(i2,k), lat, k, ttwat(i2), &
! thion(i2), tso4(i2), tso2(i2), -log10(thion(i2)), ' begin'
end do
! Aqueous chemistry calculations begin here
do i2=1,ilw
patm = tprs(i2)/101325.
khs = ek(1.23_r8,3120._r8,temp(i2))
khsk1s = khs * ek(1.3e-2_r8,2015._r8,temp(i2))
kh12s = khsk1s * ek(6.31e-8_r8,1505._r8,temp(i2))
kh2o2 = 6.2534e14*exp(-4751./temp(i2))
molo3 = to3(i2)
eko3 = ek(1.15e-2_r8,2560._r8,temp(i2))
ko3hs = 4.28e13*exp(-5533./temp(i2))
ko3so3= 7.436e16*exp(-5280./temp(i2))
! if(lat .eq. 8 .and. k .eq. 14 .and. i2 .gt. 8) then
! write(6,'(i3,7e11.4)') i2, tso2(i2), th2o2(i2), thion(i2), &
! theff(i2), temp(i2), patm, molo3
! write(6,'(3x,e11.4,f8.1)') ttwat(i2), dtchem
! endif
!cdir expand=ncyc
do n=1,ncyc !Do aqchem on a smaller dt 5/5/97
if(n .eq. 1) then
qso2 = tso2(i2)
qso4 = tso4(i2)
qh2o2 = th2o2(i2)
ahion = thion(i2)
sheff = theff(i2)
endif
fa1 = khsk1s/max(ahion,1.e-30_r8)
fa2 = kh12s/max((ahion*ahion),1.e-30_r8)
! Calculate rate of HSO3 + H2O2 reaction
molso2 = qso2 * mwa2so2
molh2o2 = qh2o2 * mwa2h2o2
ahso3 = fa1 * patm*molso2
molarh2o2= tkh2o2(i2)*patm*molh2o2
h2o2rate = kh2o2 * ahion * ahso3 * molarh2o2/(1. + 13.*ahion)
! Convert to mol/(mol-s)
molrate = h2o2rate*ttwat(i2)*28.97e-3
! check for conservation in kg/(kg-s)
tmprate = mwa2h2o2*min(molrate/mwa2h2o2, qh2o2/dtchem)
tmprate = mwa2so2*min(tmprate/mwa2so2, qso2/dtchem)
h2o2rate = tmprate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3)
! Calculate rate of o3 oxidation
siv = sheff * molso2 * patm ![S(IV)]
a1 = fa1 / sheff
a2 = fa2 / sheff
o3rate = (ko3hs*a1 + ko3so3*a2) * eko3*molo3 * patm*siv
! Convert to mol/(mol-s)
molrate = o3rate*ttwat(i2)*28.97e-3
! check for conservation in kg/(kg-s)
tmprate = mwa2so2*min(molrate/mwa2so2, qso2/dtchem)
o3rate = tmprate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3)
! check total rate for conservation
sumrate = h2o2rate + o3rate
! if(lat .eq. 8 .and. k .eq. 14 .and. i2 .gt. 8) then
! write(6,'(i3,5e11.4)') n, sumrate, h2o2rate, o3rate, &
! qso2/dtchem, qh2o2/dtchem
! endif
molrate = sumrate*ttwat(i2)*28.97e-3
molrate = mwa2so2*min(molrate/mwa2so2, qso2/dtchem)
tmprate = molrate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3)
h2o2rate = h2o2rate/max(sumrate, 1.e-30_r8) * tmprate
o3rate = o3rate/max(sumrate, 1.e-30_r8) * tmprate
! Update so2, so4, h2o2
so2snk = (h2o2rate + o3rate)*ttwat(i2)*1.e-3 * 64.
h2o2snk = h2o2rate*ttwat(i2)*1.e-3 * 34.
! if(qh2o2 .lt. omsm*h2o2snk*dtchem) then
! write(6,'(a,2e12.5,3i4)') 'AQCHEM qh2o2 < h2o2snk*dt ', qh2o2, &
! h2o2snk*dtchem, lat, ind(i2,k), k
! endif
! if(qso2 .lt. omsm*so2snk*dtchem) then
! write(6,'(a,2e12.5,3i4)') 'AQCHEM qso2 < so2snk*dt ', qso2, &
! so2snk*dtchem, lat, ind(i2,k), k
! endif
so4src = 1.5*so2snk
qso2 = max(qso2 - so2snk*dtchem, 0._r8)
qso4 = qso4 + so4src*dtchem
qh2o2 = max(qh2o2 - h2o2snk*dtchem, 0._r8)
tso4src(i2) = tso4src(i2) + so4src
th2o2snk(i2) = th2o2snk(i2) + h2o2snk
! calculate new ahion, sheff
molso2 = qso2 * mwa2so2
molarso4 = qso4/(max(ttwat(i2), 1.e-30_r8)*96.e-3)
ahso3 = khsk1s * molso2 * patm
!++pjr
! rewrote this expression so we dont take the sqrt(max(a,1.e-60))
! but rather max(sqrt(a),1.e-30)
! ahion = (molarso4 + sqrt(max(molarso4*molarso4 + 4.*ahso3, &
! 1.e-60_r8)))*0.5
!--pjr
ahion = (molarso4 + sqrt(max(molarso4*molarso4 + 4._r8*ahso3, &
0._r8)))*0.5
ahion = max(ahion,1.e-30_r8)
sheff = khs + khsk1s/ahion + kh12s/max((ahion*ahion),1.e-30_r8)
end do !ncyc
tso2n(i2) = qso2
tso4n(i2) = qso4
th2o2n(i2)= qh2o2
thion(i2) = ahion
theff(i2) = sheff
tso4src(i2) = tso4src(i2)/ncyc
th2o2snk(i2) = th2o2snk(i2)/ncyc
end do
do i2=1,ilw
! so2new is modified for just the cloudy portion of the grid cell. Now
! adjust it for the entire grid cell.
so2new(ind(i2,k),k) = cldv(ind(i2,k),k)*tso2n(i2) + &
(1.-cldv(ind(i2,k),k))*tso2(i2)
so4new(ind(i2,k),k) = cldv(ind(i2,k),k)*tso4n(i2) + &
(1.-cldv(ind(i2,k),k))*tso4(i2)
so4srca(ind(i2,k),k) = cldv(ind(i2,k),k)*tso4src(i2)
h2o2snka(ind(i2,k),k)= cldv(ind(i2,k),k)*th2o2snk(i2)
h2o2new(ind(i2,k),k) = cldv(ind(i2,k),k)*th2o2n(i2) + &
(1.-cldv(ind(i2,k),k))*th2o2(i2)
hion(ind(i2,k),k) = thion(i2)
ekso2(ind(i2,k),k) = theff(i2)
! write(6,'(3i5,4e12.5,f8.3,a)') ind(i2,k), lat, k, ttwat(i2), &
! thion(i2), tso4(i2), tso2(i2), -log10(thion(i2)), ' end'
end do
do i=1,ncol
totsb = 0.5*so2(i,k) + 1./3. * so4(i,k)
totse = 0.5*so2new(i,k) + 1./3. * so4new(i,k)
! if(abs(totsb-totse) .gt. .001*totsb) then
! write(6,*) 'Total sulfur not conserved in aqueous chemistry'
! write(6,*) i, k, totsb, totse
! stop
! endif
end do
end do
! call t_stopf ('aqchem')
return
end subroutine aqchem
!##############################################################################
subroutine cldychmini( tair, pres, qc, qr, cldf, & 1
cldv, icwmr1, icwmr2, so2, so4, &
ph, hion, ekso2, ekh2o2, ind, &
ncpts, ncol )
!-----------------------------------------------------------------------
! Purpose:
! Calculate the pH, solubility constants of SO2 and H2O2, and the
! location of the cloudy grid points
!
! Author: M. Barth
!-----------------------------------------------------------------------
implicit none
real(r8), intent(in) ::&
tair(pcols,pver), &! air temperature (K)
pres(pcols,pver), &! midpoint pressures
qc(pcols,pver), &! cloud water mixing ratio
qr(pcols,pver), &! rain mixing ratio
cldf(pcols,pver), &! fraction cloud cover
cldv(pcols,pver), &! cloudy volume which is occupied by rain or cloud water
so2(pcols,pver), &! so2 array
so4(pcols,pver), &! so4 array
icwmr1(pcols,pver), &! in cloud water mixing ration for zhang scheme
icwmr2(pcols,pver) ! in cloud water mixing ration for hack scheme
real(r8), intent(out) ::&
ph(pcols,pver), &! pH of drops
hion(pcols,pver) ! hydrogen ion concentration (mol/L)
real(DBLKIND) &
theff(pcols), &! ekso2 at cloudy grid points
ekh2o2(pcols,pver), &! H2O2 Henry's Law coefficient
ekso2(pcols,pver) ! SO2 effective Henry's Law coefficient
integer, intent(in) :: &
ncol ! number of atmospheric columns used in chunk
integer, intent(out) ::&
ind(pcols,pver), &! indices for cloudy grid points
ncpts(pver) ! number of cloudy grid points
! Local variables:
integer :: &
i, k, &
i2, &! index for cloudy grid points
ilw ! number of cloudy grid points
real(r8), dimension(pcols,pver) ::&
cwat, &! cloud water
rwat, &! rain
twat ! total liquid water = cwat + rwat
real(r8), dimension(pcols) ::&
tso2, &! so2 at cloudy grid points
tso4, &! so4 at cloudy grid points
ttwat, &! twat at cloudy grid points
temp, &! tair at cloudy grid points
tprs, &! pres at cloudy grid points
tph, &! ph at cloudy grid points
thion ! hion at cloudy grid points
real(r8) :: &
molarso4, &! total so4 in units mol SO4/L H2O
ahso3, &! (molar)**2 concentration of HSO3-
molso2, &! so2 in mol/mol units
weight, &
tc ! temperature in deg C
! Function:
!-drb real(r8) e298, dhr, tt, ek
real(r8) e298, dhr, tt
real(DBLKIND) ek
ek(e298, dhr, tt) = e298*exp(dhr*(1./tt - 1./298.))
!----------------------------------------------------------------------
do k=1,pver
do i=1,ncol
ekh2o2(i,k) = ek(7.4e4_r8, 6621._r8, tair(i,k))
ekso2(i,k) = ek(1.23_r8, 3120._r8, tair(i,k))
! Determine cloud water and rain water mixing ratios
tc = tair(i,k) - 273.16
weight = max(0._r8,min(-tc*0.05_r8,1.0_r8)) ! fraction of condensate that is ice
cwat(i,k) = (qc(i,k)/max(cldf(i,k), 0.00001_r8) + &
! add suspended water from convection and do aqchem in only the liquid phase
icwmr1(i,k) + icwmr2(i,k)) * (1.-weight)
if(tair(i,k) .gt. 273.16) then
weight = 0. ! fraction of condensate that is ice
else
weight = 1.
endif
rwat(i,k) = qr(i,k)/max(cldv(i,k), 0.00001_r8) &
*(1.-weight) ! aqchem in only the liquid phase
! Sum cwat and rwat
twat(i,k) = cwat(i,k) + rwat(i,k)
end do
end do
! Initialize arrays
do k=1,pver
do i=1,ncol
ph(i,k) = 99.
hion(i,k) = 0.
ind(i,k) = 0
end do
end do
! Find which grid points have liquid water
do k=1,pver
ncpts(k) = 0
ilw = 0
do i=1,ncol
if(cldv(i,k) .ge. 0.02 .and. twat(i,k) .ge. 1.e-12) then
ilw = ilw + 1
ind(ilw,k) = i
endif
end do
ncpts(k) = ilw !assign number of cloudy grid points to array
do i2=1,ilw
tso2(i2) = so2(ind(i2,k),k)
tso4(i2) = so4(ind(i2,k),k)
ttwat(i2) = twat(ind(i2,k),k)
temp(i2) = tair(ind(i2,k),k)
tprs(i2) = pres(ind(i2,k),k)
! Set pH as a function of HSO3- and SO4=
! Say NH4+ = 1.5 SO4= and NO3- = 0.5 SO4=
! Thus, H+ = HSO3- + SO4= as done in ECHAM
molso2 = tso2(i2) * mwa2so2
molarso4 = tso4(i2)/(max(ttwat(i2), 1.e-30_r8)*96.e-3)
ahso3 = ek(1.23_r8,3120._r8,temp(i2)) * &
ek(1.3e-2_r8,2015._r8,temp(i2)) * &
molso2 * tprs(i2)/101325.
!++pjr
! rewrote this expression so we dont take the sqrt(max(a,1.e-50))
! but rather max(sqrt(a),1.e-25)
! thion(i2) = (molarso4 + sqrt(max(molarso4*molarso4 + 4.*ahso3, &
! 1.e-50_r8)))*0.5
!--pjr
thion(i2) = (molarso4 + &
sqrt(max(molarso4*molarso4 + 4._r8*ahso3,0._r8)))*0.5
thion(i2) = max(thion(i2),1.e-25_r8)
tph(i2) = -log10(thion(i2))
theff(i2) = ek(1.23_r8, 3120._r8, temp(i2)) * (1. + &
ek(1.3e-2_r8, 2015._r8, temp(i2))/thion(i2) * (1. + &
ek(6.31e-8_r8, 1505._r8, temp(i2))/thion(i2)))
end do
do i2=1,ilw
hion(ind(i2,k),k) = thion(i2)
ph(ind(i2,k),k) = tph(i2)
ekso2(ind(i2,k),k) = theff(i2)
end do
end do
end subroutine cldychmini
!##############################################################################
subroutine coszen( calday, dodiavg, clat, lat, lon, ncol, coszrs ) 1
!-----------------------------------------------------------------------
! Purpose:
! Compute solar distance factor and cosine solar zenith angle usi
! day value where a round day (such as 213.0) refers to 0z at
! Greenwich longitude.
!
! Use formulas from Paltridge, G.W. and C.M.R. Platt 1976: Radiative
! Processes in Meterology and Climatology, Elsevier Scientific
! Publishing Company, New York p. 57, p. 62,63.
!
! Author: CCM2
!-----------------------------------------------------------------------
implicit none
real(r8), intent(in) ::&
calday, &! Calendar day, including fraction
clat(pcols) ! latitudes (radians)
integer, intent(in) ::&
lat(pcols), &! latitude indices
lon(pcols), &! longitude indices
ncol ! number of atmospheric columns used in chunk
logical, intent(in) ::&
dodiavg ! true => do diurnal averaging
real(r8), intent(out) ::&
coszrs(pcols) ! Cosine solar zenith angle
! Local variables
integer :: i ! Longitude loop index
real(r8) ::&
phi, &! Greenwich calendar day + local time + long offset
theta, &! Earth orbit seasonal angle in radians
delta, &! Solar declination angle in radians
sinc, &! Sine of latitude
cosc, &! Cosine of latitude
sind, &! Sine of declination
cosd ! Cosine of declination
real(r8) ::&
frac, &! Daylight fraction
arg, &!
tsun, &! temporary term in diurnal averaging
coszrsu ! uniform cosine zenith solar angle
!-----------------------------------------------------------------------
theta = 2.*pie*calday/dayspy
! Solar declination in radians:
delta = .006918 - .399912*cos(theta) + .070257*sin(theta) - &
.006758*cos(2.*theta) + .000907*sin(2.*theta) - &
.002697*cos(3.*theta) + .001480*sin(3.*theta)
do i=1,ncol
! Compute local cosine solar zenith angle,
sinc = sin(clat(i))
sind = sin(delta)
cosc = cos(clat(i))
cosd = cos(delta)
! If using diurnal averaging, then compute the average local cosine solar
! zenith angle using formulas from paltridge and platt 1976 p. 57, p. 62,63.
if (dodiavg) then
arg = -(sinc/cosc)*(sind/cosd)
if (arg .lt. -1.) then
frac = 1.0
else if (arg .gt. 1.) then
frac = 0.0
else
frac = (1./pie)*acos(arg)
endif
tsun = pie*frac
if (tsun .gt. 0.) then
coszrsu = sinc*sind + (cosc*cosd*sin(tsun))/tsun
else
coszrsu = 0.0
endif
coszrs(i) = coszrsu
else ! No diurnal averaging
! Calday is the calender day for Greenwich, including fraction
! of day; the fraction of the day represents a local time at
! Greenwich; to adjust this to produce a true instantaneous time
! For other longitudes, we must correct for the local time change:
! local time based on the longitude and day of year
! then compute the local cosine solar zenith angle
phi = calday + (real(lon(i)-1)/real(plon))
coszrs(i) = sinc*sind - cosc*cosd*cos(2.*pie*phi)
end if
end do
end subroutine coszen
!##############################################################################
subroutine dicor( lat, calday, clat, ncol, corr ) 1
!-----------------------------------------------------------------------
! Purpose:
! Calculate a correction factor for the ho2 reaction rate to account for
! the diurnal variation of ho2
!
! Author: Mary Barth
!-----------------------------------------------------------------------
implicit none
integer, intent(in) :: lat(pcols) ! latitude indices
real(r8), intent(in) :: calday ! Calendar day, including fraction
real(r8), intent(in) :: clat(pcols) ! latitudes (radians)
integer, intent(in) :: ncol ! number of atmospheric columns used in chunk
real(r8), intent(out) :: corr(pcols) ! correction factor
!---------------------------Local variables-----------------------------
integer ic ! index in chunk
integer i ! Longitude loop index
real(r8) phi ! Greenwich calendar day + local time + long offset
real(r8) theta ! Earth orbit seasonal angle in radians
real(r8) delta ! Solar declination angle in radians
real(r8) sinc ! Sine of latitude
real(r8) cosc ! Cosine of latitude
real(r8) sind ! Sine of declination
real(r8) cosd ! Cosine of declination
real(r8) coszrs ! Cosine solar zenith angle
real(r8) xxx ! Work space
real(r8) yyy ! Work space
integer nc ! counter
real(r8) :: dicor_corr ! h2o reaction rate correction
integer ntimes ! number of points to evaluate the diurnal average
!-----------------------------------------------------------------------
!
do ic=1,ncol
! Compute solar distance factor and cosine solar zenith angle using
! day value where a round day (such as 213.0) refers to 0z at
! Greenwich longitude.
!
! Use formulas from Paltridge, G.W. and C.M.R. Platt 1976: Radiative
! Processes in Meterology and Climatology, Elsevier Scientific
! Publishing Company, New York p. 57, p. 62,63.
theta = 2.*pie*calday/dayspy
!
! Solar declination in radians:
!
delta = .006918 - .399912*cos(theta) + .070257*sin(theta) - &
.006758*cos(2.*theta) + .000907*sin(2.*theta) - &
.002697*cos(3.*theta) + .001480*sin(3.*theta)
! now calculate the solar zenith angle
! and some useful quantities for the whole day
!
! Compute local cosine solar zenith angle,
!
sinc = sin(clat(ic))
sind = sin(delta)
cosc = cos(clat(ic))
cosd = cos(delta)
!
! Calday is the calender day for Greenwich, including fraction
! of day;
! since all we care about in this calculation is the diurnal variation
! of some quantities we just increment calday over one day.
!
xxx = 0.
yyy = 0.
nc = 0
ntimes = 48
do i=1,ntimes
phi = calday + (real(i)-1)/real(ntimes)
coszrs = sinc*sind - cosc*cosd*cos(2.*pie*phi)
if (coszrs.gt.0) then
nc = nc + 1
xxx = xxx + coszrs**2
yyy = yyy + coszrs
endif
end do
if (yyy.gt.0.) then
dicor_corr = ntimes*xxx/(yyy**2)
else
dicor_corr = 1.
endif
!
! Return correction
!
corr(ic) = dicor_corr
end do
end subroutine dicor
!##############################################################################
subroutine gaschem( so2, so4, dms, h2o2, tair, & 2,1
pres, so2new, so4new, dmsnew, h2o2new, &
h2o23d, oh, no3, ho2, h2o, &
jh2o2, so4srcg, dmssnk, so2srcg, so2srcg2, &
h2o2src, h2o2snkg, ztodt, lat, ncol )
!-----------------------------------------------------------------------
! Purpose:
! Calculate the sources and sinks from the gas chemistry:
! DMS + OH --> a*SO2 + (1-a)*MSA
! DMS + NO3 --> SO2
! SO2 + OH + M --> SO4 + M
! HO2 + HO2 --> H2O2
! H2O2+ OH --> HO2 + H2O
! H2O2+ hv --> 2OH
!
! where a is the yield of SO2 from DMS oxidation, MSA is methane sulfonic
! acid which is not carried in this model, and M is an air molecule and
! is accounted for in the rate coefficient expression.
!
! Author: M. Barth
!-----------------------------------------------------------------------
implicit none
real(r8), intent(in) :: &
so2(pcols,pver), &! so2 array
so4(pcols,pver), &! so4 from gas chemistry
dms(pcols,pver), &! dms
h2o2(pcols,pver), &! h2o2
oh(pcols,pver), &! oh
no3(pcols,pver), &! no3
ho2(pcols,pver), &! ho2
h2o(pcols,pver), &! h2o
jh2o2(pcols,pver), &! photolysis rate of H2O2
tair(pcols,pver), &! air temperature (K)
pres(pcols,pver), &! midpoint pressures
ztodt, &! time step
h2o23d(pcols,pver) ! h2o2 -- 3d field from images
integer, intent(in) :: &
lat(pcols), &! latitude indices
ncol ! number of atmospheric columns used in chunk
real(r8), intent(out) :: &
so4new(pcols,pver), &! so4 from gas chemistry
dmsnew(pcols,pver), &! dms
h2o2new(pcols,pver), &! h2o2
so4srcg(pcols,pver), &
so2srcg2(pcols,pver)
real(r8), intent(inout) :: &
so2new(pcols,pver), &! so2 array
dmssnk(pcols,pver), &
so2srcg(pcols,pver), &
h2o2src(pcols,pver), &
h2o2snkg(pcols,pver)
! Local variables:
integer ::&
i, k, ik
logical :: found
real(r8), parameter ::& ! from NASA(1992)
rl300=3.e-31, &! rate constant at low pressure
an=3.3, &! constant: (300/Tair)**an
rh300=1.5e-12, &! rate constant at high pressure
am=0. ! constant: (300/Tair)**am
real(r8) ::&
air, &! concentration of air (molec/cm3)
rlt, &! used in calculation of SO2 rate coefficient
rht, &! used in calculation of SO2 rate coefficient
ax, bx, &! used in calculation of SO2 rate coefficient
psi, &! used in calculation of SO2 rate coefficient
so2rate, &! first order rate coefficient for SO2 + OH
so2snk, &
dmsrate1, &! first order rate coefficient for DMS + OH
dmsno3, &! first order rate coefficient for DMS + NO3
dmstot ! first order rate coefficient for both rxns
real(r8) rk2 ! rate coeff. in HO2 + HO2 rxn
real(r8) rk3 ! rate coeff. in HO2 + HO2 rxn
real(r8) fh2o ! rate coeff. in HO2 + HO2 rxn
real(r8) ho2rate ! total rate coeff for HO2 + HO2 rxn
real(r8) h2o2des ! rate for H2O2*OH
real(r8) h2o2phot ! rate jh2o2*H2O2
! real(r8) tauh2o2 !relaxation time for H2O2 generation
! real(r8) xh2o2 !zonally avgd H2O2 in kg/kg
!----------------------------------------------------------------------c
! Assume a 1.5 day relaxation time for H2O2 generation
! tauh2o2 = 129600. !seconds
! call t_startf ('gaschem')
found = .false.
do k=1,pver
do i=1,ncol
! SO2 + OH rate coefficient:
if(tair(i,k) .le. 0.) then
found = .true.
exit
endif
end do
end do
if (found) then
do k=1,pver
do i=1,ncol
! SO2 + OH rate coefficient:
if(tair(i,k) .le. 0.) then
write(6,*) 'TAIR bad ', i, k, lat(i), tair(i,k)
call endrun ('GASCHEM')
endif
end do
end do
endif
do k=1,pver
do i=1,ncol
so2rate = 0.
! determine concentration of air (molec/cm3) = rho(kg/m3)*Na/MWair
air = pres(i,k)/(287.*tair(i,k)) * 1.e-3/28.966 * 6.022e23
! air = 0.8*air !N2 concentration
! SO2 + OH rate coefficient:
rlt = rl300 * (300./tair(i,k))**an
rht = rh300 * (300./tair(i,k))**am
ax = rlt*air/rht
psi = 1./(1. + (log10(ax))**2)
bx = rlt*air/(1. + ax)
so2rate = bx * 0.6**psi !rate coef. in cm3/molec-s
! rewrite so2rate as first order rate coefficient
so2rate = so2rate * oh(i,k)
so2snk = so2(i,k)*(1.-exp(-so2rate*ztodt))
! Update sulfur species
so2new(i,k) = so2(i,k) - so2snk
so4new(i,k) = so4(i,k) + 1.5*so2snk
so4srcg(i,k) = 1.5*so2snk/ztodt
! Note: 1.5 factor accounts for different molecular weights of so4 and so2
! mw(so4)=96; mw(so2)=64; mw(so4)/mw(so2)=3/2=1.5
end do
end do
do k=1,pver
do i=1,ncol
! DMS oxidation
!c DMS + OH rate coefficient and SO2 yield:
dmstot = 0.
dmsrate1 = 0.
dmsno3 = 0.
! Get dmsrate from tabulated data
ik = int(tair(i,k)-195.)
ik = min(max(ik,1), 121)
if(dms(i,k) .ge. 1.e-30) then
dmsrate1 = dmsrate(ik)*oh(i,k)
! DMS + NO3
dmsno3 = 1.e-12*exp(-500.*(0.0033557 - 1./tair(i,k))) * &
no3(i,k)
! Total first order rate coeff for DMS
dmstot = dmsrate1 + dmsno3
endif
! Update sulfur species
dmssnk(i,k) = dms(i,k)*(1.-exp(-dmstot*ztodt)) / ztodt
dmsnew(i,k) = dms(i,k) - dmssnk(i,k)*ztodt
so2srcg(i,k) = 1.032258 * (yield(ik) * dms(i,k)*(1.-exp(-dmsrate1*ztodt)) + &
dms(i,k)*(1.-exp(-dmsno3*ztodt)) ) / ztodt
so2new(i,k) = so2new(i,k) + so2srcg(i,k)*ztodt
so2srcg2(i,k) = 1.032258 / ztodt * &
yield(ik) * dms(i,k)*(1.-exp(-dmsrate1*ztodt))
! Note: 1.032258 = 64/62 accounts for different MW of so2 and dms
end do
end do
do k=1,pver
do i=1,ncol
!-----------------------------------------------------------------------
!c Generate H2O2 via gas phase processes -- parameterizing chemistry
! xh2o2 = h2o23d(i,k) * (287.*tair(i,k)/pres(i,k))*34./6.022e20
! h2o2new(i,k) = h2o2(i,k) + ztodt/tauh2o2 * (xh2o2 - h2o2(i,k))
!
! Generate H2O2 via the HO2 + HO2 reaction
! HO2 concentrations are prescribed
! water vapor concentrations taken from q(i,k,1)
!
! determine concentration of air (molec/cm3) = rho(kg/m3)*Na/MWair
air = pres(i,k)/(287.*tair(i,k)) * 1.e-3/28.966 * 6.022e23
rk2 = 1.7e-12*exp( 600.*(1./tair(i,k) - 0.0033557))
rk3 = 4.9e-32*exp(1000.*(1./tair(i,k) - 0.0033557))
fh2o = 1.4e-21*exp(2200./tair(i,k))
ho2rate = (rk2 + rk3*air)*(1. + fh2o*h2o(i,k)) * &
ho2(i,k)*ho2(i,k)
! xh2o2 = h2o23d(i,k) * (287.*tair(i,k)/pres(i,k))*34./6.022e20
! Gas-phase destruction of H2O2 occurs via:
! H2O2 + hv --> 2OH
! H2O2 + OH --> HO2 + H2O
! These reactions seem to have the same reaction rate. Therefore
! parameterize this chemistry as doubling the rate of reaction of
! H2O2 + OH.
! h2o2*air*28.966/34. = h2o2 in molec/cm3
h2o2des = 1.7e-12*exp(-200.*(1./tair(i,k) - 0.0033557)) * &
oh(i,k) * h2o2(i,k)*air*mwa2h2o2
h2o2phot = jh2o2(i,k) * h2o2(i,k)*air*mwa2h2o2
h2o2src(i,k) = ho2rate/(air*mwa2h2o2)
h2o2snkg(i,k) = (h2o2des + h2o2phot)/(air*mwa2h2o2)
! h2o2new(i,k) = h2o2(i,k) + (ho2rate - h2o2des - h2o2phot) * &
! ztodt/(air*mwa2h2o2)
h2o2new(i,k) = h2o2(i,k) + (h2o2src(i,k) - h2o2snkg(i,k)) * ztodt
!c h2o2new(i,k) = min(h2o2(i,k) +
!c _ ho2rate * ztodt*(287.*tair(i,k)/pres(i,k))*34./6.022e20,
!c _ xh2o2)
end do
end do
! call t_stopf ('gaschem')
return
end subroutine gaschem
!##############################################################################
subroutine getj( clat, lat, lon, calday, zm, ncol, jout ) 1,1
!-----------------------------------------------------------------------
! Purpose:
! Set h2o2 photolysis rates.
!
! Author: M. Barth
!-----------------------------------------------------------------------
implicit none
real(r8), intent(in) ::&
clat(pcols), &! latitudes (radians)
calday, &! Current calendar day = julian day + fraction
zm(pcols,pver) ! height of midpoint of layer
integer, intent(in) ::&
lat(pcols), &! latitude indices
lon(pcols), &! longitude indices
ncol ! number of atmospheric columns used in chunk
real(r8), intent(out) ::&
jout(pcols,pver) ! photolysis rate for location
! Local variables
logical dodiavg ! when true, diurnally avg zenith angle
real(r8) cosza(pcols) ! cosine of zenith angle
real(r8) dicosza ! diurnally avgd cosine of zenith angle
real(r8) disecza ! diurnally avgd cosine of zenith angle
integer i, is, k
integer ind(3)
integer mz
real(r8) psum
real(r8) ratza !ratio of zenith angles
real(r8) ratalb !ratio of albedos = (0.3-0.2)/(0.5-0.2)
real(r8) sumrat !sum of the ratios
real(r8) w1 !1-sumrat
!-----------------------------------------------------------------------c
dodiavg = .true.
call coszen( calday, dodiavg, clat, lat, lon, ncol, cosza )
! TBH: Loops are not in optimal order for memory access...
do i=1,ncol
dicosza = cosza(i)
if( dicosza .lt. 0.02 ) then
do k=1,pver
jout(i,k) = 0.
end do
else
disecza = min( 1._r8/max( dicosza, 1.e-7_r8 ), 50._r8 )
! Interpolate data onto current level and zenith angle for alb=0.3
! zenith angle:
do is = 1, 8 !zasec = za in table
if( zasec(is) .gt. disecza ) go to 5 !secza = za at gridpt
end do
is = is - 1
5 is = max( is-1, 1 )
ratza = (disecza - zasec(is)) * delang(is)
! albedo -- set to 0.3 (data is for 0.2 and 0.5)
ratalb = (0.3 - alb(1)) * delalb
sumrat = ratza + ratalb
w1 = 1. - sumrat
! Set indices
! note albedo index is either 1 or 2; therefore only one jh2o2 will
! have ialb=2
ind(1) = nz2*(is-1)
ind(2) = nz2*is
ind(3) = nz2*nza + nz2*(is-1)
do k = 1, pver
mz = nint( min( 50._r8, zm(i,k)*0.001_r8 + 1._r8 ) )
psum = w1*jper2(ind(1)+mz) + ratza*jper2(ind(2)+mz) + &
ratalb*jper2(ind(3)+mz)
jout(i,k) = exp( psum )
end do
endif
end do
end subroutine getj
!##############################################################################
end module sulchem