module sulchem,5
!-----------------------------------------------------------------------
! 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 abortutils, only: endrun
use physconst, only: tmelt
use perf_mod
use cam_logfile, only: iulog
implicit none
save
private
public :: &
inisulchem, &! Initialize sulfur cycle chemistry.
sulf_chemdr, &
cldychmini, dicor
! private module data
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_r8/64._r8 ! ratio molecular wt of air to so2
real(r8), parameter :: mwa2h2o2 = 28.966_r8/34._r8 ! ratio molecular wt of air to h2o2
real(r8), parameter, dimension(nsiz) :: jper = &! Tabulated j values
(/-11.63010025_r8,-11.53059959_r8,-11.47570038_r8,-11.44130039_r8,-11.40859985_r8,-11.40690041_r8,-11.42399979_r8,-11.46020031_r8,&
-11.48509979_r8,-11.49100018_r8,-11.46759987_r8,-11.41129971_r8,-11.32050037_r8,-11.20450020_r8,-10.96660042_r8,-10.65079975_r8,&
-10.31389999_r8,-10.15680027_r8,-11.93789959_r8,-11.81239986_r8,-11.73900032_r8,-11.68900013_r8,-11.62989998_r8,-11.60620022_r8,&
-11.60490036_r8,-11.61820030_r8,-11.62650013_r8,-11.62040043_r8,-11.58790016_r8,-11.52509975_r8,-11.43029976_r8,-11.31110001_r8,&
-11.06599998_r8,-10.75179958_r8,-10.39190006_r8,-10.20240021_r8,-12.20409966_r8,-12.06019974_r8,-11.97229958_r8,-11.90939999_r8,&
-11.82750034_r8,-11.78409958_r8,-11.76560020_r8,-11.75730038_r8,-11.74950027_r8,-11.73139954_r8,-11.68959999_r8,-11.62010002_r8,&
-11.52070045_r8,-11.39859962_r8,-11.14780045_r8,-10.83170033_r8,-10.46100044_r8,-10.24470043_r8,-12.51459980_r8,-12.35309982_r8,&
-12.24969959_r8,-12.17240047_r8,-12.06379986_r8,-11.99619961_r8,-11.95629978_r8,-11.92070007_r8,-11.89220047_r8,-11.85890007_r8,&
-11.80539989_r8,-11.72710037_r8,-11.62139988_r8,-11.49520016_r8,-11.23849964_r8,-10.91720009_r8,-10.54150009_r8,-10.29699993_r8,&
-13.14570045_r8,-12.95919991_r8,-12.83010006_r8,-12.72649956_r8,-12.56420040_r8,-12.44439983_r8,-12.35589981_r8,-12.25739956_r8,&
-12.18089962_r8,-12.11240005_r8,-12.03289986_r8,-11.93480015_r8,-11.81389999_r8,-11.67720032_r8,-11.40919971_r8,-11.07330036_r8,&
-10.69830036_r8,-10.41100025_r8,-14.34010029_r8,-14.15550041_r8,-14.01509953_r8,-13.88949966_r8,-13.65480042_r8,-13.43579960_r8,&
-13.23670006_r8,-12.97910023_r8,-12.77519989_r8,-12.61489964_r8,-12.47119999_r8,-12.32730007_r8,-12.17119980_r8,-12.00699997_r8,&
-11.70969963_r8,-11.35060024_r8,-10.96660042_r8,-10.65600014_r8,-15.41930008_r8,-15.26990032_r8,-15.16590023_r8,-15.07689953_r8,&
-14.89920044_r8,-14.68239975_r8,-14.41380024_r8,-13.96300030_r8,-13.55609989_r8,-13.23820019_r8,-12.98649979_r8,-12.77019978_r8,&
-12.56400013_r8,-12.36279964_r8,-12.01889992_r8,-11.63490009_r8,-11.22630024_r8,-10.92059994_r8,-17.23080063_r8,-17.09090042_r8,&
-16.99939919_r8,-16.92959976_r8,-16.82430077_r8,-16.74169922_r8,-16.66230011_r8,-16.49410057_r8,-16.09469986_r8,-15.37609959_r8,&
-14.65410042_r8,-14.07880020_r8,-13.62609959_r8,-13.26119995_r8,-12.75020027_r8,-12.25749969_r8,-11.79109955_r8,-11.45209980_r8,&
-11.17109966_r8,-11.14850044_r8,-11.14080048_r8,-11.14089966_r8,-11.15730000_r8,-11.18900013_r8,-11.23019981_r8,-11.29240036_r8,&
-11.33609962_r8,-11.35649967_r8,-11.34700012_r8,-11.30399990_r8,-11.22620010_r8,-11.12189960_r8,-10.90159988_r8,-10.60239983_r8,&
-10.27840042_r8,-10.12619972_r8,-11.52420044_r8,-11.47599983_r8,-11.44919968_r8,-11.43290043_r8,-11.42099953_r8,-11.42879963_r8,&
-11.44970036_r8,-11.48649979_r8,-11.51080036_r8,-11.51659966_r8,-11.49520016_r8,-11.44289970_r8,-11.35809994_r8,-11.24800014_r8,&
-11.01669979_r8,-10.71510029_r8,-10.36550045_r8,-10.18050003_r8,-11.81890011_r8,-11.75230026_r8,-11.71059990_r8,-11.68099976_r8,&
-11.64519978_r8,-11.63220024_r8,-11.63479996_r8,-11.64830017_r8,-11.65499973_r8,-11.64729977_r8,-11.61489964_r8,-11.55410004_r8,&
-11.46290016_r8,-11.34809971_r8,-11.10849953_r8,-10.80249977_r8,-10.44029999_r8,-10.22789955_r8,-12.15250015_r8,-12.06820011_r8,&
-12.01109982_r8,-11.96700001_r8,-11.90400028_r8,-11.86610031_r8,-11.84650040_r8,-11.83150005_r8,-11.81620026_r8,-11.79199982_r8,&
-11.74639988_r8,-11.67520046_r8,-11.57610035_r8,-11.45580006_r8,-11.20800018_r8,-10.89459991_r8,-10.52560043_r8,-10.28439999_r8,&
-12.80440044_r8,-12.69659996_r8,-12.61540031_r8,-12.54609966_r8,-12.43089962_r8,-12.34109974_r8,-12.27270031_r8,-12.19359970_r8,&
-12.12870026_r8,-12.06789970_r8,-11.99429989_r8,-11.90139961_r8,-11.78509998_r8,-11.65240002_r8,-11.39009953_r8,-11.05939960_r8,&
-10.68850040_r8,-10.40359974_r8,-13.98359966_r8,-13.87979984_r8,-13.79139996_r8,-13.70489979_r8,-13.52700043_r8,-13.34549999_r8,&
-13.17129993_r8,-12.93599987_r8,-12.74400043_r8,-12.59039974_r8,-12.45110035_r8,-12.31060028_r8,-12.15719986_r8,-11.99520016_r8,&
-11.70090008_r8,-11.34430027_r8,-10.96230030_r8,-10.65270042_r8,-15.04539967_r8,-14.97089958_r8,-14.91469955_r8,-14.86209965_r8,&
-14.74139977_r8,-14.57019997_r8,-14.33740044_r8,-13.92109966_r8,-13.53120041_r8,-13.22140026_r8,-12.97410011_r8,-12.76060009_r8,&
-12.55630016_r8,-12.35649967_r8,-12.01439953_r8,-11.63179970_r8,-11.22420025_r8,-10.91909981_r8,-16.85479927_r8,-16.78750038_r8,&
-16.74150085_r8,-16.70509911_r8,-16.64769936_r8,-16.59880066_r8,-16.54520035_r8,-16.40889931_r8,-16.04310036_r8,-15.35249996_r8,&
-14.64319992_r8,-14.07279968_r8,-13.62240028_r8,-13.25860023_r8,-12.74860001_r8,-12.25650024_r8,-11.79049969_r8,-11.45160007_r8 /)
real(r8), parameter, dimension(nz) :: zz = &! Tabulated heights
(/ 0._r8, 1._r8, 2._r8, 3._r8, 5._r8, 7._r8, 9._r8, 12._r8, 15._r8, 18._r8, 21._r8, 24._r8, &
27._r8, 30._r8, 35._r8, 40._r8, 45._r8, 50._r8 /)
real(r8), parameter, dimension(nza) :: zasec = &! Tabulated sec(zen.ang)
(/ 1._r8, 1.3_r8, 1.6_r8, 2._r8, 3._r8, 6._r8, 12._r8, 50._r8 /)
real(r8), parameter, dimension(nalb) :: alb = &! Tabulated albedo = 0.2 and 0.5
(/ 0.2_r8, 0.5_r8 /)
real(r8), parameter, dimension(121) :: yield = &! yield of DMS+OH rxn going to SO2
(/ 0.76905_r8, 0.76953_r8, 0.77002_r8, 0.77050_r8, 0.77098_r8, 0.77146_r8, &
0.77195_r8, 0.77243_r8, 0.77291_r8, 0.77340_r8, 0.77388_r8, 0.77436_r8, &
0.77485_r8, 0.77533_r8, 0.77581_r8, 0.77629_r8, 0.77678_r8, 0.77726_r8, &
0.77774_r8, 0.77822_r8, 0.77870_r8, 0.77918_r8, 0.77967_r8, 0.78015_r8, &
0.78063_r8, 0.78111_r8, 0.78159_r8, 0.78207_r8, 0.78255_r8, 0.78303_r8, &
0.78351_r8, 0.78399_r8, 0.78447_r8, 0.78495_r8, 0.78543_r8, 0.78591_r8, &
0.78640_r8, 0.78688_r8, 0.78737_r8, 0.78786_r8, 0.78835_r8, 0.78884_r8, &
0.78934_r8, 0.78984_r8, 0.79035_r8, 0.79086_r8, 0.79137_r8, 0.79190_r8, &
0.79243_r8, 0.79297_r8, 0.79353_r8, 0.79409_r8, 0.79467_r8, 0.79526_r8, &
0.79587_r8, 0.79650_r8, 0.79715_r8, 0.79783_r8, 0.79853_r8, 0.79927_r8, &
0.80004_r8, 0.80084_r8, 0.80169_r8, 0.80258_r8, 0.80352_r8, 0.80452_r8, &
0.80558_r8, 0.80671_r8, 0.80790_r8, 0.80918_r8, 0.81055_r8, 0.81200_r8, &
0.81356_r8, 0.81522_r8, 0.81700_r8, 0.81890_r8, 0.82093_r8, 0.82309_r8, &
0.82540_r8, 0.82786_r8, 0.83047_r8, 0.83325_r8, 0.83618_r8, 0.83928_r8, &
0.84255_r8, 0.84598_r8, 0.84957_r8, 0.85332_r8, 0.85722_r8, 0.86126_r8, &
0.86543_r8, 0.86973_r8, 0.87412_r8, 0.87861_r8, 0.88316_r8, 0.88777_r8, &
0.89240_r8, 0.89705_r8, 0.90169_r8, 0.90631_r8, 0.91087_r8, 0.91537_r8, &
0.91978_r8, 0.92409_r8, 0.92829_r8, 0.93236_r8, 0.93630_r8, 0.94009_r8, &
0.94373_r8, 0.94721_r8, 0.95054_r8, 0.95370_r8, 0.95670_r8, 0.95954_r8, &
0.96223_r8, 0.96476_r8, 0.96714_r8, 0.96938_r8, 0.97148_r8, 0.97344_r8, &
0.97528_r8 /)
real(r8), parameter, dimension(121) :: dmsrate = &! rate coeff of DMS+OH
(/ 0.21261E-10_r8,0.21112E-10_r8,0.20966E-10_r8,0.20823E-10_r8,0.20683E-10_r8, &
0.20546E-10_r8,0.20411E-10_r8,0.20280E-10_r8,0.20151E-10_r8,0.20024E-10_r8, &
0.19900E-10_r8,0.19779E-10_r8,0.19660E-10_r8,0.19543E-10_r8,0.19428E-10_r8, &
0.19316E-10_r8,0.19206E-10_r8,0.19098E-10_r8,0.18991E-10_r8,0.18887E-10_r8, &
0.18785E-10_r8,0.18684E-10_r8,0.18585E-10_r8,0.18488E-10_r8,0.18393E-10_r8, &
0.18299E-10_r8,0.18207E-10_r8,0.18116E-10_r8,0.18027E-10_r8,0.17939E-10_r8, &
0.17852E-10_r8,0.17766E-10_r8,0.17682E-10_r8,0.17598E-10_r8,0.17516E-10_r8, &
0.17434E-10_r8,0.17354E-10_r8,0.17274E-10_r8,0.17194E-10_r8,0.17115E-10_r8, &
0.17036E-10_r8,0.16958E-10_r8,0.16880E-10_r8,0.16801E-10_r8,0.16722E-10_r8, &
0.16643E-10_r8,0.16563E-10_r8,0.16482E-10_r8,0.16401E-10_r8,0.16317E-10_r8, &
0.16233E-10_r8,0.16146E-10_r8,0.16057E-10_r8,0.15966E-10_r8,0.15871E-10_r8, &
0.15774E-10_r8,0.15673E-10_r8,0.15567E-10_r8,0.15458E-10_r8,0.15343E-10_r8, &
0.15223E-10_r8,0.15097E-10_r8,0.14965E-10_r8,0.14826E-10_r8,0.14680E-10_r8, &
0.14526E-10_r8,0.14365E-10_r8,0.14195E-10_r8,0.14016E-10_r8,0.13828E-10_r8, &
0.13632E-10_r8,0.13426E-10_r8,0.13211E-10_r8,0.12987E-10_r8,0.12754E-10_r8, &
0.12514E-10_r8,0.12265E-10_r8,0.12009E-10_r8,0.11747E-10_r8,0.11479E-10_r8, &
0.11207E-10_r8,0.10932E-10_r8,0.10655E-10_r8,0.10376E-10_r8,0.10098E-10_r8, &
0.98217E-11_r8,0.95482E-11_r8,0.92787E-11_r8,0.90145E-11_r8,0.87565E-11_r8, &
0.85057E-11_r8,0.82630E-11_r8,0.80289E-11_r8,0.78042E-11_r8,0.75893E-11_r8, &
0.73845E-11_r8,0.71899E-11_r8,0.70058E-11_r8,0.68321E-11_r8,0.66687E-11_r8, &
0.65155E-11_r8,0.63722E-11_r8,0.62384E-11_r8,0.61140E-11_r8,0.59985E-11_r8, &
0.58915E-11_r8,0.57926E-11_r8,0.57013E-11_r8,0.56173E-11_r8,0.55402E-11_r8, &
0.54694E-11_r8,0.54047E-11_r8,0.53456E-11_r8,0.52917E-11_r8,0.52426E-11_r8, &
0.51981E-11_r8,0.51578E-11_r8,0.51214E-11_r8,0.50886E-11_r8,0.50591E-11_r8, &
0.50327E-11_r8 /)
!##############################################################################
contains
!##############################################################################
subroutine inisulchem(),19
use cam_history, only: addfld, add_default, phys_decomp
!-----------------------------------------------------------------------
! 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._r8 / (zasec(i+1) - zasec(i))
end do
delalb = 1._r8 / (alb(2) - alb(1))
! Constants from radiation routines needed for dicor()
dayspy = 365._r8
pie = 4._r8*atan(1._r8)
call addfld('DMSSNK ','kg/kg/s' ,pver, 'A', 'DMSsink ',phys_decomp)
call add_default ('DMSSNK', 1, ' ')
call addfld('SO2SRCG ','kg/kg/s' ,pver, 'A', 'SO2 Src G ',phys_decomp)
call add_default ('SO2SRCG', 1, ' ')
call addfld('SO2SRCG2','kg/kg/s' ,pver, 'A', 'SO2 Src G2 ',phys_decomp)
call add_default ('SO2SRCG2', 1, ' ')
call addfld('SO4SRC ','kg/kg/s' ,pver, 'A', 'SO4 Src Tot ',phys_decomp)
call add_default ('SO4SRC', 1, ' ')
call addfld('SO4SRCG ','kg/kg/s' ,pver, 'A', 'SO4 Src G ',phys_decomp)
call add_default ('SO4SRCG', 1, ' ')
call addfld('H2O2SRC ','kg/kg/s' ,pver, 'A', 'H2O2Src ',phys_decomp)
call add_default ('H2O2SRC', 1, ' ')
call addfld('H2O2SNKA','kg/kg/s' ,pver, 'A', 'H2O2SnkA',phys_decomp)
call add_default ('H2O2SNKA', 1, ' ')
call addfld('H2O2SNKG','kg/kg/s' ,pver, 'A', 'H2O2SnkG',phys_decomp)
call add_default ('H2O2SNKG', 1, ' ')
call addfld ('PH','pH',pver, 'A','Cloud Water pH',phys_decomp)
call add_default ('PH', 1, ' ')
end subroutine inisulchem
!##############################################################################
subroutine sulf_chemdr( t, pmid, cldwat, rain, cldf, cldv, &,13
icwmr1, icwmr2, indcp, ncldypts, hion, &
so2, so4, dms, h2o2, ekso2, ekh2o2, &
o3, h2o23d, oh, no3, ho2, q, &
zm, dtime, calday, clat, clon, ncol, lchnk )
!-----------------------------------------------------------------------
! Purpose:
! Sulfur chemistry driver separated from chemwdepdr. The prognostic species
! are DMS, SO2, SO4, and H2O2. Needed a chem driver without the wet
! deposition for MOZART chemistry
!
! Author: F Vitt
!-----------------------------------------------------------------------
use cam_history, only: outfld
real(r8), intent(in) ::&
clat(pcols), &
clon(pcols), &
calday, &
dtime, &! time step
pmid(pcols,pver), &! pressure at layer midpoints
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
cldv(pcols,pver), &! cloudy volume which is occupied by rain or cloud water
rain(pcols,pver), &! total precip mixing ratio
icwmr1(pcols,pver), &! in cloud water mixing ratio for zhang+hack scheme
icwmr2(pcols,pver) ! in cloud water mixing ratio for zhang+hack scheme
real(r8), dimension(pcols,pver), intent(inout) ::&
so2, &! so2 array
so4, &! so4 array
dms, &! dms
h2o2 ! h2o2
real(r8), dimension(pcols,pver), intent(in) ::&
o3, &
no3, &
ho2, &
oh, &
h2o23d
integer, intent(in) :: ncol
integer, intent(in) :: lchnk
real(r8), intent(in) ::&
ekh2o2(pcols,pver) ! H2O2 Henry's Law coefficient
real(r8), intent(inout) ::&
ekso2(pcols,pver) ! SO2 effective Henry's Law coefficient
integer, intent(in) ::&
indcp(pcols,pver), &! Indices of cloudy grid points
ncldypts(pver) ! Number of cloudy grid points
real(r8), dimension(pcols,pver), intent(inout) ::&
hion ! Hydrogen ion concentration in drops
! Local variables:
real(r8), dimension(pcols,pver) ::&
h2o2new, &! h2o2
h2o2tmp, &! temporary 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) ::&
jh2o2 ! photolysis rate of H2O2
integer :: i, k
!----------------------------------------------------------------------
call getj( clat, clon, calday, zm, ncol, jh2o2 )
call aqchem( so2, so4, h2o2, t, pmid, &
cldwat, rain, so2new, so4new, h2o2new, &
o3, cldf, so4srca, hion, h2o2snka, &
ekso2, ekh2o2, dtime, indcp, ncldypts, &
cldv, icwmr1, icwmr2, ncol )
so2tmp(:ncol,:) = so2new(:ncol,:)
so4tmp(:ncol,:) = so4new(:ncol,:)
h2o2tmp(:ncol,:) = h2o2new(:ncol,:)
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, ncol )
#ifdef SCYC_MASSBGT
call endrun ('sulf_chemdr: 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
so2( :ncol,:) = so2new( :ncol,:)
so4( :ncol,:) = so4new( :ncol,:)
dms( :ncol,:) = dmsnew( :ncol,:)
h2o2(:ncol,:) = h2o2new(:ncol,:)
do k = 1, pver
do i = 1, ncol
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 )
endsubroutine sulf_chemdr
subroutine aqchem( so2, so4, h2o2, tair, pres, & 1
qc, qr, so2new, so4new, h2o2new, &
o3, cldf, so4srca, hion, h2o2snka, &
ekso2, ekh2o2, ztodt, ind, ncpts, &
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(r8), intent(in) ::&
ekh2o2(pcols,pver) ! H2O2 Henry's Law coefficient
real(r8), intent(inout) ::&
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:
real(r8) &
theff(pcols), & ! ekso2 at cloudy grid points
sheff ! temp array for theff
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 (r8) ek
ek(e298, dhr, tt) = e298*exp(dhr*(1._r8/tt - 1._r8/298._r8))
!----------------------------------------------------------------------
! call t_startf ('aqchem')
dtchem = ztodt/ncyc
omsm = 1._r8 - 1.e-10_r8
do k=1,pver
do i=1,ncol
! Determine cloud water and rain water mixing ratios
tc = tair(i,k) - tmelt
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._r8-weight)
if(tair(i,k) .gt. tmelt) then
weight = 0._r8 ! fraction of condensate that is ice
else
weight = 1._r8
endif
rwat(i,k) = qr(i,k)/max(cldv(i,k), 0.01_r8) &
*(1._r8-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) - tmelt
! weight = max(0._r8,min(-tc*0.05,1.0_r8))
! write(iulog,'(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(iulog,*) ' 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._r8
h2o2snka(i,k) = 0._r8
end do
end do
! Initialize cloudy grid point arrays
! write(iulog,*) ' 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._r8
th2o2snk(i2) = 0._r8
! write(iulog,'(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._r8
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_r8*exp(-4751._r8/temp(i2))
molo3 = to3(i2)
eko3 = ek(1.15e-2_r8,2560._r8,temp(i2))
ko3hs = 4.28e13_r8*exp(-5533._r8/temp(i2))
ko3so3= 7.436e16_r8*exp(-5280._r8/temp(i2))
! if(lat .eq. 8 .and. k .eq. 14 .and. i2 .gt. 8) then
! write(iulog,'(i3,7e11.4)') i2, tso2(i2), th2o2(i2), thion(i2), &
! theff(i2), temp(i2), patm, molo3
! write(iulog,'(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._r8 + 13._r8*ahion)
! Convert to mol/(mol-s)
molrate = h2o2rate*ttwat(i2)*28.97e-3_r8
! 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_r8)
! 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_r8
! 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_r8)
! check total rate for conservation
sumrate = h2o2rate + o3rate
! if(lat .eq. 8 .and. k .eq. 14 .and. i2 .gt. 8) then
! write(iulog,'(i3,5e11.4)') n, sumrate, h2o2rate, o3rate, &
! qso2/dtchem, qh2o2/dtchem
! endif
molrate = sumrate*ttwat(i2)*28.97e-3_r8
molrate = mwa2so2*min(molrate/mwa2so2, qso2/dtchem)
tmprate = molrate/(max(ttwat(i2), 1.e-30_r8)*28.97e-3_r8)
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_r8 * 64._r8
h2o2snk = h2o2rate*ttwat(i2)*1.e-3_r8 * 34._r8
! if(qh2o2 .lt. omsm*h2o2snk*dtchem) then
! write(iulog,'(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(iulog,'(a,2e12.5,3i4)') 'AQCHEM qso2 < so2snk*dt ', qso2, &
! so2snk*dtchem, lat, ind(i2,k), k
! endif
so4src = 1.5_r8*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_r8)
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_r8
! Limit to keep pH > 0.0 (ion mixing ratio < 1)
ahion = min(max(ahion,1.e-30_r8),0.999_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._r8-cldv(ind(i2,k),k))*tso2(i2)
so4new(ind(i2,k),k) = cldv(ind(i2,k),k)*tso4n(i2) + &
(1._r8-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._r8-cldv(ind(i2,k),k))*th2o2(i2)
hion(ind(i2,k),k) = thion(i2)
ekso2(ind(i2,k),k) = theff(i2)
! write(iulog,'(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_r8*so2(i,k) + 1._r8/3._r8 * so4(i,k)
totse = 0.5_r8*so2new(i,k) + 1._r8/3._r8 * so4new(i,k)
! if(abs(totsb-totse) .gt. .001*totsb) then
! write(iulog,*) 'Total sulfur not conserved in aqueous chemistry'
! write(iulog,*) 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, &
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(r8), intent(out) :: &
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:
real(r8) :: &
theff(pcols) ! ekso2 at cloudy grid points
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(r8) ek
ek(e298, dhr, tt) = e298*exp(dhr*(1._r8/tt - 1._r8/298._r8))
!----------------------------------------------------------------------
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) - tmelt
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._r8-weight)
if(tair(i,k) .gt. tmelt) then
weight = 0._r8 ! fraction of condensate that is ice
else
weight = 1._r8
endif
rwat(i,k) = qr(i,k)/max(cldv(i,k), 0.00001_r8) &
*(1._r8-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._r8
hion(i,k) = 0._r8
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_r8 .and. twat(i,k) .ge. 1.e-12_r8) 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_r8)
ahso3 = ek(1.23_r8,3120._r8,temp(i2)) * &
ek(1.3e-2_r8,2015._r8,temp(i2)) * &
molso2 * tprs(i2)/101325._r8
!++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_r8
!++ag Limit to keep pH > 0.0 (ion mixing ratio < 1)
! thion(i2) = max(thion(i2),1.e-25_r8)
thion(i2) = min(max(thion(i2),1.e-25_r8),0.999_r8)
!--ag
tph(i2) = -log10(thion(i2))
theff(i2) = ek(1.23_r8, 3120._r8, temp(i2)) * (1._r8 + &
ek(1.3e-2_r8, 2015._r8, temp(i2))/thion(i2) * (1._r8 + &
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, clon, 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)
clon(pcols)
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._r8*pie*calday/dayspy
! Solar declination in radians:
delta = .006918_r8 - .399912_r8*cos(theta) + .070257_r8*sin(theta) - &
.006758_r8*cos(2._r8*theta) + .000907_r8*sin(2._r8*theta) - &
.002697_r8*cos(3._r8*theta) + .001480_r8*sin(3._r8*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._r8) then
frac = 1.0_r8
else if (arg .gt. 1._r8) then
frac = 0.0_r8
else
frac = (1._r8/pie)*acos(arg)
endif
tsun = pie*frac
if (tsun .gt. 0._r8) then
coszrsu = sinc*sind + (cosc*cosd*sin(tsun))/tsun
else
coszrsu = 0.0_r8
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
coszrs(i) = sinc*sind - cosc*cosd*cos(2._r8*pie*calday + clon(i))
end if
end do
end subroutine coszen
!##############################################################################
subroutine dicor( calday, clat, ncol, corr )
!-----------------------------------------------------------------------
! Purpose:
! Calculate a correction factor for the ho2 reaction rate to account for
! the diurnal variation of ho2
!
! Author: Mary Barth
!-----------------------------------------------------------------------
implicit none
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._r8*pie*calday/dayspy
!
! Solar declination in radians:
!
delta = .006918_r8 - .399912_r8*cos(theta) + .070257_r8*sin(theta) - &
.006758_r8*cos(2._r8*theta) + .000907_r8*sin(2._r8*theta) - &
.002697_r8*cos(3._r8*theta) + .001480_r8*sin(3._r8*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._r8
yyy = 0._r8
nc = 0
ntimes = 48
do i=1,ntimes
phi = calday + (real(i,r8)-1)/real(ntimes,r8)
coszrs = sinc*sind - cosc*cosd*cos(2._r8*pie*phi)
if (coszrs.gt.0) then
nc = nc + 1
xxx = xxx + coszrs**2
yyy = yyy + coszrs
endif
end do
if (yyy.gt.0._r8) then
dicor_corr = ntimes*xxx/(yyy**2)
else
dicor_corr = 1._r8
endif
!
! Return correction
!
corr(ic) = dicor_corr
end do
end subroutine dicor
!##############################################################################
subroutine gaschem( so2, so4, dms, h2o2, tair, & 1,1
pres, so2new, so4new, dmsnew, h2o2new, &
h2o23d, oh, no3, ho2, h2o, &
jh2o2, so4srcg, dmssnk, so2srcg, so2srcg2, &
h2o2src, h2o2snkg, ztodt, 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_r8, &! rate constant at low pressure
an=3.3_r8, &! constant: (300/Tair)**an
rh300=1.5e-12_r8, &! rate constant at high pressure
am=0._r8 ! 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._r8) 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._r8) then
write(iulog,*) 'TAIR bad ', i, k, tair(i,k)
call endrun ('GASCHEM')
endif
end do
end do
endif
do k=1,pver
do i=1,ncol
so2rate = 0._r8
! determine concentration of air (molec/cm3) = rho(kg/m3)*Na/MWair
air = pres(i,k)/(287._r8*tair(i,k)) * 1.e-3_r8/28.966_r8 * 6.022e23_r8
! air = 0.8*air !N2 concentration
! SO2 + OH rate coefficient:
rlt = rl300 * (300._r8/tair(i,k))**an
rht = rh300 * (300._r8/tair(i,k))**am
ax = rlt*air/rht
psi = 1._r8/(1._r8 + (log10(ax))**2)
bx = rlt*air/(1._r8 + ax)
so2rate = bx * 0.6_r8**psi !rate coef. in cm3/molec-s
! rewrite so2rate as first order rate coefficient
so2rate = so2rate * oh(i,k)
so2snk = so2(i,k)*(1._r8-exp(-so2rate*ztodt))
! Update sulfur species
so2new(i,k) = so2(i,k) - so2snk
so4new(i,k) = so4(i,k) + 1.5_r8*so2snk
so4srcg(i,k) = 1.5_r8*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._r8
dmsrate1 = 0._r8
dmsno3 = 0._r8
! Get dmsrate from tabulated data
ik = int(tair(i,k)-195._r8)
ik = min(max(ik,1), 121)
if(dms(i,k) .ge. 1.e-30_r8) then
dmsrate1 = dmsrate(ik)*oh(i,k)
! DMS + NO3
dmsno3 = 1.e-12_r8*exp(-500._r8*(0.0033557_r8 - 1._r8/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._r8-exp(-dmstot*ztodt)) / ztodt
dmsnew(i,k) = dms(i,k) - dmssnk(i,k)*ztodt
so2srcg(i,k) = 1.032258_r8 * (yield(ik) * dms(i,k)*(1._r8-exp(-dmsrate1*ztodt)) + &
dms(i,k)*(1._r8-exp(-dmsno3*ztodt)) ) / ztodt
so2new(i,k) = so2new(i,k) + so2srcg(i,k)*ztodt
so2srcg2(i,k) = 1.032258_r8 / ztodt * &
yield(ik) * dms(i,k)*(1._r8-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._r8*tair(i,k)) * 1.e-3_r8/28.966_r8 * 6.022e23_r8
rk2 = 1.7e-12_r8*exp( 600._r8*(1._r8/tair(i,k) - 0.0033557_r8))
rk3 = 4.9e-32_r8*exp(1000._r8*(1._r8/tair(i,k) - 0.0033557_r8))
fh2o = 1.4e-21_r8*exp(2200._r8/tair(i,k))
ho2rate = (rk2 + rk3*air)*(1._r8 + 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_r8*exp(-200._r8*(1._r8/tair(i,k) - 0.0033557_r8)) * &
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, clon, calday, zm, ncol, jout ) 1,1
!-----------------------------------------------------------------------
! Purpose:
! Set h2o2 photolysis rates.
!
! Author: M. Barth
!-----------------------------------------------------------------------
implicit none
real(r8), intent(in) ::&
clat(pcols), &! latitudes (radians)
clon(pcols), &! longitude (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, clon, ncol, cosza )
! TBH: Loops are not in optimal order for memory access...
do i=1,ncol
dicosza = cosza(i)
if( dicosza .lt. 0.02_r8 ) then
do k=1,pver
jout(i,k) = 0._r8
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_r8 - alb(1)) * delalb
sumrat = ratza + ratalb
w1 = 1._r8 - 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