module_mosaic_coag.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_coag
10
11
12
13 use module_peg_util
14
15
16
17 implicit none
18
19
20
21 contains
22
23
24
25 !-----------------------------------------------------------------------
26 subroutine mosaic_coag_1clm( istat_coag, &
27 it, jt, kclm_calcbgn, kclm_calcend, &
28 idiagbb_in, &
29 dtchem, dtcoag_in, &
30 id, ktau, ktauc, its, ite, jts, jte, kts, kte )
31 !
32 ! calculates aerosol particle coagulation for grid points in
33 ! the i=it, j=jt vertical column
34 ! over timestep dtcoag_in
35 !
36 ! uses a version of subr. coagsolv (see additional disclaimer below)
37 ! that was obtained from mark jacobson in june 2003,
38 ! and then modified to work with mosaic aerosol routines
39 !
40 use module_data_mosaic_asect
41 use module_data_mosaic_other
42 use module_state_description, only: param_first_scalar
43
44 ! subr arguments
45 integer, intent(inout) :: istat_coag ! =0 if no problems
46 integer, intent(in) :: &
47 it, jt, kclm_calcbgn, kclm_calcend, &
48 idiagbb_in, &
49 id, ktau, ktauc, its, ite, jts, jte, kts, kte
50 real, intent(in) :: dtchem, dtcoag_in
51
52 ! NOTE - much information is passed via arrays in
53 ! module_data_mosaic_asect and module_data_mosaic_other
54 !
55 ! rsub (inout) - aerosol mixing ratios
56 ! aqmassdry_sub, aqvoldry_sub (inout) -
57 ! aerosol dry-mass and dry-volume mixing ratios
58 ! adrydens_sub (inout) - aerosol dry density
59 ! rsub(ktemp,:,:), ptotclm, cairclm (in) -
60 ! air temperature, pressure, and molar density
61
62
63 ! local variables
64 integer, parameter :: coag_method = 1
65 integer, parameter :: ncomp_coag_maxd = maxd_acomp + 3
66
67 integer :: k, l, ll, lla, llb, llx, m
68 integer :: isize, itype, iphase
69 integer :: iconform_numb
70 integer :: idiagbb, idiag_coag_onoff, iok
71 integer :: lunout_coag
72 integer :: ncomp_coag, nsize, nsubstep_coag, ntau_coag
73 integer,save :: ncountaa(10), ncountbb(0:ncomp_coag_maxd)
74 integer :: p1st
75
76 real, parameter :: densdefault = 2.0
77 real, parameter :: factsafe = 1.00001
78 real, parameter :: onethird = 1.0/3.0
79 real, parameter :: piover6 = pi/6.0
80 real, parameter :: smallmassbb = 1.0e-30
81
82 real :: cair_box
83 real :: dtcoag
84 real :: duma, dumb, dumc
85 real :: patm_box
86 real :: temp_box
87 real :: xxdens, xxmass, xxnumb, xxvolu
88 real :: xxmasswet, xxvoluwet
89
90 real :: dpdry(maxd_asize), dpwet(maxd_asize), denswet(maxd_asize)
91 real :: num_distrib(maxd_asize)
92 real :: sumnew(ncomp_coag_maxd), sumold(ncomp_coag_maxd)
93 real :: vol_distrib(maxd_asize,ncomp_coag_maxd)
94 real :: xxvolubb(maxd_asize)
95
96 character(len=120) :: msg
97
98
99
100 istat_coag = 0
101
102 lunout_coag = 6
103 if (lunout .gt. 0) lunout_coag = lunout
104
105
106 ! when dtcoag_in > dtchem, do not perform coagulation calcs
107 ! on every chemistry time step
108 ntau_coag = nint( dtcoag_in/dtchem )
109 ntau_coag = max( 1, ntau_coag )
110 if (mod(ktau,ntau_coag) .ne. 0) return
111 dtcoag = dtchem*ntau_coag
112
113
114 ! set variables that do not change
115 p1st = PARAM_FIRST_SCALAR
116 idiagbb = idiagbb_in
117
118 ! NOTE - code currently just does 1 type
119 itype = 1
120 iphase = ai_phase
121 nsize = nsize_aer(itype)
122 ncomp_coag = ncomp_plustracer_aer(itype) + 3
123
124
125 ! loop over subareas (currently only 1) and vertical levels
126 do 2900 m = 1, nsubareas
127
128 do 2800 k = kclm_calcbgn, kclm_calcend
129
130
131 ! temporary diagnostics
132 if ((it .eq. its) .and. &
133 (jt .eq. jts) .and. (k .eq. kclm_calcbgn)) then
134 ncountaa(:) = 0
135 ncountbb(:) = 0
136 end if
137
138
139 ncountaa(1) = ncountaa(1) + 1
140 if (afracsubarea(k,m) .lt. 1.e-4) goto 2700
141
142 cair_box = cairclm(k)
143 temp_box = rsub(ktemp,k,m)
144 patm_box = ptotclm(k)/1.01325e6
145
146 nsubstep_coag = 1
147 idiag_coag_onoff = -1
148
149 ! do initial calculation of dry mass, volume, and density,
150 ! and initial number conforming (as needed)
151 vol_distrib(:,:) = 0.0
152 sumold(:) = 0.0
153 do isize = 1, nsize
154 l = numptr_aer(isize,itype,iphase)
155 xxnumb = rsub(l,k,m)
156 xxmass = aqmassdry_sub(isize,itype,k,m)
157 xxvolu = aqvoldry_sub( isize,itype,k,m)
158 xxdens = adrydens_sub( isize,itype,k,m)
159 iconform_numb = 1
160
161 duma = max( abs(xxmass), abs(xxvolu*xxdens), 0.1*smallmassbb )
162 dumb = abs(xxmass - xxvolu*xxdens)/duma
163 if ( (xxdens .lt. 0.1) .or. (xxdens .gt. 20.0) &
164 .or. (dumb .gt. 1.0e-4) ) then
165 ! (exception) case of drydensity not valid, or mass /= volu*dens
166 ! so compute dry mass and volume from rsub
167 ncountaa(2) = ncountaa(2) + 1
168 if (idiagbb .ge. 200) &
169 write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-2a', &
170 ktau, it, jt, k, isize, xxmass, xxvolu, xxdens
171 xxmass = 0.0
172 xxvolu = 0.0
173 do ll = 1, ncomp_aer(itype)
174 l = massptr_aer(ll,isize,itype,iphase)
175 if (l .ge. p1st) then
176 duma = max( 0.0, rsub(l,k,m) )*mw_aer(ll,itype)
177 xxmass = xxmass + duma
178 xxvolu = xxvolu + duma/dens_aer(ll,itype)
179 end if
180 end do
181 if (xxmass .gt. 0.99*smallmassbb) then
182 xxdens = xxmass/xxvolu
183 xxvolu = xxmass/xxdens
184 if (idiagbb .ge. 200) &
185 write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-2b', &
186 ktau, it, jt, k, isize, xxmass, xxvolu, xxdens
187 end if
188 end if
189
190 if (xxmass .le. 1.01*smallmassbb) then
191 ! when drymass extremely small, use default density and bin center size,
192 ! and zero out water
193 ncountaa(3) = ncountaa(3) + 1
194 xxdens = densdefault
195 xxvolu = xxmass/xxdens
196 xxnumb = xxmass/(volumcen_sect(isize,itype)*xxdens)
197 l = hyswptr_aer(isize,itype)
198 if (l .ge. p1st) rsub(l,k,m) = 0.0
199 l = waterptr_aer(isize,itype)
200 if (l .ge. p1st) rsub(l,k,m) = 0.0
201 iconform_numb = 0
202 if (idiagbb .ge. 200) &
203 write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-2c', &
204 ktau, it, jt, k, isize, xxmass, xxvolu, xxdens
205 else
206 xxvolu = xxmass/xxdens
207 end if
208
209 ! check for mean dry-size 'safely' within bounds, and conform number if not
210 if (iconform_numb .gt. 0) then
211 if (xxnumb .gt. &
212 xxvolu/(factsafe*volumlo_sect(isize,itype))) then
213 ncountaa(4) = ncountaa(4) + 1
214 duma = xxnumb
215 xxnumb = xxvolu/(factsafe*volumlo_sect(isize,itype))
216 if (idiagbb .ge. 200) &
217 write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-4 ', &
218 ktau, it, jt, k, isize, xxvolu, duma, xxnumb
219 else if (xxnumb .lt. &
220 xxvolu*factsafe/volumhi_sect(isize,itype)) then
221 ncountaa(5) = ncountaa(5) + 1
222 duma = xxnumb
223 xxnumb = xxvolu*factsafe/volumhi_sect(isize,itype)
224 if (idiagbb .ge. 200) &
225 write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-5 ', &
226 ktau, it, jt, k, isize, xxvolu, duma, xxnumb
227 end if
228 end if
229
230 ! load numb, mass, volu, dens back into mosaic_asect arrays
231 l = numptr_aer(isize,itype,iphase)
232 rsub(l,k,m) = xxnumb
233 adrydens_sub( isize,itype,k,m) = xxdens
234 aqmassdry_sub(isize,itype,k,m) = xxmass
235 aqvoldry_sub( isize,itype,k,m) = xxvolu
236
237 !
238 ! load coagsolv arrays with number, mass, and volume mixing ratios
239 !
240 ! *** NOTE ***
241 ! num_distrib must be (#/cm3)
242 ! vol_distrib units do not matter - they can be either masses or volumes,
243 ! and either mixing ratios or concentrations
244 ! (e.g., g/cm3-air, ug/kg-air, cm3/cm3-air, cm3/kg-air, ...)
245 !
246 num_distrib(isize) = xxnumb*cair_box
247
248 do ll = 1, ncomp_plustracer_aer(itype)
249 l = massptr_aer(ll,isize,itype,iphase)
250 duma = 0.0
251 if (l .ge. p1st) duma = max( 0.0, rsub(l,k,m) )
252 vol_distrib(isize,ll) = duma
253 sumold(ll) = sumold(ll) + duma
254 end do
255
256 do llx = 1, 3
257 ll = (ncomp_coag-3) + llx
258 duma = 0.0
259 if (llx .eq. 1) then
260 l = hyswptr_aer(isize,itype)
261 if (l .ge. p1st) duma = max( 0.0, rsub(l,k,m) )
262 else if (llx .eq. 2) then
263 l = waterptr_aer(isize,itype)
264 if (l .ge. p1st) duma = max( 0.0, rsub(l,k,m) )
265 else
266 duma = max( 0.0, aqvoldry_sub( isize,itype,k,m) )
267 end if
268 vol_distrib(isize,ll) = duma
269 sumold(ll) = sumold(ll) + duma
270 end do
271
272 ! calculate dry and wet diameters and wet density
273 if (xxmass .le. 1.01*smallmassbb) then
274 dpdry(isize) = dcen_sect(isize,itype)
275 dpwet(isize) = dpdry(isize)
276 denswet(isize) = xxdens
277 else
278 dpdry(isize) = (xxvolu/(xxnumb*piover6))**onethird
279 dpdry(isize) = max( dpdry(isize), dlo_sect(isize,itype) )
280 l = waterptr_aer(isize,itype)
281 if (l .ge. p1st) then
282 duma = max( 0.0, rsub(l,k,m) )*mw_water_aer
283 xxmasswet = xxmass + duma
284 xxvoluwet = xxvolu + duma/dens_water_aer
285 dpwet(isize) = (xxvoluwet/(xxnumb*piover6))**onethird
286 dpwet(isize) = min( dpwet(isize), 30.0*dhi_sect(isize,itype) )
287 denswet(isize) = (xxmasswet/xxvoluwet)
288 else
289 dpwet(isize) = dpdry(isize)
290 denswet(isize) = xxdens
291 end if
292 end if
293 end do
294
295
296 !
297 ! make call to coagulation routine
298 !
299 call coagsolv( &
300 nsize, maxd_asize, ncomp_coag, ncomp_coag_maxd, &
301 temp_box, patm_box, dtcoag, nsubstep_coag, &
302 lunout_coag, idiag_coag_onoff, iok, &
303 dpdry, dpwet, denswet, num_distrib, vol_distrib )
304
305 if (iok .lt. 0) then
306 msg = '*** subr mosaic_coag_1clm -- fatal error in coagsolv'
307 call peg_message( lunout, msg )
308 call peg_error_fatal( lunout, msg )
309 else if (iok .gt. 0) then
310 ncountaa(6) = ncountaa(6) + 1
311 goto 2700
312 end if
313
314
315 !
316 ! unload coagsolv arrays
317 !
318 sumnew(:) = 0.0
319 do isize = 1, nsize
320 do ll = 1, ncomp_coag
321 sumnew(ll) = sumnew(ll) + max( 0.0, vol_distrib(isize,ll) )
322 end do
323
324 l = numptr_aer(isize,itype,iphase)
325 rsub(l,k,m) = max( 0.0, num_distrib(isize)/cair_box )
326
327 ! unload mass mixing ratios into rsub; also calculate dry mass and volume
328 xxmass = 0.0
329 xxvolu = 0.0
330 do ll = 1, ncomp_aer(itype)
331 l = massptr_aer(ll,isize,itype,iphase)
332 if (l .ge. p1st) then
333 duma = max( 0.0, vol_distrib(isize,ll) )
334 rsub(l,k,m) = duma
335 duma = duma*mw_aer(ll,itype)
336 xxmass = xxmass + duma
337 xxvolu = xxvolu + duma/dens_aer(ll,itype)
338 end if
339 end do
340 aqmassdry_sub(isize,itype,k,m) = xxmass
341 xxvolubb(isize) = xxvolu
342
343 ll = (ncomp_coag-3) + 1
344 l = hyswptr_aer(isize,itype)
345 if (l .ge. p1st) rsub(l,k,m) = max( 0.0, vol_distrib(isize,ll) )
346
347 ll = (ncomp_coag-3) + 2
348 l = waterptr_aer(isize,itype)
349 if (l .ge. p1st) rsub(l,k,m) = max( 0.0, vol_distrib(isize,ll) )
350
351 ll = (ncomp_coag-3) + 3
352 aqvoldry_sub( isize,itype,k,m) = max( 0.0, vol_distrib(isize,ll) )
353 end do
354
355
356 ! check for mass and volume conservation
357 do ll = 1, ncomp_coag
358 duma = max( sumold(ll), sumnew(ll), 1.0e-35 )
359 if (abs(sumold(ll)-sumnew(ll)) .gt. 1.0e-6*duma) then
360 ncountbb(ll) = ncountbb(ll) + 1
361 ncountbb(0) = ncountbb(0) + 1
362 end if
363 end do
364
365
366 !
367 ! calculate new dry density,
368 ! and check for new mean dry-size within bounds
369 !
370 do isize = 1, nsize
371
372 xxmass = aqmassdry_sub(isize,itype,k,m)
373 xxvolu = aqvoldry_sub( isize,itype,k,m)
374 l = numptr_aer(isize,itype,iphase)
375 xxnumb = rsub(l,k,m)
376 iconform_numb = 1
377
378 ! do a cautious calculation of density, using volume from coagsolv
379 if (xxmass .le. smallmassbb) then
380 ncountaa(8) = ncountaa(8) + 1
381 xxdens = densdefault
382 xxvolu = xxmass/xxdens
383 xxnumb = xxmass/(volumcen_sect(isize,itype)*xxdens)
384 l = hyswptr_aer(isize,itype)
385 if (l .ge. p1st) rsub(l,k,m) = 0.0
386 l = waterptr_aer(isize,itype)
387 if (l .ge. p1st) rsub(l,k,m) = 0.0
388 iconform_numb = 0
389 if (idiagbb .ge. 200) &
390 write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-7a', &
391 ktau, it, jt, k, isize, xxmass, xxvolu, xxdens
392 else if (xxmass .gt. 1000.0*xxvolu) then
393 ! in this case, density is too large. setting density=1000 forces
394 ! next IF block while avoiding potential divide by zero or overflow
395 xxdens = 1000.0
396 else
397 xxdens = xxmass/xxvolu
398 end if
399
400 if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
401 ! (exception) case -- new dry density is unrealistic,
402 ! so use dry volume from rsub instead of that from coagsolv
403 ncountaa(7) = ncountaa(7) + 1
404 if (idiagbb .ge. 200) &
405 write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-7b', &
406 ktau, it, jt, k, isize, xxmass, xxvolu, xxdens
407 xxvolu = xxvolubb(isize)
408 xxdens = xxmass/xxvolu
409 if (idiagbb .ge. 200) &
410 write(*,'(a,26x,1p,4e10.2)') 'coagaa-7c', &
411 xxmass, xxvolu, xxdens
412 end if
413
414 ! check for mean size ok, and conform number if not
415 if (iconform_numb .gt. 0) then
416 if (xxnumb .gt. xxvolu/volumlo_sect(isize,itype)) then
417 ncountaa(9) = ncountaa(9) + 1
418 duma = xxnumb
419 xxnumb = xxvolu/volumlo_sect(isize,itype)
420 if (idiagbb .ge. 200) &
421 write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-9 ', &
422 ktau, it, jt, k, isize, xxvolu, duma, xxnumb
423 else if (xxnumb .lt. xxvolu/volumhi_sect(isize,itype)) then
424 ncountaa(10) = ncountaa(10) + 1
425 duma = xxnumb
426 xxnumb = xxvolu/volumhi_sect(isize,itype)
427 if (idiagbb .ge. 200) &
428 write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-10', &
429 ktau, it, jt, k, isize, xxvolu, duma, xxnumb
430 end if
431 end if
432
433 ! load number, mass, volume, dry-density back into arrays
434 l = numptr_aer(isize,itype,iphase)
435 rsub(l,k,m) = xxnumb
436 adrydens_sub( isize,itype,k,m) = xxdens
437 aqmassdry_sub(isize,itype,k,m) = xxmass
438 aqvoldry_sub( isize,itype,k,m) = xxvolu
439
440 end do
441
442
443 2700 continue
444
445 ! temporary diagnostics
446 if ((idiagbb .ge. 100) .and. &
447 (it .eq. ite) .and. &
448 (jt .eq. jte) .and. (k .eq. kclm_calcend)) then
449 write(msg,93030) 'coagbb ncntaa ', ncountaa(1:10)
450 call peg_message( lunerr, msg )
451 if (ncountbb(0) .gt. 0) then
452 do llx = 1, (ncomp_coag+9)/10
453 llb = llx*10
454 lla = llb - 9
455 llb = min( llb, ncomp_coag)
456 write(msg,93032) 'coagbb ncntbb', &
457 mod(llx,10), ncountbb(lla:llb)
458 call peg_message( lunerr, msg )
459 end do
460 end if
461 end if
462 93020 format( a, 1p, 10e10.2 )
463 93030 format( a, 1p, 10i10 )
464 93032 format( a, 1p, i1, 10i10 )
465
466
467 2800 continue ! k levels
468
469 2900 continue ! subareas
470
471
472 return
473 end subroutine mosaic_coag_1clm
474
475
476 !----------------------------------------------------------------------
477 subroutine coagsolv( &
478 nbin, nbin_maxd, ncomp, ncomp_maxd, &
479 tkelvin_inp, patm_inp, deltat_inp, nsubstep, &
480 lunout, idiag_onoff, iok, &
481 dpdry_inp, dpwet_inp, denswet_inp, &
482 num_distrib_inp, vol_distrib_inp )
483
484 implicit none
485
486 !
487 ! *********************************************************************
488 ! ***************** written by mark jacobson (1993) *******************
489 ! *** © copyright, 1993 by mark z. jacobson ***
490 ! *** this version modified december, 2001 ***
491 ! *** (650) 723-6836 ***
492 ! *********************************************************************
493 !
494 ! cccccc ooooo a gggggg ssssss ooooo l v v
495 ! c o o a a g s o o l v v
496 ! c o o a a g ggg sssss o o l v v
497 ! c o o aaaaaaa g g s o o l v v
498 ! cccccc ooooo a a gggggg ssssss ooooo lllllll v
499 !
500 ! *********************************************************************
501 ! the use of this module implies that the user agrees not to sell the
502 ! module or any portion of it, not to remove the copyright notice from it,
503 ! and not to change the name of the module, "coagsolv". users may
504 ! modify the module as needed or incorporate it into a larger model.
505 ! any special requests with regard to this module may be directed to
506 ! jacobson@stanford.edu.
507 ! *********************************************************************
508 !
509 ! *********************************************************************
510 ! * this is a box-model version of "coagsolv," a semi-implicit *
511 ! * aerosol coagulation solver. the solver is exactly volume *
512 ! * conserving, unconditionally stable (regardless of time step), *
513 ! * and noniterative. *
514 ! * *
515 ! * this program calculates brownian coagulation kernels and solves *
516 ! * self-coagulation among any number of size bins, one *
517 ! * particle type and any number of volume fractions. *
518 ! * *
519 ! * the program is set up as a box-model. *
520 ! * *
521 ! * the volume ratio of adjacent size bins can be any number > 1 *
522 ! * *
523 ! * the scheme can be used to solve coagulation with any size bin *
524 ! * structure. *
525 ! * *
526 ! * the initial size distribution here is monomer or lognormal *
527 ! * (ifsmoluc = 1 or 0) with a fixed size bin structure *
528 ! * *
529 ! *********************************************************************
530 ! * references *
531 ! *********************************************************************
532 ! * semi-implicit scheme using movable or variable bins and with *
533 ! * any volume ratio (vrat) > 1 *
534 ! * *
535 ! * jacobson m. z., turco r. p., jensen, e. j. and toon o. b. (1994) *
536 ! * modeling coagulation among particles of different compostion *
537 ! * and size. atmos. environ. 28a, 1327 - 1338. *
538 ! * *
539 ! * jacobson m. z. (1999) fundamentals of atmospheric modeling. *
540 ! * cambridge university press, new york, 656 pp. *
541 ! * *
542 ! *********************************************************************
543 ! * semi-implicit scheme using fixed bins with volume ratio >,= 2 *
544 ! * *
545 ! * toon o. b., turco r. p., westphal d., malone r. and liu m. s. *
546 ! * (1988) a multidimensional model for aerosols: description of *
547 ! * computational analogs. j. atmos. sci. 45, 2123 - 2143 *
548 ! * *
549 ! *********************************************************************
550 ! * orig semi-implicit scheme using fixed bins with volume ratio = 2 *
551 ! * *
552 ! * turco r. p., hamill p., toon o. b., whitten r. c. and kiang c. s. *
553 ! * (1979) the nasa-ames research center stratospheric aerosol *
554 ! * model: i physical processes and computational analogs. nasa *
555 ! * tech. publ. (tp) 1362, iii-94. *
556 ! *********************************************************************
557 !
558 ! modified by y. zhang for incorporation into pnnl's mirage
559 ! june-july, 2003
560 !
561 ! *********************************************************************
562 !
563 ! modified by r.c. easter for incorporation into pnnl's wrf-chem
564 ! feb 2005 (a)
565 ! added "_inp" to all subr arguments
566 ! added iradmaxd_inp & lunout
567 ! define iradmaxd, iaertymaxd, iaeromaxd locally
568 ! no include statements
569 ! changed real*16 to real*8
570 ! feb 2005 (b)
571 ! in coagsolv, distrib_inp is 2-d array that holds both number and
572 ! volume concentrations; distrib is initialized from it;
573 ! the "fracnaer" stuff is gone
574 ! feb 2005 ©
575 ! change distrib & cc2 to be 2-d arrays (which simplifies indexing!);
576 ! iaeromaxd and iaero are gone;
577 ! deleted many commented-out lines in coagsolv
578 ! feb 2005 (d)
579 ! changed cc2(i,1) to be number instead of all-species volume;
580 ! prod term for number distrib now follows jgr-2002 article
581 ! [multiply by volume(j), divide by volume(k)]
582 ! apr 2006 (a)
583 ! considerable cleanup (mostly cosmetic)
584 ! pass in the bin sizes directly, rather than vrat & dmin_um
585 ! pass in dry and wet sizes separately
586 ! pass in/out number and volume distributions in separate arrays
587 !
588 ! *********************************************************************
589
590
591 ! IN arguments
592 integer nbin ! actual number of size bins
593 integer nbin_maxd ! size-bin dimension for input (argument) arrays
594 integer ncomp ! actual number of aerosol volume components
595 integer ncomp_maxd ! volume-component dimension for input (argument) arrays
596 integer lunout ! logical unit for warning & diagnostic output
597 integer idiag_onoff ! if positive, write some mass-conservation diagnostics
598 integer nsubstep ! number of time sub-steps for the integration
599
600 real tkelvin_inp ! air temperature (K)
601 real patm_inp ! air pressure (atm)
602 real deltat_inp ! integration time (s)
603 real dpdry_inp(nbin_maxd) ! bin dry diameter (cm)
604 real dpwet_inp(nbin_maxd) ! bin wet (ambient) diameter (cm)
605 real denswet_inp(nbin_maxd) ! bin wet (ambient) density (g/cm3)
606
607 ! IN-OUT arguments
608 real num_distrib_inp(nbin_maxd)
609 ! num_distrib_inp(i) = number concentration for bin i (#/cm3)
610 real vol_distrib_inp(nbin_maxd,ncomp_maxd)
611 ! vol_distrib_inp(i,j) = volume concentration for bin i,
612 ! component j (cm3/cm3)
613
614 ! OUT arguments
615 integer iok ! status flag (0=success, anything else=failure)
616
617 !
618 ! NOTE -- The sectional representation is based on dry particle size.
619 ! Dry sizes are used to calculate the receiving bin indices and product fractions
620 ! for each (bin-i1, bin-i2) coagulation combination.
621 ! Wet sizes and densities are used to calculate the coagulation rate coefficients.
622 !
623 ! NOTE -- Units for num_distrib and vol_distrib
624 ! num_distrib units MUST BE (#/cm3)
625 ! vol_distrib units do not really matter. They can be either masses
626 ! or volumes, and either mixing ratios or concentrations,
627 ! (e.g., g/cm3-air, ug/kg-air, cm3/m3-air, cm3/kg-air, ...).
628 ! Use whatever is convenient.
629 !
630
631 !
632 ! local variables
633 !
634 integer iradmaxd_wrk, iaertymaxd_wrk
635 parameter (iradmaxd_wrk=16)
636 parameter (iaertymaxd_wrk=32)
637
638 integer irad
639 integer iaerty
640 ! irad == nbin = actual number of size bins
641 ! iradmaxd_wrk = size-bin dimension for local (working) arrays
642 !
643 ! (iaerty-1) == ncomp = actual number of aerosol volume components
644 ! iaertymaxd_wrk = volume-component dimension for local (working) arrays
645 !
646 ! iaerty = 1 --> calc number concentration only
647 ! > 1 --> calc number concentration + (iaerty-1) volume concentrations
648
649 integer i, isubstep, j, jaer, jb, k, kout
650
651 integer jbinik(iradmaxd_wrk,iradmaxd_wrk,iradmaxd_wrk)
652 integer jbins(iradmaxd_wrk,iradmaxd_wrk)
653
654 real*8 aknud, aloss, amu, avg, &
655 boltg, bpm, &
656 cbr, consmu, cpi, &
657 deltat, deltp1, deltr, &
658 divis, dti1, dti2, &
659 finkhi, finklow, fourpi, fourrsq, &
660 ggr, gmfp, &
661 onepi, &
662 patm, pmfp, prod, &
663 radi, radiust, radj, ratior, &
664 rgas2, rho3, rsuma, rsumsq, &
665 sumdc, sumnaft, sumnbef, &
666 term1, term2, third, tk, tkelvin, tworad, &
667 viscosk, vk1, voltota, vthermg, &
668 wtair
669
670 real*8 cc2(iradmaxd_wrk,iaertymaxd_wrk), &
671 conc(iradmaxd_wrk), &
672 denav(iradmaxd_wrk) , &
673 difcof(iradmaxd_wrk), &
674 distrib(iradmaxd_wrk,iaertymaxd_wrk), &
675 fink(iradmaxd_wrk,iradmaxd_wrk,iradmaxd_wrk), &
676 floss(iradmaxd_wrk,iradmaxd_wrk), &
677 radius(iradmaxd_wrk), &
678 radwet(iradmaxd_wrk), &
679 rrate(iradmaxd_wrk,iradmaxd_wrk), &
680 sumdp(iradmaxd_wrk), &
681 sumvt(iradmaxd_wrk), &
682 tvolfin(iaertymaxd_wrk), &
683 tvolinit(iaertymaxd_wrk), &
684 volume(iradmaxd_wrk), &
685 volwet(iradmaxd_wrk), &
686 vthermp(iradmaxd_wrk)
687
688
689 ! *********************************************************************
690 ! set some parameters
691 ! *********************************************************************
692 irad = nbin
693 iaerty = ncomp + 1
694
695 kout = lunout
696
697 if (irad .gt. iradmaxd_wrk) then
698 write(lunout,*) '*** coagsolv: irad > iradmaxd_wrk: ', &
699 irad, iradmaxd_wrk
700 iok = -1
701 return
702 endif
703 if (iaerty .gt. iaertymaxd_wrk) then
704 write(lunout,*) '*** coagsolv: iaerty > iaertymaxd_wrk: ', &
705 iaerty, iaertymaxd_wrk
706 iok = -2
707 return
708 endif
709 !
710 ! tkelvin = temperature (k)
711 ! denav = particle density (g cm-3)
712 ! patm = air pressure (atm)
713 !
714 tkelvin = tkelvin_inp
715 patm = patm_inp
716 do i = 1, irad
717 denav(i) = denswet_inp(i)
718 end do
719
720 !
721 ! deltat_inp = total integration time (s)
722 ! deltat = individual time step (s)
723 ! nsubstep = number of time steps
724 !
725 deltat = deltat_inp/nsubstep
726
727 !
728 ! boltg = boltzmann's 1.38054e-16 (erg k-1 = g cm**2 sec-1 k-1)
729 ! wtair = molecular weight of air (g mol-1)
730 ! avg = avogadro's number (molec. mol-1)
731 ! rgas2 = gas constant (l-atm mol-1 k-1)
732 ! amu = dynamic viscosity of air (g cm-1 sec-1)
733 ! est value at 20 c = 1.815e-04
734 ! rho3 = air density (g cm-3)
735 ! viscosk = kinematic viscosity = amu / denair = (cm**2 sec-1)
736 ! vthermg = mean thermal velocity of air molecules (cm sec-1)
737 ! gmfp = mean free path of an air molecule = 2 x viscosk /
738 ! thermal velocity of air molecules (gmfp units of cm)
739 !
740 tk = tkelvin
741 third = 1. / 3.
742 onepi = 3.14159265358979
743 fourpi = 4. * onepi
744 boltg = 1.38054e-16
745 wtair = 28.966
746 avg = 6.02252e+23
747 rgas2 = 0.08206
748 consmu = 1.8325e-04 * (296.16 + 120.0)
749 amu = (consmu / (tk + 120.)) * (tk / 296.16)**1.5
750 rho3 = patm * wtair * 0.001 / (rgas2 * tk)
751 viscosk = amu / rho3
752 vthermg = sqrt(8. * boltg * tk * avg / (onepi * wtair))
753 gmfp = 2.0 * viscosk / vthermg
754
755 !
756 ! *********************************************************************
757 ! * set volume ratio size bin grid *
758 ! *********************************************************************
759 !
760 cpi = fourpi / 3.
761
762 do 30 j = 1, irad
763 radius( j) = dpdry_inp(j) * 0.5
764 volume( j) = cpi * radius(j) * radius(j) * radius(j)
765 radwet( j) = dpwet_inp(j) * 0.5
766 volwet( j) = cpi * radwet(j) * radwet(j) * radwet(j)
767 30 continue
768
769 !
770 ! *********************************************************************
771 ! * determine bins where coagulated terms go *
772 ! *********************************************************************
773 ! finklow = volume fraction of i+j going to lower (k ) bin
774 ! finkhi = volume fraction of i+j going to higher (k+1) bin
775 ! fink(i,j,k) = volume fraction of i+j going to bin k
776 ! floss = simplified fink term used in loss rates
777 ! no self-coagulation loss out of largest bin
778 ! jbins(i,k) = number of production terms into bin k from bin i
779 ! jbinik(i,k,jb) = identifies each bin j that coagulates with bin i
780 ! to produce bin k
781 !
782 do 35 i = 1, irad
783 do 34 j = 1, irad
784 jbins(i,j) = 0
785 floss(i,j) = 0.
786 do 33 k = 1, irad
787 jbinik(i,j,k) = 0
788 fink( i,j,k) = 0.
789 33 continue
790 34 continue
791 35 continue
792
793 do 40 i = 1, irad
794 do 39 j = 1, irad
795 voltota = volume(i) + volume(j)
796 if (voltota.ge.volume(irad)) then
797
798 if (i.eq.irad) then
799 floss(i,j) = 0.0
800 else
801 floss(i,j) = 1.0
802 endif
803
804 if (j.lt.irad) then
805 jbins(i,irad) = jbins(i,irad) + 1
806 jb = jbins(i,irad)
807 jbinik(i,irad,jb) = j
808 fink( i,jb,irad) = 1.0
809 endif
810
811 else
812 do 45 k = 1, irad - 1
813 vk1 = volume(k+1)
814 if (voltota.ge.volume(k).and.voltota.lt.vk1) then
815 finklow = ((vk1 - voltota)/(vk1 - volume(k))) &
816 * volume(k) / voltota
817 finkhi = 1. - finklow
818
819 if (i.lt.k) then
820 floss(i,j) = 1.
821 elseif (i.eq.k) then
822 floss(i,j) = finkhi
823 endif
824
825 if (j.lt.k) then
826 jbins(i,k) = jbins(i,k) + 1
827 jb = jbins(i,k)
828 jbinik(i,k,jb) = j
829 fink( i,jb,k) = finklow
830 endif
831
832 jbins(i,k+1) = jbins(i,k+1) + 1
833 jb = jbins(i,k+1)
834 jbinik(i,k+1,jb) = j
835 fink( i,jb,k+1) = finkhi
836
837 endif
838 45 continue
839 endif
840
841 do 38 k = 1, irad
842 if (jbins(i,k).gt.irad) then
843 write(21,*)'coagsolv: jbins > irad: ',jbins(i,k),i,k
844 iok = -3
845 return
846 endif
847 38 continue
848
849 39 continue
850 40 continue
851
852 !
853 ! ***********************************************************************
854 ! initialize size distribution
855 !
856 ! copy initial number concentration (#/cm3-air) and volume concentrations
857 ! (cm3/cm3-air) from num/vol_distrib_inp (real*4) to distrib (real*8)
858 ! ***********************************************************************
859 do i = 1, irad
860 distrib(i,1) = num_distrib_inp(i)
861 end do
862
863 do jaer = 2, iaerty
864 do i = 1, irad
865 distrib(i,jaer) = vol_distrib_inp(i,jaer-1)
866 end do
867 end do
868
869
870 !
871 ! *********************************************************************
872 ! coagulation kernel from fuchs equations
873 ! *********************************************************************
874 ! difcof = brownian particle diffus coef (cm**2 sec-1)
875 ! = boltg * t * bpm / (6 * pi * mu * r(i))
876 ! = 5.25e-04 (diam = 0.01um); 6.23e-7 (diam = 0.5 um) seinfeld p.325.
877 ! vthermp = avg thermal vel of particle (cm sec-1)
878 ! = (8 * boltg * t / (pi * masmolec)**0.5
879 ! pmfp = mean free path of particle (cm)
880 ! = 8 * di / (pi * vthermp)
881 ! bpm = correction for particle shape and for effects of low mean
882 ! free path of air. (toon et al., 1989, physical processes in
883 ! polar stratospheric ice clouds, jgr 94, 11,359. f1 and f2 are
884 ! included in the expression below (shape effects correction)
885 ! = 1 + kn[a + bexp(-c/kn)]
886 ! deltp1 = if particles with mean free path pmfp leave the surface of
887 ! an absorbing sphere in all directions, each being probably
888 ! equal, then deltp1 is the mean distance from the surface of the
889 ! sphere reached by particles after covering a distance pmfp. (cm).
890 ! cbr = coag kernel (cm3 partic-1 s-1) due to brownian motion (fuchs, 1964)
891 ! rrate = coag kernal (cm3 partic-1 s-1) * time step (s)
892 !
893 ! *** use the wet (ambient) radius and volume here ***
894 do 145 i = 1, irad
895 radiust = radwet(i)
896 tworad = radiust + radiust
897 fourrsq = 4. * radiust * radiust
898 aknud = gmfp / radiust
899 bpm = 1. + aknud * (1.257 + 0.42*exp(-1.1/aknud))
900 difcof(i) = boltg * tk * bpm / (6. * onepi * radiust * amu)
901
902 ! use bin-varied aerosol density from mirage
903 ! vthermp(i) = sqrt(8. * boltg * tk / (onepi*denav * volume(i)))
904 vthermp(i) = sqrt(8. * boltg * tk / (onepi*denav(i) * &
905 volwet(i)))
906 sumvt(i) = vthermp(i) * vthermp(i)
907 pmfp = 8. * difcof(i) / (onepi * vthermp(i))
908 dti1 = tworad + pmfp
909 dti2 = (fourrsq + pmfp * pmfp)**1.5
910 divis = 0.166667 / (radiust * pmfp)
911 deltp1 = divis * (dti1 * dti1 * dti1 - dti2) - tworad
912 sumdp(i) = deltp1 * deltp1
913 145 continue
914
915 !yy write(kout,9165)
916 do 155 i = 1, irad
917 do 154 j = 1, irad
918 radi = radwet(i)
919 radj = radwet(j)
920 rsuma = radi + radj
921 rsumsq = rsuma * rsuma
922 ratior = radi / radj
923 sumdc = difcof(i) + difcof(j)
924 ggr = sqrt(sumvt(i) + sumvt(j))
925 deltr = sqrt(sumdp(i) + sumdp(j))
926 term1 = rsuma / (rsuma + deltr)
927 term2 = 4. * sumdc / (rsuma * ggr)
928 cbr = fourpi * rsuma * sumdc / (term1 + term2)
929 rrate(i,j) = cbr * deltat
930 !yy write(kout,9190) (2.0e4*radius(i)), (2.0e4*radius(j)), cbr
931 154 continue
932 155 continue
933
934 9165 format(16x,'coagulation coefficients (cm**3 #-1 sec-1)'/ &
935 'diam1_um diam2_um brownian ')
936 9190 format(1pe15.7,7(1x,1pe15.7))
937
938
939 !
940 ! *********************************************************************
941 ! ****************** solve coagulation equations **********************
942 ! *********************************************************************
943 !
944 ! *********************************************************************
945 ! inialize new arrays
946 ! *********************************************************************
947 ! distrib = initial conc (# cm-3 for num conc.; cm3 cm-3 for volume fractions)
948 ! cc2 = initial, final volume conc (cm3 cm-3) for all particle types.
949 ! conc = initial, final number conc (# cm-3) for number distribution
950 ! volume = volume (cm3 #-1) of particles in bin
951 !
952 !
953
954 do i = 1, irad
955 ! cc2( i,1) = distrib(i,1) * volume(i)
956 cc2( i,1) = distrib(i,1)
957 conc(i) = distrib(i,1)
958 end do
959
960 do jaer = 2, iaerty
961 do i = 1, irad
962 cc2(i,jaer) = distrib(i,jaer)
963 end do
964 end do
965
966 ! *********************************************************************
967 ! solve coagulation over several time steps
968 ! *********************************************************************
969 !
970 do 700 isubstep = 1, nsubstep
971 !
972 do 485 jaer = 1, iaerty
973 do 484 k = 1, irad
974 !
975 ! production terms
976 !
977 prod = 0.
978 if (k.gt.1) then
979 if (jaer .eq. 1) then
980 do 465 i = 1, k
981 do 464 jb = 1, jbins(i,k)
982 j = jbinik(i,k,jb)
983 prod = prod + fink(i,jb,k)*rrate(i,j)* &
984 cc2(j,jaer)*volume(j)*conc(i)
985 464 continue
986 465 continue
987 prod = prod/volume(k)
988 else
989 do 469 i = 1, k
990 do 468 jb = 1, jbins(i,k)
991 j = jbinik(i,k,jb)
992 prod = prod + &
993 fink(i,jb,k)*rrate(i,j)*cc2(j,jaer)*conc(i)
994 468 continue
995 469 continue
996 endif
997 endif
998 !
999 ! loss terms
1000 !
1001 aloss = 0.
1002 if (k.lt.irad) then
1003 do 475 j = 1, irad
1004 aloss = aloss + floss(k,j) * rrate(k,j) * conc(j)
1005 475 continue
1006 endif
1007 !
1008 ! final concentrations
1009 !
1010 cc2(k,jaer) = (cc2(k,jaer) + prod) / (1. + aloss)
1011 484 continue
1012 485 continue
1013 !
1014 ! put updated number concentration into conc array
1015 do i = 1, irad
1016 conc(i) = cc2(i,1)
1017 end do
1018
1019 700 continue
1020
1021
1022 !
1023 ! *********************************************************
1024 ! ** put the advanced concentration back on the grid **
1025 ! (copy them from the real*8 working array to the
1026 ! real*4 caller arrays)
1027 ! *********************************************************
1028 !
1029 do i = 1, irad
1030 num_distrib_inp(i) = cc2(i,1)
1031 end do
1032 do jaer = 2, iaerty
1033 do i = 1, irad
1034 vol_distrib_inp(i,jaer-1) = cc2(i,jaer)
1035 end do
1036 end do
1037
1038 !
1039 ! if no diagnostics, then return
1040 !
1041 iok = 0
1042 if (idiag_onoff .le. 0) return
1043
1044 !
1045 ! *********************************************************************
1046 ! ************* check whether mass is conserved ***********
1047 ! *********************************************************************
1048 ! tvolinit = initial volume conc (cm3 cm-3), summed over all bins
1049 ! tvolfin = final volume conc (cm3 cm-3), summed over all bins
1050 ! sumnbef = summed number conc (# cm-3) before coagulation
1051 ! sumnaft = summed number conc (# cm-3) after coagulation
1052 !
1053 do jaer = 1, iaerty
1054 tvolinit(jaer) = 0.
1055 tvolfin( jaer) = 0.
1056 end do
1057
1058 sumnbef = 0.
1059 sumnaft = 0.
1060 !
1061 ! distrib(jaer=1, i=1:nrad) = initial total number conc in # cm-3
1062 ! cc2 (jaer=1, i=1:nrad) = final total number conc in # cm-3
1063 !
1064 do i = 1, irad
1065 tvolinit(1) = tvolinit(1) + distrib(i,1) * volume(i)
1066 tvolfin( 1) = tvolfin( 1) + cc2( i,1) * volume(i)
1067 sumnbef = sumnbef + distrib(i,1)
1068 sumnaft = sumnaft + cc2( i,1)
1069 end do
1070
1071 !
1072 ! distrib(jaer=2:iaerty, i=1:nrad) = initial component volume conc in cm3 cm-3
1073 ! cc2 (jaer=2:iaerty, i=1:nrad) = initial component volume conc in cm3 cm-3
1074 !
1075 do jaer = 2, iaerty
1076 do i = 1, irad
1077 tvolinit(jaer) = tvolinit(jaer) + distrib(i,jaer)
1078 tvolfin( jaer) = tvolfin( jaer) + cc2( i,jaer)
1079 end do
1080 end do
1081
1082 write(kout,*)
1083 write(kout,9434) sumnbef, sumnaft, &
1084 tvolinit(1)*1.0e+12,tvolfin(1)*1.0e+12
1085
1086 write(kout,9435) sumnaft-sumnbef
1087 write(kout,9439) tvolinit(1) / tvolfin(1)
1088
1089 do jaer = 2, iaerty
1090 if (abs(tvolfin(jaer)) .gt. 0.0) then
1091 write(kout,9440) jaer-1, tvolinit(jaer) / tvolfin(jaer)
1092 else
1093 write(kout,9441) jaer-1, tvolinit(jaer), tvolfin(jaer)
1094 end if
1095 end do
1096
1097 9434 format('# bef, aft; vol bef, aft =',/4(1pe16.10,1x)/)
1098 9435 format('total partic change in # cm-3 = ',1pe17.10)
1099 9439 format('total partic vol bef / vol aft = ',1pe17.10, &
1100 ': if 1.0 -> exact vol conserv')
1101 9440 format('vol comp',i4,' vol bef / vol aft = ',1pe17.10, &
1102 ': if 1.0 -> exact vol conserv')
1103 9441 format('vol comp',i4,' vol bef, vol aft = ',2(1pe17.10))
1104
1105
1106 !
1107 ! ************************** print results ****************************
1108 ! distrib = initial conc in # cm-3 for total particle (jaer=1, i=1,nrad)
1109 ! conc = final conc in # cm-3 for total particle (jaer=1, i=1,nrad)
1110 !
1111
1112
1113 !
1114 ! *********************************************************************
1115 ! end of program coagsolv.f
1116 ! *********************************************************************
1117 !
1118 return
1119 end subroutine coagsolv
1120
1121
1122
1123
1124
1125 !-----------------------------------------------------------------------
1126
1127
1128
1129 end module module_mosaic_coag