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.
7 ! MOSAIC module: see module_mosaic_driver.F for references and terms of use
8 !**********************************************************************************
9 module module_mosaic_movesect
12 use module_data_mosaic_asect
13 use module_data_mosaic_other
18 ! if apply_n1_inflow = 1, then subr move_sections_apply_n1_inflow
19 ! applies an ad_hoc correction to bin 1 to account for growth of
20 ! smaller particles (which are not simulated when aernucnew_onoff=0)
22 ! if apply_n1_inflow = other, this correction is not applied
24 integer, parameter :: apply_n1_inflow = 0
30 !-----------------------------------------------------------------------
32 ! zz08movesect.f - created 24-nov-01 from scm movebin code
33 ! 04-dec-01 rce - added code to treat bins that had state="no_aerosol"
34 ! 10-dec-01 rce - diagnostic writes to fort.96 deleted
35 ! 08-aug-02 rce - this is latest version from freshair scaqs-aerosol code
36 ! 03-aug-02 rce - bypass this routine when msectional=701
37 ! output nnewsave values to lunout when msectional=702
38 ! 29-aug-03 rce - use nspec_amode_nontracer in first "do ll" loop
39 ! 12-nov-03 rce - two approaches now available
40 ! jacobson moving-center (old) - when msectional=10
41 ! tzivion mass-number advection (new) - when msectional=20
42 ! 28-jan-04 rce - eliminated the move_sections_old code
43 ! (no reason to keep it)
45 ! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
46 ! variables in module_data_mosaic_asect
47 ! rce 2005-feb-17 - fixes involving drydens_pre/aftgrow, drymass_...,
48 ! & mw_aer. All are dimensioned (isize,itype) but were being used
49 ! as (itype). In old compiler, this was OK, and treated as (itype,1).
50 ! In new compiler, this caused an error.
51 ! Also, in subr test_move_sections, set iphase based on iflag.
52 !-----------------------------------------------------------------------
55 !-----------------------------------------------------------------------
56 subroutine move_sections( iflag, iclm, jclm, k, m)
58 ! routine transfers aerosol number and mass between sections
59 ! to account for continuous aerosol growth
60 ! this routine is called after the gas condensation module (MOSAIC) or
61 ! aqueous chemistry module has increased the mass within sections
63 ! moving-center algorithm or mass-number advection algorithm is used,
64 ! depending on value of mod(msectional,100)
65 ! section mean diameter is given by
66 ! vtot = ntot * (pi/6) * (dmean**3)
67 ! where vtot and ntot are total dry-volume and number for the section
68 ! if dmean is outside the section boundaries (dlo_sect & dhi_sect), then
69 ! all the mass and number in the section are transfered to the
70 ! section with dlo_sect(nnew) < dmean < dhi_sect(nnew)
72 ! mass mixing ratios (mole/mole-air or g/mole-air) are in
73 ! rsub(massptr_aer(ll,n,itype,iphase),k,m)
74 ! number mixing ratios (particles/mole-air) are in
75 ! rsub(numptr_aer(n,itype,iphase),k,m)
76 ! these values are over-written with new values
77 ! the following are also updated:
78 ! adrydens_sub(n,itype,k,m), admeandry_sub(n,itype,k,m),
79 ! awetdens_sub(n,itype,k,m), admeanwet_sub(n,itype,k,m)
81 ! particle sizes are in cm
84 ! iflag = 1 - do transfer after trace-gas condensation/evaporation
85 ! = 2 - do transfer after aqueous chemistry
86 ! = -1/-2 - do some "first entry" tasks for the iflag=+1/+2 cases
87 ! iclm, jclm, k = current i,j,k indices
88 ! m = current subarea index
89 ! drymass_pregrow(n,itype) = dry-mass (g/mole-air) for section n
91 ! drymass_aftgrow(n,itype) = dry-mass (g/mole-air) for section n
92 ! after the growth but before inter-section transfer
93 ! drydens_pregrow(n,itype) = dry-density (g/cm3) for section n
95 ! drydens_aftgrow(n,itype) = dry-density (g/cm3) for section n
96 ! after the growth but before inter-section transfer
97 ! (drymass_pregrow and drydens_pregrow are not used by the moving-center
98 ! algorithm, but would be needed for other sectional algorithms)
100 ! this routine is the top level module for the post-october-2003 code
101 ! and will do either moving-center or mass-number advection
112 integer iflag, iclm, jclm, k, m
115 integer idiag_movesect, iphase, itype, &
116 l, ll, llhysw, llwater, lnew, lold, l3, &
117 method_movesect, n, nnew, nold
118 integer nnewsave(2,maxd_asize)
120 real densdefault, densh2o, smallmassaa, smallmassbb
121 real delta_water_conform1, delta_numb_conform1
123 real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
124 drydensxx(maxd_asize), drydensyy(maxd_asize)
125 real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
126 drymassxx(maxd_asize), drymassyy(maxd_asize)
127 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
128 real rmassxx(maxd_acomp+2,maxd_asize), &
129 rmassyy(maxd_acomp+2,maxd_asize)
130 real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
131 rnumbxx(maxd_asize), rnumbyy(maxd_asize)
132 real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
133 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
134 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
135 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
141 ! check for valid inputs
143 if (msectional .le. 0) return
144 if (ntype_aer .le. 0) return
145 if (nphase_aer .le. 0) return
149 ! run diagnostic tests
150 ! (these will only be run for certain values of idiag_movesect
151 ! rce 2004-nov-29 - to avoid recursion, test_move_sections
152 ! is now called from mosaic_driver
154 ! call test_move_sections( iflag, iclm, jclm, k, m )
157 ! get "method_movesect" from digits 1-2 of msectional (treat 1-9 as 10)
158 method_movesect = mod( msectional, 100 )
159 if (method_movesect .le. 10) method_movesect = 10
161 ! get "idiag_movesect" from digits 3-4 of msectional
162 idiag_movesect = mod( msectional, 10000 )/100
164 if ((method_movesect .eq. 10) .or. &
165 (method_movesect .eq. 11) .or. &
166 (method_movesect .eq. 20) .or. &
167 (method_movesect .eq. 21) .or. &
168 (method_movesect .eq. 30) .or. &
169 (method_movesect .eq. 31)) then
171 else if ((method_movesect .eq. 19) .or. &
172 (method_movesect .eq. 29) .or. &
173 (method_movesect .eq. 39)) then
176 msg = '*** subr move_sections error - ' // &
177 'illegal value for msectional'
178 call peg_error_fatal( lunerr, msg )
181 if (iabs(iflag) .eq. 1) then
183 else if (iabs(iflag) .eq. 2) then
185 if (nphase_aer .lt. 2) then
186 msg = '*** subr move_sections error - ' // &
187 'iflag=2 (after aqueous chemistry) but nphase_aer < 2'
188 call peg_error_fatal( lunerr, msg )
189 else if (cw_phase .ne. 2) then
190 msg = '*** subr move_sections error - ' // &
191 'iflag=2 (after aqueous chemistry) but cw_phase .ne. 2'
192 call peg_error_fatal( lunerr, msg )
195 msg = '*** subr move_sections error - ' // &
196 'iabs(iflag) must be 1 or 2'
197 call peg_error_fatal( lunerr, msg )
201 ! when iflag=-1/-2, call move_sections_checkptrs then return
202 ! if ((ncorecnt .le. 0) .and. (k .le. 1)) then
203 if (iflag .le. 0) then
204 write(msg,9040) 'method', method_movesect
205 call peg_message( lunout, msg )
206 write(msg,9040) 'idiag ', idiag_movesect
207 call peg_message( lunout, msg )
208 call move_sections_checkptrs( iflag, iclm, jclm, k, m )
211 9040 format( '*** subr move_sections - ', a, ' =', i6 )
214 if (idiag_movesect .eq. 70) then
216 call peg_message( lunout, msg )
217 write(msg,9060) iclm, jclm, k, m, msectional
218 call peg_message( lunout, msg )
220 9060 format( '*** subr move_sections diagnostics i,j,k,m,msect =', 4i4, i6 )
226 ! if bin mass mixrat < smallmassaa (1.e-20 g/mol), then assume no growth
227 ! AND no water AND conform number so that size is within bin limits
228 smallmassaa = 1.0e-20
229 ! if bin mass OR volume mixrat < smallmassbb (1.e-30 g/mol), then
230 ! assume default density to avoid divide by zero
231 smallmassbb = 1.0e-30
234 ! process each type, one at a time
235 do 1900 itype = 1, ntype_aer
237 if (nsize_aer(itype) .le. 0) goto 1900
239 call move_sections_initial_conform( &
240 iflag, iclm, jclm, k, m, iphase, itype, &
241 method_movesect, idiag_movesect, llhysw, llwater, &
242 densdefault, densh2o, smallmassaa, smallmassbb, &
243 delta_water_conform1, delta_numb_conform1, &
244 drydenspp, drymasspp, rnumbpp, &
245 drydensxx0, drymassxx0, rnumbxx0, &
246 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
247 wetmassxx, wetvolxx, &
248 specmwxx, specdensxx )
250 if (method_movesect .le. 19) then
251 call move_sections_calc_movingcenter( &
252 iflag, iclm, jclm, k, m, iphase, itype, &
253 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
254 densdefault, densh2o, smallmassaa, smallmassbb, &
255 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
256 wetmassxx, wetvolxx, &
257 xferfracvol, xferfracnum )
259 call move_sections_calc_masnumadvect( &
260 iflag, iclm, jclm, k, m, iphase, itype, &
261 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
262 densdefault, densh2o, smallmassaa, smallmassbb, &
263 drydenspp, drymasspp, rnumbpp, &
264 drydensxx0, drymassxx0, rnumbxx0, &
265 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
266 wetmassxx, wetvolxx, &
267 xferfracvol, xferfracnum )
270 call move_sections_apply_moves( &
271 iflag, iclm, jclm, k, m, iphase, itype, &
272 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
273 densdefault, densh2o, smallmassaa, smallmassbb, &
274 delta_water_conform1, delta_numb_conform1, &
275 drydenspp, drymasspp, rnumbpp, &
276 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
277 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
278 xferfracvol, xferfracnum )
280 ! rce 08-nov-2004 - this call is new (and perhaps temporary)
281 ! rce 05-jul-2005 - do n1_inflow for aerchem growth but not for cldchem
282 if ((apply_n1_inflow .eq. 1) .and. &
283 (iphase .eq. ai_phase)) then
284 call move_sections_apply_n1_inflow( &
285 iflag, iclm, jclm, k, m, iphase, itype, &
286 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
287 densdefault, densh2o, smallmassaa, smallmassbb, &
288 delta_water_conform1, delta_numb_conform1, &
289 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
290 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
291 xferfracvol, xferfracnum, &
292 specmwxx, specdensxx )
295 ! call move_sections_final_conform( &
296 ! iflag, iclm, jclm, k, m, iphase, itype )
301 end subroutine move_sections
304 !-----------------------------------------------------------------------
305 subroutine move_sections_initial_conform( &
306 iflag, iclm, jclm, k, m, iphase, itype, &
307 method_movesect, idiag_movesect, llhysw, llwater, &
308 densdefault, densh2o, smallmassaa, smallmassbb, &
309 delta_water_conform1, delta_numb_conform1, &
310 drydenspp, drymasspp, rnumbpp, &
311 drydensxx0, drymassxx0, rnumbxx0, &
312 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
313 wetmassxx, wetvolxx, &
314 specmwxx, specdensxx )
317 ! routine does some initial tasks for the section movement
318 ! load rmassxx & rnumbxx from rsub
319 ! load specmwxx & specdensxx
320 ! set drymassxx & dryvolxx from drymass_aftgrow & drydens_aftgrow,
321 ! OR compute them from rmassxx, specmwxx, specdensxx if need be
322 ! set wetmassxx & wetvolxx from dry values & water mass
323 ! conform rnumbxx so that the mean particle size of each section
324 ! (= dryvolxx/rnumbxx) is within the section limits
335 integer iflag, iclm, jclm, iphase, itype, k, &
336 m, method_movesect, idiag_movesect, llhysw, llwater
337 real densdefault, densh2o, smallmassaa, smallmassbb
338 real delta_water_conform1, delta_numb_conform1
339 real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
340 drydensxx(maxd_asize)
341 real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
342 drymassxx(maxd_asize)
343 real dryvolxx(maxd_asize)
344 real rmassxx(maxd_acomp+2,maxd_asize)
345 real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
347 real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
348 real wetvolxx(maxd_asize)
349 real wetmassxx(maxd_asize)
353 integer l, ll, lnew, lold, l3, n, nnew, nold
355 real dummass, dumnum, dumnum_at_dhi, dumnum_at_dlo, dumr, &
356 dumvol, dumvol1p, dumwatrmass
359 ! assure positive definite
361 rsub(l,k,m) = max( 0., rsub(l,k,m) )
364 ! load mixrats into working arrays and assure positive definite
365 llhysw = ncomp_plustracer_aer(itype) + 1
366 llwater = ncomp_plustracer_aer(itype) + 2
367 do n = 1, nsize_aer(itype)
368 do ll = 1, ncomp_plustracer_aer(itype)
369 l = massptr_aer(ll,n,itype,iphase)
370 rmassxx(ll,n) = rsub(l,k,m)
372 rmassxx(llhysw,n) = 0.
374 if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
375 if (l .gt. 0) rmassxx(llhysw,n) = rsub(l,k,m)
376 rmassxx(llwater,n) = 0.
378 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
379 if (l .gt. 0) rmassxx(llwater,n) = rsub(l,k,m)
381 rnumbxx(n) = rsub(numptr_aer(n,itype,iphase),k,m)
382 rnumbxx0(n) = rnumbxx(n)
383 rnumbpp(n) = rnumbxx(n)
385 drydenspp(n) = drydens_pregrow(n,itype)
386 drymasspp(n) = drymass_pregrow(n,itype)
388 drydensxx(n) = drydens_aftgrow(n,itype)
389 drymassxx(n) = drymass_aftgrow(n,itype)
390 drydensxx0(n) = drydensxx(n)
391 drymassxx0(n) = drymassxx(n)
394 ! load specmw and specdens also
395 do ll = 1, ncomp_plustracer_aer(itype)
396 specmwxx(ll) = mw_aer(ll,itype)
397 specdensxx(ll) = dens_aer(ll,itype)
400 delta_water_conform1 = 0.0
401 delta_numb_conform1 = 0.0
404 do 1390 n = 1, nsize_aer(itype)
407 ! if drydens_aftgrow < 0.1, then bin had state="no_aerosol"
408 ! compute volume using default dry-densities, set water=0,
409 ! and conform the number
410 ! also do this if mass is extremely small (below smallmassaa)
411 ! OR if drydens_aftgrow > 20 (which is unrealistic)
413 if ( (drydensxx(n) .lt. 0.1) .or. &
414 (drydensxx(n) .gt. 20.0) .or. &
415 (drymassxx(n) .le. smallmassaa) ) then
418 do ll = 1, ncomp_aer(itype)
419 dumr = rmassxx(ll,n)*specmwxx(ll)
420 dummass = dummass + dumr
421 dumvol = dumvol + dumr/specdensxx(ll)
423 drymassxx(n) = dummass
424 if (min(dummass,dumvol) .le. smallmassbb) then
425 drydensxx(n) = densdefault
426 dumvol = dummass/densdefault
427 dumnum = dummass/(volumcen_sect(n,itype)*densdefault)
429 drydensxx(n) = dummass/dumvol
431 dumnum_at_dhi = dumvol/volumhi_sect(n,itype)
432 dumnum_at_dlo = dumvol/volumlo_sect(n,itype)
433 dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
435 delta_numb_conform1 = delta_numb_conform1 + dumnum - rnumbxx(n)
437 rnumbpp(n) = rnumbxx(n)
438 delta_water_conform1 = delta_water_conform1 - rmassxx(llwater,n)
439 rmassxx(llwater,n) = 0.
442 ! load dry/wet mass and volume into "xx" arrays
443 ! which hold values before inter-mode transferring
444 dryvolxx(n) = drymassxx(n)/drydensxx(n)
445 dumwatrmass = rmassxx(llwater,n)*mw_water_aer
446 wetmassxx(n) = drymassxx(n) + dumwatrmass
447 wetvolxx(n) = dryvolxx(n) + dumwatrmass/densh2o
452 end subroutine move_sections_initial_conform
455 !-----------------------------------------------------------------------
456 subroutine move_sections_calc_movingcenter( &
457 iflag, iclm, jclm, k, m, iphase, itype, &
458 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
459 densdefault, densh2o, smallmassaa, smallmassbb, &
460 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
461 wetmassxx, wetvolxx, &
462 xferfracvol, xferfracnum )
464 ! routine calculates section movements for the moving-center approach
466 ! material in section n will be transfered to section nnewsave(1,n)
468 ! the nnewsave are calculated here
469 ! the actual transfer is done in another routine
480 integer iflag, iclm, jclm, iphase, itype, k, &
481 m, method_movesect, idiag_movesect, llhysw, llwater
482 integer nnewsave(2,maxd_asize)
483 real densdefault, densh2o, smallmassaa, smallmassbb
484 real drydensxx(maxd_asize)
485 real drymassxx(maxd_asize)
486 real dryvolxx(maxd_asize)
487 real rmassxx(maxd_acomp+2,maxd_asize)
488 real rnumbxx(maxd_asize)
489 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
490 real wetmassxx(maxd_asize)
491 real wetvolxx(maxd_asize)
494 integer ll, n, ndum, nnew, nold
495 real dumnum, dumvol, dumvol1p, sixoverpi, third
503 ! compute mean size after growth (and corresponding section)
504 ! particles in section n will be transferred to section nnewsave(1,n)
506 do 1390 n = 1, nsize_aer(itype)
510 ! don't bother to transfer bins whose mass is extremely small
511 if (drymassxx(n) .le. smallmassaa) goto 1290
516 ! check for number so small that particle volume is
517 ! above that of largest section
518 if (dumnum .le. dumvol/volumhi_sect(nsize_aer(itype),itype)) then
519 nnew = nsize_aer(itype)
521 ! or below that of smallest section
522 else if (dumnum .ge. dumvol/volumlo_sect(1,itype)) then
527 ! dumvol1p is mean particle volume (cm3) for the section
528 dumvol1p = dumvol/dumnum
529 if (dumvol1p .gt. volumhi_sect(n,itype)) then
530 do while ( (nnew .lt. nsize_aer(itype)) .and. &
531 (dumvol1p .gt. volumhi_sect(nnew,itype)) )
535 else if (dumvol1p .lt. volumlo_sect(n,itype)) then
536 do while ( (nnew .gt. 1) .and. &
537 (dumvol1p .lt. volumlo_sect(nnew,itype)) )
543 1290 nnewsave(1,n) = nnew
546 xferfracvol(1,n) = 1.0
547 xferfracvol(2,n) = 0.0
548 xferfracnum(1,n) = 1.0
549 xferfracnum(2,n) = 0.0
555 ! output nnewsave values when msectional=7xxx
556 if (idiag_movesect .eq. 70) then
558 do n = 1, nsize_aer(itype)
559 if (nnewsave(1,n) .ne. n) ndum = ndum + 1
561 if (ndum .gt. 0) then
562 write(msg,97751) 'YES', iclm, jclm, k, m, &
563 ndum, (nnewsave(1,n), n=1,nsize_aer(itype))
564 call peg_message( lunout, msg )
566 write(msg,97751) 'NO ', iclm, jclm, k, m, &
567 ndum, (nnewsave(1,n), n=1,nsize_aer(itype))
568 call peg_message( lunout, msg )
571 97751 format( 'movesect', a, 4i3, 3x, i3, 3x, 14i3 )
574 end subroutine move_sections_calc_movingcenter
577 !-----------------------------------------------------------------------
578 subroutine move_sections_calc_masnumadvect( &
579 iflag, iclm, jclm, k, m, iphase, itype, &
580 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
581 densdefault, densh2o, smallmassaa, smallmassbb, &
582 drydenspp, drymasspp, rnumbpp, &
583 drydensxx0, drymassxx0, rnumbxx0, &
584 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
585 wetmassxx, wetvolxx, &
586 xferfracvol, xferfracnum )
588 ! routine calculates section movements for the mass-number-advection approach
590 ! material in section n will be transfered to sections
591 ! nnewsave(1,n) and nnewsave(2,n)
592 ! the fractions of mass/volume transfered to each are
593 ! xferfracvol(1,n) and xferfracvol(2,n)
594 ! the fractions of number transfered to each are
595 ! xferfracnum(1,n) and xferfracnum(2,n)
597 ! the nnewsave, xferfracvol, and xferfracnum are calculated here
598 ! the actual transfer is done in another routine
609 integer iflag, iclm, jclm, iphase, itype, k, &
610 m, method_movesect, idiag_movesect, llhysw, llwater
611 integer nnewsave(2,maxd_asize)
613 real densdefault, densh2o, smallmassaa, smallmassbb
614 real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
615 drydensxx(maxd_asize)
616 real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
617 drymassxx(maxd_asize)
618 real dryvolxx(maxd_asize)
619 real rmassxx(maxd_acomp+2,maxd_asize)
620 real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
622 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
623 real wetvolxx(maxd_asize)
624 real wetmassxx(maxd_asize)
627 integer ierr, n, nnew, nnew2
628 integer iforce_movecenter(maxd_asize)
630 real dum1, dum2, dum3
631 real dumaa, dumbb, dumgamma, dumratio
632 real dumfracnum, dumfracvol
635 real dumvbar_aft, dumvbar_pre
636 real dumvcutlo_nnew_pre, dumvcuthi_nnew_pre
637 real dumvlo_pre, dumvhi_pre, dumvdel_pre
638 real dumvtot_aft, dumvtot_pre
640 real sixoverpi, third
651 ! compute mean size after growth (and corresponding section)
652 ! some of the particles in section n will be transferred to section nnewsave(1,n)
654 ! if the aftgrow mass is extremely small,
655 ! OR if the aftgrow mean size is outside of
656 ! [dlo_sect(1,itype), dhi_sect(nsize_aer(itype),itype)]
657 ! then use the moving-center method_movesect for this bin
658 ! (don't try to compute the pregrow within-bin distribution)
660 do 3900 n = 1, nsize_aer(itype)
663 iforce_movecenter(n) = 0
665 xferfracvol(1,n) = 1.0
666 xferfracvol(2,n) = 0.0
667 xferfracnum(1,n) = 1.0
668 xferfracnum(2,n) = 0.0
678 dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
679 dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
683 ! don't bother to transfer bins whose mass is extremely small
684 if (drymassxx(n) .le. smallmassaa) then
685 iforce_movecenter(n) = 1
689 dumvtot_aft = dryvolxx(n)
692 ! check for particle volume above that of largest section
693 ! or below that of smallest section
694 if (dumntot .le. dumvtot_aft/volumhi_sect(nsize_aer(itype),itype)) then
695 nnew = nsize_aer(itype)
696 iforce_movecenter(n) = 2
698 else if (dumntot .ge. dumvtot_aft/volumlo_sect(1,itype)) then
700 iforce_movecenter(n) = 3
704 ! dumvbar_aft is mean particle volume (cm3) for the section
705 ! find the section that encloses this volume
706 dumvbar_aft = dumvtot_aft/dumntot
707 if (dumvbar_aft .gt. volumhi_sect(n,itype)) then
708 do while ( (nnew .lt. nsize_aer(itype)) .and. &
709 (dumvbar_aft .gt. volumhi_sect(nnew,itype)) )
713 else if (dumvbar_aft .lt. volumlo_sect(n,itype)) then
714 do while ( (nnew .gt. 1) .and. &
715 (dumvbar_aft .lt. volumlo_sect(nnew,itype)) )
721 1290 nnewsave(1,n) = nnew
724 if (iforce_movecenter(n) .gt. 0) goto 3700
727 ! if drydenspp (pregrow) < 0.1 (because bin had state="no_aerosol" before
728 ! growth was computed, so its initial mass was very small)
729 ! then use the moving-center method_movesect for this bin
730 ! (don't try to compute the pregrow within-bin distribution)
731 ! also do this if pregrow mass is extremely small (below smallmassaa)
732 ! OR if drydenspp > 20 (unphysical)
733 if ( (drydenspp(n) .lt. 0.1) .or. &
734 (drydenspp(n) .gt. 20.0) .or. &
735 (drymasspp(n) .le. smallmassaa) ) then
736 iforce_movecenter(n) = 11
740 dumvtot_pre = drymasspp(n)/drydenspp(n)
742 dumvlo_pre = volumlo_sect(n,itype)
743 dumvhi_pre = volumhi_sect(n,itype)
744 dumvdel_pre = dumvhi_pre - dumvlo_pre
746 ! if the pregrow mean size is outside of OR very close to the bin limits,
747 ! then use moving-center approach for this bin
748 dumv = dumvhi_pre - 0.01*dumvdel_pre
749 if (dumntot .le. dumvtot_pre/dumv) then
750 iforce_movecenter(n) = 12
753 dumv = dumvlo_pre + 0.01*dumvdel_pre
754 if (dumntot .ge. dumvtot_pre/dumv) then
755 iforce_movecenter(n) = 13
759 ! calculate the pregrow within-section size distribution
760 dumvbar_pre = dumvtot_pre/dumntot
761 dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
762 dumratio = dumvbar_pre/dumvlo_pre
764 if (dumratio .le. (1.0001 + dumgamma/3.0)) then
765 dumv = dumvlo_pre + 3.0*(dumvbar_pre-dumvlo_pre)
766 dumvhi_pre = min( dumvhi_pre, dumv )
767 dumvdel_pre = dumvhi_pre - dumvlo_pre
768 dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
769 dumratio = dumvbar_pre/dumvlo_pre
770 else if (dumratio .ge. (0.9999 + dumgamma*2.0/3.0)) then
771 dumv = dumvhi_pre + 3.0*(dumvbar_pre-dumvhi_pre)
772 dumvlo_pre = max( dumvlo_pre, dumv )
773 dumvdel_pre = dumvhi_pre - dumvlo_pre
774 dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
775 dumratio = dumvbar_pre/dumvlo_pre
778 dumbb = (dumratio - 1.0 - 0.5*dumgamma)*12.0/dumgamma
779 dumaa = 1.0 - 0.5*dumbb
781 ! calculate pregrow volumes corresponding to the nnew
783 dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
784 dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
786 ! if the [dumvlo_pre, dumvhi_pre] falls completely within
787 ! the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval,
788 ! then all mass and number go to nnew
789 if (nnew .eq. 1) then
790 if (dumvhi_pre .le. dumvcuthi_nnew_pre) then
791 iforce_movecenter(n) = 21
795 else if (nnew .eq. nsize_aer(itype)) then
796 if (dumvlo_pre .ge. dumvcutlo_nnew_pre) then
797 iforce_movecenter(n) = 22
802 if ((dumvlo_pre .ge. dumvcutlo_nnew_pre) .and. &
803 (dumvhi_pre .le. dumvcuthi_nnew_pre)) then
804 iforce_movecenter(n) = 23
805 else if (dumvlo_pre .lt. dumvcutlo_nnew_pre) then
811 if (iforce_movecenter(n) .gt. 0) goto 3700
813 ! calculate the fraction of ntot and vtot that are within
814 ! the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval
815 dumzlo = (dumvcutlo_nnew_pre - dumvlo_pre)/dumvdel_pre
816 dumzhi = (dumvcuthi_nnew_pre - dumvlo_pre)/dumvdel_pre
817 dumzlo = max( dumzlo, 0.0 )
818 dumzhi = min( dumzhi, 1.0 )
819 dum1 = dumzhi - dumzlo
820 dum2 = (dumzhi**2 - dumzlo**2)*0.5
821 dum3 = (dumzhi**3 - dumzlo**3)/3.0
822 dumfracnum = dumaa*dum1 + dumbb*dum2
823 dumfracvol = (dumvlo_pre/dumvbar_pre) * (dumaa*dum1 + &
824 (dumaa*dumgamma + dumbb)*dum2 + (dumbb*dumgamma)*dum3)
826 if ((dumfracnum .le. 0.0) .or. (dumfracvol .le. 0.0)) then
827 iforce_movecenter(n) = 31
828 nnewsave(1,n) = nnew2
829 else if ((dumfracnum .ge. 1.0) .or. (dumfracvol .ge. 1.0)) then
830 iforce_movecenter(n) = 32
832 if (iforce_movecenter(n) .gt. 0) goto 3700
834 nnewsave(2,n) = nnew2
836 xferfracvol(1,n) = dumfracvol
837 xferfracvol(2,n) = 1.0 - dumfracvol
838 xferfracnum(1,n) = dumfracnum
839 xferfracnum(2,n) = 1.0 - dumfracnum
843 ! output nnewsave values when msectional=7xxx
844 if (idiag_movesect .ne. 70) goto 3800
846 if (nnewsave(2,n) .eq. 0) then
847 if (nnewsave(1,n) .eq. 0) then
849 else if (nnewsave(1,n) .eq. n) then
854 else if (nnewsave(1,n) .eq. 0) then
855 if (nnewsave(2,n) .eq. n) then
860 else if (nnewsave(2,n) .eq. n) then
861 if (nnewsave(1,n) .eq. n) then
866 else if (nnewsave(1,n) .eq. n) then
873 if (drymasspp(n) .gt. drymassxx(n)) dumch1 = '-'
876 call peg_message( lunout, msg )
877 write(msg,97010) dumch1, dumch4, iclm, jclm, k, m, &
878 n, nnewsave(1,n), nnewsave(2,n), iforce_movecenter(n)
879 call peg_message( lunout, msg )
880 write(msg,97020) 'pre mass, dens ', &
881 drymasspp(n), drydenspp(n)
882 call peg_message( lunout, msg )
883 write(msg,97020) 'aft mass, dens, numb', &
884 drymassxx(n), drydensxx(n), rnumbxx(n)
885 call peg_message( lunout, msg )
886 if ((drydensxx(n) .ne. drydensxx0(n)) .or. &
887 (drymassxx(n) .ne. drymassxx0(n)) .or. &
888 (rnumbxx(n) .ne. rnumbxx0(n) )) then
889 write(msg,97020) 'aft0 mas, dens, numb', &
890 drymassxx0(n), drydensxx0(n), rnumbxx0(n)
891 call peg_message( lunout, msg )
893 write(msg,97020) 'vlop0, vbarp, vhip0', &
894 volumlo_sect(n,itype), dumvbar_pre, volumhi_sect(n,itype)
895 call peg_message( lunout, msg )
896 write(msg,97020) 'vlop , vbarp, vhip ', &
897 dumvlo_pre, dumvbar_pre, dumvhi_pre
898 call peg_message( lunout, msg )
899 write(msg,97020) 'vloax, vbarax, vhiax', &
900 dumvcutlo_nnew_pre, dumvbar_pre, dumvcuthi_nnew_pre
901 call peg_message( lunout, msg )
902 write(msg,97020) 'vloa0, vbara, vhia0', &
903 volumlo_sect(nnew,itype), dumvbar_aft, volumhi_sect(nnew,itype)
904 call peg_message( lunout, msg )
905 write(msg,97020) 'dumfrvol, num, ratio', &
906 dumfracvol, dumfracnum, dumratio
907 call peg_message( lunout, msg )
908 write(msg,97020) 'frvol,num1; vol,num2', &
909 xferfracvol(1,n), xferfracnum(1,n), &
910 xferfracvol(2,n), xferfracnum(2,n)
911 call peg_message( lunout, msg )
913 97010 format( 'movesect', 2a, 7x, 4i3, 4x, &
914 'n,nnews', 3i3, 4x, 'iforce', i3.2 )
915 97020 format( a, 1p, 4e13.4 )
920 ! check for legal combinations of nnewsave(1,n) & nnewsave(2,n)
923 ! both are non-zero AND iabs(nnew1-nnew2) != 1
925 if (nnewsave(1,n) .eq. nnewsave(2,n)) then
927 else if (nnewsave(1,n)*nnewsave(2,n) .ne. 0) then
928 if (iabs(nnewsave(1,n)-nnewsave(2,n)) .ne. 1) ierr = 1
930 if (ierr .gt. 0) then
931 write(msg,97010) 'E', 'RROR', iclm, jclm, k, m, &
932 n, nnewsave(1,n), nnewsave(2,n), iforce_movecenter(n)
933 call peg_message( lunout, msg )
937 ! if method_movesect == 30-31 then force moving center
938 ! this is just for testing purposes
939 if ((method_movesect .ge. 30) .and. (method_movesect .le. 39)) then
942 xferfracvol(1,n) = 1.0
943 xferfracvol(2,n) = 0.0
944 xferfracnum(1,n) = 1.0
945 xferfracnum(2,n) = 0.0
951 end subroutine move_sections_calc_masnumadvect
954 !-----------------------------------------------------------------------
955 subroutine move_sections_apply_moves( &
956 iflag, iclm, jclm, k, m, iphase, itype, &
957 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
958 densdefault, densh2o, smallmassaa, smallmassbb, &
959 delta_water_conform1, delta_numb_conform1, &
960 drydenspp, drymasspp, rnumbpp, &
961 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
962 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
963 xferfracvol, xferfracnum )
965 ! routine performs the actual transfer of aerosol number and mass
977 integer iflag, iclm, jclm, iphase, itype, k, &
978 m, method_movesect, idiag_movesect, llhysw, llwater
979 integer nnewsave(2,maxd_asize)
981 real densdefault, densh2o, smallmassaa, smallmassbb
982 real delta_water_conform1, delta_numb_conform1
983 real drydenspp(maxd_asize)
984 real drymasspp(maxd_asize)
985 real drymassxx(maxd_asize), drymassyy(maxd_asize)
986 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
987 real rmassxx(maxd_acomp+2,maxd_asize), &
988 rmassyy(maxd_acomp+2,maxd_asize)
989 real rnumbpp(maxd_asize)
990 real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
991 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
992 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
993 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
996 integer jj, l, ll, n, ndum, nnew, nold
997 integer jja, jjb, jjc
999 real delta_numb_conform2, &
1000 dumbot, dumnum, dumnum_at_dhi, dumnum_at_dlo, &
1001 dumvol, dumvol1p, dumxfvol, dumxfnum, sixoverpi, third
1002 real dumpp(maxd_asize), dumxx(maxd_asize), dumyy(maxd_asize), &
1014 ! initialize "yy" arrays that hold values after inter-mode transferring
1015 ! "yy" = "xx" for sections that do not move at all
1016 ! "yy" = 0.0 for sections that do move (partially or completely)
1018 do 1900 n = 1, nsize_aer(itype)
1020 if ( (nnewsave(1,n) .eq. n) .and. &
1021 (nnewsave(2,n) .eq. 0) ) then
1022 ! if nnew == n, then material in section n will not be transferred, and
1023 ! section n will contain its initial material plus any material
1024 ! transferred from other sections
1025 ! so initialize "yy" arrays with "xx" values
1026 drymassyy(n) = drymassxx(n)
1027 dryvolyy(n) = dryvolxx(n)
1028 wetmassyy(n) = wetmassxx(n)
1029 wetvolyy(n) = wetvolxx(n)
1030 rnumbyy(n) = rnumbxx(n)
1031 do ll = 1, ncomp_plustracer_aer(itype) + 2
1032 rmassyy(ll,n) = rmassxx(ll,n)
1036 ! if nnew .ne. n, then material in section n will be transferred, and
1037 ! section n will only contain material that is transferred from
1039 ! so initialize "yy" arrays to zero
1045 do ll = 1, ncomp_plustracer_aer(itype) + 2
1054 ! do the transfer of mass and number
1056 do 2900 nold = 1, nsize_aer(itype)
1058 if ( (nnewsave(1,nold) .eq. nold) .and. &
1059 (nnewsave(2,nold) .eq. 0 ) ) goto 2900
1063 nnew = nnewsave(jj,nold)
1064 if (nnew .le. 0) goto 2800
1066 dumxfvol = xferfracvol(jj,nold)
1067 dumxfnum = xferfracnum(jj,nold)
1069 do ll = 1, ncomp_plustracer_aer(itype) + 2
1070 rmassyy(ll,nnew) = rmassyy(ll,nnew) + rmassxx(ll,nold)*dumxfvol
1072 rnumbyy(nnew) = rnumbyy(nnew) + rnumbxx(nold)*dumxfnum
1074 drymassyy(nnew) = drymassyy(nnew) + drymassxx(nold)*dumxfvol
1075 dryvolyy(nnew) = dryvolyy(nnew) + dryvolxx(nold)*dumxfvol
1076 wetmassyy(nnew) = wetmassyy(nnew) + wetmassxx(nold)*dumxfvol
1077 wetvolyy(nnew) = wetvolyy(nnew) + wetvolxx(nold)*dumxfvol
1084 ! transfer among sections is completed
1085 ! - check for conservation of mass/volume/number
1086 ! - conform number again
1087 ! - compute/store densities and mean sizes
1088 ! - if k=1, save values for use by dry deposition routine
1089 ! - copy new mixrats back to rsub array
1091 call move_sections_conserve_check( &
1092 1, iflag, iclm, jclm, k, m, iphase, itype, &
1093 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1094 densdefault, densh2o, smallmassaa, smallmassbb, &
1095 delta_water_conform1, delta_numb_conform1, &
1096 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1097 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1099 delta_numb_conform2 = 0.0
1101 do 3900 n = 1, nsize_aer(itype)
1103 dumvol = dryvolyy(n)
1104 if (min(drymassyy(n),dumvol) .le. smallmassbb) then
1105 dumvol = drymassyy(n)/densdefault
1106 dumnum = drymassyy(n)/(volumlo_sect(n,itype)*densdefault)
1107 delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n)
1109 adrydens_sub(n,itype,k,m) = densdefault
1110 awetdens_sub(n,itype,k,m) = densdefault
1111 admeandry_sub(n,itype,k,m) = dcen_sect(n,itype)
1112 admeanwet_sub(n,itype,k,m) = dcen_sect(n,itype)
1115 dumnum_at_dhi = dumvol/volumhi_sect(n,itype)
1116 dumnum_at_dlo = dumvol/volumlo_sect(n,itype)
1117 dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
1118 delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n)
1120 adrydens_sub(n,itype,k,m) = drymassyy(n)/dumvol
1121 dumvol1p = dumvol/dumnum
1122 admeandry_sub(n,itype,k,m) = (dumvol1p*sixoverpi)**third
1123 awetdens_sub(n,itype,k,m) = wetmassyy(n)/wetvolyy(n)
1124 dumvol1p = wetvolyy(n)/dumnum
1125 admeanwet_sub(n,itype,k,m) = min( 100.*dcen_sect(n,itype), &
1126 (dumvol1p*sixoverpi)**third )
1128 aqvoldry_sub(n,itype,k,m) = dumvol
1129 aqmassdry_sub(n,itype,k,m) = drymassyy(n)
1131 ! if (k .eq. 1) then
1132 ! awetdens_sfc(n,itype,iclm,jclm) = awetdens_sub(n,itype,k,m)
1133 ! admeanwet_sfc(n,itype,iclm,jclm) = admeanwet_sub(n,itype,k,m)
1136 do ll = 1, ncomp_plustracer_aer(itype)
1137 l = massptr_aer(ll,n,itype,iphase)
1138 rsub(l,k,m) = rmassyy(ll,n)
1141 if (iphase .eq. ai_phase) then
1142 l = waterptr_aer(n,itype)
1143 if (l .gt. 0) rsub(l,k,m) = rmassyy(llwater,n)
1144 l = hyswptr_aer(n,itype)
1145 if (l .gt. 0) rsub(l,k,m) = rmassyy(llhysw,n)
1147 rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
1151 delta_numb_conform1 = delta_numb_conform1 + delta_numb_conform2
1153 call move_sections_conserve_check( &
1154 2, iflag, iclm, jclm, k, m, iphase, itype, &
1155 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1156 densdefault, densh2o, smallmassaa, smallmassbb, &
1157 delta_water_conform1, delta_numb_conform1, &
1158 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1159 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1163 ! output nnewsave values when msectional=7xxx
1164 if (idiag_movesect .ne. 70) goto 4900
1167 do n = 1, nsize_aer(itype)
1168 if (nnewsave(1,n)+nnewsave(2,n) .ne. n) ndum = ndum + 1
1170 ! if (ndum .gt. 0) then
1171 ! write(msg,97010) 'SOME', iclm, jclm, k, m, &
1172 ! ndum, (nnewsave(1,n), nnewsave(2,n), n=1,nsize_aer(itype))
1173 ! call peg_message( lunout, msg )
1175 ! write(msg,97010) 'NONE', iclm, jclm, k, m, &
1176 ! ndum, (nnewsave(1,n), nnewsave(2,n), n=1,nsize_aer(itype))
1177 ! call peg_message( lunout, msg )
1181 if (ndum .gt. 0) dumch4 = 'SOME'
1183 call peg_message( lunout, msg )
1184 write(msg,97010) dumch4, iclm, jclm, k, m, ndum
1185 call peg_message( lunout, msg )
1186 do jjb = 1, nsize_aer(itype), 10
1187 jjc = min( jjb+9, nsize_aer(itype) )
1188 write(msg,97011) (nnewsave(1,n), nnewsave(2,n), n=jjb,jjc)
1189 call peg_message( lunout, msg )
1192 ! write(msg,97020) 'rnumbold', (rnumbxx(n), n=1,nsize_aer(itype))
1193 ! call peg_message( lunout, msg )
1194 ! write(msg,97020) 'rnumbnew', (rnumbyy(n), n=1,nsize_aer(itype))
1195 ! call peg_message( lunout, msg )
1197 ! write(msg,97020) 'drvolold', (dryvolxx(n), n=1,nsize_aer(itype))
1198 ! call peg_message( lunout, msg )
1199 ! write(msg,97020) 'drvolnew', (dryvolyy(n), n=1,nsize_aer(itype))
1200 ! call peg_message( lunout, msg )
1202 dumbot = log( volumhi_sect(1,itype)/volumlo_sect(1,itype) )
1203 do n = 1, nsize_aer(itype)
1207 if ( (drydenspp(n) .gt. 0.5) .and. &
1208 (drymasspp(n) .gt. smallmassaa) ) then
1209 dumvol = drymasspp(n)/drydenspp(n)
1210 if ((rnumbpp(n) .ge. 1.0e-35) .and. &
1211 (dumvol .ge. 1.0e-35)) then
1212 dumvol1p = dumvol/rnumbpp(n)
1213 dumpp(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
1216 if ((rnumbxx(n) .ge. 1.0e-35) .and. &
1217 (dryvolxx(n) .ge. 1.0e-35)) then
1218 dumvol1p = dryvolxx(n)/rnumbxx(n)
1219 dumxx(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
1221 if ((rnumbyy(n) .ge. 1.0e-35) .and. &
1222 (dryvolyy(n) .ge. 1.0e-35)) then
1223 dumvol1p = dryvolyy(n)/rnumbyy(n)
1224 dumyy(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
1228 ! write(msg,97030) 'lnvolold', (dumxx(n), n=1,nsize_aer(itype))
1229 ! call peg_message( lunout, msg )
1230 ! write(msg,97030) 'lnvolnew', (dumyy(n), n=1,nsize_aer(itype))
1231 ! call peg_message( lunout, msg )
1234 if (jja .eq. 1) then
1236 dumout(:) = rnumbxx(:)
1237 else if (jja .eq. 2) then
1239 dumout(:) = rnumbyy(:)
1240 else if (jja .eq. 3) then
1242 dumout(:) = dryvolxx(:)
1243 else if (jja .eq. 4) then
1245 dumout(:) = dryvolyy(:)
1246 else if (jja .eq. 5) then
1248 dumout(:) = dumxx(:)
1249 else if (jja .eq. 6) then
1251 dumout(:) = dumyy(:)
1252 else if (jja .eq. 7) then
1254 dumout(:) = dumpp(:)
1256 do jjb = 1, nsize_aer(itype), 10
1257 jjc = min( jjb+9, nsize_aer(itype) )
1258 if (jja .le. 4) then
1259 write(msg,97020) dumch8, (dumout(n), n=jjb,jjc)
1261 write(msg,97030) dumch8, (dumout(n), n=jjb,jjc)
1263 call peg_message( lunout, msg )
1268 !97010 format( / 'movesectapply', a, 4i3, 3x, i3 / 5x, 10(3x,2i3) )
1269 !97020 format( a, 1p, 10e9.1 / (( 8x, 1p, 10e9.1 )) )
1270 !97030 format( a, 10f9.3 / (( 8x, 10f9.3 )) )
1271 97010 format( 'movesectapply', a, 4i3, 3x, i3 )
1272 97011 format( 5x, 10(3x,2i3) )
1273 97020 format( a, 1p, 10e9.1 )
1274 97030 format( a, 10f9.3 )
1278 end subroutine move_sections_apply_moves
1281 !-----------------------------------------------------------------------
1282 ! rce 08-nov-2004 - this routine is new (and perhaps temporary)
1284 subroutine move_sections_apply_n1_inflow( &
1285 iflag, iclm, jclm, k, m, iphase, itype, &
1286 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1287 densdefault, densh2o, smallmassaa, smallmassbb, &
1288 delta_water_conform1, delta_numb_conform1, &
1289 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1290 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
1291 xferfracvol, xferfracnum, &
1292 specmwxx, specdensxx )
1294 ! routine applies an ad_hoc correction to bin 1 to account for
1295 ! growth of smaller particles (which are not simulated) growing into
1297 ! the correction to particle number balances the loss of particles from
1298 ! bin 1 to larger bins by growth
1299 ! the correction to mass assumes that the particles coming into bin 1
1300 ! are slightly larger than dlo_sect(1)
1307 ! include 'v33com9a'
1308 ! include 'v33com9b'
1311 integer iflag, iclm, jclm, iphase, itype, k, &
1312 m, method_movesect, idiag_movesect, llhysw, llwater
1313 integer nnewsave(2,maxd_asize)
1315 real densdefault, densh2o, smallmassaa, smallmassbb
1316 real delta_water_conform1, delta_numb_conform1
1317 real drymassxx(maxd_asize), drymassyy(maxd_asize)
1318 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
1319 real rmassxx(maxd_acomp+2,maxd_asize), &
1320 rmassyy(maxd_acomp+2,maxd_asize)
1321 real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
1322 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
1323 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
1324 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
1325 real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
1329 integer jj, l, ll, n, nnew, nold
1331 real deltanum, deltavol, dumvol1p, dumxfvol, dumxfnum
1335 ! compute fraction of number transferred out of bin 1 by growth
1339 if ( (nnewsave(1,nold) .eq. nold) .and. &
1340 (nnewsave(2,nold) .eq. 0 ) ) goto 3900
1344 nnew = nnewsave(jj,nold)
1345 if (nnew .le. 0) goto 2800
1346 if (nnew .eq. nold) goto 2800
1347 dumxfnum = dumxfnum + xferfracnum(jj,nold)
1351 ! compute "inflow" change to number and volume
1352 ! number change matches that lost by growth
1353 ! volume change assume inflow particles slightly bigger then dlo_sect
1355 dumvol1p = 0.95*volumlo_sect(n,itype) + 0.05*volumhi_sect(n,itype)
1356 deltanum = dumxfnum*rnumbxx(n)
1357 deltavol = deltanum*dumvol1p
1358 if (dumxfnum .le. 0.0) goto 3900
1359 if (deltanum .le. 0.0) goto 3900
1360 if (deltavol .le. 0.0) goto 3900
1363 ! increment the number and masses
1364 ! if the old dryvol (dryvolxx) > smallmassbb, then compute mass increment for
1365 ! each species from the old masses (rmassxx)
1366 ! otherwise only increment the first mass species
1368 if (dryvolxx(n) .gt. smallmassbb) then
1369 dumxfvol = deltavol/dryvolxx(n)
1370 do ll = 1, ncomp_plustracer_aer(itype) + 2
1371 rmassyy(ll,n) = rmassyy(ll,n) + dumxfvol*rmassxx(ll,n)
1375 rmassyy(ll,n) = rmassyy(ll,n) + deltavol*specdensxx(ll)/specmwxx(ll)
1377 rnumbyy(n) = rnumbyy(n) + deltanum
1380 ! transfer results into rsub
1382 do ll = 1, ncomp_plustracer_aer(itype)
1383 l = massptr_aer(ll,n,itype,iphase)
1384 rsub(l,k,m) = rmassyy(ll,n)
1386 if (iphase .eq. ai_phase) then
1387 l = waterptr_aer(n,itype)
1388 if (l .gt. 0) rsub(l,k,m) = rmassyy(llwater,n)
1389 l = hyswptr_aer(n,itype)
1390 if (l .gt. 0) rsub(l,k,m) = rmassyy(llhysw,n)
1392 rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
1397 end subroutine move_sections_apply_n1_inflow
1400 !-----------------------------------------------------------------------
1401 subroutine move_sections_conserve_check( ipass, &
1402 iflag, iclm, jclm, k, m, iphase, itype, &
1403 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1404 densdefault, densh2o, smallmassaa, smallmassbb, &
1405 delta_water_conform1, delta_numb_conform1, &
1406 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1407 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1409 ! routine checks for conservation of number, mass, and volume
1410 ! by the move_sections algorithm
1413 ! initialize all thesum(jj,ll) to zero
1414 ! computes thesum(1,ll) from rmassxx, rnumbxx, ...
1415 ! computes thesum(2,ll) from rmassyy, rnumbyy, ...
1416 ! compares thesum(1,ll) with thesum(2,ll)
1417 ! computes thesum(3,ll) from rsub before section movement
1419 ! computes thesum(4,ll) from rsub after section movement
1420 ! compares thesum(3,ll) with thesum(4,ll)
1422 ! currently only implemented for condensational growth (iflag=1)
1429 ! include 'v33com9a'
1430 ! include 'v33com9b'
1433 integer ipass, iflag, iclm, jclm, iphase, itype, k, &
1434 m, method_movesect, idiag_movesect, llhysw, llwater
1435 integer nnewsave(2,maxd_asize)
1436 real densdefault, densh2o, smallmassaa, smallmassbb
1437 real delta_water_conform1, delta_numb_conform1
1438 real drymassxx(maxd_asize), drymassyy(maxd_asize)
1439 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
1440 real rmassxx(maxd_acomp+2,maxd_asize), &
1441 rmassyy(maxd_acomp+2,maxd_asize)
1442 real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
1443 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
1444 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
1447 integer jj, l, ll, llworst, llworstb, n
1448 integer nerr, nerrmax
1450 data nerr, nerrmax / 0, 999 /
1452 real dumbot, dumtop, dumtoler, dumerr, dumworst, dumworstb
1453 real duma, dumb, dumc, dume
1454 real thesum(4,maxd_acomp+7)
1457 character*8 dumname(maxd_acomp+7)
1461 if (ipass .eq. 2) goto 2000
1463 do ll = 1, ncomp_plustracer_aer(itype)+7
1469 do n = 1, nsize_aer(itype)
1470 do ll = 1, ncomp_plustracer_aer(itype)+2
1471 thesum(1,ll) = thesum(1,ll) + rmassxx(ll,n)
1472 thesum(2,ll) = thesum(2,ll) + rmassyy(ll,n)
1474 ll = ncomp_plustracer_aer(itype)+3
1475 thesum(1,ll) = thesum(1,ll) + rnumbxx(n)
1476 thesum(2,ll) = thesum(2,ll) + rnumbyy(n)
1477 ll = ncomp_plustracer_aer(itype)+4
1478 thesum(1,ll) = thesum(1,ll) + drymassxx(n)
1479 thesum(2,ll) = thesum(2,ll) + drymassyy(n)
1480 ll = ncomp_plustracer_aer(itype)+5
1481 thesum(1,ll) = thesum(1,ll) + dryvolxx(n)
1482 thesum(2,ll) = thesum(2,ll) + dryvolyy(n)
1483 ll = ncomp_plustracer_aer(itype)+6
1484 thesum(1,ll) = thesum(1,ll) + wetmassxx(n)
1485 thesum(2,ll) = thesum(2,ll) + wetmassyy(n)
1486 ll = ncomp_plustracer_aer(itype)+7
1487 thesum(1,ll) = thesum(1,ll) + wetvolxx(n)
1488 thesum(2,ll) = thesum(2,ll) + wetvolyy(n)
1494 ! calc sum over bins for each species
1495 ! for water, account for loss in initial conform (delta_water_conform1)
1497 do n = 1, nsize_aer(itype)
1498 do ll = 1, ncomp_plustracer_aer(itype)+3
1499 if (ll .le. ncomp_plustracer_aer(itype)) then
1500 l = massptr_aer(ll,n,itype,iphase)
1501 else if (ll .eq. ncomp_plustracer_aer(itype)+1) then
1503 if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
1504 else if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1506 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1508 l = numptr_aer(n,itype,iphase)
1511 thesum(ipass+2,ll) = thesum(ipass+2,ll) + rsub(l,k,m)
1514 if (ipass .eq. 2) then
1515 ll = ncomp_plustracer_aer(itype)+2
1516 thesum(3,ll) = thesum(3,ll) + delta_water_conform1
1517 ll = ncomp_plustracer_aer(itype)+3
1518 thesum(3,ll) = thesum(3,ll) + delta_numb_conform1
1522 ! now compare either sum1-sum2 or sum3-sum4
1523 ! on ipass=1, jj=1, so compare sum1 & sum2
1524 ! on ipass=2, jj=3, so compare sum3 & sum4
1526 do ll = 1, ncomp_plustracer_aer(itype)+7
1528 write(dumname(ll),'(i4.4)') ll
1529 if (ll .le. ncomp_plustracer_aer(itype)) dumname(ll) = &
1530 name(massptr_aer(ll,1,itype,iphase))(1:4)
1531 if (ll .eq. ncomp_plustracer_aer(itype)+1) dumname(ll) = 'hysw'
1532 if (ll .eq. ncomp_plustracer_aer(itype)+2) dumname(ll) = 'watr'
1533 if (ll .eq. ncomp_plustracer_aer(itype)+3) dumname(ll) = 'numb'
1534 if (ll .eq. ncomp_plustracer_aer(itype)+4) dumname(ll) = 'drymass'
1535 if (ll .eq. ncomp_plustracer_aer(itype)+5) dumname(ll) = 'dryvol'
1536 if (ll .eq. ncomp_plustracer_aer(itype)+6) dumname(ll) = 'wetmass'
1537 if (ll .eq. ncomp_plustracer_aer(itype)+7) dumname(ll) = 'wetvol'
1545 do ll = 1, ncomp_plustracer_aer(itype)+7
1546 dumtop = thesum(jj+1,ll) - thesum(jj,ll)
1547 ! dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35 )
1548 ! change minimum dumbot to 1.0e-30 to reduce unimportant error messages
1549 dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-30 )
1550 dumerr = dumtop/dumbot
1552 ! rce 21-jul-2006 - encountered some cases when delta_*_conform1 is negative
1553 ! and large in magnitude relative to thesum, which causes the mass
1554 ! conservation to be less accurate due to roundoff
1555 ! following section recomputes relative error with delta_*_conform1
1556 ! added onto each of thesum
1557 ! also increased dumtoler slightly
1558 if (ipass .eq. 2) then
1560 if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1561 dumc = delta_water_conform1
1562 else if (ll .eq. ncomp_plustracer_aer(itype)+3) then
1563 dumc = delta_numb_conform1
1565 if (dumc .lt. 0.0) then
1566 duma = thesum(3,ll) - dumc
1567 dumb = thesum(4,ll) - dumc
1568 dumtop = dumb - duma
1569 ! dumbot = max( abs(duma), abs(dumb), 1.0e-35 )
1570 dumbot = max( abs(duma), abs(dumb), 1.0e-30 )
1571 dume = dumtop/dumbot
1572 if (abs(dume) .lt. abs(dumerr)) dumerr = dume
1576 if (abs(dumerr) .gt. abs(dumworst)) then
1578 dumworstb = dumworst
1585 if (abs(dumworst) .gt. dumtoler) then
1587 if (nerr .le. nerrmax) then
1589 call peg_message( lunout, msg )
1590 write(msg,97110) iclm, jclm, k, m, ipass, llworst
1591 call peg_message( lunout, msg )
1592 write(msg,97120) ' nnew(1,n)', &
1593 (nnewsave(1,n), n=1,nsize_aer(itype))
1594 call peg_message( lunout, msg )
1595 write(msg,97120) ' nnew(2,n)', &
1596 (nnewsave(2,n), n=1,nsize_aer(itype))
1597 call peg_message( lunout, msg )
1600 if (ll .eq. 0) ll = ncomp_plustracer_aer(itype)+2
1601 write(msg,97130) 'name/relerr/thesum', jj, '/thesum', jj+1, &
1602 dumname(ll), dumworst, thesum(jj,ll), thesum(jj+1,ll)
1603 call peg_message( lunout, msg )
1605 if ( (ll .eq. ncomp_plustracer_aer(itype)+3) .and. &
1606 (abs(dumworstb) .gt. dumtoler) ) then
1607 ll = max( 1, llworstb )
1608 dumtop = thesum(jj+1,ll) - thesum(jj,ll)
1609 ! dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35 )
1610 dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-30 )
1611 dumerr = dumtop/dumbot
1612 write(msg,97130) 'name/relerr/thesum', jj, '/thesum', jj+1, &
1613 dumname(ll), dumerr, thesum(jj,ll), thesum(jj+1,ll)
1614 call peg_message( lunout, msg )
1619 97110 format( 'movesect conserve ERROR - i/j/k/m/pass/llworst', &
1621 97120 format( a, 64i3 )
1622 97130 format( a, i1, a, i1, 2x, a, 1p, 3e16.7 )
1625 end subroutine move_sections_conserve_check
1628 !-----------------------------------------------------------------------
1629 subroutine test_move_sections( iflag_test, iclm, jclm, k, m )
1631 ! routine runs tests on move_sections, using a matrix of
1632 ! pregrow and aftgrow masses for each section
1639 ! include 'v33com9a'
1640 ! include 'v33com9b'
1643 integer iflag_test, iclm, jclm, k, m
1646 integer idiag_movesect, iflag, ii, iphase, itype, jj, &
1652 real dumnumb, dumvolpre, dumvolaft
1653 real dumsv_rsub(l2maxd)
1654 real dumsv_drymass_pregrow(maxd_asize,maxd_atype)
1655 real dumsv_drymass_aftgrow(maxd_asize,maxd_atype)
1656 real dumsv_drydens_pregrow(maxd_asize,maxd_atype)
1657 real dumsv_drydens_aftgrow(maxd_asize,maxd_atype)
1661 integer maxvolfactpre, maxvolfactaft
1662 parameter (maxvolfactpre=15, maxvolfactaft=23)
1664 real dumvolfactpre(maxvolfactpre)
1665 data dumvolfactpre / &
1666 2.0, 0.0, 1.0e-20, 0.5, 0.9, &
1667 1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0, &
1670 real dumvolfactaft(maxvolfactaft)
1671 data dumvolfactaft / &
1672 4.0, 0.0, 1.0e-20, 0.01, 0.02, 0.05, 0.1, 0.5, 0.9, &
1673 1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0, &
1674 8.1, 16.0, 32., 64., 128., 256. /
1678 ! check for valid inputs
1681 if (msectional .le. 0) return
1682 if (nsize_aer(1) .le. 0) return
1684 idiag_movesect = mod( msectional, 10000 )/100
1685 if (idiag_movesect .ne. 70) return
1687 ientryno = ientryno + 1
1688 if (ientryno .gt. 2) return
1691 ! save rsub and drymass/dens_aft/pregrow
1694 dumsv_rsub(l) = rsub(l,k,m)
1696 do itype = 1, ntype_aer
1697 do n = 1, nsize_aer(itype)
1698 dumsv_drymass_pregrow(n,itype) = drymass_pregrow(n,itype)
1699 dumsv_drymass_aftgrow(n,itype) = drymass_aftgrow(n,itype)
1700 dumsv_drydens_pregrow(n,itype) = drydens_pregrow(n,itype)
1701 dumsv_drydens_aftgrow(n,itype) = drydens_aftgrow(n,itype)
1706 ! make test calls to move_sections
1708 do 3900 iflag = 1, 1
1711 if (iabs(iflag) .eq. 2) iphase = cw_phase
1713 do 3800 itype = 1, ntype_aer
1715 do 2900 nn = 1, nsize_aer(itype)
1717 do 2800 ii = 1, maxvolfactpre
1719 do 2700 jj = 1, maxvolfactaft
1721 ! zero out rsub and dryxxx_yyygrow arrays
1722 do n = 1, nsize_aer(itype)
1723 do ll = 1, ncomp_plustracer_aer(itype)
1724 rsub(massptr_aer(ll,n,itype,iphase),k,m) = 0.0
1727 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1728 if (l .gt. 0) rsub(l,k,m) = 0.0
1729 l = numptr_aer(n,itype,iphase)
1730 if (l .gt. 0) rsub(l,k,m) = 0.0
1731 drymass_pregrow(n,itype) = 0.0
1732 drymass_aftgrow(n,itype) = 0.0
1733 drydens_pregrow(n,itype) = -1.0
1734 drydens_aftgrow(n,itype) = -1.0
1737 ! fill in values for section nn
1740 rsub(numptr_aer(n,itype,iphase),k,m) = dumnumb
1742 l = massptr_aer(ll,n,itype,iphase)
1744 dumvolpre = volumlo_sect(n,itype)*dumvolfactpre(ii)*dumnumb
1745 drydens_pregrow(n,itype) = dens_aer(ll,itype)
1746 drymass_pregrow(n,itype) = dumvolpre*drydens_pregrow(n,itype)
1747 if (ii .eq. 1) drydens_pregrow(n,itype) = -1.0
1749 dumvolaft = volumlo_sect(n,itype)*dumvolfactaft(jj)*dumnumb
1750 drydens_aftgrow(n,itype) = dens_aer(ll,itype)
1751 drymass_aftgrow(n,itype) = dumvolaft*drydens_aftgrow(n,itype)
1752 if (jj .eq. 1) drydens_aftgrow(n,itype) = -1.0
1754 rsub(l,k,m) = drymass_aftgrow(n,itype)/mw_aer(ll,itype)
1757 call peg_message( lunout, msg )
1758 write(msg,98010) nn, ii, jj
1759 call peg_message( lunout, msg )
1760 write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
1761 call peg_message( lunout, msg )
1763 ! make test call to move_sections
1764 call move_sections( iflag, iclm, jclm, k, m )
1767 call peg_message( lunout, msg )
1768 write(msg,98010) nn, ii, jj
1769 call peg_message( lunout, msg )
1770 write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
1771 call peg_message( lunout, msg )
1780 98010 format( 'test_move_sections output - nn, ii, jj =', 3i3 )
1781 98011 format( 'volfactpre, volfactaft =', 1p, 2e12.4 )
1784 ! restore rsub and drymass/dens_aft/pregrow
1787 rsub(l,k,m) = dumsv_rsub(l)
1789 do itype = 1, ntype_aer
1790 do n = 1, nsize_aer(itype)
1791 drymass_pregrow(n,itype) = dumsv_drymass_pregrow(n,itype)
1792 drymass_aftgrow(n,itype) = dumsv_drymass_aftgrow(n,itype)
1793 drydens_pregrow(n,itype) = dumsv_drydens_pregrow(n,itype)
1794 drydens_aftgrow(n,itype) = dumsv_drydens_aftgrow(n,itype)
1799 end subroutine test_move_sections
1802 !-----------------------------------------------------------------------
1803 subroutine move_sections_checkptrs( iflag, iclm, jclm, k, m )
1805 ! checks for valid number and water pointers
1812 ! include 'v33com9a'
1813 ! include 'v33com9b'
1816 integer iflag, iclm, jclm, k, m
1819 integer l, itype, iphase, n, ndum
1822 do 1900 itype = 1, ntype_aer
1823 do 1800 iphase = 1, nphase_aer
1826 do n = 1, nsize_aer(itype)
1827 l = numptr_aer(n,itype,iphase)
1828 if ((l .le. 0) .or. (l .gt. ltot2)) then
1829 msg = '*** subr move_sections error - ' // &
1830 'numptr_amode not defined'
1831 call peg_message( lunerr, msg )
1832 write(msg,9030) 'mode, numptr =', n, l
1833 call peg_message( lunerr, msg )
1834 write(msg,9030) 'iphase, itype =', iphase, itype
1835 call peg_message( lunerr, msg )
1836 call peg_error_fatal( lunerr, msg )
1838 ! checks involving nspec_amode and nspec_amode_nontracer
1839 ! being the same for all sections are no longer needed
1841 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1842 if ((l .gt. 0) .and. (l .le. ltot2)) ndum = ndum + 1
1844 if ((ndum .ne. 0) .and. (ndum .ne. nsize_aer(itype))) then
1845 msg = '*** subr move_sections error - ' // &
1846 'waterptr_aer must be on/off for all modes'
1847 call peg_message( lunerr, msg )
1848 write(msg,9030) 'iphase, itype =', iphase, itype
1849 call peg_message( lunerr, msg )
1850 call peg_error_fatal( lunerr, msg )
1852 9030 format( a, 2(1x,i6) )
1858 end subroutine move_sections_checkptrs
1862 end module module_mosaic_movesect