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 ! set ddvel
190 !~ ddvel(it,jt,:) = 0.0
191
192 ! rce 23-nov-2004 - change to using the "..._aer" species pointers
193 do iphase = 1, nphase_aer
194 do itype = 1, ntype_aer
195 do n = 1, nsize_aer(itype)
196 do ll = -2, ncomp_plustracer_aer(itype)
197 if (ll .eq. -2) then
198 l1 = numptr_aer(n,itype,iphase)
199 else if (ll .eq. -1) then
200 l1 = -1
201 if (iphase .eq. ai_phase) l1 = waterptr_aer(n,itype)
202 else if (ll .eq. 0) then
203 l1 = -1
204 if (iphase .eq. ai_phase) l1 = hyswptr_aer(n,itype)
205 else
206 l1 = massptr_aer(ll,n,itype,iphase)
207 end if
208 if (l1 .ge. param_first_scalar) then
209 ddvel(it,jt,l1) = vdep_aer(n,itype,iphase)
210 end if
211 end do
212 end do
213 end do
214 end do
215
216
217 2910 continue
218 2920 continue
219 ! print 93010, 'leaving mosaic_drydep_driver'
220
221
222 return
223 end subroutine mosaic_drydep_driver
224
225
226 !-----------------------------------------------------------------------
227 subroutine mosaic_drydep_1clm( idiagaa, it, jt, &
228 ustar_in, depresist_a_in, vdep_aer )
229
230 use module_configure, only: grid_config_rec_type
231
232 use module_data_mosaic_asect
233 use module_data_mosaic_other
234 use module_mosaic_driver, only: mapaer_tofrom_host
235 use module_peg_util, only: peg_error_fatal
236
237 implicit none
238
239 ! subr arguments
240 integer, intent(in) :: idiagaa, it, jt
241
242 ! friction velocity (m/s)
243 real, intent(in) :: ustar_in
244 ! aerodynamic resistance (s/m)
245 real, intent(in) :: depresist_a_in
246
247 ! deposition velocities (m/s)
248 real, intent(inout) :: vdep_aer(maxd_asize,maxd_atype,maxd_aphase)
249
250 ! local variables
251 real, parameter :: densdefault = 2.0
252 real, parameter :: smallmassaa = 1.0e-20
253 real, parameter :: smallmassbb = 1.0e-30
254 real, parameter :: piover6 = pi/6.0
255 real, parameter :: onethird = 1.0/3.0
256
257 integer iphase, itype, k, ll, l1, m, n
258
259 real airdens, airkinvisc
260 real depresist_a, depresist_unstabpblfact
261 real depresist_d0, depresist_d3
262 real depvel_a0, depvel_a3
263 real drydens, drydp, drymass, dryvol
264 real dum, dumalnsg, dumfact, dummass
265 real freepath
266 real rnum
267 real temp
268 real ustar
269 real vsettl_0, vsettl_3
270 real wetdgnum, wetdens, wetdp, wetmass, wetvol
271
272
273 ! if (idiagaa>0) print *, ' '
274 k = 1
275 m = 1
276
277 ! temperature (K)
278 temp = rsub(ktemp,k,m)
279 ! air density (g/cm^3)
280 ! airdens = cairclm(1)*xmwair
281 airdens = cairclm(1)*28.966
282 ! air kinematic viscosity (cm^2/s)
283 airkinvisc = ( 1.8325e-4 * (416.16/(temp+120.0)) * &
284 ((temp/296.16)**1.5) ) / airdens
285 ! air molecular freepath (cm)
286 freepath = 7.39758e-4 * airkinvisc / sqrt(temp)
287 ! friction velocity (cm/s)
288 ustar = ustar_in * 100.0
289 ! aerodynamic resistance (s/cm)
290 depresist_a = depresist_a_in * 0.01
291
292 ! enhancement factor for unstable pbl
293 depresist_unstabpblfact = 1.0
294
295
296 ! initialize vdep_aer
297 vdep_aer(:,:,:) = 0.0
298
299 ! *** for now, just calc vdep_aer for iphase = ai_phase
300 iphase = ai_phase
301
302 ! calculate vdep_aer
303 do 2900 itype = 1, ntype_aer
304 do 2800 n = 1, nsize_aer(itype)
305
306 dryvol = 0.0
307 drymass = 0.0
308 do ll = 1, ncomp_aer(itype)
309 l1 = massptr_aer(ll,n,itype,iphase)
310 dummass = rsub(l1,k,m)*mw_aer(ll,itype)
311 drymass = drymass + dummass
312 dryvol = dryvol + dummass/dens_aer(ll,itype)
313 end do
314
315 l1 = waterptr_aer(n,itype)
316 dummass = rsub(l1,k,m)*mw_water_aer
317 wetmass = drymass + dummass
318 wetvol = dryvol + dummass/dens_water_aer
319
320 l1 = numptr_aer(n,itype,iphase)
321 rnum = rsub(l1,k,m)
322
323 if (drymass .le. smallmassbb) then
324 drydp = dcen_sect(n,itype)
325 drydens = densdefault
326 wetdp = drydp
327 wetdens = drydens
328 goto 1900
329 end if
330
331 if (drymass .le. smallmassbb) then
332 wetmass = drymass
333 wetvol = dryvol
334 end if
335 drydens = drymass/dryvol
336 wetdens = wetmass/wetvol
337
338
339 if (rnum .ge. dryvol/volumlo_sect(n,itype)) then
340 drydp = dlo_sect(n,itype)
341 else if (rnum .le. dryvol/volumhi_sect(n,itype)) then
342 drydp = dhi_sect(n,itype)
343 else
344 drydp = (dryvol/(piover6*rnum))**onethird
345 end if
346
347 dumfact = (wetvol/dryvol)**onethird
348 dumfact = min( dumfact, 10.0 )
349 wetdp = drydp*dumfact
350
351 1900 continue
352
353
354 !
355 ! get surface resistance and settling velocity for mass (moment 3)
356 ! and number (moment 0)
357 !
358 dumalnsg = log( 1.0 )
359 wetdgnum = wetdp * exp( -1.5*dumalnsg*dumalnsg )
360 call aerosol_depvel_2( &
361 wetdgnum, dumalnsg, wetdens, &
362 temp, airdens, airkinvisc, freepath, &
363 ustar, depresist_unstabpblfact, &
364 depresist_d0, vsettl_0, &
365 depresist_d3, vsettl_3 )
366
367 !
368 ! compute overall deposition velocity (binkowski/shankar eqn a33)
369 !
370 dum = depresist_a + depresist_d3 + &
371 depresist_a*depresist_d3*vsettl_3
372 depvel_a3 = vsettl_3 + (1./dum)
373
374 dum = depresist_a + depresist_d0 + &
375 depresist_a*depresist_d0*vsettl_0
376 depvel_a0 = vsettl_0 + (1./dum)
377
378 ! cm/s --> m/s
379 vdep_aer(n,itype,iphase) = depvel_a3 * 0.01
380
381
382 !
383 ! diagnostic output
384 !
385 if (idiagaa>0) print 9310, it, jt, n, itype, iphase, &
386 dcen_sect(n,itype), drydp, wetdp, &
387 drydens, wetdens, vdep_aer(n,itype,iphase), &
388 vsettl_3, depresist_d3, depresist_a
389 9310 format( 'aerdep', 2i4, 3i3, 1p, 3e10.2, &
390 2x, 0p, 2f5.2, 2x, 1p, 4e10.2 )
391
392
393 2800 continue
394 2900 continue
395
396
397 return
398 end subroutine mosaic_drydep_1clm
399
400
401 !------------------------------------------------------------------------
402 !------------------------------------------------------------------------
403 subroutine aerosol_depvel_2( &
404 dgnum, alnsg, aerodens, &
405 temp, airdens, airkinvisc, freepath, &
406 ustar, depresist_unstabpblfact, &
407 depresist_d0, vsettl_0, &
408 depresist_d3, vsettl_3 )
409 !
410 ! computes the surface layer resistance term and the
411 ! gravitational settling velocity term for the 3rd moment
412 ! of a log-normal aerosol mode
413 !
414 ! input parameters
415 ! dgnum - geometric mean diameter for aerosol number (cm)
416 ! alnsg - natural logarithm of the geometric standard deviation
417 ! for aerosol number
418 ! aerodens - aerosol density (dgnum and aerodens are for the
419 ! actual wet distribution)
420 ! temp - temperature (K)
421 ! airdens - air density (g/cm^3)
422 ! airkinvisc - air kinematic viscosity (cm^2/s)
423 ! freepath - air molecular freepath (cm)
424 ! ustar - friction velocity (cm/s)
425 ! depresist_unstabpblfact = weseley et al. 1985 factor for increasing
426 ! depvel in unstable pbl -- either
427 ! 1. + (-0.3*zi/L)**0.667 OR
428 ! 1. + 0.24*((wstar/ustar)**2)
429 ! output parameters
430 ! depresist_d3 - surface layer resistance for 3rd (mass) moment (s/cm)
431 ! vsettl_3 - gravitational settling velocity for 3rd moment (cm/s)
432 ! depresist_d0 - surface layer resistance for 0th (number) moment (s/cm)
433 ! vsettl_0 - gravitational settling velocity for 0th moment (cm/s)
434 !
435
436 implicit none
437
438 real dgnum, alnsg, aerodens, &
439 temp, airdens, airkinvisc, freepath, &
440 ustar, depresist_unstabpblfact, &
441 depresist_d0, vsettl_0, &
442 depresist_d3, vsettl_3
443
444 real aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0
445 real aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
446 common / aerosol_depvel_cmn01 / &
447 aerodiffus_0, schmidt_0, stokes_0, facdepresist_d0, &
448 aerodiffus_3, schmidt_3, stokes_3, facdepresist_d3
449
450 real xknudsen, xknudsenfact, alnsg2, duma, dumb, &
451 vsettl_dgnum, aerodiffus_dgnum
452
453 real pi
454 parameter (pi = 3.1415926536)
455 ! gravity = gravitational acceleration in cm/s^2
456 real gravity
457 parameter (gravity = 980.616)
458 ! boltzmann constant in erg/deg-K
459 real boltzmann
460 parameter (boltzmann = 1.3807e-16)
461
462 xknudsen = 2.*freepath/dgnum
463 xknudsenfact = xknudsen*1.246
464 alnsg2 = alnsg*alnsg
465
466 vsettl_dgnum = (gravity*aerodens*dgnum*dgnum)/ &
467 (18.*airkinvisc*airdens)
468 vsettl_0 = vsettl_dgnum * &
469 ( exp(2.*alnsg2) + xknudsenfact*exp(0.5*alnsg2) )
470 vsettl_3 = vsettl_dgnum * &
471 ( exp(8.*alnsg2) + xknudsenfact*exp(3.5*alnsg2) )
472
473 aerodiffus_dgnum = (boltzmann*temp)/ &
474 (3.*pi*airkinvisc*airdens*dgnum)
475 aerodiffus_0 = aerodiffus_dgnum * &
476 ( exp(+0.5*alnsg2) + xknudsenfact*exp(+2.*alnsg2) )
477 aerodiffus_3 = aerodiffus_dgnum * &
478 ( exp(-2.5*alnsg2) + xknudsenfact*exp(-4.*alnsg2) )
479
480 schmidt_0 = airkinvisc/aerodiffus_0
481 schmidt_3 = airkinvisc/aerodiffus_3
482
483 stokes_0 = ustar*ustar*vsettl_0/(gravity*airkinvisc)
484 stokes_3 = ustar*ustar*vsettl_3/(gravity*airkinvisc)
485
486 duma = (schmidt_0**(-0.66666666)) + &
487 (10.**(-3./max(0.03,stokes_0)))
488 ! dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
489 dumb = duma*ustar*depresist_unstabpblfact
490 depresist_d0 = 1./max( dumb, 1.e-20 )
491 facdepresist_d0 = duma
492
493 duma = (schmidt_3**(-0.66666666)) + &
494 (10.**(-3./max(0.03,stokes_3)))
495 ! dumb = duma*ustar*(1. + 0.24*wstaroverustar*wstaroverustar)
496 dumb = duma*ustar*depresist_unstabpblfact
497 depresist_d3 = 1./max( dumb, 1.e-20 )
498 facdepresist_d3 = duma
499
500 return
501 end subroutine aerosol_depvel_2
502
503
504 !-----------------------------------------------------------------------
505 end module module_mosaic_drydep