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