chem/module_mosaic_movesect1d.F

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