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