Line data Source code
1 :
2 : module mo_sethet
3 :
4 : !
5 : ! LKE (10/11/2010): added HCN, CH3CN, HCOOH to cesm1_0_beta07_offline version
6 : ! HCN, CH3CN have new Henry's Law coefficients, HCOOH is set to CH3COOH
7 : ! LKE (10/18/2010): SO2 washout corrected based on recommendation of R.Easter, PNNL
8 : !
9 : use shr_kind_mod, only: r8 => shr_kind_r8
10 : use cam_logfile, only: iulog
11 : use gas_wetdep_opts, only: gas_wetdep_cnt, gas_wetdep_method, gas_wetdep_list
12 : use phys_control, only: phys_getopts
13 :
14 : private
15 : public :: sethet_inti, sethet
16 :
17 : save
18 :
19 : integer :: h2o2_ndx, hno3_ndx, ch2o_ndx, ch3ooh_ndx, ch3coooh_ndx, &
20 : ho2no2_ndx, ch3cocho_ndx, xooh_ndx, onitr_ndx, glyald_ndx, &
21 : ch3cho_ndx, mvk_ndx, macr_ndx, pooh_ndx, c2h5ooh_ndx, &
22 : c3h7ooh_ndx, rooh_ndx, isopno3_ndx, onit_ndx, Pb_ndx, &
23 : macrooh_ndx, isopooh_ndx, ch3oh_ndx, c2h5oh_ndx, hyac_ndx, hydrald_ndx
24 : integer :: spc_h2o2_ndx, spc_hno3_ndx
25 : integer :: spc_so2_ndx
26 : integer :: spc_sogm_ndx, spc_sogi_ndx, spc_sogt_ndx, spc_sogb_ndx, spc_sogx_ndx
27 :
28 : integer :: alkooh_ndx, mekooh_ndx, tolooh_ndx, terpooh_ndx, ch3cooh_ndx
29 : integer :: so2_ndx, soa_ndx, so4_ndx, cb2_ndx, oc2_ndx, nh3_ndx, nh4no3_ndx, &
30 : sa1_ndx, sa2_ndx, sa3_ndx, sa4_ndx, nh4_ndx, h2so4_ndx
31 : integer :: xisopno3_ndx,xho2no2_ndx,xonitr_ndx,xhno3_ndx,xonit_ndx
32 : integer :: clono2_ndx, brono2_ndx, hcl_ndx, n2o5_ndx, hocl_ndx, hobr_ndx, hbr_ndx
33 : integer :: ch3cn_ndx, hcn_ndx, hcooh_ndx
34 : integer, allocatable :: wetdep_map(:)
35 : integer :: sogm_ndx, sogi_ndx, sogt_ndx, sogb_ndx, sogx_ndx
36 : logical :: do_wetdep
37 :
38 : ! prognostic modal aerosols
39 : logical :: prog_modal_aero
40 :
41 : contains
42 :
43 1536 : subroutine sethet_inti
44 : !-----------------------------------------------------------------------
45 : ! ... intialize the wet removal rate constants routine
46 : !-----------------------------------------------------------------------
47 :
48 : use mo_chem_utls, only : get_het_ndx, get_spc_ndx
49 : use spmd_utils, only : masterproc
50 : use cam_abortutils, only : endrun
51 :
52 : integer :: k, m
53 :
54 1536 : do_wetdep = gas_wetdep_cnt>0 .and. gas_wetdep_method=='MOZ'
55 1536 : if ( .not. do_wetdep) return
56 :
57 0 : call phys_getopts( prog_modal_aero_out = prog_modal_aero )
58 :
59 0 : allocate( wetdep_map(gas_wetdep_cnt))
60 :
61 0 : do k=1,gas_wetdep_cnt
62 0 : m = get_het_ndx( trim(gas_wetdep_list(k)))
63 0 : if (m>0) then
64 0 : wetdep_map(k) = m
65 : else
66 0 : call endrun('sethet_inti: cannot map '//trim(gas_wetdep_list(k)))
67 : endif
68 : enddo
69 :
70 0 : xisopno3_ndx = get_het_ndx( 'XISOPNO3' )
71 0 : xho2no2_ndx = get_het_ndx( 'XHO2NO2' )
72 0 : xonitr_ndx = get_het_ndx( 'XONITR' )
73 0 : xhno3_ndx = get_het_ndx( 'XHNO3' )
74 0 : xonit_ndx = get_het_ndx( 'XONIT' )
75 :
76 0 : spc_h2o2_ndx = get_spc_ndx( 'H2O2' )
77 0 : spc_hno3_ndx = get_spc_ndx( 'HNO3' )
78 0 : spc_so2_ndx = get_spc_ndx( 'SO2' )
79 :
80 0 : clono2_ndx = get_het_ndx( 'CLONO2' )
81 0 : brono2_ndx = get_het_ndx( 'BRONO2' )
82 0 : hcl_ndx = get_het_ndx( 'HCL' )
83 0 : n2o5_ndx = get_het_ndx( 'N2O5' )
84 0 : hocl_ndx = get_het_ndx( 'HOCL' )
85 0 : hobr_ndx = get_het_ndx( 'HOBR' )
86 0 : hbr_ndx = get_het_ndx( 'HBR' )
87 :
88 0 : h2o2_ndx = get_het_ndx( 'H2O2' )
89 0 : hno3_ndx = get_het_ndx( 'HNO3' )
90 0 : ch2o_ndx = get_het_ndx( 'CH2O' )
91 0 : ch3ooh_ndx = get_het_ndx( 'CH3OOH' )
92 0 : ch3coooh_ndx = get_het_ndx( 'CH3COOOH' )
93 0 : ho2no2_ndx = get_het_ndx( 'HO2NO2' )
94 0 : ch3cocho_ndx = get_het_ndx( 'CH3COCHO' )
95 0 : xooh_ndx = get_het_ndx( 'XOOH' )
96 0 : onitr_ndx = get_het_ndx( 'ONITR' )
97 0 : glyald_ndx = get_het_ndx( 'GLYALD' )
98 0 : ch3cho_ndx = get_het_ndx( 'CH3CHO' )
99 0 : mvk_ndx = get_het_ndx( 'MVK' )
100 0 : macr_ndx = get_het_ndx( 'MACR' )
101 0 : pooh_ndx = get_het_ndx( 'POOH' )
102 0 : c2h5ooh_ndx = get_het_ndx( 'C2H5OOH' )
103 0 : c3h7ooh_ndx = get_het_ndx( 'C3H7OOH' )
104 0 : rooh_ndx = get_het_ndx( 'ROOH' )
105 0 : isopno3_ndx = get_het_ndx( 'ISOPNO3' )
106 0 : onit_ndx = get_het_ndx( 'ONIT' )
107 0 : Pb_ndx = get_het_ndx( 'Pb' )
108 0 : macrooh_ndx = get_het_ndx( 'MACROOH' )
109 0 : isopooh_ndx = get_het_ndx( 'ISOPOOH' )
110 0 : ch3oh_ndx = get_het_ndx( 'CH3OH' )
111 0 : c2h5oh_ndx = get_het_ndx( 'C2H5OH' )
112 0 : hyac_ndx = get_het_ndx( 'HYAC' )
113 0 : hydrald_ndx = get_het_ndx( 'HYDRALD' )
114 0 : alkooh_ndx = get_het_ndx( 'ALKOOH' )
115 0 : mekooh_ndx = get_het_ndx( 'MEKOOH' )
116 0 : tolooh_ndx = get_het_ndx( 'TOLOOH' )
117 0 : terpooh_ndx = get_het_ndx( 'TERPOOH' )
118 0 : ch3cooh_ndx = get_het_ndx( 'CH3COOH' )
119 0 : so2_ndx = get_het_ndx( 'SO2' )
120 0 : soa_ndx = get_het_ndx( 'SOA' )
121 0 : sogb_ndx = get_het_ndx( 'SOGB' )
122 0 : sogi_ndx = get_het_ndx( 'SOGI' )
123 0 : sogm_ndx = get_het_ndx( 'SOGM' )
124 0 : sogt_ndx = get_het_ndx( 'SOGT' )
125 0 : sogx_ndx = get_het_ndx( 'SOGX' )
126 0 : so4_ndx = get_het_ndx( 'SO4' )
127 0 : cb2_ndx = get_het_ndx( 'CB2' )
128 0 : oc2_ndx = get_het_ndx( 'OC2' )
129 0 : nh3_ndx = get_het_ndx( 'NH3' )
130 0 : nh4no3_ndx = get_het_ndx( 'NH4NO3' )
131 0 : nh4_ndx = get_het_ndx( 'NH4' )
132 0 : h2so4_ndx = get_het_ndx( 'H2SO4' )
133 0 : sa1_ndx = get_het_ndx( 'SA1' )
134 0 : sa2_ndx = get_het_ndx( 'SA2' )
135 0 : sa3_ndx = get_het_ndx( 'SA3' )
136 0 : sa4_ndx = get_het_ndx( 'SA4' )
137 0 : ch3cn_ndx = get_het_ndx( 'CH3CN' )
138 0 : hcn_ndx = get_het_ndx( 'HCN' )
139 0 : hcooh_ndx = get_het_ndx( 'HCOOH' )
140 :
141 0 : if (masterproc) then
142 0 : write(iulog,*) 'sethet_inti: new ndx ',so2_ndx,soa_ndx,so4_ndx,cb2_ndx,oc2_ndx, &
143 0 : nh3_ndx,nh4no3_ndx,sa1_ndx,sa2_ndx,sa3_ndx,sa4_ndx
144 0 : write(iulog,*) ' '
145 0 : write(iulog,*) 'sethet_inti: diagnotics '
146 0 : write(iulog,'(10i5)') h2o2_ndx, hno3_ndx, ch2o_ndx, ch3ooh_ndx, ch3coooh_ndx, &
147 0 : ho2no2_ndx, ch3cocho_ndx, xooh_ndx, onitr_ndx, glyald_ndx, &
148 0 : ch3cho_ndx, mvk_ndx, macr_ndx, pooh_ndx, c2h5ooh_ndx, &
149 0 : c3h7ooh_ndx, rooh_ndx, isopno3_ndx, onit_ndx, Pb_ndx, &
150 0 : macrooh_ndx, isopooh_ndx, ch3oh_ndx, c2h5oh_ndx, hyac_ndx, hydrald_ndx
151 : endif
152 :
153 : end subroutine sethet_inti
154 :
155 0 : subroutine sethet( het_rates, press, zmid, phis, tfld, &
156 0 : cmfdqr, nrain, nevapr, delt, xhnm, &
157 0 : qin, ncol, lchnk )
158 : !-----------------------------------------------------------------------
159 : ! ... compute rainout loss rates (1/s)
160 : !-----------------------------------------------------------------------
161 :
162 : use physconst, only : rga,pi
163 : use chem_mods, only : gas_pcnst
164 : use ppgrid, only : pver, pcols
165 : use phys_grid, only : get_rlat_all_p
166 : use cam_abortutils, only : endrun
167 : use mo_constants, only : avo => avogadro, boltz_cgs
168 :
169 : implicit none
170 : !-----------------------------------------------------------------------
171 : ! ... dummy arguments
172 : !-----------------------------------------------------------------------
173 : integer, intent(in) :: ncol ! columns in chunk
174 : integer, intent(in) :: lchnk ! chunk index
175 : real(r8), intent(in) :: delt ! time step ( s )
176 : real(r8), intent(in) :: press(pcols,pver) ! pressure in pascals
177 : real(r8), intent(in) :: cmfdqr(ncol,pver) ! dq/dt for convection
178 : real(r8), intent(in) :: nrain(ncol,pver) ! stratoform precip
179 : real(r8), intent(in) :: nevapr(ncol,pver) ! evaporation
180 : real(r8), intent(in) :: qin(ncol,pver,gas_pcnst) ! xported species ( vmr )
181 : real(r8), intent(in) :: zmid(ncol,pver) ! midpoint geopot (km)
182 : real(r8), intent(in) :: phis(pcols) ! surf geopot
183 : real(r8), intent(in) :: tfld(pcols,pver) ! temperature (k)
184 : real(r8), intent(in) :: xhnm(ncol,pver) ! total atms density ( /cm^3)
185 : real(r8), intent(out) :: het_rates(ncol,pver,gas_pcnst) ! rainout loss rates
186 :
187 : !-----------------------------------------------------------------------
188 : ! ... local variables
189 : !-----------------------------------------------------------------------
190 : real(r8), parameter :: xrm = .189_r8 ! mean diameter of rain drop (cm)
191 : real(r8), parameter :: xum = 748._r8 ! mean rain drop terminal velocity (cm/s)
192 : real(r8), parameter :: xvv = 6.18e-2_r8 ! kinetic viscosity (cm^2/s)
193 : real(r8), parameter :: xdg = .112_r8 ! mass transport coefficient (cm/s)
194 : real(r8), parameter :: t0 = 298._r8 ! reference temperature (k)
195 : real(r8), parameter :: xph0 = 1.e-5_r8 ! cloud [h+]
196 : real(r8), parameter :: satf_hno3 = .016_r8 ! saturation factor for hno3 in clouds
197 : real(r8), parameter :: satf_h2o2 = .016_r8 ! saturation factor for h2o2 in clouds
198 : real(r8), parameter :: satf_so2 = .016_r8 ! saturation factor for so2 in clouds
199 : real(r8), parameter :: satf_ch2o = .1_r8 ! saturation factor for ch2o in clouds
200 : real(r8), parameter :: satf_sog = .016_r8 ! saturation factor for sog in clouds
201 : real(r8), parameter :: const0 = boltz_cgs * 1.e-6_r8 ! (atmospheres/deg k/cm^3)
202 : real(r8), parameter :: hno3_diss = 15.4_r8 ! hno3 dissociation constant
203 : real(r8), parameter :: geo_fac = 6._r8 ! geometry factor (surf area/volume = geo_fac/diameter)
204 : real(r8), parameter :: mass_air = 29._r8 ! mass of background atmosphere (amu)
205 : real(r8), parameter :: mass_h2o = 18._r8 ! mass of water vapor (amu)
206 : real(r8), parameter :: h2o_mol = 1.e3_r8/mass_h2o ! (gm/mol water)
207 : real(r8), parameter :: km2cm = 1.e5_r8 ! convert km to cm
208 : real(r8), parameter :: m2km = 1.e-3_r8 ! convert m to km
209 : real(r8), parameter :: cm3_2_m3 = 1.e-6_r8 ! convert cm^3 to m^3
210 : real(r8), parameter :: m3_2_cm3 = 1.e6_r8 ! convert m^3 to cm^3
211 : real(r8), parameter :: liter_per_gram = 1.e-3_r8
212 : real(r8), parameter :: avo2 = avo * liter_per_gram * cm3_2_m3 ! (liter/gm/mol*(m/cm)^3)
213 :
214 : integer :: i, m, k, kk ! indicies
215 : real(r8) :: xkgm ! mass flux on rain drop
216 : real(r8) :: all1, all2 ! work variables
217 : real(r8) :: stay ! fraction of layer traversed by falling drop in timestep delt
218 : real(r8) :: xeqca1, xeqca2, xca1, xca2, xdtm
219 : real(r8) :: xxx1, xxx2, yhno3, yh2o2
220 : real(r8) :: all3, xeqca3, xca3, xxx3, yso2, so2_diss
221 : real(r8) :: all4, xeqca4, xca4, xxx4
222 : real(r8) :: all5, xeqca5, xca5, xxx5
223 : real(r8) :: all6, xeqca6, xca6, xxx6
224 : real(r8) :: all7, xeqca7, xca7, xxx7
225 : real(r8) :: all8, xeqca8, xca8, xxx8
226 : real(r8) :: ysogm,ysogi,ysogt,ysogb,ysogx
227 :
228 : real(r8), dimension(ncol) :: &
229 0 : xk0, work1, work2, work3, zsurf
230 : real(r8), dimension(pver) :: &
231 : xgas1, xgas2
232 : real(r8), dimension(pver) :: xgas3, xgas4, xgas5, xgas6, xgas7, xgas8
233 : real(r8), dimension(ncol) :: &
234 0 : tmp0_rates, tmp1_rates
235 : real(r8), dimension(ncol,pver) :: &
236 0 : delz, & ! layer depth about interfaces (cm)
237 0 : xhno3, & ! hno3 concentration (molecules/cm^3)
238 0 : xh2o2, & ! h2o2 concentration (molecules/cm^3)
239 0 : xso2, & ! so2 concentration (molecules/cm^3)
240 0 : xsogm, & ! sogm concentration (molecules/cm^3)
241 0 : xsogi, & ! sogi concentration (molecules/cm^3)
242 0 : xsogt, & ! sogt concentration (molecules/cm^3)
243 0 : xsogb, & ! sogb concentration (molecules/cm^3)
244 0 : xsogx, & ! sogx concentration (molecules/cm^3)
245 0 : xliq, & ! liquid rain water content in a grid cell (gm/m^3)
246 0 : rain ! conversion rate of water vapor into rain water (molecules/cm^3/s)
247 : real(r8), dimension(ncol,pver) :: &
248 0 : xhen_hno3, xhen_h2o2, xhen_ch2o, xhen_ch3ooh, xhen_ch3co3h, &
249 0 : xhen_ch3cocho, xhen_xooh, xhen_onitr, xhen_ho2no2, xhen_glyald, &
250 0 : xhen_ch3cho, xhen_mvk, xhen_macr,xhen_sog
251 : real(r8), dimension(ncol,pver) :: &
252 0 : xhen_nh3, xhen_ch3cooh
253 0 : real(r8), dimension(ncol,pver,8) :: tmp_hetrates
254 0 : real(r8), dimension(ncol,pver) :: precip
255 0 : real(r8), dimension(ncol,pver) :: xhen_hcn, xhen_ch3cn, xhen_so2
256 :
257 : integer :: ktop_all
258 0 : integer :: ktop(ncol) ! 100 mb level
259 :
260 : real(r8) :: rlat(pcols) ! latitude in radians for columns
261 : real(r8) :: p_limit
262 : real(r8), parameter :: d2r = pi/180._r8
263 : !
264 : ! jfl : new variables for rescaling sum of positive values to actual amount
265 : !
266 : real(r8) :: total_rain,total_pos
267 : character(len=3) :: hetratestrg
268 : real(r8), parameter :: MISSING = -999999._r8
269 : integer :: mm
270 :
271 : !
272 : !-----------------------------------------------------------------
273 : ! note: the press array is in pascals and must be
274 : ! mutiplied by 10 to yield dynes/cm**2.
275 : !-----------------------------------------------------------------
276 : ! ... set wet deposition for
277 : ! 1. h2o2 2. hno3
278 : ! 3. ch2o 4. ch3ooh
279 : ! 5. pooh 6. ch3coooh
280 : ! 7. ho2no2 8. onit
281 : ! 9. mvk 10. macr
282 : ! 11. c2h5ooh 12. c3h7ooh
283 : ! 13. rooh 14. ch3cocho
284 : ! 15. pb 16. macrooh
285 : ! 17. xooh 18. onitr
286 : ! 19. isopooh 20. ch3oh
287 : ! 21. c2h5oh 22. glyald
288 : ! 23. hyac 24. hydrald
289 : ! 25. ch3cho 26. isopno3
290 : !-----------------------------------------------------------------
291 :
292 0 : het_rates(:,:,:) = 0._r8
293 :
294 0 : if ( .not. do_wetdep) return
295 :
296 0 : call get_rlat_all_p(lchnk, ncol, rlat)
297 :
298 0 : do mm = 1,gas_wetdep_cnt
299 0 : m = wetdep_map(mm)
300 0 : if ( m>0 ) then
301 0 : het_rates(:,:,m) = MISSING
302 : endif
303 : end do
304 :
305 : !-----------------------------------------------------------------
306 : ! ... the 2 and .6 multipliers are from a formula by frossling (1938)
307 : !-----------------------------------------------------------------
308 : xkgm = xdg/xrm * 2._r8 + xdg/xrm * .6_r8 * sqrt( xrm*xum/xvv ) * (xvv/xdg)**(1._r8/3._r8)
309 :
310 : !-----------------------------------------------------------------
311 : ! ... Find 100 mb level
312 : !-----------------------------------------------------------------
313 0 : do i = 1,ncol
314 0 : if ( abs(rlat(i)) > 60._r8*d2r ) then
315 : p_limit = 300.e2_r8
316 : else
317 0 : p_limit = 100.e2_r8
318 : endif
319 0 : k_loop: do k = pver,1,-1
320 0 : if( press(i,k) < p_limit ) then
321 0 : ktop(i) = k
322 0 : exit k_loop
323 : end if
324 : end do k_loop
325 : end do
326 0 : ktop_all = minval( ktop(:) )
327 : !
328 : ! jfl
329 : !
330 : ! this is added to rescale the variable precip (which can only be positive)
331 : ! to the actual vertical integral of positive and negative values. This
332 : ! removes point storms
333 : !
334 0 : do i = 1,ncol
335 : total_rain = 0._r8
336 : total_pos = 0._r8
337 0 : do k = 1,pver
338 0 : precip(i,k) = cmfdqr(i,k) + nrain(i,k) - nevapr(i,k)
339 0 : total_rain = total_rain + precip(i,k)
340 0 : if ( precip(i,k) < 0._r8 ) precip(i,k) = 0._r8
341 0 : total_pos = total_pos + precip(i,k)
342 : end do
343 0 : if ( total_rain <= 0._r8 ) then
344 0 : precip(i,:) = 0._r8
345 : else
346 0 : do k = 1,pver
347 0 : precip(i,k) = precip(i,k) * total_rain/total_pos
348 : end do
349 : end if
350 : end do
351 :
352 0 : do k = 1,pver
353 : !jfl precip(:ncol,k) = cmfdqr(:ncol,k) + nrain(:ncol,k) - nevapr(:ncol,k)
354 0 : rain(:ncol,k) = mass_air*precip(:ncol,k)*xhnm(:ncol,k) / mass_h2o
355 0 : xliq(:ncol,k) = precip(:ncol,k) * delt * xhnm(:ncol,k) / avo*mass_air * m3_2_cm3
356 0 : if( spc_hno3_ndx > 0 ) then
357 0 : xhno3(:ncol,k) = qin(:ncol,k,spc_hno3_ndx) * xhnm(:ncol,k)
358 : else
359 0 : xhno3(:ncol,k) = 0._r8
360 : end if
361 0 : if( spc_h2o2_ndx > 0 ) then
362 0 : xh2o2(:ncol,k) = qin(:ncol,k,spc_h2o2_ndx) * xhnm(:ncol,k)
363 : else
364 0 : xh2o2(:ncol,k) = 0._r8
365 : end if
366 0 : if( spc_sogm_ndx > 0 ) then
367 0 : xsogm(:ncol,k) = qin(:ncol,k,spc_sogm_ndx) * xhnm(:ncol,k)
368 : else
369 0 : xsogm(:ncol,k) = 0._r8
370 : end if
371 0 : if( spc_sogi_ndx > 0 ) then
372 0 : xsogi(:ncol,k) = qin(:ncol,k,spc_sogi_ndx) * xhnm(:ncol,k)
373 : else
374 0 : xsogi(:ncol,k) = 0._r8
375 : end if
376 0 : if( spc_sogt_ndx > 0 ) then
377 0 : xsogt(:ncol,k) = qin(:ncol,k,spc_sogt_ndx) * xhnm(:ncol,k)
378 : else
379 0 : xsogt(:ncol,k) = 0._r8
380 : end if
381 0 : if( spc_sogb_ndx > 0 ) then
382 0 : xsogb(:ncol,k) = qin(:ncol,k,spc_sogb_ndx) * xhnm(:ncol,k)
383 : else
384 0 : xsogb(:ncol,k) = 0._r8
385 : end if
386 0 : if( spc_sogx_ndx > 0 ) then
387 0 : xsogx(:ncol,k) = qin(:ncol,k,spc_sogx_ndx) * xhnm(:ncol,k)
388 : else
389 0 : xsogx(:ncol,k) = 0._r8
390 : end if
391 0 : if( spc_so2_ndx > 0 ) then
392 0 : xso2(:ncol,k) = qin(:ncol,k,spc_so2_ndx) * xhnm(:ncol,k)
393 : else
394 0 : xso2(:ncol,k) = 0._r8
395 : end if
396 : end do
397 :
398 0 : zsurf(:ncol) = m2km * phis(:ncol) * rga
399 0 : do k = ktop_all,pver-1
400 0 : delz(:ncol,k) = abs( (zmid(:ncol,k) - zmid(:ncol,k+1))*km2cm )
401 : end do
402 0 : delz(:ncol,pver) = abs( (zmid(:ncol,pver) - zsurf(:ncol) )*km2cm )
403 :
404 : !-----------------------------------------------------------------
405 : ! ... part 0b, for temperature dependent of henrys
406 : ! xxhe1 = henry con for hno3
407 : ! xxhe2 = henry con for h2o2
408 : !lwh 10/00 -- take henry''s law constants from brasseur et al. [1999],
409 : ! appendix j. for hno3, also consider dissociation to
410 : ! get effective henry''s law constant; equilibrium
411 : ! constant for dissociation from brasseur et al. [1999],
412 : ! appendix k. assume ph=5 (set as xph0 above).
413 : ! heff = h*k/[h+] for hno3 (complete dissociation)
414 : ! heff = h for h2o2 (no dissociation)
415 : ! heff = h * (1 + k/[h+]) (in general)
416 : !-----------------------------------------------------------------
417 0 : do k = ktop_all,pver
418 0 : work1(:ncol) = (t0 - tfld(:ncol,k))/(t0*tfld(:ncol,k))
419 : !-----------------------------------------------------------------
420 : ! ... effective henry''s law constants:
421 : ! hno3, h2o2, ch2o, ch3ooh, ch3coooh (brasseur et al., 1999)
422 : ! xooh, onitr, macrooh (j.-f. muller; brocheton, 1999)
423 : ! isopooh (equal to hno3, as for macrooh)
424 : ! ho2no2 (mozart-1)
425 : ! ch3cocho, hoch2cho (betterton and hoffman, environ. sci. technol., 1988)
426 : ! ch3cho (staudinger and roberts, crit. rev. sci. technol., 1996)
427 : ! mvk, macr (allen et al., environ. toxicol. chem., 1998)
428 : ! soag_bg(0-4), soag_ff_bb(0-4) (Hodzic et al., 2014)
429 : !-----------------------------------------------------------------
430 0 : xk0(:) = 2.1e5_r8 *exp( 8700._r8*work1(:) )
431 0 : xhen_hno3(:,k) = xk0(:) * ( 1._r8 + hno3_diss / xph0 )
432 0 : xhen_h2o2(:,k) = 7.45e4_r8 * exp( 6620._r8 * work1(:) )
433 0 : xhen_ch2o(:,k) = 6.3e3_r8 * exp( 6460._r8 * work1(:) )
434 0 : xhen_ch3ooh(:,k) = 2.27e2_r8 * exp( 5610._r8 * work1(:) )
435 0 : xhen_ch3co3h(:,k) = 4.73e2_r8 * exp( 6170._r8 * work1(:) )
436 0 : xhen_ch3cocho(:,k) = 3.70e3_r8 * exp( 7275._r8 * work1(:) )
437 0 : xhen_xooh(:,k) = 90.5_r8 * exp( 5607._r8 * work1(:) )
438 0 : xhen_onitr(:,k) = 7.51e3_r8 * exp( 6485._r8 * work1(:) )
439 0 : xhen_ho2no2(:,k) = 2.e4_r8
440 0 : xhen_glyald(:,k) = 4.1e4_r8 * exp( 4600._r8 * work1(:) )
441 0 : xhen_ch3cho(:,k) = 1.4e1_r8 * exp( 5600._r8 * work1(:) )
442 0 : xhen_mvk(:,k) = 21._r8 * exp( 7800._r8 * work1(:) )
443 0 : xhen_macr(:,k) = 4.3_r8 * exp( 5300._r8 * work1(:) )
444 0 : xhen_ch3cooh(:,k) = 4.1e3_r8 * exp( 6300._r8 * work1(:) )
445 0 : xhen_sog(:,k) = 5.e5_r8 * exp (12._r8 * work1(:) )
446 : !
447 : ! calculation for NH3 using the parameters in drydep_tables.F90
448 : !
449 0 : xhen_nh3 (:,k) = 1.e6_r8
450 0 : xhen_ch3cn(:,k) = 50._r8 * exp( 4000._r8 * work1(:) )
451 0 : xhen_hcn(:,k) = 12._r8 * exp( 5000._r8 * work1(:) )
452 0 : do i = 1, ncol
453 0 : so2_diss = 1.23e-2_r8 * exp( 1960._r8 * work1(i) )
454 0 : xhen_so2(i,k) = 1.23_r8 * exp( 3120._r8 * work1(i) ) * ( 1._r8 + so2_diss / xph0 )
455 : end do
456 : !
457 0 : tmp_hetrates(:,k,:) = 0._r8
458 : end do
459 :
460 : !-----------------------------------------------------------------
461 : ! ... part 1, solve for high henry constant ( hno3, h2o2)
462 : !-----------------------------------------------------------------
463 0 : col_loop : do i = 1,ncol
464 0 : xgas1(:) = xhno3(i,:) ! xgas will change during
465 0 : xgas2(:) = xh2o2(i,:) ! different levels wash
466 0 : xgas3(:) = xso2 (i,:)
467 0 : xgas4(:) = xsogm(i,:)
468 0 : xgas5(:) = xsogi(i,:)
469 0 : xgas6(:) = xsogt(i,:)
470 0 : xgas7(:) = xsogb(i,:)
471 0 : xgas8(:) = xsogx(i,:)
472 0 : level_loop1 : do kk = ktop(i),pver
473 0 : stay = 1._r8
474 0 : if( rain(i,kk) /= 0._r8 ) then ! finding rain cloud
475 0 : all1 = 0._r8 ! accumulation to justisfy saturation
476 0 : all2 = 0._r8
477 0 : all3 = 0._r8
478 0 : all4 = 0._r8
479 0 : all5 = 0._r8
480 0 : all6 = 0._r8
481 0 : all7 = 0._r8
482 0 : all8 = 0._r8
483 0 : stay = ((zmid(i,kk) - zsurf(i))*km2cm)/(xum*delt)
484 0 : stay = min( stay,1._r8 )
485 : !-----------------------------------------------------------------
486 : ! ... calculate the saturation concentration eqca
487 : !-----------------------------------------------------------------
488 0 : do k = kk,pver ! cal washout below cloud
489 0 : xeqca1 = xgas1(k) &
490 0 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_hno3(i,k)*const0*tfld(i,k))) &
491 0 : * xliq(i,kk)*avo2
492 : xeqca2 = xgas2(k) &
493 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_h2o2(i,k)*const0*tfld(i,k))) &
494 0 : * xliq(i,kk)*avo2
495 : xeqca3 = xgas3(k) &
496 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_so2( i,k)*const0*tfld(i,k))) &
497 0 : * xliq(i,kk)*avo2
498 : xeqca4 = xgas4(k) &
499 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
500 0 : * xliq(i,kk)*avo2
501 : xeqca5 = xgas5(k) &
502 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
503 0 : * xliq(i,kk)*avo2
504 : xeqca6 = xgas6(k) &
505 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
506 0 : * xliq(i,kk)*avo2
507 : xeqca7 = xgas7(k) &
508 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
509 0 : * xliq(i,kk)*avo2
510 : xeqca8 = xgas8(k) &
511 : / (xliq(i,kk)*avo2 + 1._r8/(xhen_sog(i,k)*const0*tfld(i,k))) &
512 0 : * xliq(i,kk)*avo2
513 :
514 : !-----------------------------------------------------------------
515 : ! ... calculate ca; inside cloud concentration in #/cm3(air)
516 : !-----------------------------------------------------------------
517 0 : xca1 = geo_fac*xkgm*xgas1(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
518 0 : xca2 = geo_fac*xkgm*xgas2(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
519 0 : xca3 = geo_fac*xkgm*xgas3(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
520 0 : xca4 = geo_fac*xkgm*xgas4(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
521 0 : xca5 = geo_fac*xkgm*xgas5(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
522 0 : xca6 = geo_fac*xkgm*xgas6(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
523 0 : xca7 = geo_fac*xkgm*xgas7(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
524 0 : xca8 = geo_fac*xkgm*xgas8(k)/(xrm*xum)*delz(i,k) * xliq(i,kk) * cm3_2_m3
525 :
526 : !-----------------------------------------------------------------
527 : ! ... if is not saturated
528 : ! hno3(gas)_new = hno3(gas)_old - hno3(h2o)
529 : ! otherwise
530 : ! hno3(gas)_new = hno3(gas)_old
531 : !-----------------------------------------------------------------
532 0 : all1 = all1 + xca1
533 0 : all2 = all2 + xca2
534 0 : if( all1 < xeqca1 ) then
535 0 : xgas1(k) = max( xgas1(k) - xca1,0._r8 )
536 : end if
537 0 : if( all2 < xeqca2 ) then
538 0 : xgas2(k) = max( xgas2(k) - xca2,0._r8 )
539 : end if
540 0 : all3 = all3 + xca3
541 0 : if( all3 < xeqca3 ) then
542 0 : xgas3(k) = max( xgas3(k) - xca3,0._r8 )
543 : end if
544 0 : all4 = all4 + xca4
545 0 : all5 = all5 + xca5
546 0 : all6 = all6 + xca6
547 0 : all7 = all7 + xca7
548 0 : all8 = all8 + xca8
549 0 : if( all4 < xeqca4 ) then
550 0 : xgas4(k) = max( xgas4(k) - xca4,0._r8 )
551 : end if
552 0 : if( all5 < xeqca5 ) then
553 0 : xgas5(k) = max( xgas5(k) - xca5,0._r8 )
554 : end if
555 0 : if( all6 < xeqca6 ) then
556 0 : xgas6(k) = max( xgas6(k) - xca6,0._r8 )
557 : end if
558 0 : if( all7 < xeqca7 ) then
559 0 : xgas7(k) = max( xgas7(k) - xca7,0._r8 )
560 : end if
561 0 : if( all8 < xeqca8 ) then
562 0 : xgas8(k) = max( xgas8(k) - xca8,0._r8 )
563 : end if
564 : end do
565 : end if
566 : !-----------------------------------------------------------------
567 : ! ... calculate the lifetime of washout (second)
568 : ! after all layers washout
569 : ! the concentration of hno3 is reduced
570 : ! then the lifetime xtt is calculated by
571 : !
572 : ! xtt = (xhno3(ini) - xgas1(new))/(dt*xhno3(ini))
573 : ! where dt = passing time (s) in vertical
574 : ! path below the cloud
575 : ! dt = dz(cm)/um(cm/s)
576 : !-----------------------------------------------------------------
577 0 : xdtm = delz(i,kk) / xum ! the traveling time in each dz
578 0 : xxx1 = (xhno3(i,kk) - xgas1(kk))
579 0 : xxx2 = (xh2o2(i,kk) - xgas2(kk))
580 0 : if( xxx1 /= 0._r8 ) then ! if no washout lifetime = 1.e29
581 0 : yhno3 = xhno3(i,kk)/xxx1 * xdtm
582 : else
583 : yhno3 = 1.e29_r8
584 : end if
585 0 : if( xxx2 /= 0._r8 ) then ! if no washout lifetime = 1.e29
586 0 : yh2o2 = xh2o2(i,kk)/xxx2 * xdtm
587 : else
588 : yh2o2 = 1.e29_r8
589 : end if
590 0 : tmp_hetrates(i,kk,1) = max( 1._r8 / yh2o2,0._r8 ) * stay
591 0 : tmp_hetrates(i,kk,2) = max( 1._r8 / yhno3,0._r8 ) * stay
592 0 : xxx3 = (xso2( i,kk) - xgas3(kk))
593 0 : if( xxx3 /= 0._r8 ) then ! if no washout lifetime = 1.e29
594 0 : yso2 = xso2( i,kk)/xxx3 * xdtm
595 : else
596 : yso2 = 1.e29_r8
597 : end if
598 0 : tmp_hetrates(i,kk,3) = max( 1._r8 / yso2, 0._r8 ) * stay
599 0 : xxx4 = (xsogm(i,kk) - xgas4(kk))
600 0 : xxx5 = (xsogi(i,kk) - xgas5(kk))
601 0 : xxx6 = (xsogt(i,kk) - xgas6(kk))
602 0 : xxx7 = (xsogb(i,kk) - xgas7(kk))
603 0 : xxx8 = (xsogx(i,kk) - xgas8(kk))
604 0 : if( xxx4 /= 0._r8 ) then ! if no washout lifetime = 1.e29
605 0 : ysogm = xsogm(i,kk)/xxx4 * xdtm
606 : else
607 : ysogm = 1.e29_r8
608 : end if
609 0 : if( xxx5 /= 0._r8 ) then ! if no washout lifetime = 1.e29
610 0 : ysogi = xsogi(i,kk)/xxx5 * xdtm
611 : else
612 : ysogi = 1.e29_r8
613 : end if
614 0 : if( xxx6 /= 0._r8 ) then ! if no washout lifetime = 1.e29
615 0 : ysogt = xsogt(i,kk)/xxx6 * xdtm
616 : else
617 : ysogt = 1.e29_r8
618 : end if
619 0 : if( xxx7 /= 0._r8 ) then ! if no washout lifetime = 1.e29
620 0 : ysogb = xsogb(i,kk)/xxx7 * xdtm
621 : else
622 : ysogb = 1.e29_r8
623 : end if
624 0 : if( xxx8 /= 0._r8 ) then ! if no washout lifetime = 1.e29
625 0 : ysogx = xsogx(i,kk)/xxx8 * xdtm
626 : else
627 : ysogx = 1.e29_r8
628 : end if
629 0 : tmp_hetrates(i,kk,4) = max( 1._r8 / ysogm,0._r8 ) * stay
630 0 : tmp_hetrates(i,kk,5) = max( 1._r8 / ysogi,0._r8 ) * stay
631 0 : tmp_hetrates(i,kk,6) = max( 1._r8 / ysogt,0._r8 ) * stay
632 0 : tmp_hetrates(i,kk,7) = max( 1._r8 / ysogb,0._r8 ) * stay
633 0 : tmp_hetrates(i,kk,8) = max( 1._r8 / ysogx,0._r8 ) * stay
634 : end do level_loop1
635 : end do col_loop
636 :
637 : !-----------------------------------------------------------------
638 : ! ... part 2, in-cloud solve for low henry constant
639 : ! hno3 and h2o2 have both in and under cloud
640 : !-----------------------------------------------------------------
641 0 : level_loop2 : do k = ktop_all,pver
642 0 : Column_loop2 : do i=1,ncol
643 0 : if ( rain(i,k) <= 0._r8 ) then
644 0 : het_rates(i,k,:) = 0._r8
645 : cycle
646 : endif
647 :
648 0 : work1(i) = avo2 * xliq(i,k)
649 0 : work2(i) = const0 * tfld(i,k)
650 : work3(i) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch2o(i,k)*work2(i)))),0._r8 ) &
651 0 : * satf_ch2o
652 0 : if( ch2o_ndx > 0 ) then
653 0 : het_rates(i,k,ch2o_ndx) = work3(i)
654 : end if
655 0 : if( isopno3_ndx > 0 ) then
656 0 : het_rates(i,k,isopno3_ndx) = work3(i)
657 : end if
658 0 : if( xisopno3_ndx > 0 ) then
659 0 : het_rates(i,k,xisopno3_ndx) = work3(i)
660 : end if
661 0 : if( hyac_ndx > 0 ) then
662 0 : het_rates(i,k,hyac_ndx) = work3(i)
663 : end if
664 0 : if( hydrald_ndx > 0 ) then
665 0 : het_rates(i,k,hydrald_ndx) = work3(i)
666 : end if
667 :
668 0 : work3(i) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3ooh(i,k)*work2(i)))),0._r8 )
669 0 : if( ch3ooh_ndx > 0 ) then
670 0 : het_rates(i,k,ch3ooh_ndx) = work3(i)
671 : end if
672 0 : if( pooh_ndx > 0 ) then
673 0 : het_rates(i,k,pooh_ndx) = work3(i)
674 : end if
675 0 : if( c2h5ooh_ndx > 0 ) then
676 0 : het_rates(i,k,c2h5ooh_ndx) = work3(i)
677 : end if
678 0 : if( c3h7ooh_ndx > 0 ) then
679 0 : het_rates(i,k,c3h7ooh_ndx) = work3(i)
680 : end if
681 0 : if( rooh_ndx > 0 ) then
682 0 : het_rates(i,k,rooh_ndx) = work3(i)
683 : end if
684 0 : if( ch3oh_ndx > 0 ) then
685 0 : het_rates(i,k,ch3oh_ndx) = work3(i)
686 : end if
687 0 : if( c2h5oh_ndx > 0 ) then
688 0 : het_rates(i,k,c2h5oh_ndx) = work3(i)
689 : end if
690 0 : if( alkooh_ndx > 0 ) then
691 0 : het_rates(i,k,alkooh_ndx) = work3(i)
692 : end if
693 0 : if( mekooh_ndx > 0 ) then
694 0 : het_rates(i,k,mekooh_ndx) = work3(i)
695 : end if
696 0 : if( tolooh_ndx > 0 ) then
697 0 : het_rates(i,k,tolooh_ndx) = work3(i)
698 : end if
699 0 : if( terpooh_ndx > 0 ) then
700 0 : het_rates(i,k,terpooh_ndx) = work3(i)
701 : end if
702 :
703 0 : if( ch3coooh_ndx > 0 ) then
704 0 : het_rates(i,k,ch3coooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3co3h(i,k)*work2(i)))),0._r8 )
705 : end if
706 0 : if( ho2no2_ndx > 0 ) then
707 0 : het_rates(i,k,ho2no2_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ho2no2(i,k)*work2(i)))),0._r8 )
708 : end if
709 0 : if( xho2no2_ndx > 0 ) then
710 0 : het_rates(i,k,xho2no2_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ho2no2(i,k)*work2(i)))),0._r8 )
711 : end if
712 0 : if( ch3cocho_ndx > 0 ) then
713 0 : het_rates(i,k,ch3cocho_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cocho(i,k)*work2(i)))),0._r8 )
714 : end if
715 0 : if( xooh_ndx > 0 ) then
716 0 : het_rates(i,k,xooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_xooh(i,k)*work2(i)))),0._r8 )
717 : end if
718 0 : if( onitr_ndx > 0 ) then
719 0 : het_rates(i,k,onitr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_onitr(i,k)*work2(i)))),0._r8 )
720 : end if
721 0 : if( xonitr_ndx > 0 ) then
722 0 : het_rates(i,k,xonitr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_onitr(i,k)*work2(i)))),0._r8 )
723 : end if
724 0 : if( glyald_ndx > 0 ) then
725 0 : het_rates(i,k,glyald_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_glyald(i,k)*work2(i)))),0._r8 )
726 : end if
727 0 : if( ch3cho_ndx > 0 ) then
728 0 : het_rates(i,k,ch3cho_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cho(i,k)*work2(i)))),0._r8 )
729 : end if
730 0 : if( mvk_ndx > 0 ) then
731 0 : het_rates(i,k,mvk_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_mvk(i,k)*work2(i)))),0._r8 )
732 : end if
733 0 : if( macr_ndx > 0 ) then
734 0 : het_rates(i,k,macr_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_macr(i,k)*work2(i)))),0._r8 )
735 : end if
736 0 : if( h2o2_ndx > 0 ) then
737 0 : work3(i) = satf_h2o2 * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_h2o2(i,k)*work2(i)))),0._r8 )
738 0 : het_rates(i,k,h2o2_ndx) = work3(i) + tmp_hetrates(i,k,1)
739 : end if
740 0 : if ( prog_modal_aero .and. so2_ndx>0 .and. h2o2_ndx>0 ) then
741 0 : het_rates(i,k,so2_ndx) = het_rates(i,k,h2o2_ndx)
742 0 : elseif( so2_ndx > 0 ) then
743 0 : work3(i) = satf_so2 * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_so2( i,k)*work2(i)))),0._r8 )
744 0 : het_rates(i,k,so2_ndx ) = work3(i) + tmp_hetrates(i,k,3)
745 : endif
746 : !
747 0 : work3(i) = satf_sog * max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_sog(i,k)*work2(i)))),0._r8 )
748 0 : if( sogm_ndx > 0 ) then
749 0 : het_rates(i,k,sogm_ndx) = work3(i) + tmp_hetrates(i,k,4)
750 : end if
751 0 : if( sogi_ndx > 0 ) then
752 0 : het_rates(i,k,sogi_ndx) = work3(i) + tmp_hetrates(i,k,5)
753 : end if
754 0 : if( sogt_ndx > 0 ) then
755 0 : het_rates(i,k,sogt_ndx) = work3(i) + tmp_hetrates(i,k,6)
756 : end if
757 0 : if( sogb_ndx > 0 ) then
758 0 : het_rates(i,k,sogb_ndx) = work3(i) + tmp_hetrates(i,k,7)
759 : end if
760 0 : if( sogx_ndx > 0 ) then
761 0 : het_rates(i,k,sogx_ndx) = work3(i) + tmp_hetrates(i,k,8)
762 : end if
763 : !
764 : work3(i) = tmp_hetrates(i,k,2) + satf_hno3 * &
765 0 : max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_hno3(i,k)*work2(i)))),0._r8 )
766 0 : tmp0_rates(i) = work3(i)
767 0 : tmp1_rates(i) = .2_r8*work3(i)
768 0 : if( hno3_ndx > 0 ) then
769 0 : het_rates(i,k,hno3_ndx) = work3(i)
770 : end if
771 0 : if( xhno3_ndx > 0 ) then
772 0 : het_rates(i,k,xhno3_ndx) = work3(i)
773 : end if
774 0 : if( onit_ndx > 0 ) then
775 0 : het_rates(i,k,onit_ndx) = work3(i)
776 : end if
777 0 : if( xonit_ndx > 0 ) then
778 0 : het_rates(i,k,xonit_ndx) = work3(i)
779 : end if
780 0 : if( Pb_ndx > 0 ) then
781 0 : het_rates(i,k,Pb_ndx) = work3(i)
782 : end if
783 0 : if( macrooh_ndx > 0 ) then
784 0 : het_rates(i,k,macrooh_ndx) = work3(i)
785 : end if
786 0 : if( isopooh_ndx > 0 ) then
787 0 : het_rates(i,k,isopooh_ndx) = work3(i)
788 : end if
789 :
790 0 : if( clono2_ndx > 0 ) then
791 0 : het_rates(i,k, clono2_ndx) = work3(i)
792 : end if
793 0 : if( brono2_ndx > 0 ) then
794 0 : het_rates(i,k, brono2_ndx) = work3(i)
795 : end if
796 0 : if( hcl_ndx > 0 ) then
797 0 : het_rates(i,k, hcl_ndx) = work3(i)
798 : end if
799 0 : if( n2o5_ndx > 0 ) then
800 0 : het_rates(i,k, n2o5_ndx) = work3(i)
801 : end if
802 0 : if( hocl_ndx > 0 ) then
803 0 : het_rates(i,k, hocl_ndx) = work3(i)
804 : end if
805 0 : if( hobr_ndx > 0 ) then
806 0 : het_rates(i,k, hobr_ndx) = work3(i)
807 : end if
808 0 : if( hbr_ndx > 0 ) then
809 0 : het_rates(i,k, hbr_ndx) = work3(i)
810 : end if
811 :
812 0 : if( soa_ndx > 0 ) then
813 0 : het_rates(i,k,soa_ndx) = tmp1_rates(i)
814 : end if
815 0 : if( oc2_ndx > 0 ) then
816 0 : het_rates(i,k,oc2_ndx) = tmp1_rates(i)
817 : end if
818 0 : if( cb2_ndx > 0 ) then
819 0 : het_rates(i,k,cb2_ndx) = tmp1_rates(i)
820 : end if
821 0 : if( so4_ndx > 0 ) then
822 0 : het_rates(i,k,so4_ndx) = tmp1_rates(i)
823 : end if
824 0 : if( sa1_ndx > 0 ) then
825 0 : het_rates(i,k,sa1_ndx) = tmp1_rates(i)
826 : end if
827 0 : if( sa2_ndx > 0 ) then
828 0 : het_rates(i,k,sa2_ndx) = tmp1_rates(i)
829 : end if
830 0 : if( sa3_ndx > 0 ) then
831 0 : het_rates(i,k,sa3_ndx) = tmp1_rates(i)
832 : end if
833 0 : if( sa4_ndx > 0 ) then
834 0 : het_rates(i,k,sa4_ndx) = tmp1_rates(i)
835 : end if
836 :
837 0 : if( h2so4_ndx > 0 ) then
838 0 : het_rates(i,k,h2so4_ndx) = tmp0_rates(i)
839 : end if
840 0 : if( nh4_ndx > 0 ) then
841 0 : het_rates(i,k,nh4_ndx) = tmp0_rates(i)
842 : end if
843 0 : if( nh4no3_ndx > 0 ) then
844 0 : het_rates(i,k,nh4no3_ndx ) = tmp0_rates(i)
845 : end if
846 0 : if( nh3_ndx > 0 ) then
847 0 : het_rates(i,k,nh3_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_nh3(i,k)*work2(i)))),0._r8 )
848 : end if
849 :
850 0 : if( ch3cooh_ndx > 0 ) then
851 0 : het_rates(i,k,ch3cooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cooh(i,k)*work2(i)))),0._r8 )
852 : end if
853 0 : if( hcooh_ndx > 0 ) then
854 0 : het_rates(i,k,hcooh_ndx) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cooh(i,k)*work2(i)))),0._r8 )
855 : endif
856 0 : if ( hcn_ndx > 0 ) then
857 0 : het_rates(i,k,hcn_ndx ) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_hcn(i,k)*work2(i)))),0._r8 )
858 : endif
859 0 : if ( ch3cn_ndx > 0 ) then
860 0 : het_rates(i,k,ch3cn_ndx ) = max( rain(i,k) / (h2o_mol*(work1(i) + 1._r8/(xhen_ch3cn(i,k)*work2(i)))),0._r8 )
861 : endif
862 :
863 : end do Column_loop2
864 : end do level_loop2
865 :
866 : !-----------------------------------------------------------------
867 : ! ... Set rates above tropopause = 0.
868 : !-----------------------------------------------------------------
869 0 : do mm = 1,gas_wetdep_cnt
870 0 : m = wetdep_map(mm)
871 0 : do i = 1,ncol
872 0 : do k = 1,ktop(i)
873 0 : het_rates(i,k,m) = 0._r8
874 : end do
875 : end do
876 0 : if ( any( het_rates(:ncol,:,m) == MISSING) ) then
877 0 : write(hetratestrg,'(I3)') m
878 0 : call endrun('sethet: het_rates (wet dep) not set for het reaction number : '//hetratestrg)
879 : endif
880 : end do
881 :
882 0 : end subroutine sethet
883 :
884 : end module mo_sethet
|
