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 !wig, 11-Dec-2006: real*8 causing core dumps on IBM's (i.e. bluesky)
655 ! real*8 aknud, aloss, amu, avg, &
656 real aknud, aloss, amu, avg, &
657 boltg, bpm, &
658 cbr, consmu, cpi, &
659 deltat, deltp1, deltr, &
660 divis, dti1, dti2, &
661 finkhi, finklow, fourpi, fourrsq, &
662 ggr, gmfp, &
663 onepi, &
664 patm, pmfp, prod, &
665 radi, radiust, radj, ratior, &
666 rgas2, rho3, rsuma, rsumsq, &
667 sumdc, sumnaft, sumnbef, &
668 term1, term2, third, tk, tkelvin, tworad, &
669 viscosk, vk1, voltota, vthermg, &
670 wtair
671
672 !wig, 11-Dec-2006: real*8 causing core dumps on IBM's (i.e. bluesky)
673 ! real*8 cc2(iradmaxd_wrk,iaertymaxd_wrk), &
674 real cc2(iradmaxd_wrk,iaertymaxd_wrk), &
675 conc(iradmaxd_wrk), &
676 denav(iradmaxd_wrk) , &
677 difcof(iradmaxd_wrk), &
678 distrib(iradmaxd_wrk,iaertymaxd_wrk), &
679 fink(iradmaxd_wrk,iradmaxd_wrk,iradmaxd_wrk), &
680 floss(iradmaxd_wrk,iradmaxd_wrk), &
681 radius(iradmaxd_wrk), &
682 radwet(iradmaxd_wrk), &
683 rrate(iradmaxd_wrk,iradmaxd_wrk), &
684 sumdp(iradmaxd_wrk), &
685 sumvt(iradmaxd_wrk), &
686 tvolfin(iaertymaxd_wrk), &
687 tvolinit(iaertymaxd_wrk), &
688 volume(iradmaxd_wrk), &
689 volwet(iradmaxd_wrk), &
690 vthermp(iradmaxd_wrk)
691
692
693 ! *********************************************************************
694 ! set some parameters
695 ! *********************************************************************
696 irad = nbin
697 iaerty = ncomp + 1
698
699 kout = lunout
700
701 if (irad .gt. iradmaxd_wrk) then
702 write(lunout,*) '*** coagsolv: irad > iradmaxd_wrk: ', &
703 irad, iradmaxd_wrk
704 iok = -1
705 return
706 endif
707 if (iaerty .gt. iaertymaxd_wrk) then
708 write(lunout,*) '*** coagsolv: iaerty > iaertymaxd_wrk: ', &
709 iaerty, iaertymaxd_wrk
710 iok = -2
711 return
712 endif
713 !
714 ! tkelvin = temperature (k)
715 ! denav = particle density (g cm-3)
716 ! patm = air pressure (atm)
717 !
718 tkelvin = tkelvin_inp
719 patm = patm_inp
720 do i = 1, irad
721 denav(i) = denswet_inp(i)
722 end do
723
724 !
725 ! deltat_inp = total integration time (s)
726 ! deltat = individual time step (s)
727 ! nsubstep = number of time steps
728 !
729 deltat = deltat_inp/nsubstep
730
731 !
732 ! boltg = boltzmann's 1.38054e-16 (erg k-1 = g cm**2 sec-1 k-1)
733 ! wtair = molecular weight of air (g mol-1)
734 ! avg = avogadro's number (molec. mol-1)
735 ! rgas2 = gas constant (l-atm mol-1 k-1)
736 ! amu = dynamic viscosity of air (g cm-1 sec-1)
737 ! est value at 20 c = 1.815e-04
738 ! rho3 = air density (g cm-3)
739 ! viscosk = kinematic viscosity = amu / denair = (cm**2 sec-1)
740 ! vthermg = mean thermal velocity of air molecules (cm sec-1)
741 ! gmfp = mean free path of an air molecule = 2 x viscosk /
742 ! thermal velocity of air molecules (gmfp units of cm)
743 !
744 tk = tkelvin
745 third = 1. / 3.
746 onepi = 3.14159265358979
747 fourpi = 4. * onepi
748 boltg = 1.38054e-16
749 wtair = 28.966
750 avg = 6.02252e+23
751 rgas2 = 0.08206
752 consmu = 1.8325e-04 * (296.16 + 120.0)
753 amu = (consmu / (tk + 120.)) * (tk / 296.16)**1.5
754 rho3 = patm * wtair * 0.001 / (rgas2 * tk)
755 viscosk = amu / rho3
756 vthermg = sqrt(8. * boltg * tk * avg / (onepi * wtair))
757 gmfp = 2.0 * viscosk / vthermg
758
759 !
760 ! *********************************************************************
761 ! * set volume ratio size bin grid *
762 ! *********************************************************************
763 !
764 cpi = fourpi / 3.
765
766 do 30 j = 1, irad
767 radius( j) = dpdry_inp(j) * 0.5
768 volume( j) = cpi * radius(j) * radius(j) * radius(j)
769 radwet( j) = dpwet_inp(j) * 0.5
770 volwet( j) = cpi * radwet(j) * radwet(j) * radwet(j)
771 30 continue
772
773 !
774 ! *********************************************************************
775 ! * determine bins where coagulated terms go *
776 ! *********************************************************************
777 ! finklow = volume fraction of i+j going to lower (k ) bin
778 ! finkhi = volume fraction of i+j going to higher (k+1) bin
779 ! fink(i,j,k) = volume fraction of i+j going to bin k
780 ! floss = simplified fink term used in loss rates
781 ! no self-coagulation loss out of largest bin
782 ! jbins(i,k) = number of production terms into bin k from bin i
783 ! jbinik(i,k,jb) = identifies each bin j that coagulates with bin i
784 ! to produce bin k
785 !
786 do 35 i = 1, irad
787 do 34 j = 1, irad
788 jbins(i,j) = 0
789 floss(i,j) = 0.
790 do 33 k = 1, irad
791 jbinik(i,j,k) = 0
792 fink( i,j,k) = 0.
793 33 continue
794 34 continue
795 35 continue
796
797 do 40 i = 1, irad
798 do 39 j = 1, irad
799 voltota = volume(i) + volume(j)
800 if (voltota.ge.volume(irad)) then
801
802 if (i.eq.irad) then
803 floss(i,j) = 0.0
804 else
805 floss(i,j) = 1.0
806 endif
807
808 if (j.lt.irad) then
809 jbins(i,irad) = jbins(i,irad) + 1
810 jb = jbins(i,irad)
811 jbinik(i,irad,jb) = j
812 fink( i,jb,irad) = 1.0
813 endif
814
815 else
816 do 45 k = 1, irad - 1
817 vk1 = volume(k+1)
818 if (voltota.ge.volume(k).and.voltota.lt.vk1) then
819 finklow = ((vk1 - voltota)/(vk1 - volume(k))) &
820 * volume(k) / voltota
821 finkhi = 1. - finklow
822
823 if (i.lt.k) then
824 floss(i,j) = 1.
825 elseif (i.eq.k) then
826 floss(i,j) = finkhi
827 endif
828
829 if (j.lt.k) then
830 jbins(i,k) = jbins(i,k) + 1
831 jb = jbins(i,k)
832 jbinik(i,k,jb) = j
833 fink( i,jb,k) = finklow
834 endif
835
836 jbins(i,k+1) = jbins(i,k+1) + 1
837 jb = jbins(i,k+1)
838 jbinik(i,k+1,jb) = j
839 fink( i,jb,k+1) = finkhi
840
841 endif
842 45 continue
843 endif
844
845 do 38 k = 1, irad
846 if (jbins(i,k).gt.irad) then
847 write(21,*)'coagsolv: jbins > irad: ',jbins(i,k),i,k
848 iok = -3
849 return
850 endif
851 38 continue
852
853 39 continue
854 40 continue
855
856 !
857 ! ***********************************************************************
858 ! initialize size distribution
859 !
860 ! copy initial number concentration (#/cm3-air) and volume concentrations
861 ! (cm3/cm3-air) from num/vol_distrib_inp (real*4) to distrib (real*8)
862 ! ***********************************************************************
863 do i = 1, irad
864 distrib(i,1) = num_distrib_inp(i)
865 end do
866
867 do jaer = 2, iaerty
868 do i = 1, irad
869 distrib(i,jaer) = vol_distrib_inp(i,jaer-1)
870 end do
871 end do
872
873
874 !
875 ! *********************************************************************
876 ! coagulation kernel from fuchs equations
877 ! *********************************************************************
878 ! difcof = brownian particle diffus coef (cm**2 sec-1)
879 ! = boltg * t * bpm / (6 * pi * mu * r(i))
880 ! = 5.25e-04 (diam = 0.01um); 6.23e-7 (diam = 0.5 um) seinfeld p.325.
881 ! vthermp = avg thermal vel of particle (cm sec-1)
882 ! = (8 * boltg * t / (pi * masmolec)**0.5
883 ! pmfp = mean free path of particle (cm)
884 ! = 8 * di / (pi * vthermp)
885 ! bpm = correction for particle shape and for effects of low mean
886 ! free path of air. (toon et al., 1989, physical processes in
887 ! polar stratospheric ice clouds, jgr 94, 11,359. f1 and f2 are
888 ! included in the expression below (shape effects correction)
889 ! = 1 + kn[a + bexp(-c/kn)]
890 ! deltp1 = if particles with mean free path pmfp leave the surface of
891 ! an absorbing sphere in all directions, each being probably
892 ! equal, then deltp1 is the mean distance from the surface of the
893 ! sphere reached by particles after covering a distance pmfp. (cm).
894 ! cbr = coag kernel (cm3 partic-1 s-1) due to brownian motion (fuchs, 1964)
895 ! rrate = coag kernal (cm3 partic-1 s-1) * time step (s)
896 !
897 ! *** use the wet (ambient) radius and volume here ***
898 do 145 i = 1, irad
899 radiust = radwet(i)
900 tworad = radiust + radiust
901 fourrsq = 4. * radiust * radiust
902 aknud = gmfp / radiust
903 bpm = 1. + aknud * (1.257 + 0.42*exp(-1.1/aknud))
904 difcof(i) = boltg * tk * bpm / (6. * onepi * radiust * amu)
905
906 ! use bin-varied aerosol density from mirage
907 ! vthermp(i) = sqrt(8. * boltg * tk / (onepi*denav * volume(i)))
908 vthermp(i) = sqrt(8. * boltg * tk / (onepi*denav(i) * &
909 volwet(i)))
910 sumvt(i) = vthermp(i) * vthermp(i)
911 pmfp = 8. * difcof(i) / (onepi * vthermp(i))
912 dti1 = tworad + pmfp
913 dti2 = (fourrsq + pmfp * pmfp)**1.5
914 divis = 0.166667 / (radiust * pmfp)
915 deltp1 = divis * (dti1 * dti1 * dti1 - dti2) - tworad
916 sumdp(i) = deltp1 * deltp1
917 145 continue
918
919 !yy write(kout,9165)
920 do 155 i = 1, irad
921 do 154 j = 1, irad
922 radi = radwet(i)
923 radj = radwet(j)
924 rsuma = radi + radj
925 rsumsq = rsuma * rsuma
926 ratior = radi / radj
927 sumdc = difcof(i) + difcof(j)
928 ggr = sqrt(sumvt(i) + sumvt(j))
929 deltr = sqrt(sumdp(i) + sumdp(j))
930 term1 = rsuma / (rsuma + deltr)
931 term2 = 4. * sumdc / (rsuma * ggr)
932 cbr = fourpi * rsuma * sumdc / (term1 + term2)
933 rrate(i,j) = cbr * deltat
934 !yy write(kout,9190) (2.0e4*radius(i)), (2.0e4*radius(j)), cbr
935 154 continue
936 155 continue
937
938 9165 format(16x,'coagulation coefficients (cm**3 #-1 sec-1)'/ &
939 'diam1_um diam2_um brownian ')
940 9190 format(1pe15.7,7(1x,1pe15.7))
941
942
943 !
944 ! *********************************************************************
945 ! ****************** solve coagulation equations **********************
946 ! *********************************************************************
947 !
948 ! *********************************************************************
949 ! inialize new arrays
950 ! *********************************************************************
951 ! distrib = initial conc (# cm-3 for num conc.; cm3 cm-3 for volume fractions)
952 ! cc2 = initial, final volume conc (cm3 cm-3) for all particle types.
953 ! conc = initial, final number conc (# cm-3) for number distribution
954 ! volume = volume (cm3 #-1) of particles in bin
955 !
956 !
957
958 do i = 1, irad
959 ! cc2( i,1) = distrib(i,1) * volume(i)
960 cc2( i,1) = distrib(i,1)
961 conc(i) = distrib(i,1)
962 end do
963
964 do jaer = 2, iaerty
965 do i = 1, irad
966 cc2(i,jaer) = distrib(i,jaer)
967 end do
968 end do
969
970 ! *********************************************************************
971 ! solve coagulation over several time steps
972 ! *********************************************************************
973 !
974 do 700 isubstep = 1, nsubstep
975 !
976 do 485 jaer = 1, iaerty
977 do 484 k = 1, irad
978 !
979 ! production terms
980 !
981 prod = 0.
982 if (k.gt.1) then
983 if (jaer .eq. 1) then
984 do 465 i = 1, k
985 do 464 jb = 1, jbins(i,k)
986 j = jbinik(i,k,jb)
987 prod = prod + fink(i,jb,k)*rrate(i,j)* &
988 cc2(j,jaer)*volume(j)*conc(i)
989 464 continue
990 465 continue
991 prod = prod/volume(k)
992 else
993 do 469 i = 1, k
994 do 468 jb = 1, jbins(i,k)
995 j = jbinik(i,k,jb)
996 prod = prod + &
997 fink(i,jb,k)*rrate(i,j)*cc2(j,jaer)*conc(i)
998 468 continue
999 469 continue
1000 endif
1001 endif
1002 !
1003 ! loss terms
1004 !
1005 aloss = 0.
1006 if (k.lt.irad) then
1007 do 475 j = 1, irad
1008 aloss = aloss + floss(k,j) * rrate(k,j) * conc(j)
1009 475 continue
1010 endif
1011 !
1012 ! final concentrations
1013 !
1014 cc2(k,jaer) = (cc2(k,jaer) + prod) / (1. + aloss)
1015 484 continue
1016 485 continue
1017 !
1018 ! put updated number concentration into conc array
1019 do i = 1, irad
1020 conc(i) = cc2(i,1)
1021 end do
1022
1023 700 continue
1024
1025
1026 !
1027 ! *********************************************************
1028 ! ** put the advanced concentration back on the grid **
1029 ! (copy them from the real*8 working array to the
1030 ! real*4 caller arrays)
1031 ! *********************************************************
1032 !
1033 do i = 1, irad
1034 num_distrib_inp(i) = cc2(i,1)
1035 end do
1036 do jaer = 2, iaerty
1037 do i = 1, irad
1038 vol_distrib_inp(i,jaer-1) = cc2(i,jaer)
1039 end do
1040 end do
1041
1042 !
1043 ! if no diagnostics, then return
1044 !
1045 iok = 0
1046 if (idiag_onoff .le. 0) return
1047
1048 !
1049 ! *********************************************************************
1050 ! ************* check whether mass is conserved ***********
1051 ! *********************************************************************
1052 ! tvolinit = initial volume conc (cm3 cm-3), summed over all bins
1053 ! tvolfin = final volume conc (cm3 cm-3), summed over all bins
1054 ! sumnbef = summed number conc (# cm-3) before coagulation
1055 ! sumnaft = summed number conc (# cm-3) after coagulation
1056 !
1057 do jaer = 1, iaerty
1058 tvolinit(jaer) = 0.
1059 tvolfin( jaer) = 0.
1060 end do
1061
1062 sumnbef = 0.
1063 sumnaft = 0.
1064 !
1065 ! distrib(jaer=1, i=1:nrad) = initial total number conc in # cm-3
1066 ! cc2 (jaer=1, i=1:nrad) = final total number conc in # cm-3
1067 !
1068 do i = 1, irad
1069 tvolinit(1) = tvolinit(1) + distrib(i,1) * volume(i)
1070 tvolfin( 1) = tvolfin( 1) + cc2( i,1) * volume(i)
1071 sumnbef = sumnbef + distrib(i,1)
1072 sumnaft = sumnaft + cc2( i,1)
1073 end do
1074
1075 !
1076 ! distrib(jaer=2:iaerty, i=1:nrad) = initial component volume conc in cm3 cm-3
1077 ! cc2 (jaer=2:iaerty, i=1:nrad) = initial component volume conc in cm3 cm-3
1078 !
1079 do jaer = 2, iaerty
1080 do i = 1, irad
1081 tvolinit(jaer) = tvolinit(jaer) + distrib(i,jaer)
1082 tvolfin( jaer) = tvolfin( jaer) + cc2( i,jaer)
1083 end do
1084 end do
1085
1086 write(kout,*)
1087 write(kout,9434) sumnbef, sumnaft, &
1088 tvolinit(1)*1.0e+12,tvolfin(1)*1.0e+12
1089
1090 write(kout,9435) sumnaft-sumnbef
1091 write(kout,9439) tvolinit(1) / tvolfin(1)
1092
1093 do jaer = 2, iaerty
1094 if (abs(tvolfin(jaer)) .gt. 0.0) then
1095 write(kout,9440) jaer-1, tvolinit(jaer) / tvolfin(jaer)
1096 else
1097 write(kout,9441) jaer-1, tvolinit(jaer), tvolfin(jaer)
1098 end if
1099 end do
1100
1101 9434 format('# bef, aft; vol bef, aft =',/4(1pe16.10,1x)/)
1102 9435 format('total partic change in # cm-3 = ',1pe17.10)
1103 9439 format('total partic vol bef / vol aft = ',1pe17.10, &
1104 ': if 1.0 -> exact vol conserv')
1105 9440 format('vol comp',i4,' vol bef / vol aft = ',1pe17.10, &
1106 ': if 1.0 -> exact vol conserv')
1107 9441 format('vol comp',i4,' vol bef, vol aft = ',2(1pe17.10))
1108
1109
1110 !
1111 ! ************************** print results ****************************
1112 ! distrib = initial conc in # cm-3 for total particle (jaer=1, i=1,nrad)
1113 ! conc = final conc in # cm-3 for total particle (jaer=1, i=1,nrad)
1114 !
1115
1116
1117 !
1118 ! *********************************************************************
1119 ! end of program coagsolv.f
1120 ! *********************************************************************
1121 !
1122 return
1123 end subroutine coagsolv
1124
1125
1126
1127
1128
1129 !-----------------------------------------------------------------------
1130
1131
1132
1133 end module module_mosaic_coag