!**********************************************************************************
! This computer software was prepared by Battelle Memorial Institute, hereinafter
! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! MOSAIC module: see module_mosaic_driver.F for references and terms of use
!**********************************************************************************
module module_mosaic_movesect
use module_data_mosaic_asect
use module_data_mosaic_other
use module_peg_util
!
! if apply_n1_inflow = 1, then subr move_sections_apply_n1_inflow
! applies an ad_hoc correction to bin 1 to account for growth of
! smaller particles (which are not simulated when aernucnew_onoff=0)
! growing into bin 1
! if apply_n1_inflow = other, this correction is not applied
!
integer, parameter :: apply_n1_inflow = 0
contains
!-----------------------------------------------------------------------
!
! zz08movesect.f - created 24-nov-01 from scm movebin code
! 04-dec-01 rce - added code to treat bins that had state="no_aerosol"
! 10-dec-01 rce - diagnostic writes to fort.96 deleted
! 08-aug-02 rce - this is latest version from freshair scaqs-aerosol code
! 03-aug-02 rce - bypass this routine when msectional=701
! output nnewsave values to lunout when msectional=702
! 29-aug-03 rce - use nspec_amode_nontracer in first "do ll" loop
! 12-nov-03 rce - two approaches now available
! jacobson moving-center (old) - when msectional=10
! tzivion mass-number advection (new) - when msectional=20
! 28-jan-04 rce - eliminated the move_sections_old code
! (no reason to keep it)
!
! rce 2004-dec-03 - many changes associated with the new aerosol "pointer"
! variables in module_data_mosaic_asect
! rce 2005-feb-17 - fixes involving drydens_pre/aftgrow, drymass_...,
! & mw_aer. All are dimensioned (isize,itype) but were being used
! as (itype). In old compiler, this was OK, and treated as (itype,1).
! In new compiler, this caused an error.
! Also, in subr test_move_sections, set iphase based on iflag.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine move_sections( iflag, iclm, jclm, k, m)
!
! routine transfers aerosol number and mass between sections
! to account for continuous aerosol growth
! this routine is called after the gas condensation module (MOSAIC) or
! aqueous chemistry module has increased the mass within sections
!
! moving-center algorithm or mass-number advection algorithm is used,
! depending on value of mod(msectional,100)
! section mean diameter is given by
! vtot = ntot * (pi/6) * (dmean**3)
! where vtot and ntot are total dry-volume and number for the section
! if dmean is outside the section boundaries (dlo_sect & dhi_sect), then
! all the mass and number in the section are transfered to the
! section with dlo_sect(nnew) < dmean < dhi_sect(nnew)
!
! mass mixing ratios (mole/mole-air or g/mole-air) are in
! rsub(massptr_aer(ll,n,itype,iphase),k,m)
! number mixing ratios (particles/mole-air) are in
! rsub(numptr_aer(n,itype,iphase),k,m)
! these values are over-written with new values
! the following are also updated:
! adrydens_sub(n,itype,k,m), admeandry_sub(n,itype,k,m),
! awetdens_sub(n,itype,k,m), admeanwet_sub(n,itype,k,m)
!
! particle sizes are in cm
!
! input parameters
! iflag = 1 - do transfer after trace-gas condensation/evaporation
! = 2 - do transfer after aqueous chemistry
! = -1/-2 - do some "first entry" tasks for the iflag=+1/+2 cases
! iclm, jclm, k = current i,j,k indices
! m = current subarea index
! drymass_pregrow(n,itype) = dry-mass (g/mole-air) for section n
! before the growth
! drymass_aftgrow(n,itype) = dry-mass (g/mole-air) for section n
! after the growth but before inter-section transfer
! drydens_pregrow(n,itype) = dry-density (g/cm3) for section n
! before the growth
! drydens_aftgrow(n,itype) = dry-density (g/cm3) for section n
! after the growth but before inter-section transfer
! (drymass_pregrow and drydens_pregrow are not used by the moving-center
! algorithm, but would be needed for other sectional algorithms)
!
! this routine is the top level module for the post-october-2003 code
! and will do either moving-center or mass-number advection
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag, iclm, jclm, k, m
! local variables
integer idiag_movesect, iphase, itype, &
l, ll, llhysw, llwater, lnew, lold, l3, &
method_movesect, n, nnew, nold
integer nnewsave(2,maxd_asize)
real densdefault, densh2o, smallmassaa, smallmassbb
real delta_water_conform1, delta_numb_conform1
real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
drydensxx(maxd_asize), drydensyy(maxd_asize)
real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
drymassxx(maxd_asize), drymassyy(maxd_asize)
real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize), &
rmassyy(maxd_acomp+2,maxd_asize)
real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
rnumbxx(maxd_asize), rnumbyy(maxd_asize)
real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
character*160 msg
!
! check for valid inputs
!
if (msectional .le. 0) return
if (ntype_aer .le. 0) return
if (nphase_aer .le. 0) return
!
! run diagnostic tests
! (these will only be run for certain values of idiag_movesect
! rce 2004-nov-29 - to avoid recursion, test_move_sections
! is now called from mosaic_driver
!
! call test_move_sections( iflag, iclm, jclm, k, m )
! get "method_movesect" from digits 1-2 of msectional (treat 1-9 as 10)
method_movesect = mod( msectional, 100 )
if (method_movesect .le. 10) method_movesect = 10
! get "idiag_movesect" from digits 3-4 of msectional
idiag_movesect = mod( msectional, 10000 )/100
if ((method_movesect .eq. 10) .or. &
(method_movesect .eq. 11) .or. &
(method_movesect .eq. 20) .or. &
(method_movesect .eq. 21) .or. &
(method_movesect .eq. 30) .or. &
(method_movesect .eq. 31)) then
continue
else if ((method_movesect .eq. 19) .or. &
(method_movesect .eq. 29) .or. &
(method_movesect .eq. 39)) then
return
else
msg = '*** subr move_sections error - ' // &
'illegal value for msectional'
call peg_error_fatal( lunerr, msg )
end if
if (iabs(iflag) .eq. 1) then
iphase = ai_phase
else if (iabs(iflag) .eq. 2) then
iphase = cw_phase
if (nphase_aer .lt. 2) then
msg = '*** subr move_sections error - ' // &
'iflag=2 (after aqueous chemistry) but nphase_aer < 2'
call peg_error_fatal( lunerr, msg )
else if (cw_phase .ne. 2) then
msg = '*** subr move_sections error - ' // &
'iflag=2 (after aqueous chemistry) but cw_phase .ne. 2'
call peg_error_fatal( lunerr, msg )
end if
else
msg = '*** subr move_sections error - ' // &
'iabs(iflag) must be 1 or 2'
call peg_error_fatal( lunerr, msg )
end if
! when iflag=-1/-2, call move_sections_checkptrs then return
! if ((ncorecnt .le. 0) .and. (k .le. 1)) then
if (iflag .le. 0) then
write(msg,9040) 'method', method_movesect
call peg_message( lunout, msg )
write(msg,9040) 'idiag ', idiag_movesect
call peg_message( lunout, msg )
call move_sections_checkptrs( iflag, iclm, jclm, k, m )
return
end if
9040 format( '*** subr move_sections - ', a, ' =', i6 )
! diagnostics
if (idiag_movesect .eq. 70) then
msg = ' '
call peg_message( lunout, msg )
write(msg,9060) iclm, jclm, k, m, msectional
call peg_message( lunout, msg )
end if
9060 format( '*** subr move_sections diagnostics i,j,k,m,msect =', 4i4, i6 )
densdefault = 2.0
densh2o = 1.0
! if bin mass mixrat < smallmassaa (1.e-20 g/mol), then assume no growth
! AND no water AND conform number so that size is within bin limits
smallmassaa = 1.0e-20
! if bin mass OR volume mixrat < smallmassbb (1.e-30 g/mol), then
! assume default density to avoid divide by zero
smallmassbb = 1.0e-30
! process each type, one at a time
do 1900 itype = 1, ntype_aer
if (nsize_aer(itype) .le. 0) goto 1900
call move_sections_initial_conform( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drydenspp, drymasspp, rnumbpp, &
drydensxx0, drymassxx0, rnumbxx0, &
drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
wetmassxx, wetvolxx, &
specmwxx, specdensxx )
if (method_movesect .le. 19) then
call move_sections_calc_movingcenter( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
wetmassxx, wetvolxx, &
xferfracvol, xferfracnum )
else
call move_sections_calc_masnumadvect( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
drydenspp, drymasspp, rnumbpp, &
drydensxx0, drymassxx0, rnumbxx0, &
drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
wetmassxx, wetvolxx, &
xferfracvol, xferfracnum )
end if
call move_sections_apply_moves( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drydenspp, drymasspp, rnumbpp, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
xferfracvol, xferfracnum )
! rce 08-nov-2004 - this call is new (and perhaps temporary)
! rce 05-jul-2005 - do n1_inflow for aerchem growth but not for cldchem
if ((apply_n1_inflow .eq. 1) .and. &
(iphase .eq. ai_phase)) then
call move_sections_apply_n1_inflow( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
xferfracvol, xferfracnum, &
specmwxx, specdensxx )
end if
! call move_sections_final_conform( &
! iflag, iclm, jclm, k, m, iphase, itype )
1900 continue
return
end subroutine move_sections
!-----------------------------------------------------------------------
subroutine move_sections_initial_conform( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drydenspp, drymasspp, rnumbpp, &
drydensxx0, drymassxx0, rnumbxx0, &
drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
wetmassxx, wetvolxx, &
specmwxx, specdensxx )
!
! routine does some initial tasks for the section movement
! load rmassxx & rnumbxx from rsub
! load specmwxx & specdensxx
! set drymassxx & dryvolxx from drymass_aftgrow & drydens_aftgrow,
! OR compute them from rmassxx, specmwxx, specdensxx if need be
! set wetmassxx & wetvolxx from dry values & water mass
! conform rnumbxx so that the mean particle size of each section
! (= dryvolxx/rnumbxx) is within the section limits
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag, iclm, jclm, iphase, itype, k, &
m, method_movesect, idiag_movesect, llhysw, llwater
real densdefault, densh2o, smallmassaa, smallmassbb
real delta_water_conform1, delta_numb_conform1
real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
drydensxx(maxd_asize)
real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
drymassxx(maxd_asize)
real dryvolxx(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize)
real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
rnumbxx(maxd_asize)
real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
real wetvolxx(maxd_asize)
real wetmassxx(maxd_asize)
! local variables
integer l, ll, lnew, lold, l3, n, nnew, nold
real dummass, dumnum, dumnum_at_dhi, dumnum_at_dlo, dumr, &
dumvol, dumvol1p, dumwatrmass
! assure positive definite
do l = 1, ltot2
rsub(l,k,m) = max( 0., rsub(l,k,m) )
end do
! load mixrats into working arrays and assure positive definite
llhysw = ncomp_plustracer_aer(itype) + 1
llwater = ncomp_plustracer_aer(itype) + 2
do n = 1, nsize_aer(itype)
do ll = 1, ncomp_plustracer_aer(itype)
l = massptr_aer(ll,n,itype,iphase)
rmassxx(ll,n) = rsub(l,k,m)
end do
rmassxx(llhysw,n) = 0.
l = 0
if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
if (l .gt. 0) rmassxx(llhysw,n) = rsub(l,k,m)
rmassxx(llwater,n) = 0.
l = 0
if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
if (l .gt. 0) rmassxx(llwater,n) = rsub(l,k,m)
rnumbxx(n) = rsub(numptr_aer(n,itype,iphase),k,m)
rnumbxx0(n) = rnumbxx(n)
rnumbpp(n) = rnumbxx(n)
drydenspp(n) = drydens_pregrow(n,itype)
drymasspp(n) = drymass_pregrow(n,itype)
drydensxx(n) = drydens_aftgrow(n,itype)
drymassxx(n) = drymass_aftgrow(n,itype)
drydensxx0(n) = drydensxx(n)
drymassxx0(n) = drymassxx(n)
end do
! load specmw and specdens also
do ll = 1, ncomp_plustracer_aer(itype)
specmwxx(ll) = mw_aer(ll,itype)
specdensxx(ll) = dens_aer(ll,itype)
end do
delta_water_conform1 = 0.0
delta_numb_conform1 = 0.0
do 1390 n = 1, nsize_aer(itype)
!
! if drydens_aftgrow < 0.1, then bin had state="no_aerosol"
! compute volume using default dry-densities, set water=0,
! and conform the number
! also do this if mass is extremely small (below smallmassaa)
! OR if drydens_aftgrow > 20 (which is unrealistic)
!
if ( (drydensxx(n) .lt. 0.1) .or. &
(drydensxx(n) .gt. 20.0) .or. &
(drymassxx(n) .le. smallmassaa) ) then
dummass = 0.
dumvol = 0.
do ll = 1, ncomp_aer(itype)
dumr = rmassxx(ll,n)*specmwxx(ll)
dummass = dummass + dumr
dumvol = dumvol + dumr/specdensxx(ll)
end do
drymassxx(n) = dummass
if (min(dummass,dumvol) .le. smallmassbb) then
drydensxx(n) = densdefault
dumvol = dummass/densdefault
dumnum = dummass/(volumcen_sect(n,itype)*densdefault)
else
drydensxx(n) = dummass/dumvol
dumnum = rnumbxx(n)
dumnum_at_dhi = dumvol/volumhi_sect(n,itype)
dumnum_at_dlo = dumvol/volumlo_sect(n,itype)
dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
end if
delta_numb_conform1 = delta_numb_conform1 + dumnum - rnumbxx(n)
rnumbxx(n) = dumnum
rnumbpp(n) = rnumbxx(n)
delta_water_conform1 = delta_water_conform1 - rmassxx(llwater,n)
rmassxx(llwater,n) = 0.
end if
! load dry/wet mass and volume into "xx" arrays
! which hold values before inter-mode transferring
dryvolxx(n) = drymassxx(n)/drydensxx(n)
dumwatrmass = rmassxx(llwater,n)*mw_water_aer
wetmassxx(n) = drymassxx(n) + dumwatrmass
wetvolxx(n) = dryvolxx(n) + dumwatrmass/densh2o
1390 continue
return
end subroutine move_sections_initial_conform
!-----------------------------------------------------------------------
subroutine move_sections_calc_movingcenter( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
wetmassxx, wetvolxx, &
xferfracvol, xferfracnum )
!
! routine calculates section movements for the moving-center approach
!
! material in section n will be transfered to section nnewsave(1,n)
!
! the nnewsave are calculated here
! the actual transfer is done in another routine
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag, iclm, jclm, iphase, itype, k, &
m, method_movesect, idiag_movesect, llhysw, llwater
integer nnewsave(2,maxd_asize)
real densdefault, densh2o, smallmassaa, smallmassbb
real drydensxx(maxd_asize)
real drymassxx(maxd_asize)
real dryvolxx(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize)
real rnumbxx(maxd_asize)
real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
real wetmassxx(maxd_asize)
real wetvolxx(maxd_asize)
! local variables
integer ll, n, ndum, nnew, nold
real dumnum, dumvol, dumvol1p, sixoverpi, third
character*160 msg
sixoverpi = 6.0/pi
third = 1.0/3.0
!
! compute mean size after growth (and corresponding section)
! particles in section n will be transferred to section nnewsave(1,n)
!
do 1390 n = 1, nsize_aer(itype)
nnew = n
! don't bother to transfer bins whose mass is extremely small
if (drymassxx(n) .le. smallmassaa) goto 1290
dumvol = dryvolxx(n)
dumnum = rnumbxx(n)
! check for number so small that particle volume is
! above that of largest section
if (dumnum .le. dumvol/volumhi_sect(nsize_aer(itype),itype)) then
nnew = nsize_aer(itype)
goto 1290
! or below that of smallest section
else if (dumnum .ge. dumvol/volumlo_sect(1,itype)) then
nnew = 1
goto 1290
end if
! dumvol1p is mean particle volume (cm3) for the section
dumvol1p = dumvol/dumnum
if (dumvol1p .gt. volumhi_sect(n,itype)) then
do while ( (nnew .lt. nsize_aer(itype)) .and. &
(dumvol1p .gt. volumhi_sect(nnew,itype)) )
nnew = nnew + 1
end do
else if (dumvol1p .lt. volumlo_sect(n,itype)) then
do while ( (nnew .gt. 1) .and. &
(dumvol1p .lt. volumlo_sect(nnew,itype)) )
nnew = nnew - 1
end do
end if
1290 nnewsave(1,n) = nnew
nnewsave(2,n) = 0
xferfracvol(1,n) = 1.0
xferfracvol(2,n) = 0.0
xferfracnum(1,n) = 1.0
xferfracnum(2,n) = 0.0
1390 continue
! diagnostic output
! output nnewsave values when msectional=7xxx
if (idiag_movesect .eq. 70) then
ndum = 0
do n = 1, nsize_aer(itype)
if (nnewsave(1,n) .ne. n) ndum = ndum + 1
end do
if (ndum .gt. 0) then
write(msg,97751) 'YES', iclm, jclm, k, m, &
ndum, (nnewsave(1,n), n=1,nsize_aer(itype))
call peg_message( lunout, msg )
else
write(msg,97751) 'NO ', iclm, jclm, k, m, &
ndum, (nnewsave(1,n), n=1,nsize_aer(itype))
call peg_message( lunout, msg )
end if
end if
97751 format( 'movesect', a, 4i3, 3x, i3, 3x, 14i3 )
return
end subroutine move_sections_calc_movingcenter
!-----------------------------------------------------------------------
subroutine move_sections_calc_masnumadvect( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
drydenspp, drymasspp, rnumbpp, &
drydensxx0, drymassxx0, rnumbxx0, &
drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx, &
wetmassxx, wetvolxx, &
xferfracvol, xferfracnum )
!
! routine calculates section movements for the mass-number-advection approach
!
! material in section n will be transfered to sections
! nnewsave(1,n) and nnewsave(2,n)
! the fractions of mass/volume transfered to each are
! xferfracvol(1,n) and xferfracvol(2,n)
! the fractions of number transfered to each are
! xferfracnum(1,n) and xferfracnum(2,n)
!
! the nnewsave, xferfracvol, and xferfracnum are calculated here
! the actual transfer is done in another routine
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag, iclm, jclm, iphase, itype, k, &
m, method_movesect, idiag_movesect, llhysw, llwater
integer nnewsave(2,maxd_asize)
real densdefault, densh2o, smallmassaa, smallmassbb
real drydenspp(maxd_asize), drydensxx0(maxd_asize), &
drydensxx(maxd_asize)
real drymasspp(maxd_asize), drymassxx0(maxd_asize), &
drymassxx(maxd_asize)
real dryvolxx(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize)
real rnumbpp(maxd_asize), rnumbxx0(maxd_asize), &
rnumbxx(maxd_asize)
real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
real wetvolxx(maxd_asize)
real wetmassxx(maxd_asize)
! local variables
integer ierr, n, nnew, nnew2
integer iforce_movecenter(maxd_asize)
real dum1, dum2, dum3
real dumaa, dumbb, dumgamma, dumratio
real dumfracnum, dumfracvol
real dumntot
real dumv
real dumvbar_aft, dumvbar_pre
real dumvcutlo_nnew_pre, dumvcuthi_nnew_pre
real dumvlo_pre, dumvhi_pre, dumvdel_pre
real dumvtot_aft, dumvtot_pre
real dumzlo, dumzhi
real sixoverpi, third
character*4 dumch4
character*1 dumch1
character*160 msg
sixoverpi = 6.0/pi
third = 1.0/3.0
!
! compute mean size after growth (and corresponding section)
! some of the particles in section n will be transferred to section nnewsave(1,n)
!
! if the aftgrow mass is extremely small,
! OR if the aftgrow mean size is outside of
! [dlo_sect(1,itype), dhi_sect(nsize_aer(itype),itype)]
! then use the moving-center method_movesect for this bin
! (don't try to compute the pregrow within-bin distribution)
!
do 3900 n = 1, nsize_aer(itype)
nnew = n
iforce_movecenter(n) = 0
xferfracvol(1,n) = 1.0
xferfracvol(2,n) = 0.0
xferfracnum(1,n) = 1.0
xferfracnum(2,n) = 0.0
dumvtot_aft = -1.0
dumvtot_pre = -1.0
dumvbar_aft = -1.0
dumvbar_pre = -1.0
dumvlo_pre = -1.0
dumvhi_pre = -1.0
dumgamma = -1.0
dumratio = -1.0
dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
dumfracvol = -1.0
dumfracnum = -1.0
! don't bother to transfer bins whose mass is extremely small
if (drymassxx(n) .le. smallmassaa) then
iforce_movecenter(n) = 1
goto 1290
end if
dumvtot_aft = dryvolxx(n)
dumntot = rnumbxx(n)
! check for particle volume above that of largest section
! or below that of smallest section
if (dumntot .le. dumvtot_aft/volumhi_sect(nsize_aer(itype),itype)) then
nnew = nsize_aer(itype)
iforce_movecenter(n) = 2
goto 1290
else if (dumntot .ge. dumvtot_aft/volumlo_sect(1,itype)) then
nnew = 1
iforce_movecenter(n) = 3
goto 1290
end if
! dumvbar_aft is mean particle volume (cm3) for the section
! find the section that encloses this volume
dumvbar_aft = dumvtot_aft/dumntot
if (dumvbar_aft .gt. volumhi_sect(n,itype)) then
do while ( (nnew .lt. nsize_aer(itype)) .and. &
(dumvbar_aft .gt. volumhi_sect(nnew,itype)) )
nnew = nnew + 1
end do
else if (dumvbar_aft .lt. volumlo_sect(n,itype)) then
do while ( (nnew .gt. 1) .and. &
(dumvbar_aft .lt. volumlo_sect(nnew,itype)) )
nnew = nnew - 1
end do
end if
1290 nnewsave(1,n) = nnew
nnewsave(2,n) = 0
if (iforce_movecenter(n) .gt. 0) goto 3700
! if drydenspp (pregrow) < 0.1 (because bin had state="no_aerosol" before
! growth was computed, so its initial mass was very small)
! then use the moving-center method_movesect for this bin
! (don't try to compute the pregrow within-bin distribution)
! also do this if pregrow mass is extremely small (below smallmassaa)
! OR if drydenspp > 20 (unphysical)
if ( (drydenspp(n) .lt. 0.1) .or. &
(drydenspp(n) .gt. 20.0) .or. &
(drymasspp(n) .le. smallmassaa) ) then
iforce_movecenter(n) = 11
goto 3700
end if
dumvtot_pre = drymasspp(n)/drydenspp(n)
dumvlo_pre = volumlo_sect(n,itype)
dumvhi_pre = volumhi_sect(n,itype)
dumvdel_pre = dumvhi_pre - dumvlo_pre
! if the pregrow mean size is outside of OR very close to the bin limits,
! then use moving-center approach for this bin
dumv = dumvhi_pre - 0.01*dumvdel_pre
if (dumntot .le. dumvtot_pre/dumv) then
iforce_movecenter(n) = 12
goto 3700
end if
dumv = dumvlo_pre + 0.01*dumvdel_pre
if (dumntot .ge. dumvtot_pre/dumv) then
iforce_movecenter(n) = 13
goto 3700
end if
! calculate the pregrow within-section size distribution
dumvbar_pre = dumvtot_pre/dumntot
dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
dumratio = dumvbar_pre/dumvlo_pre
if (dumratio .le. (1.0001 + dumgamma/3.0)) then
dumv = dumvlo_pre + 3.0*(dumvbar_pre-dumvlo_pre)
dumvhi_pre = min( dumvhi_pre, dumv )
dumvdel_pre = dumvhi_pre - dumvlo_pre
dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
dumratio = dumvbar_pre/dumvlo_pre
else if (dumratio .ge. (0.9999 + dumgamma*2.0/3.0)) then
dumv = dumvhi_pre + 3.0*(dumvbar_pre-dumvhi_pre)
dumvlo_pre = max( dumvlo_pre, dumv )
dumvdel_pre = dumvhi_pre - dumvlo_pre
dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
dumratio = dumvbar_pre/dumvlo_pre
end if
dumbb = (dumratio - 1.0 - 0.5*dumgamma)*12.0/dumgamma
dumaa = 1.0 - 0.5*dumbb
! calculate pregrow volumes corresponding to the nnew
! section boundaries
dumvcutlo_nnew_pre = volumlo_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
dumvcuthi_nnew_pre = volumhi_sect(nnew,itype)*(dumvbar_pre/dumvbar_aft)
! if the [dumvlo_pre, dumvhi_pre] falls completely within
! the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval,
! then all mass and number go to nnew
if (nnew .eq. 1) then
if (dumvhi_pre .le. dumvcuthi_nnew_pre) then
iforce_movecenter(n) = 21
else
nnew2 = nnew + 1
end if
else if (nnew .eq. nsize_aer(itype)) then
if (dumvlo_pre .ge. dumvcutlo_nnew_pre) then
iforce_movecenter(n) = 22
else
nnew2 = nnew - 1
end if
else
if ((dumvlo_pre .ge. dumvcutlo_nnew_pre) .and. &
(dumvhi_pre .le. dumvcuthi_nnew_pre)) then
iforce_movecenter(n) = 23
else if (dumvlo_pre .lt. dumvcutlo_nnew_pre) then
nnew2 = nnew - 1
else
nnew2 = nnew + 1
end if
end if
if (iforce_movecenter(n) .gt. 0) goto 3700
! calculate the fraction of ntot and vtot that are within
! the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval
dumzlo = (dumvcutlo_nnew_pre - dumvlo_pre)/dumvdel_pre
dumzhi = (dumvcuthi_nnew_pre - dumvlo_pre)/dumvdel_pre
dumzlo = max( dumzlo, 0.0 )
dumzhi = min( dumzhi, 1.0 )
dum1 = dumzhi - dumzlo
dum2 = (dumzhi**2 - dumzlo**2)*0.5
dum3 = (dumzhi**3 - dumzlo**3)/3.0
dumfracnum = dumaa*dum1 + dumbb*dum2
dumfracvol = (dumvlo_pre/dumvbar_pre) * (dumaa*dum1 + &
(dumaa*dumgamma + dumbb)*dum2 + (dumbb*dumgamma)*dum3)
if ((dumfracnum .le. 0.0) .or. (dumfracvol .le. 0.0)) then
iforce_movecenter(n) = 31
nnewsave(1,n) = nnew2
else if ((dumfracnum .ge. 1.0) .or. (dumfracvol .ge. 1.0)) then
iforce_movecenter(n) = 32
end if
if (iforce_movecenter(n) .gt. 0) goto 3700
nnewsave(2,n) = nnew2
xferfracvol(1,n) = dumfracvol
xferfracvol(2,n) = 1.0 - dumfracvol
xferfracnum(1,n) = dumfracnum
xferfracnum(2,n) = 1.0 - dumfracnum
3700 continue
! output nnewsave values when msectional=7xxx
if (idiag_movesect .ne. 70) goto 3800
if (nnewsave(2,n) .eq. 0) then
if (nnewsave(1,n) .eq. 0) then
dumch4 = 'NO X'
else if (nnewsave(1,n) .eq. n) then
dumch4 = 'NO A'
else
dumch4 = 'YESA'
end if
else if (nnewsave(1,n) .eq. 0) then
if (nnewsave(2,n) .eq. n) then
dumch4 = 'NO B'
else
dumch4 = 'YESB'
end if
else if (nnewsave(2,n) .eq. n) then
if (nnewsave(1,n) .eq. n) then
dumch4 = 'NO Y'
else
dumch4 = 'YESC'
end if
else if (nnewsave(1,n) .eq. n) then
dumch4 = 'YESD'
else
dumch4 = 'YESE'
end if
dumch1 = '+'
if (drymasspp(n) .gt. drymassxx(n)) dumch1 = '-'
msg = ' '
call peg_message( lunout, msg )
write(msg,97010) dumch1, dumch4, iclm, jclm, k, m, &
n, nnewsave(1,n), nnewsave(2,n), iforce_movecenter(n)
call peg_message( lunout, msg )
write(msg,97020) 'pre mass, dens ', &
drymasspp(n), drydenspp(n)
call peg_message( lunout, msg )
write(msg,97020) 'aft mass, dens, numb', &
drymassxx(n), drydensxx(n), rnumbxx(n)
call peg_message( lunout, msg )
if ((drydensxx(n) .ne. drydensxx0(n)) .or. &
(drymassxx(n) .ne. drymassxx0(n)) .or. &
(rnumbxx(n) .ne. rnumbxx0(n) )) then
write(msg,97020) 'aft0 mas, dens, numb', &
drymassxx0(n), drydensxx0(n), rnumbxx0(n)
call peg_message( lunout, msg )
end if
write(msg,97020) 'vlop0, vbarp, vhip0', &
volumlo_sect(n,itype), dumvbar_pre, volumhi_sect(n,itype)
call peg_message( lunout, msg )
write(msg,97020) 'vlop , vbarp, vhip ', &
dumvlo_pre, dumvbar_pre, dumvhi_pre
call peg_message( lunout, msg )
write(msg,97020) 'vloax, vbarax, vhiax', &
dumvcutlo_nnew_pre, dumvbar_pre, dumvcuthi_nnew_pre
call peg_message( lunout, msg )
write(msg,97020) 'vloa0, vbara, vhia0', &
volumlo_sect(nnew,itype), dumvbar_aft, volumhi_sect(nnew,itype)
call peg_message( lunout, msg )
write(msg,97020) 'dumfrvol, num, ratio', &
dumfracvol, dumfracnum, dumratio
call peg_message( lunout, msg )
write(msg,97020) 'frvol,num1; vol,num2', &
xferfracvol(1,n), xferfracnum(1,n), &
xferfracvol(2,n), xferfracnum(2,n)
call peg_message( lunout, msg )
97010 format( 'movesect', 2a, 7x, 4i3, 4x, &
'n,nnews', 3i3, 4x, 'iforce', i3.2 )
97020 format( a, 1p, 4e13.4 )
3800 continue
!
! check for legal combinations of nnewsave(1,n) & nnewsave(2,n)
! error if
! nnew1 == nnew2
! both are non-zero AND iabs(nnew1-nnew2) != 1
ierr = 0
if (nnewsave(1,n) .eq. nnewsave(2,n)) then
ierr = 1
else if (nnewsave(1,n)*nnewsave(2,n) .ne. 0) then
if (iabs(nnewsave(1,n)-nnewsave(2,n)) .ne. 1) ierr = 1
end if
if (ierr .gt. 0) then
write(msg,97010) 'E', 'RROR', iclm, jclm, k, m, &
n, nnewsave(1,n), nnewsave(2,n), iforce_movecenter(n)
call peg_message( lunout, msg )
end if
! if method_movesect == 30-31 then force moving center
! this is just for testing purposes
if ((method_movesect .ge. 30) .and. (method_movesect .le. 39)) then
nnewsave(1,n) = nnew
nnewsave(2,n) = 0
xferfracvol(1,n) = 1.0
xferfracvol(2,n) = 0.0
xferfracnum(1,n) = 1.0
xferfracnum(2,n) = 0.0
end if
3900 continue
return
end subroutine move_sections_calc_masnumadvect
!-----------------------------------------------------------------------
subroutine move_sections_apply_moves( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drydenspp, drymasspp, rnumbpp, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
xferfracvol, xferfracnum )
!
! routine performs the actual transfer of aerosol number and mass
! between sections
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag, iclm, jclm, iphase, itype, k, &
m, method_movesect, idiag_movesect, llhysw, llwater
integer nnewsave(2,maxd_asize)
real densdefault, densh2o, smallmassaa, smallmassbb
real delta_water_conform1, delta_numb_conform1
real drydenspp(maxd_asize)
real drymasspp(maxd_asize)
real drymassxx(maxd_asize), drymassyy(maxd_asize)
real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize), &
rmassyy(maxd_acomp+2,maxd_asize)
real rnumbpp(maxd_asize)
real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
! local variables
integer jj, l, ll, n, ndum, nnew, nold
integer jja, jjb, jjc
real delta_numb_conform2, &
dumbot, dumnum, dumnum_at_dhi, dumnum_at_dlo, &
dumvol, dumvol1p, dumxfvol, dumxfnum, sixoverpi, third
real dumpp(maxd_asize), dumxx(maxd_asize), dumyy(maxd_asize), &
dumout(maxd_asize)
character*160 msg
character*8 dumch8
character*4 dumch4
sixoverpi = 6.0/pi
third = 1.0/3.0
!
! initialize "yy" arrays that hold values after inter-mode transferring
! "yy" = "xx" for sections that do not move at all
! "yy" = 0.0 for sections that do move (partially or completely)
!
do 1900 n = 1, nsize_aer(itype)
if ( (nnewsave(1,n) .eq. n) .and. &
(nnewsave(2,n) .eq. 0) ) then
! if nnew == n, then material in section n will not be transferred, and
! section n will contain its initial material plus any material
! transferred from other sections
! so initialize "yy" arrays with "xx" values
drymassyy(n) = drymassxx(n)
dryvolyy(n) = dryvolxx(n)
wetmassyy(n) = wetmassxx(n)
wetvolyy(n) = wetvolxx(n)
rnumbyy(n) = rnumbxx(n)
do ll = 1, ncomp_plustracer_aer(itype) + 2
rmassyy(ll,n) = rmassxx(ll,n)
end do
else
! if nnew .ne. n, then material in section n will be transferred, and
! section n will only contain material that is transferred from
! other sections
! so initialize "yy" arrays to zero
drymassyy(n) = 0.0
dryvolyy(n) = 0.0
wetmassyy(n) = 0.0
wetvolyy(n) = 0.0
rnumbyy(n) = 0.0
do ll = 1, ncomp_plustracer_aer(itype) + 2
rmassyy(ll,n) = 0.0
end do
end if
1900 continue
!
! do the transfer of mass and number
!
do 2900 nold = 1, nsize_aer(itype)
if ( (nnewsave(1,nold) .eq. nold) .and. &
(nnewsave(2,nold) .eq. 0 ) ) goto 2900
do 2800 jj = 1, 2
nnew = nnewsave(jj,nold)
if (nnew .le. 0) goto 2800
dumxfvol = xferfracvol(jj,nold)
dumxfnum = xferfracnum(jj,nold)
do ll = 1, ncomp_plustracer_aer(itype) + 2
rmassyy(ll,nnew) = rmassyy(ll,nnew) + rmassxx(ll,nold)*dumxfvol
end do
rnumbyy(nnew) = rnumbyy(nnew) + rnumbxx(nold)*dumxfnum
drymassyy(nnew) = drymassyy(nnew) + drymassxx(nold)*dumxfvol
dryvolyy(nnew) = dryvolyy(nnew) + dryvolxx(nold)*dumxfvol
wetmassyy(nnew) = wetmassyy(nnew) + wetmassxx(nold)*dumxfvol
wetvolyy(nnew) = wetvolyy(nnew) + wetvolxx(nold)*dumxfvol
2800 continue
2900 continue
!
! transfer among sections is completed
! - check for conservation of mass/volume/number
! - conform number again
! - compute/store densities and mean sizes
! - if k=1, save values for use by dry deposition routine
! - copy new mixrats back to rsub array
!
call move_sections_conserve_check( &
1, iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
delta_numb_conform2 = 0.0
do 3900 n = 1, nsize_aer(itype)
dumvol = dryvolyy(n)
if (min(drymassyy(n),dumvol) .le. smallmassbb) then
dumvol = drymassyy(n)/densdefault
dumnum = drymassyy(n)/(volumlo_sect(n,itype)*densdefault)
delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n)
rnumbyy(n) = dumnum
adrydens_sub(n,itype,k,m) = densdefault
awetdens_sub(n,itype,k,m) = densdefault
admeandry_sub(n,itype,k,m) = dcen_sect(n,itype)
admeanwet_sub(n,itype,k,m) = dcen_sect(n,itype)
else
dumnum = rnumbyy(n)
dumnum_at_dhi = dumvol/volumhi_sect(n,itype)
dumnum_at_dlo = dumvol/volumlo_sect(n,itype)
dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n)
rnumbyy(n) = dumnum
adrydens_sub(n,itype,k,m) = drymassyy(n)/dumvol
dumvol1p = dumvol/dumnum
admeandry_sub(n,itype,k,m) = (dumvol1p*sixoverpi)**third
awetdens_sub(n,itype,k,m) = wetmassyy(n)/wetvolyy(n)
dumvol1p = wetvolyy(n)/dumnum
admeanwet_sub(n,itype,k,m) = min( 100.*dcen_sect(n,itype), &
(dumvol1p*sixoverpi)**third )
end if
aqvoldry_sub(n,itype,k,m) = dumvol
aqmassdry_sub(n,itype,k,m) = drymassyy(n)
! if (k .eq. 1) then
! awetdens_sfc(n,itype,iclm,jclm) = awetdens_sub(n,itype,k,m)
! admeanwet_sfc(n,itype,iclm,jclm) = admeanwet_sub(n,itype,k,m)
! end if
do ll = 1, ncomp_plustracer_aer(itype)
l = massptr_aer(ll,n,itype,iphase)
rsub(l,k,m) = rmassyy(ll,n)
end do
l = 0
if (iphase .eq. ai_phase) then
l = waterptr_aer(n,itype)
if (l .gt. 0) rsub(l,k,m) = rmassyy(llwater,n)
l = hyswptr_aer(n,itype)
if (l .gt. 0) rsub(l,k,m) = rmassyy(llhysw,n)
end if
rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
3900 continue
delta_numb_conform1 = delta_numb_conform1 + delta_numb_conform2
call move_sections_conserve_check( &
2, iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
! diagnostic output
! output nnewsave values when msectional=7xxx
if (idiag_movesect .ne. 70) goto 4900
ndum = 0
do n = 1, nsize_aer(itype)
if (nnewsave(1,n)+nnewsave(2,n) .ne. n) ndum = ndum + 1
end do
! if (ndum .gt. 0) then
! write(msg,97010) 'SOME', iclm, jclm, k, m, &
! ndum, (nnewsave(1,n), nnewsave(2,n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
! else
! write(msg,97010) 'NONE', iclm, jclm, k, m, &
! ndum, (nnewsave(1,n), nnewsave(2,n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
! end if
dumch4 = 'NONE'
if (ndum .gt. 0) dumch4 = 'SOME'
msg = ' '
call peg_message( lunout, msg )
write(msg,97010) dumch4, iclm, jclm, k, m, ndum
call peg_message( lunout, msg )
do jjb = 1, nsize_aer(itype), 10
jjc = min( jjb+9, nsize_aer(itype) )
write(msg,97011) (nnewsave(1,n), nnewsave(2,n), n=jjb,jjc)
call peg_message( lunout, msg )
end do
! write(msg,97020) 'rnumbold', (rnumbxx(n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
! write(msg,97020) 'rnumbnew', (rnumbyy(n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
! write(msg,97020) 'drvolold', (dryvolxx(n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
! write(msg,97020) 'drvolnew', (dryvolyy(n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
dumbot = log( volumhi_sect(1,itype)/volumlo_sect(1,itype) )
do n = 1, nsize_aer(itype)
dumpp(n) = -9.99
dumxx(n) = -9.99
dumyy(n) = -9.99
if ( (drydenspp(n) .gt. 0.5) .and. &
(drymasspp(n) .gt. smallmassaa) ) then
dumvol = drymasspp(n)/drydenspp(n)
if ((rnumbpp(n) .ge. 1.0e-35) .and. &
(dumvol .ge. 1.0e-35)) then
dumvol1p = dumvol/rnumbpp(n)
dumpp(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
end if
end if
if ((rnumbxx(n) .ge. 1.0e-35) .and. &
(dryvolxx(n) .ge. 1.0e-35)) then
dumvol1p = dryvolxx(n)/rnumbxx(n)
dumxx(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
end if
if ((rnumbyy(n) .ge. 1.0e-35) .and. &
(dryvolyy(n) .ge. 1.0e-35)) then
dumvol1p = dryvolyy(n)/rnumbyy(n)
dumyy(n) = 1.0 + log(dumvol1p/volumlo_sect(1,itype))/dumbot
end if
end do
! write(msg,97030) 'lnvolold', (dumxx(n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
! write(msg,97030) 'lnvolnew', (dumyy(n), n=1,nsize_aer(itype))
! call peg_message( lunout, msg )
do jja = 1, 7
if (jja .eq. 1) then
dumch8 = 'rnumbold'
dumout(:) = rnumbxx(:)
else if (jja .eq. 2) then
dumch8 = 'rnumbnew'
dumout(:) = rnumbyy(:)
else if (jja .eq. 3) then
dumch8 = 'drvolold'
dumout(:) = dryvolxx(:)
else if (jja .eq. 4) then
dumch8 = 'drvolnew'
dumout(:) = dryvolyy(:)
else if (jja .eq. 5) then
dumch8 = 'lnvolold'
dumout(:) = dumxx(:)
else if (jja .eq. 6) then
dumch8 = 'lnvolnew'
dumout(:) = dumyy(:)
else if (jja .eq. 7) then
dumch8 = 'lnvolpre'
dumout(:) = dumpp(:)
end if
do jjb = 1, nsize_aer(itype), 10
jjc = min( jjb+9, nsize_aer(itype) )
if (jja .le. 4) then
write(msg,97020) dumch8, (dumout(n), n=jjb,jjc)
else
write(msg,97030) dumch8, (dumout(n), n=jjb,jjc)
end if
call peg_message( lunout, msg )
dumch8 = ' '
end do
end do
!97010 format( / 'movesectapply', a, 4i3, 3x, i3 / 5x, 10(3x,2i3) )
!97020 format( a, 1p, 10e9.1 / (( 8x, 1p, 10e9.1 )) )
!97030 format( a, 10f9.3 / (( 8x, 10f9.3 )) )
97010 format( 'movesectapply', a, 4i3, 3x, i3 )
97011 format( 5x, 10(3x,2i3) )
97020 format( a, 1p, 10e9.1 )
97030 format( a, 10f9.3 )
4900 continue
return
end subroutine move_sections_apply_moves
!-----------------------------------------------------------------------
! rce 08-nov-2004 - this routine is new (and perhaps temporary)
!
subroutine move_sections_apply_n1_inflow( &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy, &
xferfracvol, xferfracnum, &
specmwxx, specdensxx )
!
! routine applies an ad_hoc correction to bin 1 to account for
! growth of smaller particles (which are not simulated) growing into
! bin 1
! the correction to particle number balances the loss of particles from
! bin 1 to larger bins by growth
! the correction to mass assumes that the particles coming into bin 1
! are slightly larger than dlo_sect(1)
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag, iclm, jclm, iphase, itype, k, &
m, method_movesect, idiag_movesect, llhysw, llwater
integer nnewsave(2,maxd_asize)
real densdefault, densh2o, smallmassaa, smallmassbb
real delta_water_conform1, delta_numb_conform1
real drymassxx(maxd_asize), drymassyy(maxd_asize)
real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize), &
rmassyy(maxd_acomp+2,maxd_asize)
real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
real xferfracvol(2,maxd_asize), xferfracnum(2,maxd_asize)
real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
real specdensxx(maxd_acomp), specmwxx(maxd_acomp)
! local variables
integer jj, l, ll, n, nnew, nold
real deltanum, deltavol, dumvol1p, dumxfvol, dumxfnum
!
! compute fraction of number transferred out of bin 1 by growth
!
nold = 1
n = nold
if ( (nnewsave(1,nold) .eq. nold) .and. &
(nnewsave(2,nold) .eq. 0 ) ) goto 3900
dumxfnum = 0.0
do 2800 jj = 1, 2
nnew = nnewsave(jj,nold)
if (nnew .le. 0) goto 2800
if (nnew .eq. nold) goto 2800
dumxfnum = dumxfnum + xferfracnum(jj,nold)
2800 continue
!
! compute "inflow" change to number and volume
! number change matches that lost by growth
! volume change assume inflow particles slightly bigger then dlo_sect
!
dumvol1p = 0.95*volumlo_sect(n,itype) + 0.05*volumhi_sect(n,itype)
deltanum = dumxfnum*rnumbxx(n)
deltavol = deltanum*dumvol1p
if (dumxfnum .le. 0.0) goto 3900
if (deltanum .le. 0.0) goto 3900
if (deltavol .le. 0.0) goto 3900
!
! increment the number and masses
! if the old dryvol (dryvolxx) > smallmassbb, then compute mass increment for
! each species from the old masses (rmassxx)
! otherwise only increment the first mass species
!
if (dryvolxx(n) .gt. smallmassbb) then
dumxfvol = deltavol/dryvolxx(n)
do ll = 1, ncomp_plustracer_aer(itype) + 2
rmassyy(ll,n) = rmassyy(ll,n) + dumxfvol*rmassxx(ll,n)
end do
else
ll = 1
rmassyy(ll,n) = rmassyy(ll,n) + deltavol*specdensxx(ll)/specmwxx(ll)
end if
rnumbyy(n) = rnumbyy(n) + deltanum
!
! transfer results into rsub
!
do ll = 1, ncomp_plustracer_aer(itype)
l = massptr_aer(ll,n,itype,iphase)
rsub(l,k,m) = rmassyy(ll,n)
end do
if (iphase .eq. ai_phase) then
l = waterptr_aer(n,itype)
if (l .gt. 0) rsub(l,k,m) = rmassyy(llwater,n)
l = hyswptr_aer(n,itype)
if (l .gt. 0) rsub(l,k,m) = rmassyy(llhysw,n)
end if
rsub(numptr_aer(n,itype,iphase),k,m) = rnumbyy(n)
3900 continue
return
end subroutine move_sections_apply_n1_inflow
!-----------------------------------------------------------------------
subroutine move_sections_conserve_check( ipass, &
iflag, iclm, jclm, k, m, iphase, itype, &
method_movesect, idiag_movesect, llhysw, llwater, nnewsave, &
densdefault, densh2o, smallmassaa, smallmassbb, &
delta_water_conform1, delta_numb_conform1, &
drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx, &
drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
!
! routine checks for conservation of number, mass, and volume
! by the move_sections algorithm
!
! ipass = 1
! initialize all thesum(jj,ll) to zero
! computes thesum(1,ll) from rmassxx, rnumbxx, ...
! computes thesum(2,ll) from rmassyy, rnumbyy, ...
! compares thesum(1,ll) with thesum(2,ll)
! computes thesum(3,ll) from rsub before section movement
! ipass = 2
! computes thesum(4,ll) from rsub after section movement
! compares thesum(3,ll) with thesum(4,ll)
!
! currently only implemented for condensational growth (iflag=1)
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer ipass, iflag, iclm, jclm, iphase, itype, k, &
m, method_movesect, idiag_movesect, llhysw, llwater
integer nnewsave(2,maxd_asize)
real densdefault, densh2o, smallmassaa, smallmassbb
real delta_water_conform1, delta_numb_conform1
real drymassxx(maxd_asize), drymassyy(maxd_asize)
real dryvolxx(maxd_asize), dryvolyy(maxd_asize)
real rmassxx(maxd_acomp+2,maxd_asize), &
rmassyy(maxd_acomp+2,maxd_asize)
real rnumbxx(maxd_asize), rnumbyy(maxd_asize)
real wetvolxx(maxd_asize), wetvolyy(maxd_asize)
real wetmassxx(maxd_asize), wetmassyy(maxd_asize)
! local variables
integer jj, l, ll, llworst, llworstb, n
integer nerr, nerrmax
save nerr, nerrmax
data nerr, nerrmax / 0, 999 /
real dumbot, dumtop, dumtoler, dumerr, dumworst, dumworstb
real duma, dumb, dumc, dume
real thesum(4,maxd_acomp+7)
save thesum
character*8 dumname(maxd_acomp+7)
character*160 msg
if (ipass .eq. 2) goto 2000
do ll = 1, ncomp_plustracer_aer(itype)+7
do jj = 1, 4
thesum(jj,ll) = 0.0
end do
end do
do n = 1, nsize_aer(itype)
do ll = 1, ncomp_plustracer_aer(itype)+2
thesum(1,ll) = thesum(1,ll) + rmassxx(ll,n)
thesum(2,ll) = thesum(2,ll) + rmassyy(ll,n)
end do
ll = ncomp_plustracer_aer(itype)+3
thesum(1,ll) = thesum(1,ll) + rnumbxx(n)
thesum(2,ll) = thesum(2,ll) + rnumbyy(n)
ll = ncomp_plustracer_aer(itype)+4
thesum(1,ll) = thesum(1,ll) + drymassxx(n)
thesum(2,ll) = thesum(2,ll) + drymassyy(n)
ll = ncomp_plustracer_aer(itype)+5
thesum(1,ll) = thesum(1,ll) + dryvolxx(n)
thesum(2,ll) = thesum(2,ll) + dryvolyy(n)
ll = ncomp_plustracer_aer(itype)+6
thesum(1,ll) = thesum(1,ll) + wetmassxx(n)
thesum(2,ll) = thesum(2,ll) + wetmassyy(n)
ll = ncomp_plustracer_aer(itype)+7
thesum(1,ll) = thesum(1,ll) + wetvolxx(n)
thesum(2,ll) = thesum(2,ll) + wetvolyy(n)
end do
2000 continue
!
! calc sum over bins for each species
! for water, account for loss in initial conform (delta_water_conform1)
!
do n = 1, nsize_aer(itype)
do ll = 1, ncomp_plustracer_aer(itype)+3
if (ll .le. ncomp_plustracer_aer(itype)) then
l = massptr_aer(ll,n,itype,iphase)
else if (ll .eq. ncomp_plustracer_aer(itype)+1) then
l = 0
if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
else if (ll .eq. ncomp_plustracer_aer(itype)+2) then
l = 0
if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
else
l = numptr_aer(n,itype,iphase)
end if
if (l .gt. 0) &
thesum(ipass+2,ll) = thesum(ipass+2,ll) + rsub(l,k,m)
end do
end do
if (ipass .eq. 2) then
ll = ncomp_plustracer_aer(itype)+2
thesum(3,ll) = thesum(3,ll) + delta_water_conform1
ll = ncomp_plustracer_aer(itype)+3
thesum(3,ll) = thesum(3,ll) + delta_numb_conform1
end if
!
! now compare either sum1-sum2 or sum3-sum4
! on ipass=1, jj=1, so compare sum1 & sum2
! on ipass=2, jj=3, so compare sum3 & sum4
!
do ll = 1, ncomp_plustracer_aer(itype)+7
dumname(ll) = ' '
write(dumname(ll),'(i4.4)') ll
if (ll .le. ncomp_plustracer_aer(itype)) dumname(ll) = &
name(massptr_aer(ll,1,itype,iphase))(1:4)
if (ll .eq. ncomp_plustracer_aer(itype)+1) dumname(ll) = 'hysw'
if (ll .eq. ncomp_plustracer_aer(itype)+2) dumname(ll) = 'watr'
if (ll .eq. ncomp_plustracer_aer(itype)+3) dumname(ll) = 'numb'
if (ll .eq. ncomp_plustracer_aer(itype)+4) dumname(ll) = 'drymass'
if (ll .eq. ncomp_plustracer_aer(itype)+5) dumname(ll) = 'dryvol'
if (ll .eq. ncomp_plustracer_aer(itype)+6) dumname(ll) = 'wetmass'
if (ll .eq. ncomp_plustracer_aer(itype)+7) dumname(ll) = 'wetvol'
end do
jj = 2*ipass - 1
dumworst = 0.0
dumworstb = 0.0
llworst = 0
llworstb = 0
do ll = 1, ncomp_plustracer_aer(itype)+7
dumtop = thesum(jj+1,ll) - thesum(jj,ll)
! dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35 )
! change minimum dumbot to 1.0e-30 to reduce unimportant error messages
dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-30 )
dumerr = dumtop/dumbot
! rce 21-jul-2006 - encountered some cases when delta_*_conform1 is negative
! and large in magnitude relative to thesum, which causes the mass
! conservation to be less accurate due to roundoff
! following section recomputes relative error with delta_*_conform1
! added onto each of thesum
! also increased dumtoler slightly
if (ipass .eq. 2) then
dumc = 1.0
if (ll .eq. ncomp_plustracer_aer(itype)+2) then
dumc = delta_water_conform1
else if (ll .eq. ncomp_plustracer_aer(itype)+3) then
dumc = delta_numb_conform1
end if
if (dumc .lt. 0.0) then
duma = thesum(3,ll) - dumc
dumb = thesum(4,ll) - dumc
dumtop = dumb - duma
! dumbot = max( abs(duma), abs(dumb), 1.0e-35 )
dumbot = max( abs(duma), abs(dumb), 1.0e-30 )
dume = dumtop/dumbot
if (abs(dume) .lt. abs(dumerr)) dumerr = dume
end if
end if
if (abs(dumerr) .gt. abs(dumworst)) then
llworstb = llworst
dumworstb = dumworst
llworst = ll
dumworst = dumerr
end if
end do
dumtoler = 1.0e-6
if (abs(dumworst) .gt. dumtoler) then
nerr = nerr + 1
if (nerr .le. nerrmax) then
msg = ' '
call peg_message( lunout, msg )
write(msg,97110) iclm, jclm, k, m, ipass, llworst
call peg_message( lunout, msg )
write(msg,97120) ' nnew(1,n)', &
(nnewsave(1,n), n=1,nsize_aer(itype))
call peg_message( lunout, msg )
write(msg,97120) ' nnew(2,n)', &
(nnewsave(2,n), n=1,nsize_aer(itype))
call peg_message( lunout, msg )
ll = llworst
if (ll .eq. 0) ll = ncomp_plustracer_aer(itype)+2
write(msg,97130) 'name/relerr/thesum', jj, '/thesum', jj+1, &
dumname(ll), dumworst, thesum(jj,ll), thesum(jj+1,ll)
call peg_message( lunout, msg )
if ( (ll .eq. ncomp_plustracer_aer(itype)+3) .and. &
(abs(dumworstb) .gt. dumtoler) ) then
ll = max( 1, llworstb )
dumtop = thesum(jj+1,ll) - thesum(jj,ll)
! dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35 )
dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-30 )
dumerr = dumtop/dumbot
write(msg,97130) 'name/relerr/thesum', jj, '/thesum', jj+1, &
dumname(ll), dumerr, thesum(jj,ll), thesum(jj+1,ll)
call peg_message( lunout, msg )
end if
end if
end if
97110 format( 'movesect conserve ERROR - i/j/k/m/pass/llworst', &
4i3, 2x, 2i3 )
97120 format( a, 64i3 )
97130 format( a, i1, a, i1, 2x, a, 1p, 3e16.7 )
return
end subroutine move_sections_conserve_check
!-----------------------------------------------------------------------
subroutine test_move_sections( iflag_test, iclm, jclm, k, m )
!
! routine runs tests on move_sections, using a matrix of
! pregrow and aftgrow masses for each section
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr arguments
integer iflag_test, iclm, jclm, k, m
! local variables
integer idiag_movesect, iflag, ii, iphase, itype, jj, &
l, ll, n, nn
integer ientryno
save ientryno
data ientryno / 0 /
real dumnumb, dumvolpre, dumvolaft
real dumsv_rsub(l2maxd)
real dumsv_drymass_pregrow(maxd_asize,maxd_atype)
real dumsv_drymass_aftgrow(maxd_asize,maxd_atype)
real dumsv_drydens_pregrow(maxd_asize,maxd_atype)
real dumsv_drydens_aftgrow(maxd_asize,maxd_atype)
character*160 msg
integer maxvolfactpre, maxvolfactaft
parameter (maxvolfactpre=15, maxvolfactaft=23)
real dumvolfactpre(maxvolfactpre)
data dumvolfactpre / &
2.0, 0.0, 1.0e-20, 0.5, 0.9, &
1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0, &
8.1, 16.0 /
real dumvolfactaft(maxvolfactaft)
data dumvolfactaft / &
4.0, 0.0, 1.0e-20, 0.01, 0.02, 0.05, 0.1, 0.5, 0.9, &
1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0, &
8.1, 16.0, 32., 64., 128., 256. /
!
! check for valid inputs
! and first entry
!
if (msectional .le. 0) return
if (nsize_aer(1) .le. 0) return
idiag_movesect = mod( msectional, 10000 )/100
if (idiag_movesect .ne. 70) return
ientryno = ientryno + 1
if (ientryno .gt. 2) return
!
! save rsub and drymass/dens_aft/pregrow
!
do l = 1, ltot2
dumsv_rsub(l) = rsub(l,k,m)
end do
do itype = 1, ntype_aer
do n = 1, nsize_aer(itype)
dumsv_drymass_pregrow(n,itype) = drymass_pregrow(n,itype)
dumsv_drymass_aftgrow(n,itype) = drymass_aftgrow(n,itype)
dumsv_drydens_pregrow(n,itype) = drydens_pregrow(n,itype)
dumsv_drydens_aftgrow(n,itype) = drydens_aftgrow(n,itype)
end do
end do
!
! make test calls to move_sections
!
do 3900 iflag = 1, 1
iphase = ai_phase
if (iabs(iflag) .eq. 2) iphase = cw_phase
do 3800 itype = 1, ntype_aer
do 2900 nn = 1, nsize_aer(itype)
do 2800 ii = 1, maxvolfactpre
do 2700 jj = 1, maxvolfactaft
! zero out rsub and dryxxx_yyygrow arrays
do n = 1, nsize_aer(itype)
do ll = 1, ncomp_plustracer_aer(itype)
rsub(massptr_aer(ll,n,itype,iphase),k,m) = 0.0
end do
l = 0
if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
if (l .gt. 0) rsub(l,k,m) = 0.0
l = numptr_aer(n,itype,iphase)
if (l .gt. 0) rsub(l,k,m) = 0.0
drymass_pregrow(n,itype) = 0.0
drymass_aftgrow(n,itype) = 0.0
drydens_pregrow(n,itype) = -1.0
drydens_aftgrow(n,itype) = -1.0
end do
! fill in values for section nn
n = nn
dumnumb = 1.0e7
rsub(numptr_aer(n,itype,iphase),k,m) = dumnumb
ll = 1
l = massptr_aer(ll,n,itype,iphase)
dumvolpre = volumlo_sect(n,itype)*dumvolfactpre(ii)*dumnumb
drydens_pregrow(n,itype) = dens_aer(ll,itype)
drymass_pregrow(n,itype) = dumvolpre*drydens_pregrow(n,itype)
if (ii .eq. 1) drydens_pregrow(n,itype) = -1.0
dumvolaft = volumlo_sect(n,itype)*dumvolfactaft(jj)*dumnumb
drydens_aftgrow(n,itype) = dens_aer(ll,itype)
drymass_aftgrow(n,itype) = dumvolaft*drydens_aftgrow(n,itype)
if (jj .eq. 1) drydens_aftgrow(n,itype) = -1.0
rsub(l,k,m) = drymass_aftgrow(n,itype)/mw_aer(ll,itype)
msg = ' '
call peg_message( lunout, msg )
write(msg,98010) nn, ii, jj
call peg_message( lunout, msg )
write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
call peg_message( lunout, msg )
! make test call to move_sections
call move_sections( iflag, iclm, jclm, k, m )
msg = ' '
call peg_message( lunout, msg )
write(msg,98010) nn, ii, jj
call peg_message( lunout, msg )
write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
call peg_message( lunout, msg )
2700 continue
2800 continue
2900 continue
3800 continue
3900 continue
98010 format( 'test_move_sections output - nn, ii, jj =', 3i3 )
98011 format( 'volfactpre, volfactaft =', 1p, 2e12.4 )
!
! restore rsub and drymass/dens_aft/pregrow
!
do l = 1, ltot2
rsub(l,k,m) = dumsv_rsub(l)
end do
do itype = 1, ntype_aer
do n = 1, nsize_aer(itype)
drymass_pregrow(n,itype) = dumsv_drymass_pregrow(n,itype)
drymass_aftgrow(n,itype) = dumsv_drymass_aftgrow(n,itype)
drydens_pregrow(n,itype) = dumsv_drydens_pregrow(n,itype)
drydens_aftgrow(n,itype) = dumsv_drydens_aftgrow(n,itype)
end do
end do
return
end subroutine test_move_sections
!-----------------------------------------------------------------------
subroutine move_sections_checkptrs( iflag, iclm, jclm, k, m )
!
! checks for valid number and water pointers
!
implicit none
! include 'v33com'
! include 'v33com2'
! include 'v33com3'
! include 'v33com9a'
! include 'v33com9b'
! subr parameters
integer iflag, iclm, jclm, k, m
! local variables
integer l, itype, iphase, n, ndum
character*160 msg
do 1900 itype = 1, ntype_aer
do 1800 iphase = 1, nphase_aer
ndum = 0
do n = 1, nsize_aer(itype)
l = numptr_aer(n,itype,iphase)
if ((l .le. 0) .or. (l .gt. ltot2)) then
msg = '*** subr move_sections error - ' // &
'numptr_amode not defined'
call peg_message( lunerr, msg )
write(msg,9030) 'mode, numptr =', n, l
call peg_message( lunerr, msg )
write(msg,9030) 'iphase, itype =', iphase, itype
call peg_message( lunerr, msg )
call peg_error_fatal( lunerr, msg )
end if
! checks involving nspec_amode and nspec_amode_nontracer
! being the same for all sections are no longer needed
l = 0
if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
if ((l .gt. 0) .and. (l .le. ltot2)) ndum = ndum + 1
end do
if ((ndum .ne. 0) .and. (ndum .ne. nsize_aer(itype))) then
msg = '*** subr move_sections error - ' // &
'waterptr_aer must be on/off for all modes'
call peg_message( lunerr, msg )
write(msg,9030) 'iphase, itype =', iphase, itype
call peg_message( lunerr, msg )
call peg_error_fatal( lunerr, msg )
end if
9030 format( a, 2(1x,i6) )
1800 continue
1900 continue
return
end subroutine move_sections_checkptrs
end module module_mosaic_movesect