module_mosaic_drydep.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
9 module module_mosaic_drydep
10
11
12 ! rce 2005-feb-18 - several fixes for indices of dlo/hi_sect, volumlo/hi_sect,
13 ! which are now (isize,itype)
14
15 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
16 ! variables in module_data_mosaic_asect
17
18
19 contains
20
21
22 !-----------------------------------------------------------------------
23 subroutine mosaic_drydep_driver( &
24 id, ktau, dtstep, config_flags, &
25 gmt, julday, &
26 t_phy, rho_phy, p_phy, &
27 ust, aer_res, &
28 moist, chem, ddvel, &
29 ids,ide, jds,jde, kds,kde, &
30 ims,ime, jms,jme, kms,kme, &
31 its,ite, jts,jte, kts,kte )
32
33 use module_configure, only: grid_config_rec_type, num_moist, num_chem
34 use module_state_description, only: param_first_scalar
35
36 use module_data_mosaic_asect
37 use module_data_mosaic_other
38 use module_mosaic_driver, only: mapaer_tofrom_host
39 use module_peg_util, only: peg_error_fatal
40
41 implicit none
42
43 ! subr arguments
44 integer, intent(in) :: &
45 id, ktau, julday
46
47 integer, intent(in) :: &
48 ids, ide, jds, jde, kds, kde, &
49 ims, ime, jms, jme, kms, kme, &
50 its, ite, jts, jte, kts, kte
51
52 real, intent(in) :: dtstep, gmt
53
54 real, intent(in), &
55 dimension( ims:ime, kms:kme, jms:jme ) :: &
56 t_phy, rho_phy, p_phy
57
58 real, intent(in), &
59 dimension( ims:ime, jms:jme ) :: &
60 ust
61
62 real, intent(in), &
63 dimension( its:ite, jts:jte ) :: &
64 aer_res
65
66 real, intent(in), &
67 dimension( ims:ime, kms:kme, jms:jme, 1:num_moist ) :: &
68 moist
69
70 real, intent(inout), &
71 dimension( ims:ime, kms:kme, jms:jme, 1:num_chem ) :: &
72 chem
73
74 real, intent(inout), &
75 dimension( its:ite, jts:jte, 1:num_chem ) :: &
76 ddvel
77
78 type(grid_config_rec_type), intent(in) :: config_flags
79
80
81 ! local variables
82 integer idum, jdum
83 integer it, jt, kt
84 integer iphase, itype
85 integer ktmaps, ktmape
86 integer ll, l1, n
87 integer levdbg_err, levdbg_info
88
89 integer idiagaa_dum, ijcount_dum
90
91 real dum
92 real vdep_aer(maxd_asize,maxd_atype,maxd_aphase)
93
94 character*100 msg
95
96
97 !!rcetestdd diagnostics --------------------------------------------------
98 ! print 93010, ' '
99 ! print 93010, 'rcetestdd diagnostics from mosaic_drydep_driver'
100 ! print 93010, 'id, chem_opt, ktau, julday ', &
101 ! id, config_flags%chem_opt, ktau, julday
102 ! print 93010, 'ims/e, j, k', ims, ime, jms, jme, kms, kme
103 ! print 93010, 'its/e, j, k', its, ite, jts, jte, kts, kte
104 !
105 ! do jdum = 0, 2
106 ! do idum = 0, 2
107 ! jt = jts + ((jte-jts)/2)*jdum
108 ! it = its + ((ite-its)/2)*idum
109 !if (idum .eq. 2) it = ite
110 !if (jdum .eq. 2) jt = jte
111 ! kt = kts
112 ! print 93050, 'it, jt, t, p, ust, aer_res', it, jt, &
113 ! t_phy(it,kt,jt), p_phy(it,kt,jt), ust(it,jt), aer_res(it,jt)
114 ! end do
115 ! end do
116 !
117 !93010 format( a, 8(1x,i6) )
118 !93020 format( a, 8(1p,e14.6) )
119 !93050 format( a, 2(1x,i4), 8(1p,e14.6) )
120 !!rcetestdd diagnostics --------------------------------------------------
121
122
123 ! ktmaps,ktmape = first/last wrf kt for which depvel is calculated
124 ! ktot = number of levels at which depvel is calculated
125 ktmaps = kts
126 ktmape = kts
127 ktot = 1
128
129
130 ! set some variables to their wrf-chem "standard" values
131 lunerr = -1
132 lunout = -1
133 levdbg_err = 0
134 levdbg_info = 15
135
136 iymdcur = 20030822
137 ihmscur = 0
138 dum = gmt*3600.0 + dtstep*(ktau-1)
139 dum = mod( dum, 86400.0 )
140 ihmscur = nint( dum )
141
142 t = dtstep*(ktau-1)
143 ncorecnt = ktau - 1
144
145 ! reset some variables to "box test" values
146 ! (*** aboxtest_get_extra_args is for "box testing" only
147 ! and should be not be called from wrf-chem ***)
148 ! call aboxtest_get_extra_args( 20, &
149 ! iymdcur, ihmscur, &
150 ! aboxtest_units_convert, aboxtest_rh_method, &
151 ! aboxtest_map_method, aboxtest_gases_fixed, &
152 ! lunerr, lunout, t, dum )
153
154
155 ! set "pegasus" grid size variables
156 itot = ite
157 jtot = jte
158 nsubareas = 1
159
160 ijcount_dum = 0
161
162
163 do 2920 jt = jts, jte
164 do 2910 it = its, ite
165
166 ijcount_dum = ijcount_dum + 1
167
168
169 call mapaer_tofrom_host( 0, &
170 ims,ime, jms,jme, kms,kme, &
171 its,ite, jts,jte, kts,kte, &
172 it, jt, ktmaps,ktmape, &
173 num_moist, num_chem, moist, chem, &
174 t_phy, p_phy, rho_phy )
175
176
177 ! compute deposition velocities
178 idiagaa_dum = 1
179 idiagaa_dum = 0
180 if ((jt.ne.jts) .and. (jt.ne.jte) .and. &
181 (jt.ne.(jts+jte)/2)) idiagaa_dum = 0
182 if ((it.ne.its) .and. (it.ne.ite) .and. &
183 (it.ne.(its+ite)/2)) idiagaa_dum = 0
184
185 call mosaic_drydep_1clm( idiagaa_dum, it, jt, &
186 ust(it,jt), aer_res(it,jt), vdep_aer )
187
188
189 ! rce 23-nov-2004 - change to using the "..._aer" species pointers
190 do iphase = 1, nphase_aer
191 do itype = 1, ntype_aer
192 do n = 1, nsize_aer(itype)
193 do ll = -2, ncomp_plustracer_aer(itype)
194 if (ll .eq. -2) then
195 l1 = numptr_aer(n,itype,iphase)
196 else if (ll .eq. -1) then
197 l1 = -1
198 if (iphase .eq. ai_phase) l1 = waterptr_aer(n,itype)
199 else if (ll .eq. 0) then
200 l1 = -1
201 if (iphase .eq. ai_phase) l1 = hyswptr_aer(n,itype)
202 else
203 l1 = massptr_aer(ll,n,itype,iphase)
204 end if
205 if (l1 .ge. param_first_scalar) then
206 ddvel(it,jt,l1) = vdep_aer(n,itype,iphase)
207 end if
208 end do
209 end do
210 end do
211 end do
212
213
214 2910 continue
215 2920 continue
216 ! print 93010, 'leaving mosaic_drydep_driver'
217
218
219 return
220 end subroutine mosaic_drydep_driver
221
222
223 !-----------------------------------------------------------------------
224 subroutine mosaic_drydep_1clm( idiagaa, it, jt, &
225 ustar_in, depresist_a_in, vdep_aer )
226
227 use module_configure, only: grid_config_rec_type
228
229 use module_data_mosaic_asect
230 use module_data_mosaic_other
231 use module_mosaic_driver, only: mapaer_tofrom_host
232 use module_peg_util, only: peg_error_fatal
233
234 implicit none
235
236 ! subr arguments
237 integer, intent(in) :: idiagaa, it, jt
238
239 ! friction velocity (m/s)
240 real, intent(in) :: ustar_in
241 ! aerodynamic resistance (s/m)
242 real, intent(in) :: depresist_a_in
243
244 ! deposition velocities (m/s)
245 real, intent(inout) :: vdep_aer(maxd_asize,maxd_atype,maxd_aphase)
246
247 ! local variables
248 real, parameter :: densdefault = 2.0
249 real, parameter :: smallmassaa = 1.0e-20
250 real, parameter :: smallmassbb = 1.0e-30
251 real, parameter :: piover6 = pi/6.0
252 real, parameter :: onethird = 1.0/3.0
253
254 integer iphase, itype, k, ll, l1, m, n
255
256 real airdens, airkinvisc
257 real depresist_a, depresist_unstabpblfact
258 real depresist_d0, depresist_d3
259 real depvel_a0, depvel_a3
260 real drydens, drydp, drymass, dryvol
261 real dum, dumalnsg, dumfact, dummass
262 real freepath
263 real rnum
264 real temp
265 real ustar
266 real vsettl_0, vsettl_3
267 real wetdgnum, wetdens, wetdp, wetmass, wetvol
268
269
270 ! if (idiagaa>0) print *, ' '
271 k = 1
272 m = 1
273
274 ! temperature (K)
275 temp = rsub(ktemp,k,m)
276 ! air density (g/cm^3)
277 ! airdens = cairclm(1)*xmwair
278 airdens = cairclm(1)*28.966
279 ! air kinematic viscosity (cm^2/s)
280 airkinvisc = ( 1.8325e-4 * (416.16/(temp+120.0)) * &
281 ((temp/296.16)**1.5) ) / airdens
282 ! air molecular freepath (cm)
283 freepath = 7.39758e-4 * airkinvisc / sqrt(temp)
284 ! friction velocity (cm/s)
285 ustar = ustar_in * 100.0
286 ! aerodynamic resistance (s/cm)
287 depresist_a = depresist_a_in * 0.01
288
289 ! enhancement factor for unstable pbl
290 depresist_unstabpblfact = 1.0
291
292
293 ! initialize vdep_aer
294 vdep_aer(:,:,:) = 0.0
295
296 ! *** for now, just calc vdep_aer for iphase = ai_phase
297 iphase = ai_phase
298
299 ! calculate vdep_aer
300 do 2900 itype = 1, ntype_aer
301 do 2800 n = 1, nsize_aer(itype)
302
303 dryvol = 0.0
304 drymass = 0.0
305 do ll = 1, ncomp_aer(itype)
306 l1 = massptr_aer(ll,n,itype,iphase)
307 dummass = rsub(l1,k,m)*mw_aer(ll,itype)
308 drymass = drymass + dummass
309 dryvol = dryvol + dummass/dens_aer(ll,itype)
310 end do
311
312 l1 = waterptr_aer(n,itype)
313 dummass = rsub(l1,k,m)*mw_water_aer
314 wetmass = drymass + dummass
315 wetvol = dryvol + dummass/dens_water_aer
316
317 l1 = numptr_aer(n,itype,iphase)
318 rnum = rsub(l1,k,m)
319
320 if (drymass .le. smallmassbb) then
321 drydp = dcen_sect(n,itype)
322 drydens = densdefault
323 wetdp = drydp
324 wetdens = drydens
325 goto 1900
326 end if
327
328 !jdf if (drymass .le. smallmassbb) then
329 if (drymass .le. smallmassaa) then
330 wetmass = drymass
331 wetvol = dryvol
332 end if
333 drydens = drymass/dryvol
334 wetdens = wetmass/wetvol
335
336
337 if (rnum .ge. dryvol/volumlo_sect(n,itype)) then
338 drydp = dlo_sect(n,itype)
339 else if (rnum .le. dryvol/volumhi_sect(n,itype)) then
340 drydp = dhi_sect(n,itype)
341 else
342 drydp = (dryvol/(piover6*rnum))**onethird
343 end if
344
345 !jdf dumfact = (wetvol/dryvol)**onethird
346 !jdf dumfact = min( dumfact, 10.0 )
347 if(abs(wetvol).gt.(1000.*abs(dryvol))) then
348 dumfact=10.0
349 else
350 dumfact=abs(wetvol/dryvol)**onethird
351 dumfact=max(1.0,min(dumfact,10.0))
352 endif
353 !jdf
354 wetdp = drydp*dumfact
355
356 1900 continue
357
358
359 !
360 ! get surface resistance and settling velocity for mass (moment 3)
361 ! and number (moment 0)
362 !
363 dumalnsg = log( 1.0 )
364 wetdgnum = wetdp * exp( -1.5*dumalnsg*dumalnsg )
365 call aerosol_depvel_2( &
366 wetdgnum, dumalnsg, wetdens, &
367 temp, airdens, airkinvisc, freepath, &
368 ustar, depresist_unstabpblfact, &
369 depresist_d0, vsettl_0, &
370 depresist_d3, vsettl_3 )
371
372 !
373 ! compute overall deposition velocity (binkowski/shankar eqn a33)
374 !
375 dum = depresist_a + depresist_d3 + &
376 depresist_a*depresist_d3*vsettl_3
377 depvel_a3 = vsettl_3 + (1./dum)
378
379 dum = depresist_a + depresist_d0 + &
380 depresist_a*depresist_d0*vsettl_0
381 depvel_a0 = vsettl_0 + (1./dum)
382
383 ! cm/s --> m/s
384 vdep_aer(n,itype,iphase) = depvel_a3 * 0.01
385
386
387 !
388 ! diagnostic output
389 !
390 if (idiagaa>0) print 9310, it, jt, n, itype, iphase, &
391 dcen_sect(n,itype), drydp, wetdp, &
392 drydens, wetdens, vdep_aer(n,itype,iphase), &
393 vsettl_3, depresist_d3, depresist_a
394 9310 format( 'aerdep', 2i4, 3i3, 1p, 3e10.2, &
395 2x, 0p, 2f5.2, 2x, 1p, 4e10.2 )
396
397
398 2800 continue
399 2900 continue
400
401
402 return
403 end subroutine mosaic_drydep_1clm
404
405
406 !------------------------------------------------------------------------
407 !------------------------------------------------------------------------
408 subroutine aerosol_depvel_2( &
409 dgnum, alnsg, aerodens, &
410 temp, airdens, airkinvisc, freepath, &
411 ustar, depresist_unstabpblfact, &
412 depresist_d0, vsettl_0, &
413 depresist_d3, vsettl_3 )
414 !
415 ! computes the surface layer resistance term and the
416 ! gravitational settling velocity term for the 3rd moment
417 ! of a log-normal aerosol mode
418 !
419 ! input parameters
420 ! dgnum - geometric mean diameter for aerosol number (cm)
421 ! alnsg - natural logarithm of the geometric standard deviation
422 ! for aerosol number
423 ! aerodens - aerosol density (dgnum and aerodens are for the
424 ! actual wet distribution)
425 ! temp - temperature (K)
426 ! airdens - air density (g/cm^3)
427 ! airkinvisc - air kinematic viscosity (cm^2/s)
428 ! freepath - air molecular freepath (cm)
429 ! ustar - friction velocity (cm/s)
430 ! depresist_unstabpblfact = weseley et al. 1985 factor for increasing
431 ! depvel in unstable pbl -- either
432 ! 1. + (-0.3*zi/L)**0.667 OR
433 ! 1. + 0.24*((wstar/ustar)**2)
434 ! output parameters
435 ! depresist_d3 - surface layer resistance for 3rd (mass) moment (s/cm)
436 ! vsettl_3 - gravitational settling velocity for 3rd moment (cm/s)
437 ! depresist_d0 - surface layer resistance for 0th (number) moment (s/cm)
438 ! vsettl_0 - gravitational settling velocity for 0th moment (cm/s)
439 !
440
441 implicit none
442
443 real dgnum, alnsg, aerodens, &
444 temp, airdens, airkinvisc, freepath, &
445 ustar, depresist_unstabpblfact, &
446 depresist_d0, vsettl_0, &
447 depresist_d3, vsettl_3
448
449 real aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0
450 real aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
451 common / aerosol_depvel_cmn01 / &
452 aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0, &
453 aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
454
455 real xknudsen, xknudsenfact, alnsg2, duma, dumb, &
456 vsettl_dgnum, aerodiffus_dgnum
457
458 real pi
459 parameter (pi = 3.1415926536)
460 ! gravity = gravitational acceleration in cm/s^2
461 real gravity
462 parameter (gravity = 980.616)
463 ! boltzmann constant in erg/deg-K
464 real boltzmann
465 parameter (boltzmann = 1.3807e-16)
466
467 xknudsen = 2.*freepath/dgnum
468 xknudsenfact = xknudsen*1.246
469 alnsg2 = alnsg*alnsg
470
471 vsettl_dgnum = (gravity*aerodens*dgnum*dgnum)/ &
472 (18.*airkinvisc*airdens)
473 vsettl_0 = vsettl_dgnum * &
474 ( exp(2.*alnsg2) + xknudsenfact*exp(0.5*alnsg2) )
475 vsettl_3 = vsettl_dgnum * &
476 ( exp(8.*alnsg2) + xknudsenfact*exp(3.5*alnsg2) )
477
478 aerodiffus_dgnum = (boltzmann*temp)/ &
479 (3.*pi*airkinvisc*airdens*dgnum)
480 aerodiffus_0 = aerodiffus_dgnum * &
481 ( exp(+0.5*alnsg2) + xknudsenfact*exp(+2.*alnsg2) )
482 aerodiffus_3 = aerodiffus_dgnum * &
483 ( exp(-2.5*alnsg2) + xknudsenfact*exp(-4.*alnsg2) )
484
485 schmidt_0 = airkinvisc/aerodiffus_0
486 schmidt_3 = airkinvisc/aerodiffus_3
487
488 stokes_0 = ustar*ustar*vsettl_0/(gravity*airkinvisc)
489 stokes_3 = ustar*ustar*vsettl_3/(gravity*airkinvisc)
490
491 duma = (schmidt_0**(-0.66666666)) + &
492 (10.**(-3./max(0.03,stokes_0)))
493 ! dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
494 dumb = duma*ustar*depresist_unstabpblfact
495 depresist_d0 = 1./max( dumb, 1.e-20 )
496 facdepresist_d0 = duma
497
498 duma = (schmidt_3**(-0.66666666)) + &
499 (10.**(-3./max(0.03,stokes_3)))
500 ! dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
501 dumb = duma*ustar*depresist_unstabpblfact
502 depresist_d3 = 1./max( dumb, 1.e-20 )
503 facdepresist_d3 = duma
504
505 return
506 end subroutine aerosol_depvel_2
507
508
509 !-----------------------------------------------------------------------
510 end module module_mosaic_drydep