added README_changes.txt
[wrffire.git] / wrfv2_fire / chem / module_mosaic_movesect.F
blobfd6c9ce3323be4e294abf159209e37db7dda08ae
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
14         use module_peg_util
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
24         integer, parameter :: apply_n1_inflow = 0
26         contains
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
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)
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
103         implicit none
105 !       include 'v33com'
106 !       include 'v33com2'
107 !       include 'v33com3'
108 !       include 'v33com9a'
109 !       include 'v33com9b'
111 !   subr arguments
112         integer iflag, iclm, jclm, k, m
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)
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)
137         character*160 msg
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
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
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
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 )
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 )
222         densdefault = 2.0
223         densh2o = 1.0
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 )
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
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 )
292         end if
294 !       call move_sections_final_conform(   &
295 !         iflag, iclm, jclm, k, m, iphase, itype )
297 1900    continue
299         return
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
325         implicit none
327 !       include 'v33com'
328 !       include 'v33com2'
329 !       include 'v33com3'
330 !       include 'v33com9a'
331 !       include 'v33com9b'
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)
351 !   local variables
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
359         do l = 1, ltot2
360             rsub(l,k,m) = max( 0., rsub(l,k,m) )
361         end do
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)
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)
391         end do
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
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
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
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
448 1390    continue
450         return
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
470         implicit none
472 !       include 'v33com'
473 !       include 'v33com2'
474 !       include 'v33com3'
475 !       include 'v33com9a'
476 !       include 'v33com9b'
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)
492 !   local variables
493         integer ll, n, ndum, nnew, nold
494         real dumnum, dumvol, dumvol1p, sixoverpi, third
495         character*160 msg
498         sixoverpi = 6.0/pi
499         third = 1.0/3.0
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)
507         nnew = n
509 !   don't bother to transfer bins whose mass is extremely small
510         if (drymassxx(n) .le. smallmassaa) goto 1290
512         dumvol = dryvolxx(n)
513         dumnum = rnumbxx(n)
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
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
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
540         end if
542 1290    nnewsave(1,n) = nnew
543         nnewsave(2,n) = 0
545         xferfracvol(1,n) = 1.0
546         xferfracvol(2,n) = 0.0
547         xferfracnum(1,n) = 1.0
548         xferfracnum(2,n) = 0.0
550 1390    continue
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 )
572         return
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
599         implicit none
601 !       include 'v33com'
602 !       include 'v33com2'
603 !       include 'v33com3'
604 !       include 'v33com9a'
605 !       include 'v33com9b'
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)
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)
625 !   local variables
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
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
641         character*4 dumch4
642         character*1 dumch1
643         character*160 msg
646         sixoverpi = 6.0/pi
647         third = 1.0/3.0
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)
661         nnew = n
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
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
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
688         dumvtot_aft = dryvolxx(n)
689         dumntot = rnumbxx(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
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
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
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
718         end if
720 1290    nnewsave(1,n) = nnew
721         nnewsave(2,n) = 0
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
736             goto 3700
737         end if
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
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
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
775         end if
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
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)
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
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
830         end if
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
840 3700    continue
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
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
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 )
912 97010   format( 'movesect', 2a, 7x, 4i3, 4x,   &
913                 'n,nnews', 3i3, 4x, 'iforce', i3.2 )
914 97020   format( a, 1p, 4e13.4 )
916 3800    continue
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
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
947 3900    continue
949         return
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
965 !       between sections
967         implicit none
969 !       include 'v33com'
970 !       include 'v33com2'
971 !       include 'v33com3'
972 !       include 'v33com9a'
973 !       include 'v33com9b'
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)
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)
994 !   local variables
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),   &
1002           dumout(maxd_asize)
1004         character*160 msg
1005         character*8 dumch8
1006         character*4 dumch4
1009         sixoverpi = 6.0/pi
1010         third = 1.0/3.0
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)
1032             end do
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
1048         end if
1050 1900    continue
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
1060         do 2800 jj = 1, 2
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
1070         end do
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
1078 2800    continue
1080 2900    continue
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)
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)
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
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)
1148 3900    continue
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 )
1161 !   diagnostic output
1162 !   output nnewsave values when msectional=7xxx
1163         if (idiag_movesect .ne. 70) goto 4900
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
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
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)
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
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 )
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
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 )
1275 4900    continue
1276         return
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
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)
1301         implicit none
1303 !       include 'v33com'
1304 !       include 'v33com2'
1305 !       include 'v33com3'
1306 !       include 'v33com9a'
1307 !       include 'v33com9b'
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)
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)
1327 !   local variables
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
1336         nold = 1
1337         n = nold
1338         if ( (nnewsave(1,nold) .eq. nold) .and.   &
1339              (nnewsave(2,nold) .eq. 0   ) ) goto 3900
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
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)
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
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)
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)
1394 3900    continue
1395         return
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
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)
1421 !   currently only implemented for condensational growth (iflag=1)
1423         implicit none
1425 !       include 'v33com'
1426 !       include 'v33com2'
1427 !       include 'v33com3'
1428 !       include 'v33com9a'
1429 !       include 'v33com9b'
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)
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 /
1451         real dumbot, dumtop, dumtoler, dumerr, dumworst, dumworstb
1452         real duma, dumb, dumc, dume
1453         real thesum(4,maxd_acomp+7)
1454         save thesum
1456         character*8 dumname(maxd_acomp+7)
1457         character*160 msg
1460         if (ipass .eq. 2) goto 2000
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
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
1491 2000    continue
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
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
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
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
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
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
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
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 )
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 )
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
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 )
1619         return
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
1629         implicit none
1631 !       include 'v33com'
1632 !       include 'v33com2'
1633 !       include 'v33com3'
1634 !       include 'v33com9a'
1635 !       include 'v33com9b'
1637 !   subr arguments
1638         integer iflag_test, iclm, jclm, k, m
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 /
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)
1654         character*160 msg
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,   &
1663         8.1, 16.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
1674 !   and first entry
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
1681         
1682         ientryno = ientryno + 1
1683         if (ientryno .gt. 2) return
1686 !   save rsub and drymass/dens_aft/pregrow
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
1701 !   make test calls to move_sections
1703         do 3900 iflag = 1, 1
1705         iphase = ai_phase
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
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
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)
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 )
1758 !   make test call to move_sections
1759         call move_sections( iflag, iclm, jclm, k, m )
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 )
1768 2700    continue
1769 2800    continue
1770 2900    continue
1772 3800    continue
1773 3900    continue
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
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
1793         return
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
1802         implicit none
1804 !       include 'v33com'
1805 !       include 'v33com2'
1806 !       include 'v33com3'
1807 !       include 'v33com9a'
1808 !       include 'v33com9b'
1810 !   subr parameters
1811         integer iflag, iclm, jclm, k, m
1813 !   local variables
1814         integer l, itype, iphase, n, ndum
1815         character*160 msg
1817         do 1900 itype = 1, ntype_aer
1818         do 1800 iphase = 1, nphase_aer
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) )
1849 1800    continue
1850 1900    continue
1852         return
1853         end subroutine move_sections_checkptrs                           
1857         end module module_mosaic_movesect