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 information 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
148 ! run diagnostic tests
149 ! (these will only be run for certain values of idiag_movesect
150 ! rce 2004-nov-29 - to avoid recursion, test_move_sections
151 ! is now called from mosaic_driver
153 ! call test_move_sections( iflag, iclm, jclm, k, m )
156 ! get "method_movesect" from digits 1-2 of msectional (treat 1-9 as 10)
157 method_movesect = mod( msectional, 100 )
158 if (method_movesect .le. 10) method_movesect = 10
160 ! get "idiag_movesect" from digits 3-4 of msectional
161 idiag_movesect = mod( msectional, 10000 )/100
163 if ((method_movesect .eq. 10) .or. &
164 (method_movesect .eq. 11) .or. &
165 (method_movesect .eq. 20) .or. &
166 (method_movesect .eq. 21) .or. &
167 (method_movesect .eq. 30) .or. &
168 (method_movesect .eq. 31)) then
170 else if ((method_movesect .eq. 19) .or. &
171 (method_movesect .eq. 29) .or. &
172 (method_movesect .eq. 39)) then
175 msg = '*** subr move_sections error - ' // &
176 'illegal value for msectional'
177 call peg_error_fatal( lunerr, msg )
180 if (iabs(iflag) .eq. 1) then
182 else if (iabs(iflag) .eq. 2) then
184 if (nphase_aer .lt. 2) then
185 msg = '*** subr move_sections error - ' // &
186 'iflag=2 (after aqueous chemistry) but nphase_aer < 2'
187 call peg_error_fatal( lunerr, msg )
188 else if (cw_phase .ne. 2) then
189 msg = '*** subr move_sections error - ' // &
190 'iflag=2 (after aqueous chemistry) but cw_phase .ne. 2'
191 call peg_error_fatal( lunerr, msg )
194 msg = '*** subr move_sections error - ' // &
195 'iabs(iflag) must be 1 or 2'
196 call peg_error_fatal( lunerr, msg )
200 ! when iflag=-1/-2, call move_sections_checkptrs then return
201 ! if ((ncorecnt .le. 0) .and. (k .le. 1)) then
202 if (iflag .le. 0) then
203 write(msg,9040) 'method', method_movesect
204 call peg_message( lunout, msg )
205 write(msg,9040) 'idiag ', idiag_movesect
206 call peg_message( lunout, msg )
207 call move_sections_checkptrs( iflag, iclm, jclm, k, m )
210 9040 format( '*** subr move_sections - ', a, ' =', i6 )
213 if (idiag_movesect .eq. 70) then
215 call peg_message( lunout, msg )
216 write(msg,9060) iclm, jclm, k, m, msectional
217 call peg_message( lunout, msg )
219 9060 format( '*** subr move_sections diagnostics i,j,k,m,msect =', 4i4, i6 )
225 ! if bin mass mixrat < smallmassaa (1.e-20 g/mol), then assume no growth
226 ! AND no water AND conform number so that size is within bin limits
227 smallmassaa = 1.0e-20
228 ! if bin mass OR volume mixrat < smallmassbb (1.e-30 g/mol), then
229 ! assume default density to avoid divide by zero
230 smallmassbb = 1.0e-30
233 ! process each type, one at a time
234 do 1900 itype = 1, ntype_aer
236 if (nsize_aer(itype) .le. 0) goto 1900
238 call move_sections_initial_conform( &
239 iflag, iclm, jclm, k, m, iphase, itype, &
240 method_movesect, idiag_movesect, llhysw, llwater, &
241 densdefault, densh2o, smallmassaa, smallmassbb, &
242 delta_water_conform1, delta_numb_conform1, &
243 drydenspp, drymasspp, rnumbpp, &
244 drydensxx0, drymassxx0, rnumbxx0, &
245 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
246 wetmassxx, wetvolxx, &
247 specmwxx, specdensxx )
249 if (method_movesect .le. 19) then
250 call move_sections_calc_movingcenter( &
251 iflag, iclm, jclm, k, m, iphase, itype, &
252 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
253 densdefault, densh2o, smallmassaa, smallmassbb, &
254 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
255 wetmassxx, wetvolxx, &
256 xferfracvol, xferfracnum )
258 call move_sections_calc_masnumadvect( &
259 iflag, iclm, jclm, k, m, iphase, itype, &
260 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
261 densdefault, densh2o, smallmassaa, smallmassbb, &
262 drydenspp, drymasspp, rnumbpp, &
263 drydensxx0, drymassxx0, rnumbxx0, &
264 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
265 wetmassxx, wetvolxx, &
266 xferfracvol, xferfracnum )
269 call move_sections_apply_moves( &
270 iflag, iclm, jclm, k, m, iphase, itype, &
271 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
272 densdefault, densh2o, smallmassaa, smallmassbb, &
273 delta_water_conform1, delta_numb_conform1, &
274 drydenspp, drymasspp, rnumbpp, &
275 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
276 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
277 xferfracvol, xferfracnum )
279 ! rce 08-nov-2004 - this call is new (and perhaps temporary)
280 ! rce 05-jul-2005 - do n1_inflow for aerchem growth but not for cldchem
281 if ((apply_n1_inflow .eq. 1) .and. &
282 (iphase .eq. ai_phase)) then
283 call move_sections_apply_n1_inflow( &
284 iflag, iclm, jclm, k, m, iphase, itype, &
285 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
286 densdefault, densh2o, smallmassaa, smallmassbb, &
287 delta_water_conform1, delta_numb_conform1, &
288 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
289 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
290 xferfracvol, xferfracnum, &
291 specmwxx, specdensxx )
294 ! call move_sections_final_conform( &
295 ! iflag, iclm, jclm, k, m, iphase, itype )
300 end subroutine move_sections
303 !-----------------------------------------------------------------------
304 subroutine move_sections_initial_conform( &
305 iflag, iclm, jclm, k, m, iphase, itype, &
306 method_movesect, idiag_movesect, llhysw, llwater, &
307 densdefault, densh2o, smallmassaa, smallmassbb, &
308 delta_water_conform1, delta_numb_conform1, &
309 drydenspp, drymasspp, rnumbpp, &
310 drydensxx0, drymassxx0, rnumbxx0, &
311 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
312 wetmassxx, wetvolxx, &
313 specmwxx, specdensxx )
316 ! routine does some initial tasks for the section movement
317 ! load rmassxx & rnumbxx from rsub
318 ! load specmwxx & specdensxx
319 ! set drymassxx & dryvolxx from drymass_aftgrow & drydens_aftgrow,
320 ! OR compute them from rmassxx, specmwxx, specdensxx if need be
321 ! set wetmassxx & wetvolxx from dry values & water mass
322 ! conform rnumbxx so that the mean particle size of each section
323 ! (= dryvolxx/rnumbxx) is within the section limits
334 integer iflag, iclm, jclm, iphase, itype, k, &
335 m, method_movesect, idiag_movesect, llhysw, llwater
336 real densdefault, densh2o, smallmassaa, smallmassbb
337 real delta_water_conform1, delta_numb_conform1
338 real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
339 drydensxx(maxd_asize)
340 real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
341 drymassxx(maxd_asize)
342 real dryvolxx(maxd_asize)
343 real rmassxx(maxd_acomp+2,maxd_asize)
344 real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
346 real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
347 real wetvolxx(maxd_asize)
348 real wetmassxx(maxd_asize)
352 integer l, ll, lnew, lold, l3, n, nnew, nold
354 real dummass, dumnum, dumnum_at_dhi, dumnum_at_dlo, dumr, &
355 dumvol, dumvol1p, dumwatrmass
358 ! assure positive definite
360 rsub(l,k,m) = max( 0., rsub(l,k,m) )
363 ! load mixrats into working arrays and assure positive definite
364 llhysw = ncomp_plustracer_aer(itype) + 1
365 llwater = ncomp_plustracer_aer(itype) + 2
366 do n = 1, nsize_aer(itype)
367 do ll = 1, ncomp_plustracer_aer(itype)
368 l = massptr_aer(ll,n,itype,iphase)
369 rmassxx(ll,n) = rsub(l,k,m)
371 rmassxx(llhysw,n) = 0.
373 if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
374 if (l .gt. 0) rmassxx(llhysw,n) = rsub(l,k,m)
375 rmassxx(llwater,n) = 0.
377 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
378 if (l .gt. 0) rmassxx(llwater,n) = rsub(l,k,m)
380 rnumbxx(n) = rsub(numptr_aer(n,itype,iphase),k,m)
381 rnumbxx0(n) = rnumbxx(n)
382 rnumbpp(n) = rnumbxx(n)
384 drydenspp(n) = drydens_pregrow(n,itype)
385 drymasspp(n) = drymass_pregrow(n,itype)
387 drydensxx(n) = drydens_aftgrow(n,itype)
388 drymassxx(n) = drymass_aftgrow(n,itype)
389 drydensxx0(n) = drydensxx(n)
390 drymassxx0(n) = drymassxx(n)
393 ! load specmw and specdens also
394 do ll = 1, ncomp_plustracer_aer(itype)
395 specmwxx(ll) = mw_aer(ll,itype)
396 specdensxx(ll) = dens_aer(ll,itype)
399 delta_water_conform1 = 0.0
400 delta_numb_conform1 = 0.0
403 do 1390 n = 1, nsize_aer(itype)
406 ! if drydens_aftgrow < 0.1, then bin had state="no_aerosol"
407 ! compute volume using default dry-densities, set water=0,
408 ! and conform the number
409 ! also do this if mass is extremely small (below smallmassaa)
410 ! OR if drydens_aftgrow > 20 (which is unrealistic)
412 if ( (drydensxx(n) .lt. 0.1) .or. &
413 (drydensxx(n) .gt. 20.0) .or. &
414 (drymassxx(n) .le. smallmassaa) ) then
417 do ll = 1, ncomp_aer(itype)
418 dumr = rmassxx(ll,n)*specmwxx(ll)
419 dummass = dummass + dumr
420 dumvol = dumvol + dumr/specdensxx(ll)
422 drymassxx(n) = dummass
423 if (min(dummass,dumvol) .le. smallmassbb) then
424 drydensxx(n) = densdefault
425 dumvol = dummass/densdefault
426 dumnum = dummass/(volumcen_sect(n,itype)*densdefault)
428 drydensxx(n) = dummass/dumvol
430 dumnum_at_dhi = dumvol/volumhi_sect(n,itype)
431 dumnum_at_dlo = dumvol/volumlo_sect(n,itype)
432 dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
434 delta_numb_conform1 = delta_numb_conform1 + dumnum - rnumbxx(n)
436 rnumbpp(n) = rnumbxx(n)
437 delta_water_conform1 = delta_water_conform1 - rmassxx(llwater,n)
438 rmassxx(llwater,n) = 0.
441 ! load dry/wet mass and volume into "xx" arrays
442 ! which hold values before inter-mode transferring
443 dryvolxx(n) = drymassxx(n)/drydensxx(n)
444 dumwatrmass = rmassxx(llwater,n)*mw_water_aer
445 wetmassxx(n) = drymassxx(n) + dumwatrmass
446 wetvolxx(n) = dryvolxx(n) + dumwatrmass/densh2o
451 end subroutine move_sections_initial_conform
454 !-----------------------------------------------------------------------
455 subroutine move_sections_calc_movingcenter( &
456 iflag, iclm, jclm, k, m, iphase, itype, &
457 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
458 densdefault, densh2o, smallmassaa, smallmassbb, &
459 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
460 wetmassxx, wetvolxx, &
461 xferfracvol, xferfracnum )
463 ! routine calculates section movements for the moving-center approach
465 ! material in section n will be transfered to section nnewsave(1,n)
467 ! the nnewsave are calculated here
468 ! the actual transfer is done in another routine
479 integer iflag, iclm, jclm, iphase, itype, k, &
480 m, method_movesect, idiag_movesect, llhysw, llwater
481 integer nnewsave(2,maxd_asize)
482 real densdefault, densh2o, smallmassaa, smallmassbb
483 real drydensxx(maxd_asize)
484 real drymassxx(maxd_asize)
485 real dryvolxx(maxd_asize)
486 real rmassxx(maxd_acomp+2,maxd_asize)
487 real rnumbxx(maxd_asize)
488 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
489 real wetmassxx(maxd_asize)
490 real wetvolxx(maxd_asize)
493 integer ll, n, ndum, nnew, nold
494 real dumnum, dumvol, dumvol1p, sixoverpi, third
502 ! compute mean size after growth (and corresponding section)
503 ! particles in section n will be transferred to section nnewsave(1,n)
505 do 1390 n = 1, nsize_aer(itype)
509 ! don't bother to transfer bins whose mass is extremely small
510 if (drymassxx(n) .le. smallmassaa) goto 1290
515 ! check for number so small that particle volume is
516 ! above that of largest section
517 if (dumnum .le. dumvol/volumhi_sect(nsize_aer(itype),itype)) then
518 nnew = nsize_aer(itype)
520 ! or below that of smallest section
521 else if (dumnum .ge. dumvol/volumlo_sect(1,itype)) then
526 ! dumvol1p is mean particle volume (cm3) for the section
527 dumvol1p = dumvol/dumnum
528 if (dumvol1p .gt. volumhi_sect(n,itype)) then
529 do while ( (nnew .lt. nsize_aer(itype)) .and. &
530 (dumvol1p .gt. volumhi_sect(nnew,itype)) )
534 else if (dumvol1p .lt. volumlo_sect(n,itype)) then
535 do while ( (nnew .gt. 1) .and. &
536 (dumvol1p .lt. volumlo_sect(nnew,itype)) )
542 1290 nnewsave(1,n) = nnew
545 xferfracvol(1,n) = 1.0
546 xferfracvol(2,n) = 0.0
547 xferfracnum(1,n) = 1.0
548 xferfracnum(2,n) = 0.0
554 ! output nnewsave values when msectional=7xxx
555 if (idiag_movesect .eq. 70) then
557 do n = 1, nsize_aer(itype)
558 if (nnewsave(1,n) .ne. n) ndum = ndum + 1
560 if (ndum .gt. 0) then
561 write(msg,97751) 'YES', iclm, jclm, k, m, &
562 ndum, (nnewsave(1,n), n=1,nsize_aer(itype))
563 call peg_message( lunout, msg )
565 write(msg,97751) 'NO ', iclm, jclm, k, m, &
566 ndum, (nnewsave(1,n), n=1,nsize_aer(itype))
567 call peg_message( lunout, msg )
570 97751 format( 'movesect', a, 4i3, 3x, i3, 3x, 14i3 )
573 end subroutine move_sections_calc_movingcenter
576 !-----------------------------------------------------------------------
577 subroutine move_sections_calc_masnumadvect( &
578 iflag, iclm, jclm, k, m, iphase, itype, &
579 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
580 densdefault, densh2o, smallmassaa, smallmassbb, &
581 drydenspp, drymasspp, rnumbpp, &
582 drydensxx0, drymassxx0, rnumbxx0, &
583 drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
584 wetmassxx, wetvolxx, &
585 xferfracvol, xferfracnum )
587 ! routine calculates section movements for the mass-number-advection approach
589 ! material in section n will be transfered to sections
590 ! nnewsave(1,n) and nnewsave(2,n)
591 ! the fractions of mass/volume transfered to each are
592 ! xferfracvol(1,n) and xferfracvol(2,n)
593 ! the fractions of number transfered to each are
594 ! xferfracnum(1,n) and xferfracnum(2,n)
596 ! the nnewsave, xferfracvol, and xferfracnum are calculated here
597 ! the actual transfer is done in another routine
608 integer iflag, iclm, jclm, iphase, itype, k, &
609 m, method_movesect, idiag_movesect, llhysw, llwater
610 integer nnewsave(2,maxd_asize)
612 real densdefault, densh2o, smallmassaa, smallmassbb
613 real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
614 drydensxx(maxd_asize)
615 real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
616 drymassxx(maxd_asize)
617 real dryvolxx(maxd_asize)
618 real rmassxx(maxd_acomp+2,maxd_asize)
619 real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
621 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
622 real wetvolxx(maxd_asize)
623 real wetmassxx(maxd_asize)
626 integer ierr, n, nnew, nnew2
627 integer iforce_movecenter(maxd_asize)
629 real dum1, dum2, dum3
630 real dumaa, dumbb, dumgamma, dumratio
631 real dumfracnum, dumfracvol
634 real dumvbar_aft, dumvbar_pre
635 real dumvcutlo_nnew_pre, dumvcuthi_nnew_pre
636 real dumvlo_pre, dumvhi_pre, dumvdel_pre
637 real dumvtot_aft, dumvtot_pre
639 real sixoverpi, third
650 ! compute mean size after growth (and corresponding section)
651 ! some of the particles in section n will be transferred to section nnewsave(1,n)
653 ! if the aftgrow mass is extremely small,
654 ! OR if the aftgrow mean size is outside of
655 ! [dlo_sect(1,itype), dhi_sect(nsize_aer(itype),itype)]
656 ! then use the moving-center method_movesect for this bin
657 ! (don't try to compute the pregrow within-bin distribution)
659 do 3900 n = 1, nsize_aer(itype)
662 iforce_movecenter(n) = 0
664 xferfracvol(1,n) = 1.0
665 xferfracvol(2,n) = 0.0
666 xferfracnum(1,n) = 1.0
667 xferfracnum(2,n) = 0.0
677 dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
678 dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
682 ! don't bother to transfer bins whose mass is extremely small
683 if (drymassxx(n) .le. smallmassaa) then
684 iforce_movecenter(n) = 1
688 dumvtot_aft = dryvolxx(n)
691 ! check for particle volume above that of largest section
692 ! or below that of smallest section
693 if (dumntot .le. dumvtot_aft/volumhi_sect(nsize_aer(itype),itype)) then
694 nnew = nsize_aer(itype)
695 iforce_movecenter(n) = 2
697 else if (dumntot .ge. dumvtot_aft/volumlo_sect(1,itype)) then
699 iforce_movecenter(n) = 3
703 ! dumvbar_aft is mean particle volume (cm3) for the section
704 ! find the section that encloses this volume
705 dumvbar_aft = dumvtot_aft/dumntot
706 if (dumvbar_aft .gt. volumhi_sect(n,itype)) then
707 do while ( (nnew .lt. nsize_aer(itype)) .and. &
708 (dumvbar_aft .gt. volumhi_sect(nnew,itype)) )
712 else if (dumvbar_aft .lt. volumlo_sect(n,itype)) then
713 do while ( (nnew .gt. 1) .and. &
714 (dumvbar_aft .lt. volumlo_sect(nnew,itype)) )
720 1290 nnewsave(1,n) = nnew
723 if (iforce_movecenter(n) .gt. 0) goto 3700
726 ! if drydenspp (pregrow) < 0.1 (because bin had state="no_aerosol" before
727 ! growth was computed, so its initial mass was very small)
728 ! then use the moving-center method_movesect for this bin
729 ! (don't try to compute the pregrow within-bin distribution)
730 ! also do this if pregrow mass is extremely small (below smallmassaa)
731 ! OR if drydenspp > 20 (unphysical)
732 if ( (drydenspp(n) .lt. 0.1) .or. &
733 (drydenspp(n) .gt. 20.0) .or. &
734 (drymasspp(n) .le. smallmassaa) ) then
735 iforce_movecenter(n) = 11
739 dumvtot_pre = drymasspp(n)/drydenspp(n)
741 dumvlo_pre = volumlo_sect(n,itype)
742 dumvhi_pre = volumhi_sect(n,itype)
743 dumvdel_pre = dumvhi_pre - dumvlo_pre
745 ! if the pregrow mean size is outside of OR very close to the bin limits,
746 ! then use moving-center approach for this bin
747 dumv = dumvhi_pre - 0.01*dumvdel_pre
748 if (dumntot .le. dumvtot_pre/dumv) then
749 iforce_movecenter(n) = 12
752 dumv = dumvlo_pre + 0.01*dumvdel_pre
753 if (dumntot .ge. dumvtot_pre/dumv) then
754 iforce_movecenter(n) = 13
758 ! calculate the pregrow within-section size distribution
759 dumvbar_pre = dumvtot_pre/dumntot
760 dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
761 dumratio = dumvbar_pre/dumvlo_pre
763 if (dumratio .le. (1.0001 + dumgamma/3.0)) then
764 dumv = dumvlo_pre + 3.0*(dumvbar_pre-dumvlo_pre)
765 dumvhi_pre = min( dumvhi_pre, dumv )
766 dumvdel_pre = dumvhi_pre - dumvlo_pre
767 dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
768 dumratio = dumvbar_pre/dumvlo_pre
769 else if (dumratio .ge. (0.9999 + dumgamma*2.0/3.0)) then
770 dumv = dumvhi_pre + 3.0*(dumvbar_pre-dumvhi_pre)
771 dumvlo_pre = max( dumvlo_pre, dumv )
772 dumvdel_pre = dumvhi_pre - dumvlo_pre
773 dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
774 dumratio = dumvbar_pre/dumvlo_pre
777 dumbb = (dumratio - 1.0 - 0.5*dumgamma)*12.0/dumgamma
778 dumaa = 1.0 - 0.5*dumbb
780 ! calculate pregrow volumes corresponding to the nnew
782 dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
783 dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
785 ! if the [dumvlo_pre, dumvhi_pre] falls completely within
786 ! the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval,
787 ! then all mass and number go to nnew
788 if (nnew .eq. 1) then
789 if (dumvhi_pre .le. dumvcuthi_nnew_pre) then
790 iforce_movecenter(n) = 21
794 else if (nnew .eq. nsize_aer(itype)) then
795 if (dumvlo_pre .ge. dumvcutlo_nnew_pre) then
796 iforce_movecenter(n) = 22
801 if ((dumvlo_pre .ge. dumvcutlo_nnew_pre) .and. &
802 (dumvhi_pre .le. dumvcuthi_nnew_pre)) then
803 iforce_movecenter(n) = 23
804 else if (dumvlo_pre .lt. dumvcutlo_nnew_pre) then
810 if (iforce_movecenter(n) .gt. 0) goto 3700
812 ! calculate the fraction of ntot and vtot that are within
813 ! the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval
814 dumzlo = (dumvcutlo_nnew_pre - dumvlo_pre)/dumvdel_pre
815 dumzhi = (dumvcuthi_nnew_pre - dumvlo_pre)/dumvdel_pre
816 dumzlo = max( dumzlo, 0.0 )
817 dumzhi = min( dumzhi, 1.0 )
818 dum1 = dumzhi - dumzlo
819 dum2 = (dumzhi**2 - dumzlo**2)*0.5
820 dum3 = (dumzhi**3 - dumzlo**3)/3.0
821 dumfracnum = dumaa*dum1 + dumbb*dum2
822 dumfracvol = (dumvlo_pre/dumvbar_pre) * (dumaa*dum1 + &
823 (dumaa*dumgamma + dumbb)*dum2 + (dumbb*dumgamma)*dum3)
825 if ((dumfracnum .le. 0.0) .or. (dumfracvol .le. 0.0)) then
826 iforce_movecenter(n) = 31
827 nnewsave(1,n) = nnew2
828 else if ((dumfracnum .ge. 1.0) .or. (dumfracvol .ge. 1.0)) then
829 iforce_movecenter(n) = 32
831 if (iforce_movecenter(n) .gt. 0) goto 3700
833 nnewsave(2,n) = nnew2
835 xferfracvol(1,n) = dumfracvol
836 xferfracvol(2,n) = 1.0 - dumfracvol
837 xferfracnum(1,n) = dumfracnum
838 xferfracnum(2,n) = 1.0 - dumfracnum
842 ! output nnewsave values when msectional=7xxx
843 if (idiag_movesect .ne. 70) goto 3800
845 if (nnewsave(2,n) .eq. 0) then
846 if (nnewsave(1,n) .eq. 0) then
848 else if (nnewsave(1,n) .eq. n) then
853 else if (nnewsave(1,n) .eq. 0) then
854 if (nnewsave(2,n) .eq. n) then
859 else if (nnewsave(2,n) .eq. n) then
860 if (nnewsave(1,n) .eq. n) then
865 else if (nnewsave(1,n) .eq. n) then
872 if (drymasspp(n) .gt. drymassxx(n)) dumch1 = '-'
875 call peg_message( lunout, msg )
876 write(msg,97010) dumch1, dumch4, iclm, jclm, k, m, &
877 n, nnewsave(1,n), nnewsave(2,n), iforce_movecenter(n)
878 call peg_message( lunout, msg )
879 write(msg,97020) 'pre mass, dens ', &
880 drymasspp(n), drydenspp(n)
881 call peg_message( lunout, msg )
882 write(msg,97020) 'aft mass, dens, numb', &
883 drymassxx(n), drydensxx(n), rnumbxx(n)
884 call peg_message( lunout, msg )
885 if ((drydensxx(n) .ne. drydensxx0(n)) .or. &
886 (drymassxx(n) .ne. drymassxx0(n)) .or. &
887 (rnumbxx(n) .ne. rnumbxx0(n) )) then
888 write(msg,97020) 'aft0 mas, dens, numb', &
889 drymassxx0(n), drydensxx0(n), rnumbxx0(n)
890 call peg_message( lunout, msg )
892 write(msg,97020) 'vlop0, vbarp, vhip0', &
893 volumlo_sect(n,itype), dumvbar_pre, volumhi_sect(n,itype)
894 call peg_message( lunout, msg )
895 write(msg,97020) 'vlop , vbarp, vhip ', &
896 dumvlo_pre, dumvbar_pre, dumvhi_pre
897 call peg_message( lunout, msg )
898 write(msg,97020) 'vloax, vbarax, vhiax', &
899 dumvcutlo_nnew_pre, dumvbar_pre, dumvcuthi_nnew_pre
900 call peg_message( lunout, msg )
901 write(msg,97020) 'vloa0, vbara, vhia0', &
902 volumlo_sect(nnew,itype), dumvbar_aft, volumhi_sect(nnew,itype)
903 call peg_message( lunout, msg )
904 write(msg,97020) 'dumfrvol, num, ratio', &
905 dumfracvol, dumfracnum, dumratio
906 call peg_message( lunout, msg )
907 write(msg,97020) 'frvol,num1; vol,num2', &
908 xferfracvol(1,n), xferfracnum(1,n), &
909 xferfracvol(2,n), xferfracnum(2,n)
910 call peg_message( lunout, msg )
912 97010 format( 'movesect', 2a, 7x, 4i3, 4x, &
913 'n,nnews', 3i3, 4x, 'iforce', i3.2 )
914 97020 format( a, 1p, 4e13.4 )
919 ! check for legal combinations of nnewsave(1,n) & nnewsave(2,n)
922 ! both are non-zero AND iabs(nnew1-nnew2) != 1
924 if (nnewsave(1,n) .eq. nnewsave(2,n)) then
926 else if (nnewsave(1,n)*nnewsave(2,n) .ne. 0) then
927 if (iabs(nnewsave(1,n)-nnewsave(2,n)) .ne. 1) ierr = 1
929 if (ierr .gt. 0) then
930 write(msg,97010) 'E', 'RROR', iclm, jclm, k, m, &
931 n, nnewsave(1,n), nnewsave(2,n), iforce_movecenter(n)
932 call peg_message( lunout, msg )
936 ! if method_movesect == 30-31 then force moving center
937 ! this is just for testing purposes
938 if ((method_movesect .ge. 30) .and. (method_movesect .le. 39)) then
941 xferfracvol(1,n) = 1.0
942 xferfracvol(2,n) = 0.0
943 xferfracnum(1,n) = 1.0
944 xferfracnum(2,n) = 0.0
950 end subroutine move_sections_calc_masnumadvect
953 !-----------------------------------------------------------------------
954 subroutine move_sections_apply_moves( &
955 iflag, iclm, jclm, k, m, iphase, itype, &
956 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
957 densdefault, densh2o, smallmassaa, smallmassbb, &
958 delta_water_conform1, delta_numb_conform1, &
959 drydenspp, drymasspp, rnumbpp, &
960 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
961 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
962 xferfracvol, xferfracnum )
964 ! routine performs the actual transfer of aerosol number and mass
976 integer iflag, iclm, jclm, iphase, itype, k, &
977 m, method_movesect, idiag_movesect, llhysw, llwater
978 integer nnewsave(2,maxd_asize)
980 real densdefault, densh2o, smallmassaa, smallmassbb
981 real delta_water_conform1, delta_numb_conform1
982 real drydenspp(maxd_asize)
983 real drymasspp(maxd_asize)
984 real drymassxx(maxd_asize), drymassyy(maxd_asize)
985 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
986 real rmassxx(maxd_acomp+2,maxd_asize), &
987 rmassyy(maxd_acomp+2,maxd_asize)
988 real rnumbpp(maxd_asize)
989 real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
990 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
991 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
992 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
995 integer jj, l, ll, n, ndum, nnew, nold
996 integer jja, jjb, jjc
998 real delta_numb_conform2, &
999 dumbot, dumnum, dumnum_at_dhi, dumnum_at_dlo, &
1000 dumvol, dumvol1p, dumxfvol, dumxfnum, sixoverpi, third
1001 real dumpp(maxd_asize), dumxx(maxd_asize), dumyy(maxd_asize), &
1013 ! initialize "yy" arrays that hold values after inter-mode transferring
1014 ! "yy" = "xx" for sections that do not move at all
1015 ! "yy" = 0.0 for sections that do move (partially or completely)
1017 do 1900 n = 1, nsize_aer(itype)
1019 if ( (nnewsave(1,n) .eq. n) .and. &
1020 (nnewsave(2,n) .eq. 0) ) then
1021 ! if nnew == n, then material in section n will not be transferred, and
1022 ! section n will contain its initial material plus any material
1023 ! transferred from other sections
1024 ! so initialize "yy" arrays with "xx" values
1025 drymassyy(n) = drymassxx(n)
1026 dryvolyy(n) = dryvolxx(n)
1027 wetmassyy(n) = wetmassxx(n)
1028 wetvolyy(n) = wetvolxx(n)
1029 rnumbyy(n) = rnumbxx(n)
1030 do ll = 1, ncomp_plustracer_aer(itype) + 2
1031 rmassyy(ll,n) = rmassxx(ll,n)
1035 ! if nnew .ne. n, then material in section n will be transferred, and
1036 ! section n will only contain material that is transferred from
1038 ! so initialize "yy" arrays to zero
1044 do ll = 1, ncomp_plustracer_aer(itype) + 2
1053 ! do the transfer of mass and number
1055 do 2900 nold = 1, nsize_aer(itype)
1057 if ( (nnewsave(1,nold) .eq. nold) .and. &
1058 (nnewsave(2,nold) .eq. 0 ) ) goto 2900
1062 nnew = nnewsave(jj,nold)
1063 if (nnew .le. 0) goto 2800
1065 dumxfvol = xferfracvol(jj,nold)
1066 dumxfnum = xferfracnum(jj,nold)
1068 do ll = 1, ncomp_plustracer_aer(itype) + 2
1069 rmassyy(ll,nnew) = rmassyy(ll,nnew) + rmassxx(ll,nold)*dumxfvol
1071 rnumbyy(nnew) = rnumbyy(nnew) + rnumbxx(nold)*dumxfnum
1073 drymassyy(nnew) = drymassyy(nnew) + drymassxx(nold)*dumxfvol
1074 dryvolyy(nnew) = dryvolyy(nnew) + dryvolxx(nold)*dumxfvol
1075 wetmassyy(nnew) = wetmassyy(nnew) + wetmassxx(nold)*dumxfvol
1076 wetvolyy(nnew) = wetvolyy(nnew) + wetvolxx(nold)*dumxfvol
1083 ! transfer among sections is completed
1084 ! - check for conservation of mass/volume/number
1085 ! - conform number again
1086 ! - compute/store densities and mean sizes
1087 ! - if k=1, save values for use by dry deposition routine
1088 ! - copy new mixrats back to rsub array
1090 call move_sections_conserve_check( &
1091 1, iflag, iclm, jclm, k, m, iphase, itype, &
1092 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1093 densdefault, densh2o, smallmassaa, smallmassbb, &
1094 delta_water_conform1, delta_numb_conform1, &
1095 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1096 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1098 delta_numb_conform2 = 0.0
1100 do 3900 n = 1, nsize_aer(itype)
1102 dumvol = dryvolyy(n)
1103 if (min(drymassyy(n),dumvol) .le. smallmassbb) then
1104 dumvol = drymassyy(n)/densdefault
1105 dumnum = drymassyy(n)/(volumlo_sect(n,itype)*densdefault)
1106 delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n)
1108 adrydens_sub(n,itype,k,m) = densdefault
1109 awetdens_sub(n,itype,k,m) = densdefault
1110 admeandry_sub(n,itype,k,m) = dcen_sect(n,itype)
1111 admeanwet_sub(n,itype,k,m) = dcen_sect(n,itype)
1114 dumnum_at_dhi = dumvol/volumhi_sect(n,itype)
1115 dumnum_at_dlo = dumvol/volumlo_sect(n,itype)
1116 dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
1117 delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n)
1119 adrydens_sub(n,itype,k,m) = drymassyy(n)/dumvol
1120 dumvol1p = dumvol/dumnum
1121 admeandry_sub(n,itype,k,m) = (dumvol1p*sixoverpi)**third
1122 awetdens_sub(n,itype,k,m) = wetmassyy(n)/wetvolyy(n)
1123 dumvol1p = wetvolyy(n)/dumnum
1124 admeanwet_sub(n,itype,k,m) = min( 100.*dcen_sect(n,itype), &
1125 (dumvol1p*sixoverpi)**third )
1127 aqvoldry_sub(n,itype,k,m) = dumvol
1128 aqmassdry_sub(n,itype,k,m) = drymassyy(n)
1130 ! if (k .eq. 1) then
1131 ! awetdens_sfc(n,itype,iclm,jclm) = awetdens_sub(n,itype,k,m)
1132 ! admeanwet_sfc(n,itype,iclm,jclm) = admeanwet_sub(n,itype,k,m)
1135 do ll = 1, ncomp_plustracer_aer(itype)
1136 l = massptr_aer(ll,n,itype,iphase)
1137 rsub(l,k,m) = rmassyy(ll,n)
1140 if (iphase .eq. ai_phase) then
1141 l = waterptr_aer(n,itype)
1142 if (l .gt. 0) rsub(l,k,m) = rmassyy(llwater,n)
1143 l = hyswptr_aer(n,itype)
1144 if (l .gt. 0) rsub(l,k,m) = rmassyy(llhysw,n)
1146 rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
1150 delta_numb_conform1 = delta_numb_conform1 + delta_numb_conform2
1152 call move_sections_conserve_check( &
1153 2, iflag, iclm, jclm, k, m, iphase, itype, &
1154 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1155 densdefault, densh2o, smallmassaa, smallmassbb, &
1156 delta_water_conform1, delta_numb_conform1, &
1157 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1158 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1162 ! output nnewsave values when msectional=7xxx
1163 if (idiag_movesect .ne. 70) goto 4900
1166 do n = 1, nsize_aer(itype)
1167 if (nnewsave(1,n)+nnewsave(2,n) .ne. n) ndum = ndum + 1
1169 ! if (ndum .gt. 0) then
1170 ! write(msg,97010) 'SOME', iclm, jclm, k, m, &
1171 ! ndum, (nnewsave(1,n), nnewsave(2,n), n=1,nsize_aer(itype))
1172 ! call peg_message( lunout, msg )
1174 ! write(msg,97010) 'NONE', iclm, jclm, k, m, &
1175 ! ndum, (nnewsave(1,n), nnewsave(2,n), n=1,nsize_aer(itype))
1176 ! call peg_message( lunout, msg )
1180 if (ndum .gt. 0) dumch4 = 'SOME'
1182 call peg_message( lunout, msg )
1183 write(msg,97010) dumch4, iclm, jclm, k, m, ndum
1184 call peg_message( lunout, msg )
1185 do jjb = 1, nsize_aer(itype), 10
1186 jjc = min( jjb+9, nsize_aer(itype) )
1187 write(msg,97011) (nnewsave(1,n), nnewsave(2,n), n=jjb,jjc)
1188 call peg_message( lunout, msg )
1191 ! write(msg,97020) 'rnumbold', (rnumbxx(n), n=1,nsize_aer(itype))
1192 ! call peg_message( lunout, msg )
1193 ! write(msg,97020) 'rnumbnew', (rnumbyy(n), n=1,nsize_aer(itype))
1194 ! call peg_message( lunout, msg )
1196 ! write(msg,97020) 'drvolold', (dryvolxx(n), n=1,nsize_aer(itype))
1197 ! call peg_message( lunout, msg )
1198 ! write(msg,97020) 'drvolnew', (dryvolyy(n), n=1,nsize_aer(itype))
1199 ! call peg_message( lunout, msg )
1201 dumbot = log( volumhi_sect(1,itype)/volumlo_sect(1,itype) )
1202 do n = 1, nsize_aer(itype)
1206 if ( (drydenspp(n) .gt. 0.5) .and. &
1207 (drymasspp(n) .gt. smallmassaa) ) then
1208 dumvol = drymasspp(n)/drydenspp(n)
1209 if ((rnumbpp(n) .ge. 1.0e-35) .and. &
1210 (dumvol .ge. 1.0e-35)) then
1211 dumvol1p = dumvol/rnumbpp(n)
1212 dumpp(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
1215 if ((rnumbxx(n) .ge. 1.0e-35) .and. &
1216 (dryvolxx(n) .ge. 1.0e-35)) then
1217 dumvol1p = dryvolxx(n)/rnumbxx(n)
1218 dumxx(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
1220 if ((rnumbyy(n) .ge. 1.0e-35) .and. &
1221 (dryvolyy(n) .ge. 1.0e-35)) then
1222 dumvol1p = dryvolyy(n)/rnumbyy(n)
1223 dumyy(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
1227 ! write(msg,97030) 'lnvolold', (dumxx(n), n=1,nsize_aer(itype))
1228 ! call peg_message( lunout, msg )
1229 ! write(msg,97030) 'lnvolnew', (dumyy(n), n=1,nsize_aer(itype))
1230 ! call peg_message( lunout, msg )
1233 if (jja .eq. 1) then
1235 dumout(:) = rnumbxx(:)
1236 else if (jja .eq. 2) then
1238 dumout(:) = rnumbyy(:)
1239 else if (jja .eq. 3) then
1241 dumout(:) = dryvolxx(:)
1242 else if (jja .eq. 4) then
1244 dumout(:) = dryvolyy(:)
1245 else if (jja .eq. 5) then
1247 dumout(:) = dumxx(:)
1248 else if (jja .eq. 6) then
1250 dumout(:) = dumyy(:)
1251 else if (jja .eq. 7) then
1253 dumout(:) = dumpp(:)
1255 do jjb = 1, nsize_aer(itype), 10
1256 jjc = min( jjb+9, nsize_aer(itype) )
1257 if (jja .le. 4) then
1258 write(msg,97020) dumch8, (dumout(n), n=jjb,jjc)
1260 write(msg,97030) dumch8, (dumout(n), n=jjb,jjc)
1262 call peg_message( lunout, msg )
1267 !97010 format( / 'movesectapply', a, 4i3, 3x, i3 / 5x, 10(3x,2i3) )
1268 !97020 format( a, 1p, 10e9.1 / (( 8x, 1p, 10e9.1 )) )
1269 !97030 format( a, 10f9.3 / (( 8x, 10f9.3 )) )
1270 97010 format( 'movesectapply', a, 4i3, 3x, i3 )
1271 97011 format( 5x, 10(3x,2i3) )
1272 97020 format( a, 1p, 10e9.1 )
1273 97030 format( a, 10f9.3 )
1277 end subroutine move_sections_apply_moves
1280 !-----------------------------------------------------------------------
1281 ! rce 08-nov-2004 - this routine is new (and perhaps temporary)
1283 subroutine move_sections_apply_n1_inflow( &
1284 iflag, iclm, jclm, k, m, iphase, itype, &
1285 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1286 densdefault, densh2o, smallmassaa, smallmassbb, &
1287 delta_water_conform1, delta_numb_conform1, &
1288 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1289 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
1290 xferfracvol, xferfracnum, &
1291 specmwxx, specdensxx )
1293 ! routine applies an ad_hoc correction to bin 1 to account for
1294 ! growth of smaller particles (which are not simulated) growing into
1296 ! the correction to particle number balances the loss of particles from
1297 ! bin 1 to larger bins by growth
1298 ! the correction to mass assumes that the particles coming into bin 1
1299 ! are slightly larger than dlo_sect(1)
1306 ! include 'v33com9a'
1307 ! include 'v33com9b'
1310 integer iflag, iclm, jclm, iphase, itype, k, &
1311 m, method_movesect, idiag_movesect, llhysw, llwater
1312 integer nnewsave(2,maxd_asize)
1314 real densdefault, densh2o, smallmassaa, smallmassbb
1315 real delta_water_conform1, delta_numb_conform1
1316 real drymassxx(maxd_asize), drymassyy(maxd_asize)
1317 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
1318 real rmassxx(maxd_acomp+2,maxd_asize), &
1319 rmassyy(maxd_acomp+2,maxd_asize)
1320 real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
1321 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
1322 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
1323 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
1324 real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
1328 integer jj, l, ll, n, nnew, nold
1330 real deltanum, deltavol, dumvol1p, dumxfvol, dumxfnum
1334 ! compute fraction of number transferred out of bin 1 by growth
1338 if ( (nnewsave(1,nold) .eq. nold) .and. &
1339 (nnewsave(2,nold) .eq. 0 ) ) goto 3900
1343 nnew = nnewsave(jj,nold)
1344 if (nnew .le. 0) goto 2800
1345 if (nnew .eq. nold) goto 2800
1346 dumxfnum = dumxfnum + xferfracnum(jj,nold)
1350 ! compute "inflow" change to number and volume
1351 ! number change matches that lost by growth
1352 ! volume change assume inflow particles slightly bigger then dlo_sect
1354 dumvol1p = 0.95*volumlo_sect(n,itype) + 0.05*volumhi_sect(n,itype)
1355 deltanum = dumxfnum*rnumbxx(n)
1356 deltavol = deltanum*dumvol1p
1357 if (dumxfnum .le. 0.0) goto 3900
1358 if (deltanum .le. 0.0) goto 3900
1359 if (deltavol .le. 0.0) goto 3900
1362 ! increment the number and masses
1363 ! if the old dryvol (dryvolxx) > smallmassbb, then compute mass increment for
1364 ! each species from the old masses (rmassxx)
1365 ! otherwise only increment the first mass species
1367 if (dryvolxx(n) .gt. smallmassbb) then
1368 dumxfvol = deltavol/dryvolxx(n)
1369 do ll = 1, ncomp_plustracer_aer(itype) + 2
1370 rmassyy(ll,n) = rmassyy(ll,n) + dumxfvol*rmassxx(ll,n)
1374 rmassyy(ll,n) = rmassyy(ll,n) + deltavol*specdensxx(ll)/specmwxx(ll)
1376 rnumbyy(n) = rnumbyy(n) + deltanum
1379 ! transfer results into rsub
1381 do ll = 1, ncomp_plustracer_aer(itype)
1382 l = massptr_aer(ll,n,itype,iphase)
1383 rsub(l,k,m) = rmassyy(ll,n)
1385 if (iphase .eq. ai_phase) then
1386 l = waterptr_aer(n,itype)
1387 if (l .gt. 0) rsub(l,k,m) = rmassyy(llwater,n)
1388 l = hyswptr_aer(n,itype)
1389 if (l .gt. 0) rsub(l,k,m) = rmassyy(llhysw,n)
1391 rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
1396 end subroutine move_sections_apply_n1_inflow
1399 !-----------------------------------------------------------------------
1400 subroutine move_sections_conserve_check( ipass, &
1401 iflag, iclm, jclm, k, m, iphase, itype, &
1402 method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
1403 densdefault, densh2o, smallmassaa, smallmassbb, &
1404 delta_water_conform1, delta_numb_conform1, &
1405 drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
1406 drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1408 ! routine checks for conservation of number, mass, and volume
1409 ! by the move_sections algorithm
1412 ! initialize all thesum(jj,ll) to zero
1413 ! computes thesum(1,ll) from rmassxx, rnumbxx, ...
1414 ! computes thesum(2,ll) from rmassyy, rnumbyy, ...
1415 ! compares thesum(1,ll) with thesum(2,ll)
1416 ! computes thesum(3,ll) from rsub before section movement
1418 ! computes thesum(4,ll) from rsub after section movement
1419 ! compares thesum(3,ll) with thesum(4,ll)
1421 ! currently only implemented for condensational growth (iflag=1)
1428 ! include 'v33com9a'
1429 ! include 'v33com9b'
1432 integer ipass, iflag, iclm, jclm, iphase, itype, k, &
1433 m, method_movesect, idiag_movesect, llhysw, llwater
1434 integer nnewsave(2,maxd_asize)
1435 real densdefault, densh2o, smallmassaa, smallmassbb
1436 real delta_water_conform1, delta_numb_conform1
1437 real drymassxx(maxd_asize), drymassyy(maxd_asize)
1438 real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
1439 real rmassxx(maxd_acomp+2,maxd_asize), &
1440 rmassyy(maxd_acomp+2,maxd_asize)
1441 real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
1442 real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
1443 real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
1446 integer jj, l, ll, llworst, llworstb, n
1447 integer nerr, nerrmax
1449 data nerr, nerrmax / 0, 999 /
1451 real dumbot, dumtop, dumtoler, dumerr, dumworst, dumworstb
1452 real duma, dumb, dumc, dume
1453 real thesum(4,maxd_acomp+7)
1456 character*8 dumname(maxd_acomp+7)
1460 if (ipass .eq. 2) goto 2000
1462 do ll = 1, ncomp_plustracer_aer(itype)+7
1468 do n = 1, nsize_aer(itype)
1469 do ll = 1, ncomp_plustracer_aer(itype)+2
1470 thesum(1,ll) = thesum(1,ll) + rmassxx(ll,n)
1471 thesum(2,ll) = thesum(2,ll) + rmassyy(ll,n)
1473 ll = ncomp_plustracer_aer(itype)+3
1474 thesum(1,ll) = thesum(1,ll) + rnumbxx(n)
1475 thesum(2,ll) = thesum(2,ll) + rnumbyy(n)
1476 ll = ncomp_plustracer_aer(itype)+4
1477 thesum(1,ll) = thesum(1,ll) + drymassxx(n)
1478 thesum(2,ll) = thesum(2,ll) + drymassyy(n)
1479 ll = ncomp_plustracer_aer(itype)+5
1480 thesum(1,ll) = thesum(1,ll) + dryvolxx(n)
1481 thesum(2,ll) = thesum(2,ll) + dryvolyy(n)
1482 ll = ncomp_plustracer_aer(itype)+6
1483 thesum(1,ll) = thesum(1,ll) + wetmassxx(n)
1484 thesum(2,ll) = thesum(2,ll) + wetmassyy(n)
1485 ll = ncomp_plustracer_aer(itype)+7
1486 thesum(1,ll) = thesum(1,ll) + wetvolxx(n)
1487 thesum(2,ll) = thesum(2,ll) + wetvolyy(n)
1493 ! calc sum over bins for each species
1494 ! for water, account for loss in initial conform (delta_water_conform1)
1496 do n = 1, nsize_aer(itype)
1497 do ll = 1, ncomp_plustracer_aer(itype)+3
1498 if (ll .le. ncomp_plustracer_aer(itype)) then
1499 l = massptr_aer(ll,n,itype,iphase)
1500 else if (ll .eq. ncomp_plustracer_aer(itype)+1) then
1502 if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
1503 else if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1505 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1507 l = numptr_aer(n,itype,iphase)
1510 thesum(ipass+2,ll) = thesum(ipass+2,ll) + rsub(l,k,m)
1513 if (ipass .eq. 2) then
1514 ll = ncomp_plustracer_aer(itype)+2
1515 thesum(3,ll) = thesum(3,ll) + delta_water_conform1
1516 ll = ncomp_plustracer_aer(itype)+3
1517 thesum(3,ll) = thesum(3,ll) + delta_numb_conform1
1521 ! now compare either sum1-sum2 or sum3-sum4
1522 ! on ipass=1, jj=1, so compare sum1 & sum2
1523 ! on ipass=2, jj=3, so compare sum3 & sum4
1525 do ll = 1, ncomp_plustracer_aer(itype)+7
1527 write(dumname(ll),'(i4.4)') ll
1528 if (ll .le. ncomp_plustracer_aer(itype)) dumname(ll) = &
1529 name(massptr_aer(ll,1,itype,iphase))(1:4)
1530 if (ll .eq. ncomp_plustracer_aer(itype)+1) dumname(ll) = 'hysw'
1531 if (ll .eq. ncomp_plustracer_aer(itype)+2) dumname(ll) = 'watr'
1532 if (ll .eq. ncomp_plustracer_aer(itype)+3) dumname(ll) = 'numb'
1533 if (ll .eq. ncomp_plustracer_aer(itype)+4) dumname(ll) = 'drymass'
1534 if (ll .eq. ncomp_plustracer_aer(itype)+5) dumname(ll) = 'dryvol'
1535 if (ll .eq. ncomp_plustracer_aer(itype)+6) dumname(ll) = 'wetmass'
1536 if (ll .eq. ncomp_plustracer_aer(itype)+7) dumname(ll) = 'wetvol'
1544 do ll = 1, ncomp_plustracer_aer(itype)+7
1545 dumtop = thesum(jj+1,ll) - thesum(jj,ll)
1546 dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35 )
1547 dumerr = dumtop/dumbot
1549 ! rce 21-jul-2006 - encountered some cases when delta_*_conform1 is negative
1550 ! and large in magnitude relative to thesum, which causes the mass
1551 ! conservation to be less accurate due to roundoff
1552 ! following section recomputes relative error with delta_*_conform1
1553 ! added onto each of thesum
1554 ! also increased dumtoler slightly
1555 if (ipass .eq. 2) then
1557 if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1558 dumc = delta_water_conform1
1559 else if (ll .eq. ncomp_plustracer_aer(itype)+3) then
1560 dumc = delta_numb_conform1
1562 if (dumc .lt. 0.0) then
1563 duma = thesum(3,ll) - dumc
1564 dumb = thesum(4,ll) - dumc
1565 dumtop = dumb - duma
1566 dumbot = max( abs(duma), abs(dumb), 1.0e-35 )
1567 dume = dumtop/dumbot
1568 if (abs(dume) .lt. abs(dumerr)) dumerr = dume
1572 if (abs(dumerr) .gt. abs(dumworst)) then
1574 dumworstb = dumworst
1581 if (abs(dumworst) .gt. dumtoler) then
1583 if (nerr .le. nerrmax) then
1585 call peg_message( lunout, msg )
1586 write(msg,97110) iclm, jclm, k, m, ipass, llworst
1587 call peg_message( lunout, msg )
1588 write(msg,97120) ' nnew(1,n)', &
1589 (nnewsave(1,n), n=1,nsize_aer(itype))
1590 call peg_message( lunout, msg )
1591 write(msg,97120) ' nnew(2,n)', &
1592 (nnewsave(2,n), n=1,nsize_aer(itype))
1593 call peg_message( lunout, msg )
1596 if (ll .eq. 0) ll = ncomp_plustracer_aer(itype)+2
1597 write(msg,97130) 'name/relerr/thesum', jj, '/thesum', jj+1, &
1598 dumname(ll), dumworst, thesum(jj,ll), thesum(jj+1,ll)
1599 call peg_message( lunout, msg )
1601 if ( (ll .eq. ncomp_plustracer_aer(itype)+3) .and. &
1602 (abs(dumworstb) .gt. dumtoler) ) then
1603 ll = max( 1, llworstb )
1604 dumtop = thesum(jj+1,ll) - thesum(jj,ll)
1605 dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35 )
1606 dumerr = dumtop/dumbot
1607 write(msg,97130) 'name/relerr/thesum', jj, '/thesum', jj+1, &
1608 dumname(ll), dumerr, thesum(jj,ll), thesum(jj+1,ll)
1609 call peg_message( lunout, msg )
1614 97110 format( 'movesect conserve ERROR - i/j/k/m/pass/llworst', &
1616 97120 format( a, 64i3 )
1617 97130 format( a, i1, a, i1, 2x, a, 1p, 3e16.7 )
1620 end subroutine move_sections_conserve_check
1623 !-----------------------------------------------------------------------
1624 subroutine test_move_sections( iflag_test, iclm, jclm, k, m )
1626 ! routine runs tests on move_sections, using a matrix of
1627 ! pregrow and aftgrow masses for each section
1634 ! include 'v33com9a'
1635 ! include 'v33com9b'
1638 integer iflag_test, iclm, jclm, k, m
1641 integer idiag_movesect, iflag, ii, iphase, itype, jj, &
1647 real dumnumb, dumvolpre, dumvolaft
1648 real dumsv_rsub(l2maxd)
1649 real dumsv_drymass_pregrow(maxd_asize,maxd_atype)
1650 real dumsv_drymass_aftgrow(maxd_asize,maxd_atype)
1651 real dumsv_drydens_pregrow(maxd_asize,maxd_atype)
1652 real dumsv_drydens_aftgrow(maxd_asize,maxd_atype)
1656 integer maxvolfactpre, maxvolfactaft
1657 parameter (maxvolfactpre=15, maxvolfactaft=23)
1659 real dumvolfactpre(maxvolfactpre)
1660 data dumvolfactpre / &
1661 2.0, 0.0, 1.0e-20, 0.5, 0.9, &
1662 1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0, &
1665 real dumvolfactaft(maxvolfactaft)
1666 data dumvolfactaft / &
1667 4.0, 0.0, 1.0e-20, 0.01, 0.02, 0.05, 0.1, 0.5, 0.9, &
1668 1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0, &
1669 8.1, 16.0, 32., 64., 128., 256. /
1673 ! check for valid inputs
1676 if (msectional .le. 0) return
1677 if (nsize_aer(1) .le. 0) return
1679 idiag_movesect = mod( msectional, 10000 )/100
1680 if (idiag_movesect .ne. 70) return
1682 ientryno = ientryno + 1
1683 if (ientryno .gt. 2) return
1686 ! save rsub and drymass/dens_aft/pregrow
1689 dumsv_rsub(l) = rsub(l,k,m)
1691 do itype = 1, ntype_aer
1692 do n = 1, nsize_aer(itype)
1693 dumsv_drymass_pregrow(n,itype) = drymass_pregrow(n,itype)
1694 dumsv_drymass_aftgrow(n,itype) = drymass_aftgrow(n,itype)
1695 dumsv_drydens_pregrow(n,itype) = drydens_pregrow(n,itype)
1696 dumsv_drydens_aftgrow(n,itype) = drydens_aftgrow(n,itype)
1701 ! make test calls to move_sections
1703 do 3900 iflag = 1, 1
1706 if (iabs(iflag) .eq. 2) iphase = cw_phase
1708 do 3800 itype = 1, ntype_aer
1710 do 2900 nn = 1, nsize_aer(itype)
1712 do 2800 ii = 1, maxvolfactpre
1714 do 2700 jj = 1, maxvolfactaft
1716 ! zero out rsub and dryxxx_yyygrow arrays
1717 do n = 1, nsize_aer(itype)
1718 do ll = 1, ncomp_plustracer_aer(itype)
1719 rsub(massptr_aer(ll,n,itype,iphase),k,m) = 0.0
1722 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1723 if (l .gt. 0) rsub(l,k,m) = 0.0
1724 l = numptr_aer(n,itype,iphase)
1725 if (l .gt. 0) rsub(l,k,m) = 0.0
1726 drymass_pregrow(n,itype) = 0.0
1727 drymass_aftgrow(n,itype) = 0.0
1728 drydens_pregrow(n,itype) = -1.0
1729 drydens_aftgrow(n,itype) = -1.0
1732 ! fill in values for section nn
1735 rsub(numptr_aer(n,itype,iphase),k,m) = dumnumb
1737 l = massptr_aer(ll,n,itype,iphase)
1739 dumvolpre = volumlo_sect(n,itype)*dumvolfactpre(ii)*dumnumb
1740 drydens_pregrow(n,itype) = dens_aer(ll,itype)
1741 drymass_pregrow(n,itype) = dumvolpre*drydens_pregrow(n,itype)
1742 if (ii .eq. 1) drydens_pregrow(n,itype) = -1.0
1744 dumvolaft = volumlo_sect(n,itype)*dumvolfactaft(jj)*dumnumb
1745 drydens_aftgrow(n,itype) = dens_aer(ll,itype)
1746 drymass_aftgrow(n,itype) = dumvolaft*drydens_aftgrow(n,itype)
1747 if (jj .eq. 1) drydens_aftgrow(n,itype) = -1.0
1749 rsub(l,k,m) = drymass_aftgrow(n,itype)/mw_aer(ll,itype)
1752 call peg_message( lunout, msg )
1753 write(msg,98010) nn, ii, jj
1754 call peg_message( lunout, msg )
1755 write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
1756 call peg_message( lunout, msg )
1758 ! make test call to move_sections
1759 call move_sections( iflag, iclm, jclm, k, m )
1762 call peg_message( lunout, msg )
1763 write(msg,98010) nn, ii, jj
1764 call peg_message( lunout, msg )
1765 write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
1766 call peg_message( lunout, msg )
1775 98010 format( 'test_move_sections output - nn, ii, jj =', 3i3 )
1776 98011 format( 'volfactpre, volfactaft =', 1p, 2e12.4 )
1779 ! restore rsub and drymass/dens_aft/pregrow
1782 rsub(l,k,m) = dumsv_rsub(l)
1784 do itype = 1, ntype_aer
1785 do n = 1, nsize_aer(itype)
1786 drymass_pregrow(n,itype) = dumsv_drymass_pregrow(n,itype)
1787 drymass_aftgrow(n,itype) = dumsv_drymass_aftgrow(n,itype)
1788 drydens_pregrow(n,itype) = dumsv_drydens_pregrow(n,itype)
1789 drydens_aftgrow(n,itype) = dumsv_drydens_aftgrow(n,itype)
1794 end subroutine test_move_sections
1797 !-----------------------------------------------------------------------
1798 subroutine move_sections_checkptrs( iflag, iclm, jclm, k, m )
1800 ! checks for valid number and water pointers
1807 ! include 'v33com9a'
1808 ! include 'v33com9b'
1811 integer iflag, iclm, jclm, k, m
1814 integer l, itype, iphase, n, ndum
1817 do 1900 itype = 1, ntype_aer
1818 do 1800 iphase = 1, nphase_aer
1821 do n = 1, nsize_aer(itype)
1822 l = numptr_aer(n,itype,iphase)
1823 if ((l .le. 0) .or. (l .gt. ltot2)) then
1824 msg = '*** subr move_sections error - ' // &
1825 'numptr_amode not defined'
1826 call peg_message( lunerr, msg )
1827 write(msg,9030) 'mode, numptr =', n, l
1828 call peg_message( lunerr, msg )
1829 write(msg,9030) 'iphase, itype =', iphase, itype
1830 call peg_message( lunerr, msg )
1831 call peg_error_fatal( lunerr, msg )
1833 ! checks involving nspec_amode and nspec_amode_nontracer
1834 ! being the same for all sections are no longer needed
1836 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1837 if ((l .gt. 0) .and. (l .le. ltot2)) ndum = ndum + 1
1839 if ((ndum .ne. 0) .and. (ndum .ne. nsize_aer(itype))) then
1840 msg = '*** subr move_sections error - ' // &
1841 'waterptr_aer must be on/off for all modes'
1842 call peg_message( lunerr, msg )
1843 write(msg,9030) 'iphase, itype =', iphase, itype
1844 call peg_message( lunerr, msg )
1845 call peg_error_fatal( lunerr, msg )
1847 9030 format( a, 2(1x,i6) )
1853 end subroutine move_sections_checkptrs
1857 end module module_mosaic_movesect