chem/module_mosaic_movesect3d.F

   1         module module_mosaic_movesect3d
   2
   3
   4         use module_data_mosaic_kind, only:  r8
   5         use module_data_mosaic_main, only:  mw_air, ntot_used, pi
   6         use module_data_mosaic_asecthp, only:  lunerr, lunout
   7         use module_peg_util
   8
   9
  10         implicit none
  11
  12
  13         contains
  14
  15
  16 !-----------------------------------------------------------------------
  17 !**********************************************************************************
  18 ! following routines adapted from
  19 !    fjaersky:/home/d37080/box/aqchem/pandis/ccboxwrf6/module_mosaic_movesect.F.saveaa
  20 !**********************************************************************************
  21
  22
  23 !-----------------------------------------------------------------------
  24         subroutine move_sect_3d_x1( iphase_flag, iclm, jclm, k, m, rbox,   &
  25           fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam,   &
  26           drydens_aftgrow, drydens_pregrow,   &
  27           drymass_aftgrow, drymass_pregrow,   &
  28           adrydens_tmp, awetdens_tmp, adrydpav_tmp, awetdpav_tmp,   &
  29           adryqmas_tmp, it_mosaic, mmovesect_flag1, idiag_movesect )
  30 !
  31 !   routine transfers aerosol number and mass between sections
  32 !       to account for continuous aerosol growth
  33 !   this routine is called after the gas condensation module (MOSAIC) or
  34 !       aqueous chemistry module has increased the mass within sections
  35 !
  36 !   moving-center algorithm or mass-number advection algorithm is used,
  37 !   depending on value of mod(mmovesect_flag1,100)
  38 !       section mean diameter is given by
  39 !           vtot = ntot * (pi/6) * (dmean**3)
  40 !       where vtot and ntot are total dry-volume and number for the section
  41 !       if dmean is outside the section boundaries (dlo_sect & dhi_sect), then
  42 !           all the mass and number in the section are transfered to the
  43 !           section with dlo_sect(nnew) < dmean < dhi_sect(nnew)
  44 !
  45 !   mass mixing ratios are in rbox(massptr_aer(ll,n,itype,iphase))
  46 !       units such that rbox*fact_apmassmr gives (g-AP/g-air)
  47 !   number mixing ratios are in rbox(numptr_aer(n,itype,iphase))
  48 !       units such that rbox*fact_apnumbmr gives (#/g-air)
  49 !   these values are over-written with new values
  50
  51 !   the following are also updated:
  52 !       adrydens_tmp(n,itype), adrydpav_tmp(n,itype),
  53 !       awetdens_tmp(n,itype), awetdpav_tmp(n,itype)
  54 !       adryqmas_tmp(n,itype)
  55 !
  56 !   the module_data_mosaic_asecthp densities are such that
  57 !       dens_aer*fact_apdens give (g/cm^3)
  58 !   the module_data_mosaic_asecthp diameters are such that
  59 !       dlo/hi_sect*fact_apdiam give (cm)
  60 !   the module_data_mosaic_asecthp volumes are such that
  61 !       volumlo/hi_sect**fact_apdiam**3) give (cm^3)
  62 !
  63 !   input parameters
  64 !       iphase_flag = 1 - do transfer after trace-gas condensation/evaporation
  65 !             = 2 - do transfer after aqueous chemistry
  66 !             = -1/-2 - do some "first entry" tasks for the iphase_flag=+1/+2 cases
  67 !
  68 !       iclm, jclm, k = current i,j,k indices
  69 !       m = current subarea index
  70 !
  71 !       drymass_pregrow(n,itype) = dry-mass for section n before the growth
  72 !       drymass_aftgrow(n,itype) = dry-mass for section n after the growth
  73 !                                  but before inter-section transfer
  74 !           (units for both are same as rbox mass entries)
  75 !       drydens_pregrow(n,itype) = dry-density for section n before the growth
  76 !       drydens_aftgrow(n,itype) = dry-density for section n after the growth
  77 !                                  but before inter-section transfer
  78 !           (units for both are same as dens_aer)
  79 !
  80 !       (drymass_pregrow and drydens_pregrow are used by the linear-discrete
  81 !           algorithm but not the moving-center algorithm)
  82 !
  83 !       method_movesect (from first two digits of movesect_flag1, which
  84 !           is a "data_module" variable)
  85 !           50-59 - do 3d moving-center algorithm
  86 !           60-69 - do 3d linear-discrete algorithm
  87 !
  88         use module_data_mosaic_asecthp, only:  &
  89             maxd_acomp, maxd_asize, maxd_atype,   &
  90             ntype_aer, nsize_aer, nphase_aer, ai_phase, cw_phase,   &
  91             dens_aer, dcen_sect,   &
  92             smallmassaa_asecthp => smallmassaa, smallmassbb_asecthp => smallmassbb,   &
  93             volumcen_sect, volumlo_sect, volumhi_sect
  94
  95         implicit none
  96
  97 !   subr arguments
  98         integer, intent(in) :: iphase_flag, iclm, jclm, k, m, it_mosaic, &
  99                                mmovesect_flag1, idiag_movesect
 100
 101         real(r8), intent(inout) :: rbox(ntot_used)
 102
 103         real(r8), intent(in) :: fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam
 104
 105         real(r8), intent(in) :: drymass_pregrow(maxd_asize,maxd_atype)
 106         real(r8), intent(in) :: drydens_pregrow(maxd_asize,maxd_atype)
 107         real(r8), intent(in) :: drymass_aftgrow(maxd_asize,maxd_atype)
 108         real(r8), intent(in) :: drydens_aftgrow(maxd_asize,maxd_atype)
 109
 110         real(r8), intent(inout) :: adrydens_tmp(maxd_asize,maxd_atype), awetdens_tmp(maxd_asize,maxd_atype)
 111         real(r8), intent(inout) :: adrydpav_tmp(maxd_asize,maxd_atype), awetdpav_tmp(maxd_asize,maxd_atype)
 112         real(r8), intent(inout) :: adryqmas_tmp(maxd_asize,maxd_atype)
 113 ! adrydens_tmp = aerosol dry density (units same as dens_aer)
 114 ! awetdens_tmp = aerosol wet density (units same as dens_aer)
 115 ! adrydpav_tmp = aerosol mean dry diameter (units same as dlo_sect)
 116 ! awetdpav_tmp = aerosol mean wet diameter (units same as dlo_sect)
 117 ! adryqmas_tmp = aerosol total-dry-mass mixing ratio (units same as rbox)
 118
 119 !   local variables
 120         integer iphase, isize, itype,   &
 121           l, ll, llhysw, llwater, lnew, lold, l3,   &
 122           method_movesect, n, nnew, nold
 123         integer nnewsave(2,maxd_asize,maxd_atype)
 124         integer itypenewsave(4,maxd_asize,maxd_atype)
 125
 126         real(r8) densh2o
 127         real(r8) delta_water_conform1, delta_numb_conform1
 128         real(r8) fact_apvolu
 129         real(r8) smallmassaa, smallmassbb
 130
 131         real(r8) dcen_stmp(maxd_asize,maxd_atype)
 132         real(r8) densdefault(maxd_atype)
 133         real(r8) drydenspp(maxd_asize,maxd_atype), drydensxx0(maxd_asize,maxd_atype),   &
 134             drydensxx(maxd_asize,maxd_atype), drydensyy(maxd_asize,maxd_atype)
 135         real(r8) drymasspp(maxd_asize,maxd_atype), drymassxx0(maxd_asize,maxd_atype),   &
 136             drymassxx(maxd_asize,maxd_atype), drymassyy(maxd_asize,maxd_atype)
 137         real(r8) dryvolxx(maxd_asize,maxd_atype), dryvolyy(maxd_asize,maxd_atype)
 138         real(r8) rmassxx(maxd_acomp+2,maxd_asize,maxd_atype),   &
 139             rmassyy(maxd_acomp+2,maxd_asize,maxd_atype)
 140         real(r8) rnumbpp(maxd_asize,maxd_atype), rnumbxx0(maxd_asize,maxd_atype),   &
 141             rnumbxx(maxd_asize,maxd_atype), rnumbyy(maxd_asize,maxd_atype)
 142         real(r8) specdensxx(maxd_acomp)
 143         real(r8) xferfracvol(2,4,maxd_asize,maxd_atype), xferfracnum(2,4,maxd_asize,maxd_atype)
 144         real(r8) volumcen_stmp(maxd_asize,maxd_atype), &
 145             volumlo_stmp(maxd_asize,maxd_atype), volumhi_stmp(maxd_asize,maxd_atype)
 146         real(r8) wetvolxx(maxd_asize,maxd_atype), wetvolyy(maxd_asize,maxd_atype)
 147         real(r8) wetmassxx(maxd_asize,maxd_atype), wetmassyy(maxd_asize,maxd_atype)
 148
 149         character*160 msg
 150
 151
 152 !
 153 !   check for valid inputs
 154 !
 155         if (mmovesect_flag1 <= 0) return
 156         if (ntype_aer       <= 0) return
 157         if (nphase_aer      <= 0) return
 158
 159
 160 !   get "method_movesect" from digits 1-2 of mmovesect_flag1
 161         method_movesect = mod( max(0,mmovesect_flag1), 100 )
 162
 163         if ((method_movesect/10 .eq. 5) .or.   &
 164             (method_movesect/10 .eq. 6)) then
 165             continue
 166         else
 167             msg = '*** subr move_sect_3d error - ' //   &
 168                 'illegal value for mmovesect_flag1'
 169             call peg_error_fatal( lunerr, msg )
 170         end if
 171
 172
 173 !   check iphase_flag
 174         if (iabs(iphase_flag) .eq. 1) then
 175             iphase = ai_phase
 176         else if (iabs(iphase_flag) .eq. 2) then
 177 !           iphase = cw_phase
 178 !           if (nphase_aer .lt. 2) then
 179 !               msg = '*** subr move_sect_3d error - ' //   &
 180 !                   'iphase_flag=2 (after aqueous chemistry) but nphase_aer < 2'
 181 !               call peg_error_fatal( lunerr, msg )
 182 !           else if (cw_phase .ne. 2) then
 183 !               msg = '*** subr move_sect_3d error - ' //   &
 184 !                   'iphase_flag=2 (after aqueous chemistry) but cw_phase .ne. 2'
 185 !               call peg_error_fatal( lunerr, msg )
 186 !           end if
 187             msg = '*** subr move_sect_3d error - ' //   &
 188                 'iphase_flag=2 (after aqueous chemistry) is not implemented'
 189             call peg_error_fatal( lunerr, msg )
 190         else
 191             msg = '*** subr move_sect_3d error - ' //   &
 192                 'iabs(iphase_flag) must be 1 or 2'
 193             call peg_error_fatal( lunerr, msg )
 194         end if
 195
 196
 197 !   when iphase_flag=-1/-2, call move_sect_3d_checkptrs then return
 198 !       if ((ncorecnt .le. 0) .and. (k .le. 1)) then
 199         if (iphase_flag .le. 0) then
 200             write(msg,9040) 'method, idiag', method_movesect, idiag_movesect
 201             call peg_message( lunout, msg )
 202             call move_sect_3d_checkptrs( iphase_flag, iclm, jclm, k, m )
 203             return
 204         else if (it_mosaic .le. 1) then
 205 !   otherwise call in on first time step
 206 !   with movesect_3d, all the types must have the same sizes, species, ...
 207             call move_sect_3d_checkptrs( -iphase_flag, iclm, jclm, k, m )
 208         end if
 209 9040    format( '*** subr move_sect_3d - ', a, ' =', 2(1x,i6) )
 210
 211
 212 !   diagnostics
 213         if (idiag_movesect .ge. 70) then
 214             msg = ' '
 215             call peg_message( lunout, msg )
 216             write(msg,9060) mmovesect_flag1, iclm, jclm, k, m, it_mosaic
 217             call peg_message( lunout, msg )
 218         end if
 219 9060    format( '*** move_sect_3d diags - ', &
 220             'msflag, ijkm, itime =', i7, 4i4, i7 )
 221
 222
 223         densh2o = 1.0
 224         smallmassaa = smallmassaa_asecthp
 225         smallmassbb = smallmassbb_asecthp
 226
 227 !   smallmassaa & smallmassbb are now set in module_data_mosaic_asecthp
 228 !
 229 !   if bin mass mixrat < smallmassaa (1.e-22 g/g), then assume no growth
 230 !   AND no water AND conform number so that size is within bin limits
 231 !
 232 !   if bin mass OR volume mixrat < smallmassbb (1.e-30 g/g), then
 233 !   assume default density to avoid divide by zero
 234
 235
 236 !   calc single particle sizes and volumes in (cm) and (cm^3)
 237         fact_apvolu = fact_apdiam*fact_apdiam*fact_apdiam
 238 !       write(*,'(/a,1p,5e12.4/)') 'fact_ap...', &
 239 !          fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu
 240         do itype = 1, ntype_aer
 241         densdefault(itype) = dens_aer(1,itype)*fact_apdens   ! use density of first component as default
 242         do isize = 1, nsize_aer(itype)
 243             dcen_stmp(    isize,itype) = dcen_sect(    isize,itype)*fact_apdiam
 244             volumcen_stmp(isize,itype) = volumcen_sect(isize,itype)*fact_apvolu
 245             volumlo_stmp( isize,itype) = volumlo_sect( isize,itype)*fact_apvolu
 246             volumhi_stmp( isize,itype) = volumhi_sect( isize,itype)*fact_apvolu
 247         end do
 248         end do
 249
 250
 251 !   process all types simultaneously
 252         do itype = 1, ntype_aer
 253         do isize = 1, nsize_aer(itype)
 254             drydenspp(isize,itype) = drydens_pregrow(isize,itype)*fact_apdens
 255             drydensxx(isize,itype) = drydens_aftgrow(isize,itype)*fact_apdens
 256
 257             drymasspp(isize,itype) = drymass_pregrow(isize,itype)*fact_apmassmr
 258             drymassxx(isize,itype) = drymass_aftgrow(isize,itype)*fact_apmassmr
 259         end do
 260         end do
 261
 262 !       write(*,*) 'move_sect_3d_x1 before initial_conform'
 263         call move_sect_3d_initial_conform(   &
 264           iphase_flag, iclm, jclm, k, m, iphase,   &
 265           method_movesect, idiag_movesect, llhysw, llwater,   &
 266           densdefault, densh2o, smallmassaa, smallmassbb,   &
 267           fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu,   &
 268           delta_water_conform1, delta_numb_conform1,   &
 269           rbox,   &
 270           drydenspp, drymasspp, rnumbpp,   &
 271           drydensxx0, drymassxx0, rnumbxx0,   &
 272           drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx,   &
 273           wetmassxx, wetvolxx, specdensxx,   &
 274           volumcen_stmp, volumlo_stmp, volumhi_stmp )
 275
 276         if (method_movesect/10 .eq. 5) then
 277 !       write(*,*) 'move_sect_3d_x1 before calc_movingcenter'
 278         call move_sect_3d_calc_movingcenter(   &
 279           iphase_flag, iclm, jclm, k, m, iphase,   &
 280           method_movesect, idiag_movesect, llhysw, llwater,   &
 281           nnewsave, itypenewsave,    &
 282           densdefault, densh2o, smallmassaa, smallmassbb,   &
 283           drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx,   &
 284           wetmassxx, wetvolxx,   &
 285           xferfracvol, xferfracnum,   &
 286           volumcen_stmp, volumlo_stmp, volumhi_stmp )
 287         else
 288         call move_sect_3d_calc_masnumadvect(   &
 289           iphase_flag, iclm, jclm, k, m, iphase,   &
 290           method_movesect, idiag_movesect, llhysw, llwater,   &
 291           nnewsave, itypenewsave,    &
 292           densdefault, densh2o, smallmassaa, smallmassbb,   &
 293           drydenspp, drymasspp, rnumbpp,   &
 294           drydensxx0, drymassxx0, rnumbxx0,   &
 295           drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx,   &
 296           wetmassxx, wetvolxx,   &
 297           xferfracvol, xferfracnum,   &
 298           volumcen_stmp, volumlo_stmp, volumhi_stmp )
 299         end if
 300
 301 !       write(*,*) 'move_sect_3d_x1 before apply_moves'
 302         call move_sect_3d_apply_moves(   &
 303           iphase_flag, iclm, jclm, k, m, iphase,   &
 304           method_movesect, idiag_movesect, llhysw, llwater,   &
 305           nnewsave, itypenewsave,    &
 306           densdefault, densh2o, smallmassaa, smallmassbb,   &
 307           fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu,   &
 308           delta_water_conform1, delta_numb_conform1,   &
 309           rbox,   &
 310           drydenspp, drymasspp, rnumbpp,   &
 311           drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx,   &
 312           drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy,   &
 313           xferfracvol, xferfracnum,   &
 314           dcen_stmp, volumcen_stmp, volumlo_stmp, volumhi_stmp,   &
 315           adrydens_tmp, awetdens_tmp, adrydpav_tmp, awetdpav_tmp,   &
 316           adryqmas_tmp )
 317
 318
 319         return
 320         end subroutine move_sect_3d_x1
 321
 322
 323 !-----------------------------------------------------------------------
 324         subroutine move_sect_3d_initial_conform(   &
 325           iphase_flag, iclm, jclm, k, m, iphase,   &
 326           method_movesect, idiag_movesect, llhysw, llwater,   &
 327           densdefault, densh2o, smallmassaa, smallmassbb,   &
 328           fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu,   &
 329           delta_water_conform1, delta_numb_conform1,   &
 330           rbox,   &
 331           drydenspp, drymasspp, rnumbpp,   &
 332           drydensxx0, drymassxx0, rnumbxx0,   &
 333           drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx,   &
 334           wetmassxx, wetvolxx, specdensxx,   &
 335           volumcen_stmp, volumlo_stmp, volumhi_stmp )
 336
 337 !
 338 !   routine does some initial tasks for the section movement
 339 !       load rmassxx & rnumbxx from rbox
 340 !       load specdensxx
 341 !       set drymassxx & dryvolxx from drymass_aftgrow & drydens_aftgrow,
 342 !           OR compute them from rmassxx, specdensxx if need be
 343 !       set wetmassxx & wetvolxx from dry values & water mass
 344 !       conform rnumbxx so that the mean particle size of each section
 345 !           (= dryvolxx/rnumbxx) is within the section limits
 346 !
 347         use module_data_mosaic_asecthp, only:  &
 348             maxd_acomp, maxd_asize, maxd_atype,   &
 349             nsize_aer, ntype_aer, ncomp_aer, ncomp_plustracer_aer, ai_phase,   &
 350             massptr_aer, numptr_aer, waterptr_aer, hyswptr_aer,   &
 351             dens_aer
 352
 353         implicit none
 354
 355 !   subr arguments
 356         integer iphase_flag, iclm, jclm, iphase, k,   &
 357           m, method_movesect, idiag_movesect, llhysw, llwater
 358         real(r8) densh2o, smallmassaa, smallmassbb
 359         real(r8) fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu
 360         real(r8) delta_water_conform1, delta_numb_conform1
 361         real(r8) densdefault(maxd_atype)
 362         real(r8) drydenspp(maxd_asize,maxd_atype), drydensxx0(maxd_asize,maxd_atype),   &
 363              drydensxx(maxd_asize,maxd_atype)
 364         real(r8) drymasspp(maxd_asize,maxd_atype), drymassxx0(maxd_asize,maxd_atype),   &
 365              drymassxx(maxd_asize,maxd_atype)
 366         real(r8) dryvolxx(maxd_asize,maxd_atype)
 367         real(r8) rbox(ntot_used)
 368         real(r8) rmassxx(maxd_acomp+2,maxd_asize,maxd_atype)
 369         real(r8) rnumbpp(maxd_asize,maxd_atype), rnumbxx0(maxd_asize,maxd_atype),   &
 370              rnumbxx(maxd_asize,maxd_atype)
 371         real(r8) specdensxx(maxd_acomp)
 372         real(r8) volumcen_stmp(maxd_asize,maxd_atype), &
 373             volumlo_stmp(maxd_asize,maxd_atype), volumhi_stmp(maxd_asize,maxd_atype)
 374         real(r8) wetvolxx(maxd_asize,maxd_atype)
 375         real(r8) wetmassxx(maxd_asize,maxd_atype)
 376
 377
 378 !   local variables
 379         integer itype, l, ll, lnew, lold, l3, n, nnew, nold
 380
 381         real(r8) dummass, dumnum, dumnum_at_dhi, dumnum_at_dlo, dumr,   &
 382           dumvol, dumvol1p, dumwatrmass
 383
 384
 385 !   load specdens
 386 !   note - all types must have same species in same order
 387 !          and subr move_sect_3d_checkptrs checks for this
 388         itype = 1
 389         do ll = 1, ncomp_plustracer_aer(itype)
 390             specdensxx(ll) = dens_aer(ll,itype)*fact_apdens
 391         end do
 392
 393 !   assure positive definite
 394         do l = 1, ntot_used
 395             rbox(l) = max( 0.0_r8, rbox(l) )
 396         end do
 397
 398 !   initialize these
 399         delta_water_conform1 = 0.0
 400         delta_numb_conform1 = 0.0
 401
 402 !
 403 !   main loop over types
 404 !
 405         do 1900 itype = 1, ntype_aer
 406
 407 !   load mixrats into working arrays and assure positive definite
 408         llhysw = ncomp_plustracer_aer(itype) + 1
 409         llwater = ncomp_plustracer_aer(itype) + 2
 410         do n = 1, nsize_aer(itype)
 411             do ll = 1, ncomp_plustracer_aer(itype)
 412                 l = massptr_aer(ll,n,itype,iphase)
 413                 rmassxx(ll,n,itype) = rbox(l)*fact_apmassmr
 414             end do
 415             rmassxx(llhysw,n,itype) = 0.
 416             l = 0
 417             if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
 418             if (l .gt. 0) rmassxx(llhysw,n,itype) = rbox(l)*fact_apmassmr
 419             rmassxx(llwater,n,itype) = 0.
 420             l = 0
 421             if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
 422             if (l .gt. 0) rmassxx(llwater,n,itype) = rbox(l)*fact_apmassmr
 423
 424             rnumbxx(n,itype)  = rbox(numptr_aer(n,itype,iphase))*fact_apnumbmr
 425             rnumbxx0(n,itype) = rnumbxx(n,itype)
 426             rnumbpp(n,itype)  = rnumbxx(n,itype)
 427
 428             drydensxx0(n,itype) = drydensxx(n,itype)
 429             drymassxx0(n,itype) = drymassxx(n,itype)
 430         end do
 431
 432         do 1390 n = 1, nsize_aer(itype)
 433
 434 !
 435 !   if drydens_aftgrow < 0.1 g/cm^3, then bin had state="no_aerosol"
 436 !   compute volume using default dry-densities, set water=0,
 437 !       and conform the number
 438 !   also do this if mass is extremely small (below smallmassaa)
 439 !       OR if drydens_aftgrow > 20 g/cm^3 (which is unreal(r8))
 440 !
 441         if ( (drydensxx(n,itype) .lt.  0.1) .or.   &
 442              (drydensxx(n,itype) .gt. 20.0) .or.   &
 443              (drymassxx(n,itype) .le. smallmassaa) ) then
 444             dummass = 0.
 445             dumvol = 0.
 446             do ll = 1, ncomp_aer(itype)
 447                 dumr = rmassxx(ll,n,itype)
 448                 dummass = dummass + dumr
 449                 dumvol  = dumvol  + dumr/specdensxx(ll)
 450             end do
 451             drymassxx(n,itype) = dummass
 452             if (min(dummass,dumvol) .le. smallmassbb) then
 453                 drydensxx(n,itype) = densdefault(itype)
 454                 dumvol = dummass/densdefault(itype)
 455                 dumnum = dummass/(volumcen_stmp(n,itype)*densdefault(itype))
 456             else
 457                 drydensxx(n,itype) = dummass/dumvol
 458                 dumnum = rnumbxx(n,itype)
 459                 dumnum_at_dhi = dumvol/volumhi_stmp(n,itype)
 460                 dumnum_at_dlo = dumvol/volumlo_stmp(n,itype)
 461                 dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
 462             end if
 463             delta_numb_conform1 = delta_numb_conform1 + dumnum - rnumbxx(n,itype)
 464             rnumbxx(n,itype) = dumnum
 465             rnumbpp(n,itype) = rnumbxx(n,itype)
 466             delta_water_conform1 = delta_water_conform1 - rmassxx(llwater,n,itype)
 467             rmassxx(llwater,n,itype) = 0.
 468         end if
 469
 470 !   load dry/wet mass and volume into "xx" arrays
 471 !   which hold values before inter-mode transferring
 472         dryvolxx(n,itype) = drymassxx(n,itype)/drydensxx(n,itype)
 473         dumwatrmass = rmassxx(llwater,n,itype)
 474         wetmassxx(n,itype) = drymassxx(n,itype) + dumwatrmass
 475         wetvolxx(n,itype) = dryvolxx(n,itype) + dumwatrmass/densh2o
 476
 477 1390    continue
 478
 479 1900    continue
 480
 481         return
 482         end subroutine move_sect_3d_initial_conform
 483
 484
 485 !-----------------------------------------------------------------------
 486         subroutine move_sect_3d_calc_movingcenter(   &
 487           iphase_flag, iclm, jclm, k, m, iphase,   &
 488           method_movesect, idiag_movesect, llhysw, llwater,   &
 489           nnewsave, itypenewsave,    &
 490           densdefault, densh2o, smallmassaa, smallmassbb,   &
 491           drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx,   &
 492           wetmassxx, wetvolxx,   &
 493           xferfracvol, xferfracnum,   &
 494           volumcen_stmp, volumlo_stmp, volumhi_stmp )
 495 !
 496 !   routine calculates section movements for the moving-center approach
 497 !
 498 !   material in section n will be transfered to section nnewsave(1,n)
 499 !
 500 !   the nnewsave are calculated here
 501 !   the actual transfer is done in another routine
 502 !
 503         use module_data_mosaic_aero, only:  &
 504             method_bcfrac, method_kappa, msectional_flag2
 505         use module_data_mosaic_asecthp, only:  &
 506             dens_aer, hygro_aer,   &
 507             itype_md1_of_itype, itype_md2_of_itype, itype_of_itype_md1md2,   &
 508             lptr_bc_aer, massptr_aer, &
 509             maxd_acomp, maxd_asize, maxd_atype,   &
 510             ncomp_aer, nsize_aer, ntype_aer, ntype_md1_aer, ntype_md2_aer,   &
 511             xcut_atype_md1, xcut_atype_md2
 512
 513         implicit none
 514
 515 !   subr arguments
 516         integer iphase_flag, iclm, jclm, iphase, k,   &
 517           m, method_movesect, idiag_movesect, llhysw, llwater
 518         integer nnewsave(2,maxd_asize,maxd_atype)
 519         integer itypenewsave(4,maxd_asize,maxd_atype)
 520         real(r8) densh2o, smallmassaa, smallmassbb
 521         real(r8) densdefault(maxd_atype)
 522         real(r8) drydensxx(maxd_asize,maxd_atype)
 523         real(r8) drymassxx(maxd_asize,maxd_atype)
 524         real(r8) dryvolxx(maxd_asize,maxd_atype)
 525         real(r8) rmassxx(maxd_acomp+2,maxd_asize,maxd_atype)
 526         real(r8) rnumbxx(maxd_asize,maxd_atype)
 527         real(r8) xferfracvol(2,4,maxd_asize,maxd_atype), xferfracnum(2,4,maxd_asize,maxd_atype)
 528         real(r8) volumcen_stmp(maxd_asize,maxd_atype), &
 529             volumlo_stmp(maxd_asize,maxd_atype), volumhi_stmp(maxd_asize,maxd_atype)
 530         real(r8) wetmassxx(maxd_asize,maxd_atype)
 531         real(r8) wetvolxx(maxd_asize,maxd_atype)
 532
 533 !   local variables
 534         integer isize, itype, itype1, itype2, itypenew, itype1new, itype2new
 535         integer itmpa, itmpb
 536         integer jtype1, jtype2
 537         integer ll, n, ndum, ndumb, nnew, nold
 538         real(r8) dumnum, dumvol, dumvol1p, sixoverpi, third
 539         real(r8) tmpa, tmpb, tmpc, tmpd, tmpe, tmpf
 540         real(r8) tmp_kappa, tmp_wbc
 541         character*160 fmtaa, msg
 542         character*11 txt11
 543
 544
 545         sixoverpi = 6.0/pi
 546         third = 1.0/3.0
 547
 548         if (idiag_movesect .ge. 70) then
 549             call peg_message( lunout, ' ' )
 550             call peg_message( lunout, &
 551                 'move_sect_3d_calc_movingcenter diagnostics' )
 552             call peg_message( lunout, ' ' )
 553         end if
 554
 555 !
 556 !   main loop over types
 557 !
 558         do 1900 itype = 1, ntype_aer
 559
 560 !
 561 !   compute mean size after growth (and corresponding section)
 562 !   particles in section n will be transferred to section nnewsave(1,n)
 563 !
 564         do 1390 n = 1, nsize_aer(itype)
 565
 566         nnew = n
 567         itypenew = itype
 568
 569 !   don't bother to transfer bins whose mass is extremely small
 570         if (drymassxx(n,itype) .le. smallmassaa) goto 1290
 571
 572         dumvol = dryvolxx(n,itype)
 573         dumnum = rnumbxx(n,itype)
 574
 575 !   check for number so small that particle volume is
 576 !   above that of largest section
 577         isize = nsize_aer(itype)
 578         if ( dumnum .le. dumvol/volumhi_stmp(isize,itype) ) then
 579             nnew = nsize_aer(itype)
 580             goto 1250
 581 !   or below that of smallest section
 582         else if ( dumnum .ge. dumvol/volumlo_stmp(1,itype) ) then
 583             nnew = 1
 584             goto 1250
 585         end if
 586
 587 !   dumvol1p is mean particle volume (cm3) for the section
 588         dumvol1p = dumvol/dumnum
 589         if ( dumvol1p .gt. volumhi_stmp(n,itype) ) then
 590             do while ( ( nnew .lt. nsize_aer(itype) ) .and.   &
 591                        ( dumvol1p .gt. volumhi_stmp(nnew,itype) ) )
 592                 nnew = nnew + 1
 593             end do
 594
 595         else if ( dumvol1p .lt. volumlo_stmp(n,itype) ) then
 596             do while ( ( nnew .gt. 1 ) .and.   &
 597                        ( dumvol1p .lt. volumlo_stmp(nnew,itype) ) )
 598                 nnew = nnew - 1
 599             end do
 600
 601         end if
 602
 603 1250    continue
 604         if ((ntype_aer <= 1) .or. (msectional_flag2 <=0)) then
 605             goto 1290   ! here itypenew = itype
 606         end if
 607
 608 !   determine new type
 609         tmpa = 0.0   ! bc mass
 610         tmpb = 0.0   ! sum of mass
 611         tmpc = 0.0   ! sum of volume*kappa
 612         tmpd = 0.0   ! sum of volume
 613         do ll = 1, ncomp_aer(itype)
 614             tmpe = rmassxx(ll,n,itype)
 615             tmpb = tmpb + tmpe
 616 !           if (ll == ibc_a) then  ! this only works in mosaic box model, when the ordering
 617 !                                  ! of species in aer(...) and massptr_aer(...) are the same
 618             if (massptr_aer(ll,n,itype,iphase) == lptr_bc_aer(n,itype,iphase)) then
 619                 tmpa = tmpa + tmpe
 620                 if (method_kappa == 12) cycle
 621                 ! in this case, tmpc & tmpd include non-bc species only
 622             end if
 623             tmpf = tmpe/dens_aer(ll,itype)
 624             tmpc = tmpc + tmpf*hygro_aer(ll,itype)
 625             tmpd = tmpd + tmpf
 626         end do
 627         tmp_wbc   = tmpa/max(tmpb,1.0e-35_r8)
 628         tmp_kappa = tmpc/max(tmpd,1.0e-35_r8)
 629
 630         itype1new = ntype_md1_aer
 631         do jtype1 = 1, ntype_md1_aer - 1
 632             if (tmp_wbc <= xcut_atype_md1(jtype1)) then
 633                 itype1new = jtype1
 634                 exit
 635             end if
 636         end do
 637
 638         itype2new = ntype_md2_aer
 639         do jtype2 = 1, ntype_md2_aer - 1
 640             if (tmp_kappa <= xcut_atype_md2(jtype2)) then
 641                 itype2new = jtype2
 642                 exit
 643             end if
 644         end do
 645
 646         itypenew = itype_of_itype_md1md2(itype1new,itype2new)
 647         itype1 = itype_md1_of_itype(itype)
 648         itype2 = itype_md2_of_itype(itype)
 649
 650         if (idiag_movesect .ge. 70) then
 651             write(msg,'(a,3i5,1p,2e10.2)') 'isz,ity,itynew,wbc,kap', n, itype, itypenew, tmp_wbc, tmp_kappa
 652             call peg_message( lunout, msg )
 653         end if
 654
 655 !   what was this for ???
 656 !       if ((itype1 == ntype_md1_aer) .and. (itype2 == 1)) then
 657 !           if ((n == 11) .or. (n == 17)) then
 658 !               write(*,'(a,2(2x,3i3),2f9.4)') 'is,it1/2, ...new, wbc, kap', &
 659 !                   n, itype1, itype2, &
 660 !                   nnew, itype1new, itype2new, tmp_wbc, tmp_kappa
 661 !           end if
 662 !       end if
 663
 664 1290    continue
 665
 666         nnewsave(:,n,itype) = 0
 667         nnewsave(1,n,itype) = nnew
 668         itypenewsave(:,n,itype) = 0
 669         itypenewsave(1,n,itype) = itypenew
 670         if (method_movesect == 51) itypenewsave(1,n,itype) = itype
 671
 672         xferfracvol(:,:,n,itype) = 0.0
 673         xferfracvol(1,1,n,itype) = 1.0
 674         xferfracnum(:,:,n,itype) = 0.0
 675         xferfracnum(1,1,n,itype) = 1.0
 676
 677 1390    continue
 678
 679
 680 !   diagnostic output
 681         if (idiag_movesect .ge. 70) then
 682             ndum = 0 ; ndumb = 0
 683             do n = 1, nsize_aer(itype)
 684                 if (nnewsave(1,n,itype) .ne. n) ndum = ndum + 1
 685                 if (itypenewsave(1,n,itype) .ne. itype) ndumb = ndumb + 1
 686             end do
 687             if (ndum .gt. 0 .and. ndumb .gt. 0) then
 688                 txt11 = 'movesectYS3'
 689             else if (ndum .gt. 0) then
 690                 txt11 = 'movesectYS1'
 691             else if (ndumb .gt. 0) then
 692                 txt11 = 'movesectYS2'
 693             else
 694                 txt11 = 'movesectNO '
 695             end if
 696             do itmpa = 1, nsize_aer(itype), 24
 697                 itmpb = min( itmpa+23, nsize_aer(itype) )
 698                 if (itmpa == 1) then
 699                     fmtaa = '( a, 4i3, i5, i6, 3x, 24i3 )'
 700                     write(msg,fmtaa) txt11, iclm, jclm, k, m, itype,   &
 701                         ndum, (nnewsave(1,n,itype), n=itmpa,itmpb)
 702                 else
 703                     fmtaa = '( 37x, 24i3 )'
 704                     write(msg,fmtaa) (nnewsave(1,n,itype), n=itmpa,itmpb)
 705                 end if
 706                 call peg_message( lunout, msg )
 707             end do
 708
 709             txt11 = 'itypenew   '
 710             do itmpa = 1, nsize_aer(itype), 24
 711                 itmpb = min( itmpa+23, nsize_aer(itype) )
 712                 if (itmpa == 1) then
 713                     fmtaa = '( 26x, a, 24i3 )'
 714                     write(msg,fmtaa) txt11, (itypenewsave(1,n,itype), n=itmpa,itmpb)
 715                 else
 716                     fmtaa = '( 37x, 24i3 )'
 717                     write(msg,fmtaa) (itypenewsave(1,n,itype), n=itmpa,itmpb)
 718                 end if
 719                 call peg_message( lunout, msg )
 720             end do
 721
 722         end if
 723
 724
 725 1900    continue
 726
 727         return
 728         end subroutine move_sect_3d_calc_movingcenter
 729
 730
 731 !-----------------------------------------------------------------------
 732         subroutine move_sect_3d_calc_masnumadvect(   &
 733           iphase_flag, iclm, jclm, k, m, iphase,   &
 734           method_movesect, idiag_movesect, llhysw, llwater,   &
 735           nnewsave, itypenewsave,    &
 736           densdefault, densh2o, smallmassaa, smallmassbb,   &
 737           drydenspp, drymasspp, rnumbpp,   &
 738           drydensxx0, drymassxx0, rnumbxx0,   &
 739           drydensxx, drymassxx, dryvolxx, rmassxx, rnumbxx,   &
 740           wetmassxx, wetvolxx,   &
 741           xferfracvol, xferfracnum,   &
 742           volumcen_stmp, volumlo_stmp, volumhi_stmp )
 743 !
 744 !   routine calculates section movements for the mass-number-advection approach
 745 !
 746 !   material in section n will be transfered to sections
 747 !       nnewsave(1,n) and nnewsave(2,n)
 748 !   the fractions of mass/volume transfered to each are
 749 !       xferfracvol(1,n) and xferfracvol(2,n)
 750 !   the fractions of number transfered to each are
 751 !       xferfracnum(1,n) and xferfracnum(2,n)
 752 !
 753 !   the nnewsave, xferfracvol, and xferfracnum are calculated here
 754 !   the actual transfer is done in another routine
 755 !
 756         use module_data_mosaic_aero, only:  &
 757             method_bcfrac, method_kappa, msectional_flag2
 758         use module_data_mosaic_asecthp, only:  &
 759             dens_aer, hygro_aer,   &
 760             itype_md1_of_itype, itype_md2_of_itype, itype_of_itype_md1md2,   &
 761             lptr_bc_aer, massptr_aer, &
 762             maxd_acomp, maxd_asize, maxd_atype,   &
 763             ncomp_aer, nsize_aer, ntype_aer, ntype_md1_aer, ntype_md2_aer,   &
 764             xcut_atype_md1, xcut_atype_md2
 765
 766         implicit none
 767
 768 !   subr arguments
 769         integer iphase_flag, iclm, jclm, iphase, k,   &
 770           m, method_movesect, idiag_movesect, llhysw, llwater
 771         integer nnewsave(2,maxd_asize,maxd_atype)
 772         integer itypenewsave(4,maxd_asize,maxd_atype)
 773
 774         real(r8) densh2o, smallmassaa, smallmassbb
 775         real(r8) densdefault(maxd_atype)
 776         real(r8) drydenspp(maxd_asize,maxd_atype), drydensxx0(maxd_asize,maxd_atype),   &
 777              drydensxx(maxd_asize,maxd_atype)
 778         real(r8) drymasspp(maxd_asize,maxd_atype), drymassxx0(maxd_asize,maxd_atype),   &
 779              drymassxx(maxd_asize,maxd_atype)
 780         real(r8) dryvolxx(maxd_asize,maxd_atype)
 781         real(r8) rmassxx(maxd_acomp+2,maxd_asize,maxd_atype)
 782         real(r8) rnumbpp(maxd_asize,maxd_atype), rnumbxx0(maxd_asize,maxd_atype),   &
 783              rnumbxx(maxd_asize,maxd_atype)
 784         real(r8) xferfracvol(2,4,maxd_asize,maxd_atype), xferfracnum(2,4,maxd_asize,maxd_atype)
 785         real(r8) volumcen_stmp(maxd_asize,maxd_atype), &
 786             volumlo_stmp(maxd_asize,maxd_atype), volumhi_stmp(maxd_asize,maxd_atype)
 787         real(r8) wetvolxx(maxd_asize,maxd_atype)
 788         real(r8) wetmassxx(maxd_asize,maxd_atype)
 789
 790 !   local variables
 791         integer ierr, itype, itype1, itype2, itypenew, itype1new, itype2new
 792         integer iforce_movecenter(maxd_asize,maxd_atype)
 793         integer jtype1, jtype2
 794         integer kktypenew
 795         integer ll
 796         integer n, nnew, nnew2
 797
 798         real(r8) dum1, dum2, dum3
 799         real(r8) dumaa, dumbb, dumgamma, dumratio
 800         real(r8) dumfracnum, dumfracvol
 801         real(r8) dumntot
 802         real(r8) dumv
 803         real(r8) dumvbar_aft, dumvbar_pre
 804         real(r8) dumvcutlo_nnew_pre, dumvcuthi_nnew_pre
 805         real(r8) dumvlo_pre, dumvhi_pre, dumvdel_pre
 806         real(r8) dumvtot_aft, dumvtot_pre
 807         real(r8) dumzlo, dumzhi
 808         real(r8) sixoverpi, third, tmpa
 809         real(r8) tmpb, tmpc, tmpd, tmpe, tmpf, tmp_kappa, tmp_wbc
 810
 811         character*4 dumch4
 812         character*1 dumch1
 813         character*160 msg
 814
 815
 816         sixoverpi = 6.0/pi
 817         third = 1.0/3.0
 818         itypenewsave(:,:,:) = 0
 819
 820 !
 821 !   compute mean size after growth (and corresponding section)
 822 !   some of the particles in section n will be transferred to section nnewsave(1,n,itype)
 823 !
 824 !   if the aftgrow mass is extremely small,
 825 !   OR if the aftgrow mean size is outside of
 826 !       [dlo_sect(1,itype), dhi_sect(nsize_aer(itype),itype)]
 827 !   then use the moving-center method_movesect for this bin
 828 !   (don't try to compute the pregrow within-bin distribution)
 829 !
 830         do 4900 itype = 1, ntype_aer
 831
 832         if (idiag_movesect .ge. 70) then
 833             call peg_message( lunout, ' ' )
 834             write( msg, '(a,i5)' ) &
 835                 'move_sect_3d_calc_massnumadvect diagnostics - itype ', &
 836                 itype
 837             call peg_message( lunout, msg )
 838         end if
 839
 840         do 3900 n = 1, nsize_aer(itype)
 841
 842         nnew = n
 843         iforce_movecenter(n,itype) = 0
 844
 845         kktypenew = 1
 846         itypenewsave(1,n,itype) = itype
 847
 848 !   the following lines set xferfrac values to be appropriate for
 849 !      either no transfer or moving-center transfer
 850 !   in these cases, nnewsave(1,n,itype)=n or nnew, and nnewsave(2,n,itype)=0
 851         xferfracvol(:,:,n,itype) = 0.0
 852         xferfracnum(:,:,n,itype) = 0.0
 853         xferfracvol(1,kktypenew,n,itype) = 1.0
 854         xferfracnum(1,kktypenew,n,itype) = 1.0
 855
 856         dumvtot_aft = -1.0
 857         dumvtot_pre = -1.0
 858         dumvbar_aft = -1.0
 859         dumvbar_pre = -1.0
 860         dumvlo_pre = -1.0
 861         dumvhi_pre = -1.0
 862         dumgamma = -1.0
 863         dumratio = -1.0
 864         dumvcutlo_nnew_pre = volumlo_stmp(nnew,itype)*(dumvbar_pre/dumvbar_aft)
 865         dumvcuthi_nnew_pre = volumhi_stmp(nnew,itype)*(dumvbar_pre/dumvbar_aft)
 866         dumfracvol = -1.0
 867         dumfracnum = -1.0
 868
 869 !   don't bother to transfer bins whose mass is extremely small
 870         if (drymassxx(n,itype) .le. smallmassaa) then
 871             iforce_movecenter(n,itype) = 1
 872             goto 1290
 873         end if
 874
 875         dumvtot_aft = dryvolxx(n,itype)
 876         dumntot = rnumbxx(n,itype)
 877
 878 !   check for particle volume above that of largest section
 879 !   or below that of smallest section
 880         if (dumntot .le. dumvtot_aft/volumhi_stmp(nsize_aer(itype),itype)) then
 881             nnew = nsize_aer(itype)
 882             iforce_movecenter(n,itype) = 2
 883             goto 1290
 884         else if (dumntot .ge. dumvtot_aft/volumlo_stmp(1,itype)) then
 885             nnew = 1
 886             iforce_movecenter(n,itype) = 3
 887             goto 1290
 888         end if
 889
 890 !   dumvbar_aft is mean particle volume (cm3) for the section
 891 !   find the section that encloses this volume
 892         dumvbar_aft = dumvtot_aft/dumntot
 893         if (dumvbar_aft .gt. volumhi_stmp(n,itype)) then
 894             do while ( (nnew .lt. nsize_aer(itype)) .and.   &
 895                        (dumvbar_aft .gt. volumhi_stmp(nnew,itype)) )
 896                 nnew = nnew + 1
 897             end do
 898
 899         else if (dumvbar_aft .lt. volumlo_stmp(n,itype)) then
 900             do while ( (nnew .gt. 1) .and.   &
 901                        (dumvbar_aft .lt. volumlo_stmp(nnew,itype)) )
 902                 nnew = nnew - 1
 903             end do
 904
 905         end if
 906
 907 1290    nnewsave(:,n,itype) = 0
 908         nnewsave(1,n,itype) = nnew
 909
 910         if (iforce_movecenter(n,itype) .gt. 0) goto 3700
 911
 912
 913 !   if drydenspp (pregrow) < 0.1 (because bin had state="no_aerosol" before
 914 !       growth was computed, so its initial mass was very small)
 915 !   then use the moving-center method_movesect for this bin
 916 !   (don't try to compute the pregrow within-bin distribution)
 917 !   also do this if pregrow mass is extremely small (below smallmassaa)
 918 !       OR if drydenspp > 20 g/cm^3 (unphysical)
 919         if ( (drydenspp(n,itype) .lt.  0.1) .or.   &
 920              (drydenspp(n,itype) .gt. 20.0) .or.   &
 921              (drymasspp(n,itype) .le. smallmassaa) ) then
 922             iforce_movecenter(n,itype) = 11
 923             goto 3700
 924         end if
 925
 926         dumvtot_pre = drymasspp(n,itype)/drydenspp(n,itype)
 927
 928         dumvlo_pre = volumlo_stmp(n,itype)
 929         dumvhi_pre = volumhi_stmp(n,itype)
 930         dumvdel_pre = dumvhi_pre - dumvlo_pre
 931
 932 !   if the pregrow mean size is outside of OR very close to the bin limits,
 933 !   then use moving-center approach for this bin
 934         dumv = dumvhi_pre - 0.01*dumvdel_pre
 935         if (dumntot .le. dumvtot_pre/dumv) then
 936             iforce_movecenter(n,itype) = 12
 937             goto 3700
 938         end if
 939         dumv = dumvlo_pre + 0.01*dumvdel_pre
 940         if (dumntot .ge. dumvtot_pre/dumv) then
 941             iforce_movecenter(n,itype) = 13
 942             goto 3700
 943         end if
 944
 945 !   calculate the pregrow within-section size distribution
 946         dumvbar_pre = dumvtot_pre/dumntot
 947         dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
 948         dumratio = dumvbar_pre/dumvlo_pre
 949
 950         if (dumratio .le. (1.0001 + dumgamma/3.0)) then
 951             dumv = dumvlo_pre + 3.0*(dumvbar_pre-dumvlo_pre)
 952             dumvhi_pre = min( dumvhi_pre, dumv )
 953             dumvdel_pre = dumvhi_pre - dumvlo_pre
 954             dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
 955             dumratio = dumvbar_pre/dumvlo_pre
 956         else if (dumratio .ge. (0.9999 + dumgamma*2.0/3.0)) then
 957             dumv = dumvhi_pre + 3.0*(dumvbar_pre-dumvhi_pre)
 958             dumvlo_pre = max( dumvlo_pre, dumv )
 959             dumvdel_pre = dumvhi_pre - dumvlo_pre
 960             dumgamma = (dumvhi_pre/dumvlo_pre) - 1.0
 961             dumratio = dumvbar_pre/dumvlo_pre
 962         end if
 963
 964         dumbb = (dumratio - 1.0 - 0.5*dumgamma)*12.0/dumgamma
 965         dumaa = 1.0 - 0.5*dumbb
 966
 967 !   calculate pregrow volumes corresponding to the nnew
 968 !   section boundaries
 969         dumvcutlo_nnew_pre = volumlo_stmp(nnew,itype)*(dumvbar_pre/dumvbar_aft)
 970         dumvcuthi_nnew_pre = volumhi_stmp(nnew,itype)*(dumvbar_pre/dumvbar_aft)
 971
 972 !   if the [dumvlo_pre, dumvhi_pre] falls completely within
 973 !   the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval,
 974 !   then all mass and number go to nnew
 975         if (nnew .eq. 1) then
 976             if (dumvhi_pre .le. dumvcuthi_nnew_pre) then
 977                 iforce_movecenter(n,itype) = 21
 978             else
 979                 nnew2 = nnew + 1
 980             end if
 981         else if (nnew .eq. nsize_aer(itype)) then
 982             if (dumvlo_pre .ge. dumvcutlo_nnew_pre) then
 983                 iforce_movecenter(n,itype) = 22
 984             else
 985                 nnew2 = nnew - 1
 986             end if
 987         else
 988             if ((dumvlo_pre .ge. dumvcutlo_nnew_pre) .and.   &
 989                 (dumvhi_pre .le. dumvcuthi_nnew_pre)) then
 990                 iforce_movecenter(n,itype) = 23
 991             else if (dumvlo_pre .lt. dumvcutlo_nnew_pre) then
 992                 nnew2 = nnew - 1
 993             else
 994                 nnew2 = nnew + 1
 995             end if
 996         end if
 997         if (iforce_movecenter(n,itype) .gt. 0) goto 3700
 998
 999 !   calculate the fraction of ntot and vtot that are within
1000 !   the [dumvcutlo_nnew_pre, dumvcuthi_nnew_pre] interval
1001         dumzlo = (dumvcutlo_nnew_pre - dumvlo_pre)/dumvdel_pre
1002         dumzhi = (dumvcuthi_nnew_pre - dumvlo_pre)/dumvdel_pre
1003         dumzlo = max( dumzlo, 0.0_r8 )
1004         dumzhi = min( dumzhi, 1.0_r8 )
1005         dum1 =  dumzhi    - dumzlo
1006         dum2 = (dumzhi**2 - dumzlo**2)*0.5
1007         dum3 = (dumzhi**3 - dumzlo**3)/3.0
1008         dumfracnum = dumaa*dum1 + dumbb*dum2
1009         dumfracvol = (dumvlo_pre/dumvbar_pre) * (dumaa*dum1 +   &
1010                 (dumaa*dumgamma + dumbb)*dum2 + (dumbb*dumgamma)*dum3)
1011
1012         if ((dumfracnum .le. 0.0) .or. (dumfracvol .le. 0.0)) then
1013             iforce_movecenter(n,itype) = 31
1014             nnewsave(1,n,itype) = nnew2
1015         else if ((dumfracnum .ge. 1.0) .or. (dumfracvol .ge. 1.0)) then
1016             iforce_movecenter(n,itype) = 32
1017         end if
1018         if (iforce_movecenter(n,itype) .gt. 0) goto 3700
1019
1020         nnewsave(2,n,itype) = nnew2
1021
1022 !   at this point, iforce_movecenter=0, so bin will be transferred
1023 !      to 2 size bins
1024 !   currently we do linear-discrete for size and moving-center for type
1025 !   thus for type, everything will go to a single
1026 !      itype = itypenewsave(1,n,itype)
1027 !
1028 !   so the xferfrac are as follows
1029 !      xferfracxxx(1,ktypenew,n,itype) is amount going to
1030 !         isize = nnewsave(1,n,itype) and itype=itypenewsave(1,n,itype)
1031 !      xferfracxxx(2,ktypenew,n,itype) is amount going to
1032 !         isize = nnewsave(2,n,itype) and itype=itypenewsave(1,n,itype)
1033         xferfracvol(1,kktypenew,n,itype) = dumfracvol
1034         xferfracvol(2,kktypenew,n,itype) = 1.0 - dumfracvol
1035         xferfracnum(1,kktypenew,n,itype) = dumfracnum
1036         xferfracnum(2,kktypenew,n,itype) = 1.0 - dumfracnum
1037
1038 3700    continue
1039
1040 !   determine new type
1041         tmpa = 0.0   ! bc mass
1042         tmpb = 0.0   ! sum of mass
1043         tmpc = 0.0   ! sum of volume*kappa
1044         tmpd = 0.0   ! sum of volume
1045         do ll = 1, ncomp_aer(itype)
1046             tmpe = rmassxx(ll,n,itype)
1047             tmpb = tmpb + tmpe
1048 !           if (ll == ibc_a) then  ! this only works in mosaic box model, when the ordering
1049 !                                  ! of species in aer(...) and massptr_aer(...) are the same
1050             if (massptr_aer(ll,n,itype,iphase) == lptr_bc_aer(n,itype,iphase)) then
1051                 tmpa = tmpa + tmpe
1052                 if (method_kappa == 12) cycle
1053                 ! in this case, tmpc & tmpd include non-bc species only
1054             end if
1055             tmpf = tmpe/dens_aer(ll,itype)
1056             tmpc = tmpc + tmpf*hygro_aer(ll,itype)
1057             tmpd = tmpd + tmpf
1058         end do
1059         tmp_wbc   = tmpa/max(tmpb,1.0e-35_r8)
1060         tmp_kappa = tmpc/max(tmpd,1.0e-35_r8)
1061
1062         itype1new = ntype_md1_aer
1063         do jtype1 = 1, ntype_md1_aer - 1
1064             if (tmp_wbc <= xcut_atype_md1(jtype1)) then
1065                 itype1new = jtype1
1066                 exit
1067             end if
1068         end do
1069
1070         itype2new = ntype_md2_aer
1071         do jtype2 = 1, ntype_md2_aer - 1
1072             if (tmp_kappa <= xcut_atype_md2(jtype2)) then
1073                 itype2new = jtype2
1074                 exit
1075             end if
1076         end do
1077
1078         itypenew = itype_of_itype_md1md2(itype1new,itype2new)
1079         itype1 = itype_md1_of_itype(itype)
1080         itype2 = itype_md2_of_itype(itype)
1081
1082         itypenewsave(:,n,itype) = 0
1083         itypenewsave(1,n,itype) = itypenew
1084         if (method_movesect == 61) itypenewsave(1,n,itype) = itype
1085
1086
1087 !   diagnostic output
1088         if (idiag_movesect .lt. 70) goto 3800
1089
1090         if (nnewsave(2,n,itype) .eq. 0) then
1091             if (nnewsave(1,n,itype) .eq. 0) then
1092                 dumch4 = 'NO X'
1093             else if (nnewsave(1,n,itype) .eq. n) then
1094                 dumch4 = 'NO A'
1095             else
1096                 dumch4 = 'YESA'
1097             end if
1098         else if (nnewsave(1,n,itype) .eq. 0) then
1099             if (nnewsave(2,n,itype) .eq. n) then
1100                 dumch4 = 'NO B'
1101             else
1102                 dumch4 = 'YESB'
1103             end if
1104         else if (nnewsave(2,n,itype) .eq. n) then
1105             if (nnewsave(1,n,itype) .eq. n) then
1106                 dumch4 = 'NO Y'
1107             else
1108                 dumch4 = 'YESC'
1109             end if
1110         else if (nnewsave(1,n,itype) .eq. n) then
1111             dumch4 = 'YESD'
1112         else
1113             dumch4 = 'YESE'
1114         end if
1115
1116         dumch1 = '+'
1117         if (drymasspp(n,itype) .gt. drymassxx(n,itype)) dumch1 = '-'
1118
1119         msg = ' '
1120         call peg_message( lunout, msg )
1121         write(msg,97010) dumch1, dumch4, iclm, jclm, k, m,   &
1122                 n, nnewsave(1,n,itype), nnewsave(2,n,itype), &
1123                 iforce_movecenter(n,itype), itypenewsave(1,n,itype)
1124         call peg_message( lunout, msg )
1125         write(msg,97020) 'pre mass, dens      ',   &
1126                 drymasspp(n,itype), drydenspp(n,itype)
1127         call peg_message( lunout, msg )
1128         write(msg,97020) 'aft mass, dens, numb',   &
1129                 drymassxx(n,itype), drydensxx(n,itype), rnumbxx(n,itype)
1130         call peg_message( lunout, msg )
1131         if ((drydensxx(n,itype) .ne. drydensxx0(n,itype)) .or.   &
1132             (drymassxx(n,itype) .ne. drymassxx0(n,itype)) .or.   &
1133             (rnumbxx(n,itype)   .ne. rnumbxx0(n,itype)  )) then
1134             write(msg,97020) 'aft0 mas, dens, numb',   &
1135                 drymassxx0(n,itype), drydensxx0(n,itype), rnumbxx0(n,itype)
1136             call peg_message( lunout, msg )
1137         end if
1138         write(msg,97020) 'vlop0, vbarp,  vhip0',   &
1139                 volumlo_stmp(n,itype), dumvbar_pre, volumhi_stmp(n,itype)
1140         call peg_message( lunout, msg )
1141         write(msg,97020) 'vlop , vbarp,  vhip ',   &
1142                 dumvlo_pre, dumvbar_pre, dumvhi_pre
1143         call peg_message( lunout, msg )
1144         write(msg,97020) 'vloax, vbarax, vhiax',   &
1145                 dumvcutlo_nnew_pre, dumvbar_pre, dumvcuthi_nnew_pre
1146         call peg_message( lunout, msg )
1147         write(msg,97020) 'vloa0, vbara,  vhia0',   &
1148                 volumlo_stmp(nnew,itype), dumvbar_aft, volumhi_stmp(nnew,itype)
1149         call peg_message( lunout, msg )
1150         write(msg,97020) 'dumfrvol, num, ratio',   &
1151                 dumfracvol, dumfracnum, dumratio
1152         call peg_message( lunout, msg )
1153         write(msg,97020) 'frvol,num1; vol,num2',   &
1154                 xferfracvol(1,kktypenew,n,itype), xferfracnum(1,kktypenew,n,itype),   &
1155                 xferfracvol(2,kktypenew,n,itype), xferfracnum(2,kktypenew,n,itype)
1156         call peg_message( lunout, msg )
1157
1158 97010   format( 'movesect', 2a, 7x, 4i3, 4x,   &
1159                 'n,nnews', 3i3, 4x, 'iforce', i3.2, 4x, 'itypenew', i5 )
1160 97020   format( a, 1p, 4e13.4 )
1161
1162 3800    continue
1163
1164 !
1165 !   check for legal combinations of nnewsave(1,:,n,:) & nnewsave(2,:,n,:)
1166 !   error if
1167 !     nnew1 == nnew2
1168 !     both are non-zero AND iabs(nnew1-nnew2) != 1
1169         ierr = 0
1170         if (nnewsave(1,n,itype) .eq. nnewsave(2,n,itype)) then
1171             ierr = 1
1172         else if (nnewsave(1,n,itype)*nnewsave(2,n,itype) .ne. 0) then
1173             if (iabs(nnewsave(1,n,itype)-nnewsave(2,n,itype)) .ne. 1) ierr = 1
1174         end if
1175         if (ierr .gt. 0) then
1176             write(msg,97010) 'E', 'RROR', iclm, jclm, k, m,   &
1177                 n, nnewsave(1,n,itype), nnewsave(2,n,itype), iforce_movecenter(n,itype)
1178             call peg_message( lunout, msg )
1179         end if
1180
1181
1182 !   if method_movesect == 30-31 then force moving center
1183 !   this is just for testing purposes
1184         if ((method_movesect .ge. 30) .and. (method_movesect .le. 39)) then
1185             nnewsave(1,n,itype) = nnew
1186             nnewsave(2,n,itype) = 0
1187             xferfracvol(:,:,n,itype) = 0.0
1188             xferfracnum(:,:,n,itype) = 0.0
1189             xferfracvol(1,kktypenew,n,itype) = 1.0
1190             xferfracnum(1,kktypenew,n,itype) = 1.0
1191         end if
1192
1193 3900    continue   ! do isize = ...
1194
1195 4900    continue   ! do itype = ...
1196
1197         return
1198         end subroutine move_sect_3d_calc_masnumadvect
1199
1200
1201 !-----------------------------------------------------------------------
1202         subroutine move_sect_3d_apply_moves(   &
1203           iphase_flag, iclm, jclm, k, m, iphase,   &
1204           method_movesect, idiag_movesect, llhysw, llwater,   &
1205           nnewsave, itypenewsave,    &
1206           densdefault, densh2o, smallmassaa, smallmassbb,   &
1207           fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu,   &
1208           delta_water_conform1, delta_numb_conform1,   &
1209           rbox,   &
1210           drydenspp, drymasspp, rnumbpp,   &
1211           drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx,   &
1212           drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy,   &
1213           xferfracvol, xferfracnum,   &
1214           dcen_stmp, volumcen_stmp, volumlo_stmp, volumhi_stmp,   &
1215           adrydens_tmp, awetdens_tmp, adrydpav_tmp, awetdpav_tmp,   &
1216           adryqmas_tmp )
1217 !
1218 !   routine performs the actual transfer of aerosol number and mass
1219 !       between sections
1220 !
1221         use module_data_mosaic_asecthp, only:  &
1222             maxd_acomp, maxd_asize, maxd_atype,   &
1223             nsize_aer, ntype_aer, ncomp_plustracer_aer, ai_phase,   &
1224             massptr_aer, numptr_aer, waterptr_aer, hyswptr_aer
1225
1226         implicit none
1227
1228 !   subr arguments
1229         integer iphase_flag, iclm, jclm, iphase, k,   &
1230           m, method_movesect, idiag_movesect, llhysw, llwater
1231         integer nnewsave(2,maxd_asize,maxd_atype)
1232         integer itypenewsave(4,maxd_asize,maxd_atype)
1233
1234         real(r8) densh2o, smallmassaa, smallmassbb
1235         real(r8) fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, fact_apvolu
1236         real(r8) delta_water_conform1, delta_numb_conform1
1237         real(r8) densdefault(maxd_atype)
1238         real(r8) drydenspp(maxd_asize,maxd_atype)
1239         real(r8) drymasspp(maxd_asize,maxd_atype)
1240         real(r8) drymassxx(maxd_asize,maxd_atype), drymassyy(maxd_asize,maxd_atype)
1241         real(r8) dryvolxx(maxd_asize,maxd_atype), dryvolyy(maxd_asize,maxd_atype)
1242         real(r8) rbox(ntot_used)
1243         real(r8) rmassxx(maxd_acomp+2,maxd_asize,maxd_atype),   &
1244              rmassyy(maxd_acomp+2,maxd_asize,maxd_atype)
1245         real(r8) rnumbpp(maxd_asize,maxd_atype)
1246         real(r8) rnumbxx(maxd_asize,maxd_atype), rnumbyy(maxd_asize,maxd_atype)
1247         real(r8) xferfracvol(2,4,maxd_asize,maxd_atype), xferfracnum(2,4,maxd_asize,maxd_atype)
1248         real(r8) dcen_stmp(maxd_asize,maxd_atype), volumcen_stmp(maxd_asize,maxd_atype), &
1249             volumlo_stmp(maxd_asize,maxd_atype), volumhi_stmp(maxd_asize,maxd_atype)
1250         real(r8) wetvolxx(maxd_asize,maxd_atype), wetvolyy(maxd_asize,maxd_atype)
1251         real(r8) wetmassxx(maxd_asize,maxd_atype), wetmassyy(maxd_asize,maxd_atype)
1252         real(r8) adrydens_tmp(maxd_asize,maxd_atype),  awetdens_tmp(maxd_asize,maxd_atype)
1253         real(r8) adrydpav_tmp(maxd_asize,maxd_atype), awetdpav_tmp(maxd_asize,maxd_atype)
1254         real(r8) adryqmas_tmp(maxd_asize,maxd_atype)
1255
1256 !   local variables
1257         integer itmpa, itype, itypenew, jj, kk, l, ll, n, ndum, nnew, nold
1258         integer jja, jjb, jjc
1259
1260         real(r8) delta_numb_conform2,   &
1261           dumbot, dumnum, dumnum_at_dhi, dumnum_at_dlo,   &
1262           dumvol, dumvol1p, dumxfvol, dumxfnum, sixoverpi, third
1263         real(r8) dumpp(maxd_asize), dumxx(maxd_asize), dumyy(maxd_asize),   &
1264           dumout(maxd_asize)
1265
1266         character*160 msg
1267         character*8 dumch8
1268         character*4 dumch4
1269
1270
1271         sixoverpi = 6.0/pi
1272         third = 1.0/3.0
1273
1274 !
1275 !   initialize "yy" arrays that hold values after inter-mode transferring
1276 !       "yy" = "xx" for sections that do not move at all
1277 !       "yy" = 0.0  for sections that do move (partially or completely)
1278 !
1279         do 1950 itype = 1, ntype_aer
1280         do 1940 n = 1, nsize_aer(itype)
1281
1282         itmpa = sum( itypenewsave(2:4,n,itype) )
1283         if ( (nnewsave(1,n,itype) .eq. n) .and.   &
1284              (nnewsave(2,n,itype) .eq. 0) .and.   &
1285              (itypenewsave(1,n,itype) .eq. itype) .and.   &
1286              (itmpa                   .eq.     0) ) then
1287 !   if nnew == n AND itypenew == itype,
1288 !       then no material will be transferred out of section n,itype,
1289 !       and section can be intialized to its "old" values
1290 !   so initialize "yy" arrays with "xx" values
1291             drymassyy(n,itype) = drymassxx(n,itype)
1292             dryvolyy(n,itype) = dryvolxx(n,itype)
1293             wetmassyy(n,itype) = wetmassxx(n,itype)
1294             wetvolyy(n,itype) = wetvolxx(n,itype)
1295             rnumbyy(n,itype) = rnumbxx(n,itype)
1296             do ll = 1, ncomp_plustracer_aer(itype) + 2
1297                 rmassyy(ll,n,itype) = rmassxx(ll,n,itype)
1298             end do
1299
1300         else
1301 !   if nnew /= n OR itypenew /= itype,
1302 !       then some material will be transferred out of section n,itype,
1303 !       and section must be initialized to zeroa
1304 !   so initialize "yy" arrays to zero
1305             drymassyy(n,itype) = 0.0
1306             dryvolyy(n,itype) = 0.0
1307             wetmassyy(n,itype) = 0.0
1308             wetvolyy(n,itype) = 0.0
1309             rnumbyy(n,itype) = 0.0
1310             do ll = 1, ncomp_plustracer_aer(itype) + 2
1311                 rmassyy(ll,n,itype) = 0.0
1312             end do
1313
1314         end if
1315
1316
1317 1940    continue   ! n=isize loop
1318 1950    continue   ! itype loop
1319
1320 !
1321 !   do the transfer of mass and number
1322 !
1323         do 2950 itype = 1, ntype_aer
1324         do 2900 n = 1, nsize_aer(itype)
1325
1326 !   check for no transfer out of n,itype
1327         itmpa = sum( itypenewsave(2:4,n,itype) )
1328         if ( (nnewsave(1,n,itype) .eq. n) .and.   &
1329              (nnewsave(2,n,itype) .eq. 0) .and.   &
1330              (itypenewsave(1,n,itype) .eq. itype) .and.   &
1331              (itmpa                   .eq.     0) ) goto 2900
1332
1333         do 2820 jj = 1, 2
1334         nnew = nnewsave(jj,n,itype)
1335         if (nnew .le. 0) goto 2820
1336
1337         do 2810 kk = 1, 4
1338         itypenew = itypenewsave(kk,n,itype)
1339         if (itypenew .le. 0) goto 2810
1340
1341         dumxfvol = xferfracvol(jj,kk,n,itype)
1342         dumxfnum = xferfracnum(jj,kk,n,itype)
1343         if ((dumxfvol .le. 0.0) .and. (dumxfnum .le. 0.0)) goto 2810
1344
1345         do ll = 1, ncomp_plustracer_aer(itype) + 2
1346             rmassyy(ll,nnew,itypenew) = rmassyy(ll,nnew,itypenew) + rmassxx(ll,n,itype)*dumxfvol
1347         end do
1348         rnumbyy(nnew,itypenew) = rnumbyy(nnew,itypenew) + rnumbxx(n,itype)*dumxfnum
1349
1350         drymassyy(nnew,itypenew) = drymassyy(nnew,itypenew) + drymassxx(n,itype)*dumxfvol
1351         dryvolyy( nnew,itypenew) = dryvolyy( nnew,itypenew) + dryvolxx( n,itype)*dumxfvol
1352         wetmassyy(nnew,itypenew) = wetmassyy(nnew,itypenew) + wetmassxx(n,itype)*dumxfvol
1353         wetvolyy( nnew,itypenew) = wetvolyy( nnew,itypenew) + wetvolxx( n,itype)*dumxfvol
1354
1355 2810    continue
1356 2820    continue
1357
1358 2900    continue   ! n=isize loop
1359 2950    continue   ! itype loop
1360
1361 !
1362 !   transfer among sections is completed
1363 !   - check for conservation of mass/volume/number
1364 !   - conform number again
1365 !   - compute/store densities and mean sizes
1366 !   - if k=1, save values for use by dry deposition routine
1367 !   - copy new mixrats back to rbox array
1368 !
1369         call move_sect_3d_conserve_check(   &
1370           1, iphase_flag, iclm, jclm, k, m, iphase,   &
1371           method_movesect, idiag_movesect, llhysw, llwater,   &
1372           nnewsave, itypenewsave,    &
1373           densdefault, densh2o, smallmassaa, smallmassbb,   &
1374           fact_apmassmr, fact_apnumbmr,   &
1375           delta_water_conform1, delta_numb_conform1,   &
1376           rbox,   &
1377           drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx,   &
1378           drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1379
1380         delta_numb_conform2 = 0.0
1381
1382         do 3950 itype = 1, ntype_aer
1383         do 3900 n = 1, nsize_aer(itype)
1384
1385         dumvol = dryvolyy(n,itype)
1386         if (min(drymassyy(n,itype),dumvol) .le. smallmassbb) then
1387             dumvol = drymassyy(n,itype)/densdefault(itype)
1388             dumnum = drymassyy(n,itype)/(volumcen_stmp(n,itype)*densdefault(itype))
1389             delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n,itype)
1390             rnumbyy(n,itype) = dumnum
1391             adrydens_tmp(n,itype) = densdefault(itype)/fact_apdens
1392             awetdens_tmp(n,itype) = densdefault(itype)/fact_apdens
1393             adrydpav_tmp(n,itype) = dcen_stmp(n,itype)/fact_apdiam
1394             awetdpav_tmp(n,itype) = dcen_stmp(n,itype)/fact_apdiam
1395         else
1396             dumnum = rnumbyy(n,itype)
1397             dumnum_at_dhi = dumvol/volumhi_stmp(n,itype)
1398             dumnum_at_dlo = dumvol/volumlo_stmp(n,itype)
1399             dumnum = max( dumnum_at_dhi, min( dumnum_at_dlo, dumnum ) )
1400             delta_numb_conform2 = delta_numb_conform2 + dumnum - rnumbyy(n,itype)
1401             rnumbyy(n,itype) = dumnum
1402             adrydens_tmp(n,itype) = drymassyy(n,itype)/dumvol/fact_apdens
1403             dumvol1p = dumvol/dumnum
1404             adrydpav_tmp(n,itype) = ((dumvol1p*sixoverpi)**third) / fact_apdiam
1405             awetdens_tmp(n,itype) = wetmassyy(n,itype)/wetvolyy(n,itype)/fact_apdens
1406             dumvol1p = wetvolyy(n,itype)/dumnum
1407             awetdpav_tmp(n,itype) = min( 100.0_r8*dcen_stmp(n,itype),   &
1408                         (dumvol1p*sixoverpi)**third ) / fact_apdiam
1409         end if
1410         adryqmas_tmp(n,itype) = drymassyy(n,itype)/fact_apmassmr
1411
1412         do ll = 1, ncomp_plustracer_aer(itype)
1413             l = massptr_aer(ll,n,itype,iphase)
1414             rbox(l) = rmassyy(ll,n,itype)/fact_apmassmr
1415         end do
1416         l = 0
1417         if (iphase .eq. ai_phase) then
1418             l = waterptr_aer(n,itype)
1419             if (l .gt. 0) rbox(l) = rmassyy(llwater,n,itype)/fact_apmassmr
1420             l = hyswptr_aer(n,itype)
1421             if (l .gt. 0) rbox(l) = rmassyy(llhysw,n,itype)/fact_apmassmr
1422         end if
1423         rbox(numptr_aer(n,itype,iphase)) = rnumbyy(n,itype)/fact_apnumbmr
1424
1425 3900    continue   ! n=isize loop
1426 3950    continue   ! itype loop
1427
1428         delta_numb_conform1 = delta_numb_conform1 + delta_numb_conform2
1429
1430         call move_sect_3d_conserve_check(   &
1431           2, iphase_flag, iclm, jclm, k, m, iphase,   &
1432           method_movesect, idiag_movesect, llhysw, llwater,   &
1433           nnewsave, itypenewsave,    &
1434           densdefault, densh2o, smallmassaa, smallmassbb,   &
1435           fact_apmassmr, fact_apnumbmr,   &
1436           delta_water_conform1, delta_numb_conform1,   &
1437           rbox,   &
1438           drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx,   &
1439           drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1440
1441
1442 !   diagnostic output
1443         if (idiag_movesect .lt. 70) goto 4900
1444
1445         do 4800 itype = 1, ntype_aer
1446
1447         msg = ' '
1448         call peg_message( lunout, msg )
1449         write(msg,97005) itype
1450         call peg_message( lunout, msg )
1451
1452         ndum = 0
1453         do n = 1, nsize_aer(itype)
1454             if (nnewsave(1,n,itype)+nnewsave(2,n,itype) .ne. n) ndum = ndum + 1
1455         end do
1456
1457         dumch4 = 'NONE'
1458         if (ndum .gt. 0) dumch4 = 'SOME'
1459         msg = ' '
1460         call peg_message( lunout, msg )
1461         write(msg,97010) dumch4, iclm, jclm, k, m, ndum
1462         call peg_message( lunout, msg )
1463         do jjb = 1, nsize_aer(itype), 10
1464             jjc = min( jjb+9, nsize_aer(itype) )
1465             write(msg,97011) (nnewsave(1,n,itype), nnewsave(2,n,itype), n=jjb,jjc)
1466             call peg_message( lunout, msg )
1467         end do
1468
1469 !       write(msg,97020) 'rnumbold', (rnumbxx(n,itype), n=1,nsize_aer(itype))
1470 !       call peg_message( lunout, msg )
1471 !       write(msg,97020) 'rnumbnew', (rnumbyy(n,itype), n=1,nsize_aer(itype))
1472 !       call peg_message( lunout, msg )
1473
1474 !       write(msg,97020) 'drvolold', (dryvolxx(n,itype), n=1,nsize_aer(itype))
1475 !       call peg_message( lunout, msg )
1476 !       write(msg,97020) 'drvolnew', (dryvolyy(n,itype), n=1,nsize_aer(itype))
1477 !       call peg_message( lunout, msg )
1478
1479         dumbot = log( volumhi_stmp(1,itype)/volumlo_stmp(1,itype) )
1480         do n = 1, nsize_aer(itype)
1481             dumpp(n) = -9.99
1482             dumxx(n) = -9.99
1483             dumyy(n) = -9.99
1484             if ( (drydenspp(n,itype) .gt. 0.5) .and.   &
1485                  (drymasspp(n,itype) .gt. smallmassaa) ) then
1486                 dumvol = drymasspp(n,itype)/drydenspp(n,itype)
1487                 if ((rnumbpp(n,itype) .ge. 1.0e-35) .and.   &
1488                     (dumvol .ge. 1.0e-35)) then
1489                     dumvol1p = dumvol/rnumbpp(n,itype)
1490                     dumpp(n) = 1.0 + log(dumvol1p/volumlo_stmp(1,itype))/dumbot
1491                 end if
1492             end if
1493             if ((rnumbxx(n,itype) .ge. 1.0e-35) .and.   &
1494                 (dryvolxx(n,itype) .ge. 1.0e-35)) then
1495                 dumvol1p = dryvolxx(n,itype)/rnumbxx(n,itype)
1496                 dumxx(n) = 1.0 + log(dumvol1p/volumlo_stmp(1,itype))/dumbot
1497             end if
1498             if ((rnumbyy(n,itype) .ge. 1.0e-35) .and.   &
1499                 (dryvolyy(n,itype) .ge. 1.0e-35)) then
1500                 dumvol1p = dryvolyy(n,itype)/rnumbyy(n,itype)
1501                 dumyy(n) = 1.0 + log(dumvol1p/volumlo_stmp(1,itype))/dumbot
1502             end if
1503         end do
1504
1505 !       write(msg,97030) 'lnvolold', (dumxx(n,itype), n=1,nsize_aer(itype))
1506 !       call peg_message( lunout, msg )
1507 !       write(msg,97030) 'lnvolnew', (dumyy(n,itype), n=1,nsize_aer(itype))
1508 !       call peg_message( lunout, msg )
1509
1510         do jja = 1, 7
1511 !       do jja = 1, 13
1512             if      (jja .eq. 1) then
1513                 dumch8 = 'rnumbold'
1514                 dumout(:) = rnumbxx(:,itype)
1515             else if (jja .eq. 2) then
1516                 dumch8 = 'rnumbnew'
1517                 dumout(:) = rnumbyy(:,itype)
1518             else if (jja .eq. 3) then
1519                 dumch8 = 'drvolold'
1520                 dumout(:) = dryvolxx(:,itype)
1521             else if (jja .eq. 4) then
1522                 dumch8 = 'drvolnew'
1523                 dumout(:) = dryvolyy(:,itype)
1524             else if (jja .eq. 5) then
1525                 dumch8 = 'lnvolold'
1526                 dumout(:) = dumxx(:)
1527             else if (jja .eq. 6) then
1528                 dumch8 = 'lnvolnew'
1529                 dumout(:) = dumyy(:)
1530             else if (jja .eq. 7) then
1531                 dumch8 = 'lnvolpre'
1532                 dumout(:) = dumpp(:)
1533             else if (jja .eq. 8) then
1534                 dumch8 = 'wtvolold'
1535                 dumout(:) = wetvolxx(:,itype)
1536             else if (jja .eq. 9) then
1537                 dumch8 = 'wtvolnew'
1538                 dumout(:) = wetvolyy(:,itype)
1539             else if (jja .eq. 10) then
1540                 dumch8 = 'wtmasold'
1541                 dumout(:) = wetmassxx(:,itype)
1542             else if (jja .eq. 11) then
1543                 dumch8 = 'wtmasnew'
1544                 dumout(:) = wetmassyy(:,itype)
1545             else if (jja .eq. 12) then
1546                 dumch8 = 'drmasold'
1547                 dumout(:) = drymassxx(:,itype)
1548             else if (jja .eq. 13) then
1549                 dumch8 = 'drmasnew'
1550                 dumout(:) = drymassyy(:,itype)
1551             end if
1552             do jjb = 1, nsize_aer(itype), 10
1553                 jjc = min( jjb+9, nsize_aer(itype) )
1554                 if ((jja .le. 4) .or. (jja .ge. 8)) then
1555                     if (jjc == nsize_aer(itype)) then
1556                         write(msg,97020) dumch8, (dumout(n), n=jjb,jjc), &
1557                             sum( dumout(1:jjc) )
1558                     else
1559                         write(msg,97020) dumch8, (dumout(n), n=jjb,jjc)
1560                     end if
1561                 else
1562                     write(msg,97030) dumch8, (dumout(n), n=jjb,jjc)
1563                 end if
1564                 call peg_message( lunout, msg )
1565                 dumch8 = ' '
1566             end do
1567         end do
1568
1569 97005   format( 'movesectapply - itype =', i6 )
1570 97010   format( 'movesectapply', a, 4i3, i6 )
1571 !97011  format( 5x, 10(3x,2i3) )
1572 !97020  format( a, 1p, 10e9.1 )
1573 !97030  format( a,     10f9.3 )
1574 97011   format( 5x, 10(5x,2i3) )
1575 97020   format( a, 1p, 11e11.3 )
1576 97030   format( a,     10f11.5 )
1577
1578 4800    continue   ! itype loop
1579
1580 4900    continue
1581         return
1582         end subroutine move_sect_3d_apply_moves
1583
1584
1585 !-----------------------------------------------------------------------
1586         subroutine move_sect_3d_conserve_check( ipass,   &
1587           iphase_flag, iclm, jclm, k, m, iphase,   &
1588           method_movesect, idiag_movesect, llhysw, llwater,   &
1589           nnewsave, itypenewsave,    &
1590           densdefault, densh2o, smallmassaa, smallmassbb,   &
1591           fact_apmassmr, fact_apnumbmr,   &
1592           delta_water_conform1, delta_numb_conform1,   &
1593           rbox,   &
1594           drymassxx, dryvolxx, rmassxx, rnumbxx, wetmassxx, wetvolxx,   &
1595           drymassyy, dryvolyy, rmassyy, rnumbyy, wetmassyy, wetvolyy )
1596 !
1597 !   routine checks for conservation of number, mass, and volume
1598 !       by the move_sect_3d algorithm
1599 !   note - move_sections does separate conserve check for each type
1600 !        - move_sect_3d  does a single conserve check for all types,
1601 !          because number & mass are transferred among types
1602 !
1603 !   ipass = 1
1604 !       initialize all thesum(jj,ll) to zero
1605 !       computes thesum(1,ll) from rmassxx, rnumbxx, ...
1606 !       computes thesum(2,ll) from rmassyy, rnumbyy, ...
1607 !       compares thesum(1,ll) with thesum(2,ll)
1608 !       computes thesum(3,ll) from rbox before section movement
1609 !   ipass = 2
1610 !       computes thesum(4,ll) from rbox after  section movement
1611 !       compares thesum(3,ll) with thesum(4,ll)
1612 !
1613 !   currently only implemented for condensational growth (iphase_flag=1)
1614 !
1615         use module_data_mosaic_asecthp, only:  &
1616             maxd_acomp, maxd_asize, maxd_atype,   &
1617             nsize_aer, ntype_aer, ncomp_plustracer_aer, ai_phase,   &
1618             massptr_aer, numptr_aer, waterptr_aer, hyswptr_aer,   &
1619             name_aer
1620
1621         implicit none
1622
1623 !   subr arguments
1624         integer ipass, iphase_flag, iclm, jclm, iphase, k,   &
1625           m, method_movesect, idiag_movesect, llhysw, llwater
1626         integer nnewsave(2,maxd_asize,maxd_atype)
1627         integer itypenewsave(2,maxd_asize,maxd_atype)
1628         real(r8) densh2o, smallmassaa, smallmassbb
1629         real(r8) fact_apmassmr, fact_apnumbmr
1630         real(r8) delta_water_conform1, delta_numb_conform1
1631         real(r8) densdefault(maxd_atype)
1632         real(r8) drymassxx(maxd_asize,maxd_atype), drymassyy(maxd_asize,maxd_atype)
1633         real(r8) dryvolxx(maxd_asize,maxd_atype), dryvolyy(maxd_asize,maxd_atype)
1634         real(r8) rbox(ntot_used)
1635         real(r8) rmassxx(maxd_acomp+2,maxd_asize,maxd_atype),   &
1636              rmassyy(maxd_acomp+2,maxd_asize,maxd_atype)
1637         real(r8) rnumbxx(maxd_asize,maxd_atype), rnumbyy(maxd_asize,maxd_atype)
1638         real(r8) wetvolxx(maxd_asize,maxd_atype), wetvolyy(maxd_asize,maxd_atype)
1639         real(r8) wetmassxx(maxd_asize,maxd_atype), wetmassyy(maxd_asize,maxd_atype)
1640
1641 !   local variables
1642         integer ii, itype, itmpa, itmpb, jj, jtype, l, ll, llworst, llworstb, n
1643         integer nerr, nerrmax
1644         save nerr, nerrmax
1645         data nerr, nerrmax / 0, 999 /
1646
1647         real(r8) dumbot, dumtop, dumtoler, dumerr, dumworst, dumworstb
1648         real(r8) duma, dumb, dumc, dume
1649         real(r8) dumerrsv(maxd_acomp+7)
1650         real(r8) thesum(4,maxd_acomp+7)
1651
1652         character*8 dumname(maxd_acomp+7)
1653         character*160 fmtaa, msg
1654
1655
1656         if (ipass .eq. 2) goto 2000
1657
1658         itype = 1
1659         do ll = 1, ncomp_plustracer_aer(itype)+7
1660         do jj = 1, 4
1661             thesum(jj,ll) = 0.0
1662         end do
1663         end do
1664
1665         do itype = 1, ntype_aer
1666         do n = 1, nsize_aer(itype)
1667             do ll = 1, ncomp_plustracer_aer(itype)+2
1668                 thesum(1,ll) = thesum(1,ll) + rmassxx(ll,n,itype)
1669                 thesum(2,ll) = thesum(2,ll) + rmassyy(ll,n,itype)
1670             end do
1671             ll = ncomp_plustracer_aer(itype)+3
1672             thesum(1,ll) = thesum(1,ll) + rnumbxx(n,itype)
1673             thesum(2,ll) = thesum(2,ll) + rnumbyy(n,itype)
1674             ll = ncomp_plustracer_aer(itype)+4
1675             thesum(1,ll) = thesum(1,ll) + drymassxx(n,itype)
1676             thesum(2,ll) = thesum(2,ll) + drymassyy(n,itype)
1677             ll = ncomp_plustracer_aer(itype)+5
1678             thesum(1,ll) = thesum(1,ll) + dryvolxx(n,itype)
1679             thesum(2,ll) = thesum(2,ll) + dryvolyy(n,itype)
1680             ll = ncomp_plustracer_aer(itype)+6
1681             thesum(1,ll) = thesum(1,ll) + wetmassxx(n,itype)
1682             thesum(2,ll) = thesum(2,ll) + wetmassyy(n,itype)
1683             ll = ncomp_plustracer_aer(itype)+7
1684             thesum(1,ll) = thesum(1,ll) + wetvolxx(n,itype)
1685             thesum(2,ll) = thesum(2,ll) + wetvolyy(n,itype)
1686         end do
1687         end do
1688
1689
1690 !
1691 !   calc sum over bins for each species
1692 !   for water, account for loss in initial conform (delta_water_conform1)
1693 !
1694 2000    continue
1695         do itype = 1, ntype_aer
1696         do n = 1, nsize_aer(itype)
1697             do ll = 1, ncomp_plustracer_aer(itype)+3
1698                 duma = fact_apmassmr
1699                 if (ll .le. ncomp_plustracer_aer(itype)) then
1700                     l = massptr_aer(ll,n,itype,iphase)
1701                 else if (ll .eq. ncomp_plustracer_aer(itype)+1) then
1702                     l = 0
1703                     if (iphase .eq. ai_phase) l = hyswptr_aer(n,itype)
1704                 else if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1705                     l = 0
1706                     if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1707                 else
1708                     l = numptr_aer(n,itype,iphase)
1709                     duma = fact_apnumbmr
1710                 end if
1711                 if (l .gt. 0)   &
1712                     thesum(ipass+2,ll) = thesum(ipass+2,ll) + rbox(l)*duma
1713             end do
1714         end do
1715         end do
1716
1717         itype = 1
1718         if (ipass .eq. 2) then
1719             ll = ncomp_plustracer_aer(itype)+2
1720             thesum(3,ll) = thesum(3,ll) + delta_water_conform1
1721             ll = ncomp_plustracer_aer(itype)+3
1722             thesum(3,ll) = thesum(3,ll) + delta_numb_conform1
1723         end if
1724
1725 !
1726 !   now compare either sum1-sum2 or sum3-sum4
1727 !       on ipass=1, jj=1, so compare sum1 & sum2
1728 !       on ipass=2, jj=3, so compare sum3 & sum4
1729 !
1730         itype = 1
1731         do ll = 1, ncomp_plustracer_aer(itype)+7
1732             dumname(ll) = ' '
1733             write(dumname(ll),'(i4.4)') ll
1734 !           if (ll .le. ncomp_plustracer_aer(itype)) dumname(ll) =   &
1735 !                       species(massptr_aer(ll,1,itype,iphase))(1:4)
1736             if (ll .le. ncomp_plustracer_aer(itype))   dumname(ll) = name_aer(ll,1)(1:4)
1737             if (ll .eq. ncomp_plustracer_aer(itype)+1) dumname(ll) = 'hysw'
1738             if (ll .eq. ncomp_plustracer_aer(itype)+2) dumname(ll) = 'watr'
1739             if (ll .eq. ncomp_plustracer_aer(itype)+3) dumname(ll) = 'numb'
1740             if (ll .eq. ncomp_plustracer_aer(itype)+4) dumname(ll) = 'drymass'
1741             if (ll .eq. ncomp_plustracer_aer(itype)+5) dumname(ll) = 'dryvol'
1742             if (ll .eq. ncomp_plustracer_aer(itype)+6) dumname(ll) = 'wetmass'
1743             if (ll .eq. ncomp_plustracer_aer(itype)+7) dumname(ll) = 'wetvol'
1744         end do
1745
1746         jj = 2*ipass - 1
1747         dumworst = 0.0
1748         dumworstb = 0.0
1749         dumerrsv(:) = 0.0
1750         llworst = 0
1751         llworstb = 0
1752         do ll = 1, ncomp_plustracer_aer(itype)+7
1753             dumtop = thesum(jj+1,ll) - thesum(jj,ll)
1754             dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35_r8 )
1755             dumerr = dumtop/dumbot
1756
1757 ! rce 21-jul-2006 - encountered some cases when delta_*_conform1 is negative
1758 !   and large in magnitude relative to thesum, which causes the mass
1759 !   conservation to be less accurate due to roundoff
1760 ! following section recomputes relative error with delta_*_conform1
1761 !   added onto each of thesum
1762 ! also increased dumtoler slightly
1763             if (ipass .eq. 2) then
1764                 dumc = 1.0
1765                 if (ll .eq. ncomp_plustracer_aer(itype)+2) then
1766                     dumc = delta_water_conform1
1767                 else if (ll .eq. ncomp_plustracer_aer(itype)+3) then
1768                     dumc = delta_numb_conform1
1769                 end if
1770                 if (dumc .lt. 0.0) then
1771                     duma = thesum(3,ll) - dumc
1772                     dumb = thesum(4,ll) - dumc
1773                     dumtop = dumb - duma
1774                     dumbot = max( abs(duma), abs(dumb), 1.0e-35_r8 )
1775                     dume = dumtop/dumbot
1776                     if (abs(dume) .lt. abs(dumerr)) dumerr = dume
1777                 end if
1778             end if
1779             dumerrsv(ll) = dumerr
1780
1781             if (abs(dumerr) .gt. abs(dumworst)) then
1782                 llworstb = llworst
1783                 dumworstb = dumworst
1784                 llworst = ll
1785                 dumworst = dumerr
1786             end if
1787         end do
1788
1789         dumtoler = 1.0e-6
1790         if (abs(dumworst) .gt. dumtoler) then
1791             nerr = nerr + 1
1792             if (nerr .le. nerrmax) then
1793                 msg = ' '
1794                 call peg_message( lunout, msg )
1795                 write(msg,97110) iclm, jclm, k, m, ipass, llworst
1796                 call peg_message( lunout, msg )
1797                 write(msg,'(a)') '    jtype, ii, nnew(ii,1:n,jtype)'
1798                 call peg_message( lunout, msg )
1799
1800                 do jtype = 1, ntype_aer
1801                 do ii = 1, 2
1802                 do itmpa = 1, nsize_aer(jtype), 24
1803                     itmpb = min( itmpa+23, nsize_aer(jtype) )
1804                     if (itmpa == 1) then
1805                         fmtaa = '(i4,i2,2x,24i4)'
1806                         write(msg,fmtaa) jtype, ii, &
1807                             (nnewsave(ii,n,jtype), n=itmpa,itmpb)
1808                     else
1809                         fmtaa = '(8x,24i4)'
1810                         write(msg,fmtaa) (nnewsave(ii,n,jtype), n=itmpa,itmpb)
1811                     end if
1812                     call peg_message( lunout, msg )
1813                 end do
1814                 end do
1815                 end do
1816
1817                 ll = llworst
1818                 if (ll .eq. 0) ll = ncomp_plustracer_aer(itype)+2
1819                 write(msg,97130) 'name/relerrA/thesum', jj, '/thesum', jj+1,   &
1820                         dumname(ll), dumworst, thesum(jj,ll), thesum(jj+1,ll)
1821                 call peg_message( lunout, msg )
1822
1823                 if ( (ll .eq. ncomp_plustracer_aer(itype)+3) .and.   &
1824                      (abs(dumworstb) .gt. dumtoler) ) then
1825                     ll = max( 1, llworstb )
1826                     dumtop = thesum(jj+1,ll) - thesum(jj,ll)
1827                     dumbot = max( abs(thesum(jj,ll)), abs(thesum(jj+1,ll)), 1.0e-35_r8 )
1828                     dumerr = dumtop/dumbot
1829                     write(msg,97130) 'name/relerrB/thesum', jj, '/thesum', jj+1,   &
1830                         dumname(ll), dumerr, thesum(jj,ll), thesum(jj+1,ll)
1831                     call peg_message( lunout, msg )
1832                 end if
1833             end if
1834         end if
1835
1836 !   following added on 10-mar-2010
1837         do ll = 1, ncomp_plustracer_aer(itype)+7
1838 !           exit   ! this deactivates it
1839             if (abs(dumerrsv(ll)) .le. dumtoler) cycle
1840             nerr = nerr + 1
1841             if (nerr .le. nerrmax) then
1842                 write(msg,97130) 'name/relerrC/thesum', jj, '/thesum', jj+1,   &
1843                         dumname(ll), dumworst, thesum(jj,ll), thesum(jj+1,ll)
1844                 call peg_message( lunout, msg )
1845             end if
1846         end do ! ll
1847
1848 97110   format( 'movesect conserve ERROR - i/j/k/m/pass/llworst',   &
1849                 4i3, 2x, 2i3 )
1850 97120   format( 64i3 )
1851 97130   format( a, i1, a, i1, 2x, a, 1p, 3e16.7 )
1852
1853
1854         return
1855         end subroutine move_sect_3d_conserve_check
1856
1857
1858 !-----------------------------------------------------------------------
1859         subroutine move_sect_3d_checkptrs( iphase_flag, iclm, jclm, k, m )
1860 !
1861 !   checks for valid number and water pointers
1862 !
1863         use module_data_mosaic_asecthp, only:  &
1864             ai_phase, ntype_aer, nsize_aer, nphase_aer,   &
1865             ncomp_aer, ncomp_plustracer_aer,   &
1866             numptr_aer, waterptr_aer, mastercompptr_aer
1867
1868         implicit none
1869
1870 !   subr parameters
1871         integer iphase_flag, iclm, jclm, k, ll, m
1872
1873 !   local variables
1874         integer l, itype, iphase, n, ndum
1875         character*160 msg
1876
1877         do 1900 itype = 1, ntype_aer
1878         do 1800 iphase = 1, nphase_aer
1879
1880         ndum = 0
1881         do n = 1, nsize_aer(itype)
1882             l = numptr_aer(n,itype,iphase)
1883             if ((l .lt. 1) .or. (l .gt. ntot_used)) then
1884                 msg = '*** subr move_sect_3d error - ' //   &
1885                         'numptr_amode not defined'
1886                 call peg_message( lunerr, msg )
1887                 write(msg,9030) 'mode, numptr =', n, l
1888                 call peg_message( lunerr, msg )
1889                 write(msg,9030) 'iphase, itype =', iphase, itype
1890                 call peg_message( lunerr, msg )
1891                 call peg_error_fatal( lunerr, msg )
1892             end if
1893 !       checks involving nspec_amode and nspec_amode_nontracer
1894 !       being the same for all sections are no longer needed
1895             l = 0
1896             if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
1897             if ((l .ge. 1) .and. (l .le. ntot_used)) ndum = ndum + 1
1898         end do
1899         if ((ndum .ne. 0) .and. (ndum .ne. nsize_aer(itype))) then
1900             msg = '*** subr move_sect_3d error - ' //   &
1901                 'waterptr_aer must be on/off for all modes'
1902             call peg_message( lunerr, msg )
1903             write(msg,9030) 'iphase, itype =', iphase, itype
1904             call peg_message( lunerr, msg )
1905             call peg_error_fatal( lunerr, msg )
1906         end if
1907 9030    format( a, 6(1x,i6) )
1908
1909         msg = ' '
1910         if (nsize_aer(itype) /= nsize_aer(1)) then
1911             write(msg,9030) 'nsize diffs', &
1912                 1, itype, nsize_aer(1), nsize_aer(itype)
1913         else if (ncomp_aer(itype) /= ncomp_aer(1)) then
1914             write(msg,9030) 'ncomp diffs', &
1915                 1, itype, ncomp_aer(1), ncomp_aer(itype)
1916         else if (ncomp_plustracer_aer(itype) /= ncomp_plustracer_aer(1)) then
1917             write(msg,9030) 'ncomp_plustracer diffs', &
1918                 1, itype, ncomp_plustracer_aer(1), ncomp_plustracer_aer(itype)
1919         end if
1920         do ll = 1, ncomp_plustracer_aer(1)
1921             if (mastercompptr_aer(ll,itype) /= mastercompptr_aer(ll,1)) then
1922                 write(msg,9030) 'mastercompptrcomp diffs', &
1923                     1, itype, ll, mastercompptr_aer(ll,1), mastercompptr_aer(ll,itype)
1924             end if
1925         end do
1926         if (msg /= ' ') then
1927             call peg_message( lunerr, msg )
1928             call peg_error_fatal( lunerr, msg )
1929         end if
1930
1931 1800    continue
1932 1900    continue
1933
1934         return
1935         end subroutine move_sect_3d_checkptrs
1936
1937
1938 !-----------------------------------------------------------------------
1939         subroutine test_move_sect_3d( iphase_flag_inp, iclm, jclm, k, m, &
1940              it_mosaic, mmovesect_flag1, idiag_movesect )
1941 !
1942 !   routine runs tests on move_sect_3d, using a matrix of
1943 !       pregrow and aftgrow masses for each section
1944 !
1945         use module_data_mosaic_asecthp, only:  &
1946             maxd_atype, maxd_asize,   &
1947             ntype_aer, nsize_aer, ncomp_plustracer_aer,   &
1948             ai_phase, cw_phase,   &
1949             massptr_aer, numptr_aer, waterptr_aer,   &
1950             dens_aer, volumlo_sect
1951
1952         implicit none
1953
1954 !   subr arguments
1955         integer, intent(in) :: iphase_flag_inp, iclm, jclm, k, m, it_mosaic, &
1956                                mmovesect_flag1, idiag_movesect
1957
1958 !   local variables
1959         integer iphase_flag, ii, iphase, itype, jj,   &
1960           l, ll, n, nn
1961         integer, save :: ientryno = 0
1962
1963         real(r8) dumnumb, dumvolpre, dumvolaft
1964         real(r8) adrydens_tmp(maxd_asize,maxd_atype),  awetdens_tmp(maxd_asize,maxd_atype)
1965         real(r8) adrydpav_tmp(maxd_asize,maxd_atype), awetdpav_tmp(maxd_asize,maxd_atype)
1966         real(r8) adryqmas_tmp(maxd_asize,maxd_atype)
1967         real(r8) fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam
1968         real(r8) tmp_drymass_pregrow(maxd_asize,maxd_atype)
1969         real(r8) tmp_drymass_aftgrow(maxd_asize,maxd_atype)
1970         real(r8) tmp_drydens_pregrow(maxd_asize,maxd_atype)
1971         real(r8) tmp_drydens_aftgrow(maxd_asize,maxd_atype)
1972         real(r8) rbox(ntot_used)
1973
1974         character*160 msg
1975
1976         integer maxvolfactpre, maxvolfactaft
1977         parameter (maxvolfactpre=15, maxvolfactaft=23)
1978
1979         real(r8) dumvolfactpre(maxvolfactpre)
1980         data dumvolfactpre /   &
1981         2.0, 0.0, 1.0e-20, 0.5, 0.9,   &
1982         1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0,   &
1983         8.1, 16.0 /
1984
1985         real(r8) dumvolfactaft(maxvolfactaft)
1986         data dumvolfactaft /   &
1987         4.0, 0.0, 1.0e-20, 0.01, 0.02, 0.05, 0.1, 0.5, 0.9,   &
1988         1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0,   &
1989         8.1, 16.0, 32., 64., 128., 256. /
1990 !       7.9, 0.0, 1.0e-20, 0.01, 0.02, 0.05, 0.1, 0.5, 0.9,   &
1991 !       1.0, 1.01, 1.1, 2.0, 4.0, 7.9, 7.99, 8.0,   &
1992 !       8.1, 16.0, 32., 64., 128., 256. /
1993
1994
1995 !
1996 !   check for valid inputs
1997 !   and first entry
1998 !
1999         if (mmovesect_flag1 .le. 0) return
2000         if (nsize_aer(1) .le. 0) return
2001
2002         ientryno = ientryno + 1
2003         if (ientryno .gt. 1) return
2004
2005         write(*,*) '*** doing move_sect_3d_checkptrs ***'
2006         call move_sect_3d_checkptrs( -1, iclm, jclm, k, m )
2007
2008         if (idiag_movesect .ne. 80) return
2009
2010         write(*,*) '*** executing test_move_sect_3d ***'
2011
2012
2013 !
2014 !   make test calls to move_sect_3d
2015 !
2016         do 3900 iphase_flag = 1, 1
2017
2018         iphase = ai_phase
2019         if (iabs(iphase_flag) .eq. 2) iphase = cw_phase
2020
2021         do 3800 itype = 1, ntype_aer
2022
2023         do 2900 nn = 1, nsize_aer(itype)
2024
2025         do 2800 ii = 1, maxvolfactpre
2026
2027         do 2700 jj = 1, maxvolfactaft
2028
2029 !----------------------------------------------------------
2030 !   following lines produce a small subset of the ii/jj/nn tests
2031 !
2032 !   following limits ii/jj to 1 and "matching dumvolfact"
2033 !       if (dumvolfactpre(ii) /= dumvolfactpre(1)) goto 2700
2034 !       if (dumvolfactaft(jj) /= dumvolfactaft(1)) goto 2700
2035 !   following limits dumvolfactpre to 4.0
2036 !       if (abs(dumvolfactpre(ii) - 4.0) > 1.0e-4) goto 2800
2037 !   following limits nn
2038 !       if (nn /= 4) goto 2900
2039 !----------------------------------------------------------
2040
2041 !   zero out rbox and dryxxx_yyygrow arrays
2042         do n = 1, nsize_aer(itype)
2043             do ll = 1, ncomp_plustracer_aer(itype)
2044                 rbox(massptr_aer(ll,n,itype,iphase)) = 0.0
2045             end do
2046             l = 0
2047             if (iphase .eq. ai_phase) l = waterptr_aer(n,itype)
2048             if (l .gt. 0) rbox(l) = 0.0
2049             l = numptr_aer(n,itype,iphase)
2050             if (l .gt. 0) rbox(l) = 0.0
2051             tmp_drymass_pregrow(n,itype) = 0.0
2052             tmp_drymass_aftgrow(n,itype) = 0.0
2053             tmp_drydens_pregrow(n,itype) = -1.0
2054             tmp_drydens_aftgrow(n,itype) = -1.0
2055         end do
2056
2057 !   fill in values for section nn
2058         n = nn
2059         dumnumb = 1.0e7/mw_air   ! 1.0e7 #/mol-air = 3.45e5 #/g-air = 345 #/mg-air
2060         rbox(numptr_aer(n,itype,iphase)) = dumnumb
2061         ll = 1
2062         l = massptr_aer(ll,n,itype,iphase)
2063
2064         dumvolpre = volumlo_sect(n,itype)*dumvolfactpre(ii)*dumnumb
2065         tmp_drydens_pregrow(n,itype) = dens_aer(ll,itype)
2066         tmp_drymass_pregrow(n,itype) = dumvolpre*tmp_drydens_pregrow(n,itype)
2067         if (ii .eq. 1) tmp_drydens_pregrow(n,itype) = -1.0
2068
2069         dumvolaft = volumlo_sect(n,itype)*dumvolfactaft(jj)*dumnumb
2070         tmp_drydens_aftgrow(n,itype) = dens_aer(ll,itype)
2071         tmp_drymass_aftgrow(n,itype) = dumvolaft*tmp_drydens_aftgrow(n,itype)
2072         if (jj .eq. 1) tmp_drydens_aftgrow(n,itype) = -1.0
2073
2074         rbox(l) = tmp_drymass_aftgrow(n,itype)   ! g/g-air
2075
2076         msg = ' '
2077         call peg_message( lunout, msg )
2078         write(msg,98010) nn, ii, jj
2079         call peg_message( lunout, msg )
2080         write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
2081         call peg_message( lunout, msg )
2082
2083         fact_apmassmr = 1.0   ! units are g-AP/g-air already
2084         fact_apnumbmr = 1.0   ! units are    #/g-air already
2085         fact_apdens = 1.0     ! units are    g/cm3   already
2086         fact_apdiam = 1.0     ! units are   cm       already
2087
2088 !   make test call to move_sect_3d
2089         call move_sect_3d_x1( iphase_flag, iclm, jclm, k, m, rbox,   &
2090           fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam,   &
2091           tmp_drydens_aftgrow, tmp_drydens_pregrow,   &
2092           tmp_drymass_aftgrow, tmp_drymass_pregrow,   &
2093           adrydens_tmp, awetdens_tmp, adrydpav_tmp, awetdpav_tmp,   &
2094           adryqmas_tmp, it_mosaic, mmovesect_flag1, idiag_movesect )
2095
2096         msg = ' '
2097         call peg_message( lunout, msg )
2098         write(msg,98010) nn, ii, jj
2099         call peg_message( lunout, msg )
2100         write(msg,98011) dumvolfactpre(ii), dumvolfactaft(jj)
2101         call peg_message( lunout, msg )
2102
2103 2700    continue
2104 2800    continue
2105 2900    continue
2106
2107 3800    continue
2108 3900    continue
2109
2110 98010   format( 'test_move_sect_3d output - nn, ii, jj =', 3i3 )
2111 98011   format( 'volfactpre, volfactaft =', 1p, 2e12.4 )
2112
2113
2114         write(*,*) '*** leaving   test_move_sect_3d at end ***'
2115         return
2116         end subroutine test_move_sect_3d
2117
2118
2119 !-----------------------------------------------------------------------
2120
2121
2122         end module module_mosaic_movesect3d
2123