module_mosaic_movesect.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
9 module module_mosaic_movesect
10
11
12 use module_data_mosaic_asect
13 use module_data_mosaic_other
14 use module_peg_util
15
16
17 !
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)
21 ! growing into bin 1
22 ! if apply_n1_inflow = other, this correction is not applied
23 !
24 integer, parameter :: apply_n1_inflow = 0
25
26 contains
27
28
29
30 !-----------------------------------------------------------------------
31 !
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)
44 !
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 !-----------------------------------------------------------------------
53
54
55 !-----------------------------------------------------------------------
56 subroutine move_sections( iflag, iclm, jclm, k, m)
57 !
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
62 !
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)
71 !
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)
80 !
81 ! particle sizes are in cm
82 !
83 ! input parameters
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
90 ! before the growth
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
94 ! before the growth
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)
99 !
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
102 !
103 implicit none
104
105 ! include 'v33com'
106 ! include 'v33com2'
107 ! include 'v33com3'
108 ! include 'v33com9a'
109 ! include 'v33com9b'
110
111 ! subr arguments
112 integer iflag, iclm, jclm, k, m
113
114 ! local variables
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)
119
120 real densdefault, densh2o, smallmassaa, smallmassbb
121 real delta_water_conform1, delta_numb_conform1
122
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)
136
137 character*160 msg
138
139
140 !
141 ! check for valid inputs
142 !
143 if (msectional .le. 0) return
144 if (ntype_aer .le. 0) return
145 if (nphase_aer .le. 0) return
146
147 !
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
152 !
153 ! call test_move_sections( iflag, iclm, jclm, k, m )
154
155
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
159
160 ! get "idiag_movesect" from digits 3-4 of msectional
161 idiag_movesect = mod( msectional, 10000 )/100
162
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
169 continue
170 else if ((method_movesect .eq. 19) .or. &
171 (method_movesect .eq. 29) .or. &
172 (method_movesect .eq. 39)) then
173 return
174 else
175 msg = '*** subr move_sections error - ' // &
176 'illegal value for msectional'
177 call peg_error_fatal( lunerr, msg )
178 end if
179
180 if (iabs(iflag) .eq. 1) then
181 iphase = ai_phase
182 else if (iabs(iflag) .eq. 2) then
183 iphase = cw_phase
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 )
192 end if
193 else
194 msg = '*** subr move_sections error - ' // &
195 'iabs(iflag) must be 1 or 2'
196 call peg_error_fatal( lunerr, msg )
197 end if
198
199
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 )
208 return
209 end if
210 9040 format( '*** subr move_sections - ', a, ' =', i6 )
211
212 ! diagnostics
213 if (idiag_movesect .eq. 70) then
214 msg = ' '
215 call peg_message( lunout, msg )
216 write(msg,9060) iclm, jclm, k, m, msectional
217 call peg_message( lunout, msg )
218 end if
219 9060 format( '*** subr move_sections diagnostics i,j,k,m,msect =', 4i4, i6 )
220
221
222 densdefault = 2.0
223 densh2o = 1.0
224
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
231
232
233 ! process each type, one at a time
234 do 1900 itype = 1, ntype_aer
235
236 if (nsize_aer(itype) .le. 0) goto 1900
237
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 )
248
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 )
257 else
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 )
267 end if
268
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 )
278
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 )
292 end if
293
294 ! call move_sections_final_conform( &
295 ! iflag, iclm, jclm, k, m, iphase, itype )
296
297 1900 continue
298
299 return
300 end subroutine move_sections
301
302
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 )
314
315 !
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
324 !
325 implicit none
326
327 ! include 'v33com'
328 ! include 'v33com2'
329 ! include 'v33com3'
330 ! include 'v33com9a'
331 ! include 'v33com9b'
332
333 ! subr arguments
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), &
345 rnumbxx(maxd_asize)
346 real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
347 real wetvolxx(maxd_asize)
348 real wetmassxx(maxd_asize)
349
350
351 ! local variables
352 integer l, ll, lnew, lold, l3, n, nnew, nold
353
354 real dummass, dumnum, dumnum_at_dhi, dumnum_at_dlo, dumr, &
355 dumvol, dumvol1p, dumwatrmass
356
357
358 ! assure positive definite
359 do l = 1, ltot2
360 rsub(l,k,m) = max( 0., rsub(l,k,m) )
361 end do
362
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)
370 end do
371 rmassxx(llhysw,n) = 0.
372 l = 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.
376 l = 0
377 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
378 if (l .gt. 0) rmassxx(llwater,n) = rsub(l,k,m)
379
380 rnumbxx(n) = rsub(numptr_aer(n,itype,iphase),k,m)
381 rnumbxx0(n) = rnumbxx(n)
382 rnumbpp(n) = rnumbxx(n)
383
384 drydenspp(n) = drydens_pregrow(n,itype)
385 drymasspp(n) = drymass_pregrow(n,itype)
386
387 drydensxx(n) = drydens_aftgrow(n,itype)
388 drymassxx(n) = drymass_aftgrow(n,itype)
389 drydensxx0(n) = drydensxx(n)
390 drymassxx0(n) = drymassxx(n)
391 end do
392
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)
397 end do
398
399 delta_water_conform1 = 0.0
400 delta_numb_conform1 = 0.0
401
402
403 do 1390 n = 1, nsize_aer(itype)
404
405 !
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)
411 !
412 if ( (drydensxx(n) .lt. 0.1) .or. &
413 (drydensxx(n) .gt. 20.0) .or. &
414 (drymassxx(n) .le. smallmassaa) ) then
415 dummass = 0.
416 dumvol = 0.
417 do ll = 1, ncomp_aer(itype)
418 dumr = rmassxx(ll,n)*specmwxx(ll)
419 dummass = dummass + dumr
420 dumvol = dumvol + dumr/specdensxx(ll)
421 end do
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)
427 else
428 drydensxx(n) = dummass/dumvol
429 dumnum = rnumbxx(n)
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 ) )
433 end if
434 delta_numb_conform1 = delta_numb_conform1 + dumnum - rnumbxx(n)
435 rnumbxx(n) = dumnum
436 rnumbpp(n) = rnumbxx(n)
437 delta_water_conform1 = delta_water_conform1 - rmassxx(llwater,n)
438 rmassxx(llwater,n) = 0.
439 end if
440
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
447
448 1390 continue
449
450 return
451 end subroutine move_sections_initial_conform
452
453
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 )
462 !
463 ! routine calculates section movements for the moving-center approach
464 !
465 ! material in section n will be transfered to section nnewsave(1,n)
466 !
467 ! the nnewsave are calculated here
468 ! the actual transfer is done in another routine
469 !
470 implicit none
471
472 ! include 'v33com'
473 ! include 'v33com2'
474 ! include 'v33com3'
475 ! include 'v33com9a'
476 ! include 'v33com9b'
477
478 ! subr arguments
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)
491
492 ! local variables
493 integer ll, n, ndum, nnew, nold
494 real dumnum, dumvol, dumvol1p, sixoverpi, third
495 character*160 msg
496
497
498 sixoverpi = 6.0/pi
499 third = 1.0/3.0
500
501 !
502 ! compute mean size after growth (and corresponding section)
503 ! particles in section n will be transferred to section nnewsave(1,n)
504 !
505 do 1390 n = 1, nsize_aer(itype)
506
507 nnew = n
508
509 ! don't bother to transfer bins whose mass is extremely small
510 if (drymassxx(n) .le. smallmassaa) goto 1290
511
512 dumvol = dryvolxx(n)
513 dumnum = rnumbxx(n)
514
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)
519 goto 1290
520 ! or below that of smallest section
521 else if (dumnum .ge. dumvol/volumlo_sect(1,itype)) then
522 nnew = 1
523 goto 1290
524 end if
525
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)) )
531 nnew = nnew + 1
532 end do
533
534 else if (dumvol1p .lt. volumlo_sect(n,itype)) then
535 do while ( (nnew .gt. 1) .and. &
536 (dumvol1p .lt. volumlo_sect(nnew,itype)) )
537 nnew = nnew - 1
538 end do
539
540 end if
541
542 1290 nnewsave(1,n) = nnew
543 nnewsave(2,n) = 0
544
545 xferfracvol(1,n) = 1.0
546 xferfracvol(2,n) = 0.0
547 xferfracnum(1,n) = 1.0
548 xferfracnum(2,n) = 0.0
549
550 1390 continue
551
552
553 ! diagnostic output
554 ! output nnewsave values when msectional=7xxx
555 if (idiag_movesect .eq. 70) then
556 ndum = 0
557 do n = 1, nsize_aer(itype)
558 if (nnewsave(1,n) .ne. n) ndum = ndum + 1
559 end do
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 )
564 else
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 )
568 end if
569 end if
570 97751 format( 'movesect', a, 4i3, 3x, i3, 3x, 14i3 )
571
572 return
573 end subroutine move_sections_calc_movingcenter
574
575
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 )
586 !
587 ! routine calculates section movements for the mass-number-advection approach
588 !
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)
595 !
596 ! the nnewsave, xferfracvol, and xferfracnum are calculated here
597 ! the actual transfer is done in another routine
598 !
599 implicit none
600
601 ! include 'v33com'
602 ! include 'v33com2'
603 ! include 'v33com3'
604 ! include 'v33com9a'
605 ! include 'v33com9b'
606
607 ! subr arguments
608 integer iflag, iclm, jclm, iphase, itype, k, &
609 m, method_movesect, idiag_movesect, llhysw, llwater
610 integer nnewsave(2,maxd_asize)
611
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), &
620 rnumbxx(maxd_asize)
621 real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
622 real wetvolxx(maxd_asize)
623 real wetmassxx(maxd_asize)
624
625 ! local variables
626 integer ierr, n, nnew, nnew2
627 integer iforce_movecenter(maxd_asize)
628
629 real dum1, dum2, dum3
630 real dumaa, dumbb, dumgamma, dumratio
631 real dumfracnum, dumfracvol
632 real dumntot
633 real dumv
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
638 real dumzlo, dumzhi
639 real sixoverpi, third
640
641 character*4 dumch4
642 character*1 dumch1
643 character*160 msg
644
645
646 sixoverpi = 6.0/pi
647 third = 1.0/3.0
648
649 !
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)
652 !
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)
658 !
659 do 3900 n = 1, nsize_aer(itype)
660
661 nnew = n
662 iforce_movecenter(n) = 0
663
664 xferfracvol(1,n) = 1.0
665 xferfracvol(2,n) = 0.0
666 xferfracnum(1,n) = 1.0
667 xferfracnum(2,n) = 0.0
668
669 dumvtot_aft = -1.0
670 dumvtot_pre = -1.0
671 dumvbar_aft = -1.0
672 dumvbar_pre = -1.0
673 dumvlo_pre = -1.0
674 dumvhi_pre = -1.0
675 dumgamma = -1.0
676 dumratio = -1.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)
679 dumfracvol = -1.0
680 dumfracnum = -1.0
681
682 ! don't bother to transfer bins whose mass is extremely small
683 if (drymassxx(n) .le. smallmassaa) then
684 iforce_movecenter(n) = 1
685 goto 1290
686 end if
687
688 dumvtot_aft = dryvolxx(n)
689 dumntot = rnumbxx(n)
690
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
696 goto 1290
697 else if (dumntot .ge. dumvtot_aft/volumlo_sect(1,itype)) then
698 nnew = 1
699 iforce_movecenter(n) = 3
700 goto 1290
701 end if
702
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)) )
709 nnew = nnew + 1
710 end do
711
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)) )
715 nnew = nnew - 1
716 end do
717
718 end if
719
720 1290 nnewsave(1,n) = nnew
721 nnewsave(2,n) = 0
722
723 if (iforce_movecenter(n) .gt. 0) goto 3700
724
725
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
736 goto 3700
737 end if
738
739 dumvtot_pre = drymasspp(n)/drydenspp(n)
740
741 dumvlo_pre = volumlo_sect(n,itype)
742 dumvhi_pre = volumhi_sect(n,itype)
743 dumvdel_pre = dumvhi_pre - dumvlo_pre
744
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
750 goto 3700
751 end if
752 dumv = dumvlo_pre + 0.01*dumvdel_pre
753 if (dumntot .ge. dumvtot_pre/dumv) then
754 iforce_movecenter(n) = 13
755 goto 3700
756 end if
757
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
762
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
775 end if
776
777 dumbb = (dumratio - 1.0 - 0.5*dumgamma)*12.0/dumgamma
778 dumaa = 1.0 - 0.5*dumbb
779
780 ! calculate pregrow volumes corresponding to the nnew
781 ! section boundaries
782 dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
783 dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
784
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
791 else
792 nnew2 = nnew + 1
793 end if
794 else if (nnew .eq. nsize_aer(itype)) then
795 if (dumvlo_pre .ge. dumvcutlo_nnew_pre) then
796 iforce_movecenter(n) = 22
797 else
798 nnew2 = nnew - 1
799 end if
800 else
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
805 nnew2 = nnew - 1
806 else
807 nnew2 = nnew + 1
808 end if
809 end if
810 if (iforce_movecenter(n) .gt. 0) goto 3700
811
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)
824
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
830 end if
831 if (iforce_movecenter(n) .gt. 0) goto 3700
832
833 nnewsave(2,n) = nnew2
834
835 xferfracvol(1,n) = dumfracvol
836 xferfracvol(2,n) = 1.0 - dumfracvol
837 xferfracnum(1,n) = dumfracnum
838 xferfracnum(2,n) = 1.0 - dumfracnum
839
840 3700 continue
841
842 ! output nnewsave values when msectional=7xxx
843 if (idiag_movesect .ne. 70) goto 3800
844
845 if (nnewsave(2,n) .eq. 0) then
846 if (nnewsave(1,n) .eq. 0) then
847 dumch4 = 'NO X'
848 else if (nnewsave(1,n) .eq. n) then
849 dumch4 = 'NO A'
850 else
851 dumch4 = 'YESA'
852 end if
853 else if (nnewsave(1,n) .eq. 0) then
854 if (nnewsave(2,n) .eq. n) then
855 dumch4 = 'NO B'
856 else
857 dumch4 = 'YESB'
858 end if
859 else if (nnewsave(2,n) .eq. n) then
860 if (nnewsave(1,n) .eq. n) then
861 dumch4 = 'NO Y'
862 else
863 dumch4 = 'YESC'
864 end if
865 else if (nnewsave(1,n) .eq. n) then
866 dumch4 = 'YESD'
867 else
868 dumch4 = 'YESE'
869 end if
870
871 dumch1 = '+'
872 if (drymasspp(n) .gt. drymassxx(n)) dumch1 = '-'
873
874 msg = ' '
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 )
891 end if
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 )
911
912 97010 format( 'movesect', 2a, 7x, 4i3, 4x, &
913 'n,nnews', 3i3, 4x, 'iforce', i3.2 )
914 97020 format( a, 1p, 4e13.4 )
915
916 3800 continue
917
918 !
919 ! check for legal combinations of nnewsave(1,n) & nnewsave(2,n)
920 ! error if
921 ! nnew1 == nnew2
922 ! both are non-zero AND iabs(nnew1-nnew2) != 1
923 ierr = 0
924 if (nnewsave(1,n) .eq. nnewsave(2,n)) then
925 ierr = 1
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
928 end if
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 )
933 end if
934
935
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
939 nnewsave(1,n) = nnew
940 nnewsave(2,n) = 0
941 xferfracvol(1,n) = 1.0
942 xferfracvol(2,n) = 0.0
943 xferfracnum(1,n) = 1.0
944 xferfracnum(2,n) = 0.0
945 end if
946
947 3900 continue
948
949 return
950 end subroutine move_sections_calc_masnumadvect
951
952
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 )
963 !
964 ! routine performs the actual transfer of aerosol number and mass
965 ! between sections
966 !
967 implicit none
968
969 ! include 'v33com'
970 ! include 'v33com2'
971 ! include 'v33com3'
972 ! include 'v33com9a'
973 ! include 'v33com9b'
974
975 ! subr arguments
976 integer iflag, iclm, jclm, iphase, itype, k, &
977 m, method_movesect, idiag_movesect, llhysw, llwater
978 integer nnewsave(2,maxd_asize)
979
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)
993
994 ! local variables
995 integer jj, l, ll, n, ndum, nnew, nold
996 integer jja, jjb, jjc
997
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), &
1002 dumout(maxd_asize)
1003
1004 character*160 msg
1005 character*8 dumch8
1006 character*4 dumch4
1007
1008
1009 sixoverpi = 6.0/pi
1010 third = 1.0/3.0
1011
1012 !
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)
1016 !
1017 do 1900 n = 1, nsize_aer(itype)
1018
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)
1032 end do
1033
1034 else
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
1037 ! other sections
1038 ! so initialize "yy" arrays to zero
1039 drymassyy(n) = 0.0
1040 dryvolyy(n) = 0.0
1041 wetmassyy(n) = 0.0
1042 wetvolyy(n) = 0.0
1043 rnumbyy(n) = 0.0
1044 do ll = 1, ncomp_plustracer_aer(itype) + 2
1045 rmassyy(ll,n) = 0.0
1046 end do
1047
1048 end if
1049
1050 1900 continue
1051
1052 !
1053 ! do the transfer of mass and number
1054 !
1055 do 2900 nold = 1, nsize_aer(itype)
1056
1057 if ( (nnewsave(1,nold) .eq. nold) .and. &
1058 (nnewsave(2,nold) .eq. 0 ) ) goto 2900
1059
1060 do 2800 jj = 1, 2
1061
1062 nnew = nnewsave(jj,nold)
1063 if (nnew .le. 0) goto 2800
1064
1065 dumxfvol = xferfracvol(jj,nold)
1066 dumxfnum = xferfracnum(jj,nold)
1067
1068 do ll = 1, ncomp_plustracer_aer(itype) + 2
1069 rmassyy(ll,nnew) = rmassyy(ll,nnew) + rmassxx(ll,nold)*dumxfvol
1070 end do
1071 rnumbyy(nnew) = rnumbyy(nnew) + rnumbxx(nold)*dumxfnum
1072
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
1077
1078 2800 continue
1079
1080 2900 continue
1081
1082 !
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
1089 !
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 )
1097
1098 delta_numb_conform2 = 0.0
1099
1100 do 3900 n = 1, nsize_aer(itype)
1101
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)
1107 rnumbyy(n) = dumnum
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)
1112 else
1113 dumnum = rnumbyy(n)
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)
1118 rnumbyy(n) = dumnum
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 )
1126 end if
1127 aqvoldry_sub(n,itype,k,m) = dumvol
1128 aqmassdry_sub(n,itype,k,m) = drymassyy(n)
1129
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)
1133 ! end if
1134
1135 do ll = 1, ncomp_plustracer_aer(itype)
1136 l = massptr_aer(ll,n,itype,iphase)
1137 rsub(l,k,m) = rmassyy(ll,n)
1138 end do
1139 l = 0
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)
1145 end if
1146 rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
1147
1148 3900 continue
1149
1150 delta_numb_conform1 = delta_numb_conform1 + delta_numb_conform2
1151
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 )
1159
1160
1161 ! diagnostic output
1162 ! output nnewsave values when msectional=7xxx
1163 if (idiag_movesect .ne. 70) goto 4900
1164
1165 ndum = 0
1166 do n = 1, nsize_aer(itype)
1167 if (nnewsave(1,n)+nnewsave(2,n) .ne. n) ndum = ndum + 1
1168 end do
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 )
1173 ! else
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 )
1177 ! end if
1178
1179 dumch4 = 'NONE'
1180 if (ndum .gt. 0) dumch4 = 'SOME'
1181 msg = ' '
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 )
1189 end do
1190
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 )
1195
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 )
1200
1201 dumbot = log( volumhi_sect(1,itype)/volumlo_sect(1,itype) )
1202 do n = 1, nsize_aer(itype)
1203 dumpp(n) = -9.99
1204 dumxx(n) = -9.99
1205 dumyy(n) = -9.99
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
1213 end if
1214 end if
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
1219 end if
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
1224 end if
1225 end do
1226
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 )
1231
1232 do jja = 1, 7
1233 if (jja .eq. 1) then
1234 dumch8 = 'rnumbold'
1235 dumout(:) = rnumbxx(:)
1236 else if (jja .eq. 2) then
1237 dumch8 = 'rnumbnew'
1238 dumout(:) = rnumbyy(:)
1239 else if (jja .eq. 3) then
1240 dumch8 = 'drvolold'
1241 dumout(:) = dryvolxx(:)
1242 else if (jja .eq. 4) then
1243 dumch8 = 'drvolnew'
1244 dumout(:) = dryvolyy(:)
1245 else if (jja .eq. 5) then
1246 dumch8 = 'lnvolold'
1247 dumout(:) = dumxx(:)
1248 else if (jja .eq. 6) then
1249 dumch8 = 'lnvolnew'
1250 dumout(:) = dumyy(:)
1251 else if (jja .eq. 7) then
1252 dumch8 = 'lnvolpre'
1253 dumout(:) = dumpp(:)
1254 end if
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)
1259 else
1260 write(msg,97030) dumch8, (dumout(n), n=jjb,jjc)
1261 end if
1262 call peg_message( lunout, msg )
1263 dumch8 = ' '
1264 end do
1265 end do
1266
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 )
1274
1275 4900 continue
1276 return
1277 end subroutine move_sections_apply_moves
1278
1279
1280 !-----------------------------------------------------------------------
1281 ! rce 08-nov-2004 - this routine is new (and perhaps temporary)
1282 !
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 )
1292 !
1293 ! routine applies an ad_hoc correction to bin 1 to account for
1294 ! growth of smaller particles (which are not simulated) growing into
1295 ! bin 1
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)
1300 !
1301 implicit none
1302
1303 ! include 'v33com'
1304 ! include 'v33com2'
1305 ! include 'v33com3'
1306 ! include 'v33com9a'
1307 ! include 'v33com9b'
1308
1309 ! subr arguments
1310 integer iflag, iclm, jclm, iphase, itype, k, &
1311 m, method_movesect, idiag_movesect, llhysw, llwater
1312 integer nnewsave(2,maxd_asize)
1313
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)
1325
1326
1327 ! local variables
1328 integer jj, l, ll, n, nnew, nold
1329
1330 real deltanum, deltavol, dumvol1p, dumxfvol, dumxfnum
1331
1332
1333 !
1334 ! compute fraction of number transferred out of bin 1 by growth
1335 !
1336 nold = 1
1337 n = nold
1338 if ( (nnewsave(1,nold) .eq. nold) .and. &
1339 (nnewsave(2,nold) .eq. 0 ) ) goto 3900
1340
1341 dumxfnum = 0.0
1342 do 2800 jj = 1, 2
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)
1347 2800 continue
1348
1349 !
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
1353 !
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
1360
1361 !
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
1366 !
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)
1371 end do
1372 else
1373 ll = 1
1374 rmassyy(ll,n) = rmassyy(ll,n) + deltavol*specdensxx(ll)/specmwxx(ll)
1375 end if
1376 rnumbyy(n) = rnumbyy(n) + deltanum
1377
1378 !
1379 ! transfer results into rsub
1380 !
1381 do ll = 1, ncomp_plustracer_aer(itype)
1382 l = massptr_aer(ll,n,itype,iphase)
1383 rsub(l,k,m) = rmassyy(ll,n)
1384 end do
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)
1390 end if
1391 rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
1392
1393
1394 3900 continue
1395 return
1396 end subroutine move_sections_apply_n1_inflow
1397
1398
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 )
1407 !
1408 ! routine checks for conservation of number, mass, and volume
1409 ! by the move_sections algorithm
1410 !
1411 ! ipass = 1
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
1417 ! ipass = 2
1418 ! computes thesum(4,ll) from rsub after section movement
1419 ! compares thesum(3,ll) with thesum(4,ll)
1420 !
1421 ! currently only implemented for condensational growth (iflag=1)
1422 !
1423 implicit none
1424
1425 ! include 'v33com'
1426 ! include 'v33com2'
1427 ! include 'v33com3'
1428 ! include 'v33com9a'
1429 ! include 'v33com9b'
1430
1431 ! subr arguments
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)
1444
1445 ! local variables
1446 integer jj, l, ll, llworst, llworstb, n
1447 integer nerr, nerrmax
1448 save nerr, nerrmax
1449 data nerr, nerrmax / 0, 999 /
1450
1451 real dumbot, dumtop, dumtoler, dumerr, dumworst, dumworstb
1452 real duma, dumb, dumc, dume
1453 real thesum(4,maxd_acomp+7)
1454 save thesum
1455
1456 character*8 dumname(maxd_acomp+7)
1457 character*160 msg
1458
1459
1460 if (ipass .eq. 2) goto 2000
1461
1462 do ll = 1, ncomp_plustracer_aer(itype)+7
1463 do jj = 1, 4
1464 thesum(jj,ll) = 0.0
1465 end do
1466 end do
1467
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)
1472 end do
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)
1488 end do
1489
1490
1491 2000 continue
1492 !
1493 ! calc sum over bins for each species
1494 ! for water, account for loss in initial conform (delta_water_conform1)
1495 !
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
1501 l = 0
1502 if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
1503 else if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1504 l = 0
1505 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1506 else
1507 l = numptr_aer(n,itype,iphase)
1508 end if
1509 if (l .gt. 0) &
1510 thesum(ipass+2,ll) = thesum(ipass+2,ll) + rsub(l,k,m)
1511 end do
1512 end do
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
1518 end if
1519
1520 !
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
1524 !
1525 do ll = 1, ncomp_plustracer_aer(itype)+7
1526 dumname(ll) = ' '
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'
1537 end do
1538
1539 jj = 2*ipass - 1
1540 dumworst = 0.0
1541 dumworstb = 0.0
1542 llworst = 0
1543 llworstb = 0
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
1548
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
1556 dumc = 1.0
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
1561 end if
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
1569 end if
1570 end if
1571
1572 if (abs(dumerr) .gt. abs(dumworst)) then
1573 llworstb = llworst
1574 dumworstb = dumworst
1575 llworst = ll
1576 dumworst = dumerr
1577 end if
1578 end do
1579
1580 dumtoler = 1.0e-6
1581 if (abs(dumworst) .gt. dumtoler) then
1582 nerr = nerr + 1
1583 if (nerr .le. nerrmax) then
1584 msg = ' '
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 )
1594
1595 ll = llworst
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 )
1600
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 )
1610 end if
1611 end if
1612 end if
1613
1614 97110 format( 'movesect conserve ERROR - i/j/k/m/pass/llworst', &
1615 4i3, 2x, 2i3 )
1616 97120 format( a, 64i3 )
1617 97130 format( a, i1, a, i1, 2x, a, 1p, 3e16.7 )
1618
1619 return
1620 end subroutine move_sections_conserve_check
1621
1622
1623 !-----------------------------------------------------------------------
1624 subroutine test_move_sections( iflag_test, iclm, jclm, k, m )
1625 !
1626 ! routine runs tests on move_sections, using a matrix of
1627 ! pregrow and aftgrow masses for each section
1628 !
1629 implicit none
1630
1631 ! include 'v33com'
1632 ! include 'v33com2'
1633 ! include 'v33com3'
1634 ! include 'v33com9a'
1635 ! include 'v33com9b'
1636
1637 ! subr arguments
1638 integer iflag_test, iclm, jclm, k, m
1639
1640 ! local variables
1641 integer idiag_movesect, iflag, ii, iphase, itype, jj, &
1642 l, ll, n, nn
1643 integer ientryno
1644 save ientryno
1645 data ientryno / 0 /
1646
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)
1653
1654 character*160 msg
1655
1656 integer maxvolfactpre, maxvolfactaft
1657 parameter (maxvolfactpre=15, maxvolfactaft=23)
1658
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, &
1663 8.1, 16.0 /
1664
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. /
1670
1671
1672 !
1673 ! check for valid inputs
1674 ! and first entry
1675 !
1676 if (msectional .le. 0) return
1677 if (nsize_aer(1) .le. 0) return
1678
1679 idiag_movesect = mod( msectional, 10000 )/100
1680 if (idiag_movesect .ne. 70) return
1681
1682 ientryno = ientryno + 1
1683 if (ientryno .gt. 2) return
1684
1685 !
1686 ! save rsub and drymass/dens_aft/pregrow
1687 !
1688 do l = 1, ltot2
1689 dumsv_rsub(l) = rsub(l,k,m)
1690 end do
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)
1697 end do
1698 end do
1699
1700 !
1701 ! make test calls to move_sections
1702 !
1703 do 3900 iflag = 1, 1
1704
1705 iphase = ai_phase
1706 if (iabs(iflag) .eq. 2) iphase = cw_phase
1707
1708 do 3800 itype = 1, ntype_aer
1709
1710 do 2900 nn = 1, nsize_aer(itype)
1711
1712 do 2800 ii = 1, maxvolfactpre
1713
1714 do 2700 jj = 1, maxvolfactaft
1715
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
1720 end do
1721 l = 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
1730 end do
1731
1732 ! fill in values for section nn
1733 n = nn
1734 dumnumb = 1.0e7
1735 rsub(numptr_aer(n,itype,iphase),k,m) = dumnumb
1736 ll = 1
1737 l = massptr_aer(ll,n,itype,iphase)
1738
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
1743
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
1748
1749 rsub(l,k,m) = drymass_aftgrow(n,itype)/mw_aer(ll,itype)
1750
1751 msg = ' '
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 )
1757
1758 ! make test call to move_sections
1759 call move_sections( iflag, iclm, jclm, k, m )
1760
1761 msg = ' '
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 )
1767
1768 2700 continue
1769 2800 continue
1770 2900 continue
1771
1772 3800 continue
1773 3900 continue
1774
1775 98010 format( 'test_move_sections output - nn, ii, jj =', 3i3 )
1776 98011 format( 'volfactpre, volfactaft =', 1p, 2e12.4 )
1777
1778 !
1779 ! restore rsub and drymass/dens_aft/pregrow
1780 !
1781 do l = 1, ltot2
1782 rsub(l,k,m) = dumsv_rsub(l)
1783 end do
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)
1790 end do
1791 end do
1792
1793 return
1794 end subroutine test_move_sections
1795
1796
1797 !-----------------------------------------------------------------------
1798 subroutine move_sections_checkptrs( iflag, iclm, jclm, k, m )
1799 !
1800 ! checks for valid number and water pointers
1801 !
1802 implicit none
1803
1804 ! include 'v33com'
1805 ! include 'v33com2'
1806 ! include 'v33com3'
1807 ! include 'v33com9a'
1808 ! include 'v33com9b'
1809
1810 ! subr parameters
1811 integer iflag, iclm, jclm, k, m
1812
1813 ! local variables
1814 integer l, itype, iphase, n, ndum
1815 character*160 msg
1816
1817 do 1900 itype = 1, ntype_aer
1818 do 1800 iphase = 1, nphase_aer
1819
1820 ndum = 0
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 )
1832 end if
1833 ! checks involving nspec_amode and nspec_amode_nontracer
1834 ! being the same for all sections are no longer needed
1835 l = 0
1836 if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1837 if ((l .gt. 0) .and. (l .le. ltot2)) ndum = ndum + 1
1838 end do
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 )
1846 end if
1847 9030 format( a, 2(1x,i6) )
1848
1849 1800 continue
1850 1900 continue
1851
1852 return
1853 end subroutine move_sections_checkptrs
1854
1855
1856
1857 end module module_mosaic_movesect