chem/module_mosaic_coag3d.F

   1         module module_mosaic_coag3d
   2
   3
   4
   5         use module_data_mosaic_kind
   6         use module_peg_util
   7
   8
   9
  10         implicit none
  11
  12
  13
  14         contains
  15
  16
  17
  18 !-----------------------------------------------------------------------
  19         subroutine mosaic_coag_3d_1box( istat_coag,   &
  20             idiagbb_in, itstep,   &
  21             dtcoag_in,   &
  22             temp_box, patm_box, rhoair_box, rbox,   &
  23             fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam,   &
  24             adrydens_box, awetdens_box,   &
  25             adrydpav_box, awetdpav_box, adryqmas_box, mcoag_flag1 )
  26 !
  27 !   calculates aerosol particle coagulation for a single grid box
  28 !       over timestep dtcoag_in
  29 !
  30 !   uses the single-type coagulation algorithm of jacobson (2002)
  31 !
  32         use module_data_mosaic_main, only:  ntot_used, pi, piover6, third
  33         !use module_data_mosaic_aero, only:  mcoag_flag1, ifreq_coag !BSINGH
  34         use module_data_mosaic_aero, only:  ifreq_coag !BSINGH
  35         use module_data_mosaic_asecthp, only: &
  36             ai_phase, hyswptr_aer, &
  37             lunerr, lunout, &
  38             maxd_acomp, maxd_asize, maxd_atype, &
  39             ncomp_aer, ncomp_plustracer_aer, nsize_aer, ntype_aer, &
  40             massptr_aer, numptr_aer, waterptr_aer, &
  41             dens_aer, dens_water_aer, &
  42             dcen_sect, dhi_sect, dlo_sect, &
  43             smallmassbb, &
  44             volumcen_sect, volumcut_sect, volumhi_sect, volumlo_sect
  45
  46
  47 !   subr arguments
  48         integer, intent(inout)  :: istat_coag  ! =0 if no problems
  49         integer, intent(in)     :: idiagbb_in  ! controls diagnostics
  50         integer, intent(in)     :: itstep      ! host code time step index
  51         integer, intent(in)     :: mcoag_flag1
  52         real(r8), intent(in)    :: dtcoag_in   ! time step (s)
  53         real(r8), intent(in)    :: temp_box    ! air temp (K)
  54         real(r8), intent(in)    :: patm_box    ! air pressure (atm)
  55         real(r8), intent(in)    :: rhoair_box    ! air density (g/cm^3)
  56
  57         real(r8), intent(inout) :: rbox(ntot_used)
  58 ! rbox - holds aerosol number and component mass mixing ratios
  59 !       units for number mr are such that (rbox*fact_apnumbmr) gives (#/g-air)
  60 !       units for mass   mr are such that (rbox*fact_apmassmr) gives (g/g-air)
  61
  62         real(r8), intent(in) :: fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam
  63 ! these are units conversion factors
  64 !   the module_data_mosaic_asecthp densities are such that
  65 !       dens_aer*fact_apdens gives (g/cm^3)
  66 !   the module_data_mosaic_asecthp diameters are such that
  67 !       d[lo,cen,hi]_sect*fact_apdiam gives (cm)
  68 !   the module_data_mosaic_asecthp volumes are such that
  69 !       volum[lo,cen,hi]_sect**fact_apdiam**3) gives (cm^3)
  70
  71         real(r8), intent(inout) :: adrydens_box(maxd_asize,maxd_atype)
  72         real(r8), intent(inout) :: awetdens_box(maxd_asize,maxd_atype)
  73         real(r8), intent(inout) :: adrydpav_box(maxd_asize,maxd_atype)
  74         real(r8), intent(inout) :: awetdpav_box(maxd_asize,maxd_atype)
  75         real(r8), intent(inout) :: adryqmas_box(maxd_asize,maxd_atype)
  76 ! adryqmas_box = aerosol total-dry-mass mixing ratio (g/mol-air)
  77 ! adrydens_box = aerosol dry density (units same as dens_aer)
  78 ! awetdens_box = aerosol wet density (units same as dens_aer)
  79 ! adrydpav_box = aerosol mean dry diameter (units same as dlo_sect)
  80 ! awetdpav_box = aerosol mean wet diameter (units same as dlo_sect)
  81 ! adryqmas_box = aerosol total-dry-mass mixing ratio (units same as rbox)
  82
  83
  84 !   local variables
  85         integer, parameter :: coag_method = 1
  86         integer, parameter :: ncomp_coag_maxd = maxd_acomp + 3
  87
  88         integer :: l, ll, lla, llb, llx
  89         integer :: isize, itype, iphase
  90         integer :: iconform_numb
  91         integer :: idiagbb, idiag_coag_onoff
  92         integer :: lunout_coag
  93         integer :: ncomp_coag, nsize, nsubstep_coag, ntau_coag
  94         integer,save :: ncountaa(10), ncountbb(0:ncomp_coag_maxd)
  95
  96         real(r8), parameter :: factsafe = 1.00001
  97
  98 ! units for local working varibles:
  99 !    size = cm
 100 !    density = g/cm^3
 101 !    number mixing ratio = #/g-air
 102 !       EXCEPT num_distrib = #/cm^3-air
 103 !    mass   mixing ratio = g/g-air
 104 !    volume mixing ratio = cm^3/g-air
 105 !
 106         real(r8) :: densdefault
 107         real(r8) :: dtcoag
 108         real(r8) :: fact_apdia3, fact_apvolumr
 109         real(r8) :: press_cgs   ! air pressure (dyne/cm2)
 110         real(r8) :: tmpa, tmpb
 111         real(r8) :: wwdens, wwdpav, wwmass, wwvolu
 112         real(r8) :: xxdens, xxdpav, xxmass, xxnumb, xxvolu
 113         real(r8) :: xxmasswet, xxvoluwet
 114
 115         real(r8) :: adryqvol_box(maxd_asize,maxd_atype) ! aerosol total-dry-volume mixing ratio (cm3/mol-air)
 116         real(r8) :: dpdry(maxd_asize,maxd_atype), dpwet(maxd_asize,maxd_atype)
 117         real(r8) :: densdry(maxd_asize,maxd_atype), denswet(maxd_asize,maxd_atype)
 118         real(r8) :: num_distrib(maxd_asize,maxd_atype)
 119         real(r8) :: sumnew(ncomp_coag_maxd), sumold(ncomp_coag_maxd)
 120         real(r8) :: sum_num(2), sum_vol(0:ncomp_coag_maxd,2)
 121         real(r8) :: vol_distrib(maxd_asize,maxd_atype,ncomp_coag_maxd)
 122         real(r8) :: vpcut(0:maxd_asize), vpdry(maxd_asize,maxd_atype)
 123         real(r8) :: xxvolubb(maxd_asize,maxd_atype)
 124
 125         character(len=120) :: msg
 126
 127
 128
 129         if (mcoag_flag1 <= 0) return
 130         if (ifreq_coag > 1) then
 131            if (mod(itstep,ifreq_coag) /= 0) return
 132            dtcoag = dtcoag_in*ifreq_coag
 133         else
 134            dtcoag = dtcoag_in
 135         end if
 136
 137         istat_coag = 0
 138
 139         lunout_coag = 6
 140         if (lunout .gt. 0) lunout_coag = lunout
 141         write(lunout_coag,'(/a,i10)') 'mosaic_coag_3d_1box - itstep ', itstep
 142
 143
 144 !   set variables that do not change
 145         idiagbb = idiagbb_in
 146         idiag_coag_onoff = +1
 147         if (idiagbb <= 0) idiag_coag_onoff = 0
 148
 149         iphase = ai_phase
 150         fact_apdia3 = fact_apdiam*fact_apdiam*fact_apdiam
 151         fact_apvolumr = fact_apmassmr/fact_apdens
 152
 153         nsubstep_coag = 1
 154         press_cgs = patm_box*1.01325e6_r8
 155
 156
 157         ncountaa(:) = 0   ! used for temporary diagnostics
 158         ncountbb(:) = 0
 159
 160 !   all types have same size cuts and same component species
 161 !   so nsize, ncomp_coag, densdefault and vpcut are same for all types
 162         itype = 1
 163         nsize = nsize_aer(itype)
 164         ncomp_coag = ncomp_plustracer_aer(itype) + 3
 165         densdefault = dens_aer(1,itype)*fact_apdens
 166         vpcut(0:nsize) = volumcut_sect(0:nsize,itype)*fact_apdia3
 167
 168         ncountaa(1) = ncountaa(1) + 1
 169
 170 !   do initial calculation of dry mass, volume, and density,
 171 !   and initial number conforming (as needed)
 172         vol_distrib(:,:,:) = 0.0
 173         sumold(:) = 0.0
 174
 175 itype_loop_aa: &
 176         do itype = 1, ntype_aer
 177
 178 isize_loop_aa: &
 179         do isize = 1, nsize
 180             l = numptr_aer(isize,itype,iphase)
 181             xxnumb = rbox(l)*fact_apnumbmr
 182             xxmass = adryqmas_box(isize,itype)*fact_apmassmr
 183             xxdens = adrydens_box(isize,itype)*fact_apdens
 184             iconform_numb = 1
 185
 186             if ( (xxdens .lt. 0.1) .or. (xxdens .gt. 20.0) ) then
 187 !   (exception) case of drydensity not valid
 188 !   so compute dry mass and volume from rbox
 189                 ncountaa(2) = ncountaa(2) + 1
 190                 xxvolu = 0.0
 191                 if (idiagbb .ge. 200)   &
 192                     write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-2a',   &
 193                     itstep, itype, isize, xxmass, xxvolu, xxdens
 194                 xxmass = 0.0
 195                 do ll = 1, ncomp_aer(itype)
 196                     l = massptr_aer(ll,isize,itype,iphase)
 197                     if (l .ge. 1) then
 198                         tmpa = max( 0.0_r8, rbox(l) )
 199                         xxmass = xxmass + tmpa
 200                         xxvolu = xxvolu + tmpa/dens_aer(ll,itype)
 201                     end if
 202                 end do
 203                 xxmass = xxmass*fact_apmassmr
 204                 xxvolu = xxvolu*fact_apvolumr
 205                 if (xxmass .gt. 0.99*smallmassbb) then
 206                     xxdens = xxmass/xxvolu
 207                     xxvolu = xxmass/xxdens
 208                     if (idiagbb .ge. 200)   &
 209                         write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-2b',   &
 210                         itstep, itype, isize, xxmass, xxvolu, xxdens
 211                 end if
 212             else
 213                 xxvolu = xxmass/xxdens
 214             end if
 215
 216             if (xxmass .le. 1.01*smallmassbb) then
 217 !   when drymass extremely small, use default density and bin center size,
 218 !   and zero out water
 219                 ncountaa(3) = ncountaa(3) + 1
 220                 xxdens = densdefault
 221                 xxvolu = xxmass/xxdens
 222                 xxnumb = xxmass/(volumcen_sect(isize,itype)*fact_apdia3*xxdens)
 223                 l = hyswptr_aer(isize,itype)
 224                 if (l .ge. 1) rbox(l) = 0.0
 225                 l = waterptr_aer(isize,itype)
 226                 if (l .ge. 1) rbox(l) = 0.0
 227                 iconform_numb = 0
 228                 if (idiagbb .ge. 210)   &
 229                     write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-2c',   &
 230                     itstep, itype, isize, xxmass, xxvolu, xxdens
 231             else
 232                 xxvolu = xxmass/xxdens
 233             end if
 234
 235 !   check for mean dry-size 'safely' within bounds, and conform number if not
 236             if (iconform_numb .gt. 0) then
 237                 if (xxnumb .gt.   &
 238                         xxvolu/(factsafe*volumlo_sect(isize,itype)*fact_apdia3)) then
 239                     ncountaa(4) = ncountaa(4) + 1
 240                     tmpa = xxnumb
 241                     xxnumb = xxvolu/(factsafe*volumlo_sect(isize,itype)*fact_apdia3)
 242                     if (idiagbb .ge. 200)   &
 243                         write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-4 ',   &
 244                         itstep, itype, isize, xxvolu, tmpa, xxnumb
 245                 else if (xxnumb .lt.   &
 246                         xxvolu*factsafe/(volumhi_sect(isize,itype)*fact_apdia3)) then
 247                     ncountaa(5) = ncountaa(5) + 1
 248                     tmpa = xxnumb
 249                     xxnumb = xxvolu*factsafe/(volumhi_sect(isize,itype)*fact_apdia3)
 250                     if (idiagbb .ge. 200)   &
 251                         write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-5 ',   &
 252                         itstep, itype, isize, xxvolu, tmpa, xxnumb
 253                     if (idiagbb .ge. 200)   &
 254                         write(lunout_coag,'(a,10x, 8x,1p,4e12.4)') 'coagcc-5 ',   &
 255                         dlo_sect(isize,itype), dlo_sect(isize,itype)*fact_apdiam, &
 256                         volumlo_sect(isize,itype), volumlo_sect(isize,itype)*fact_apdia3
 257                     if (idiagbb .ge. 200)   &
 258                         write(lunout_coag,'(a,10x, 8x,1p,4e12.4)') 'coagcc-5 ',   &
 259                         dhi_sect(isize,itype), dhi_sect(isize,itype)*fact_apdiam, &
 260                         volumhi_sect(isize,itype), volumhi_sect(isize,itype)*fact_apdia3
 261                 end if
 262             end if
 263
 264 !   load numb, mass, volu, dens back into mosaic_asecthp arrays
 265             l = numptr_aer(isize,itype,iphase)
 266             rbox(l) = xxnumb/fact_apnumbmr
 267             adrydens_box(isize,itype) = xxdens/fact_apdens
 268             adryqmas_box(isize,itype) = xxmass/fact_apmassmr
 269             adryqvol_box(isize,itype) = xxvolu/fact_apvolumr
 270
 271 !
 272 !   load number/mass/volume into working arrays
 273 !
 274 !   *** NOTE ***
 275 !     num_distrib must be (#/cm3)
 276 !     vol_distrib units do not matter - they can be either masses or volumes,
 277 !        and either mixing ratios or concentrations
 278 !        (e.g., g/cm3-air, ug/kg-air, cm3/cm3-air, cm3/kg-air, ...)
 279 !
 280             num_distrib(isize,itype) = xxnumb*rhoair_box
 281
 282             do ll = 1, ncomp_plustracer_aer(itype)
 283                 l = massptr_aer(ll,isize,itype,iphase)
 284                 tmpa = 0.0
 285                 if (l .ge. 1) tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr
 286                 vol_distrib(isize,itype,ll)   = tmpa
 287                 sumold(ll) = sumold(ll) + tmpa
 288             end do
 289
 290             do llx = 1, 3
 291                 ll = (ncomp_coag-3) + llx
 292                 tmpa = 0.0
 293                 if (llx .eq. 1) then
 294                     l = hyswptr_aer(isize,itype)
 295                     if (l .ge. 1) tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr
 296                 else if (llx .eq. 2) then
 297                     l = waterptr_aer(isize,itype)
 298                     if (l .ge. 1) tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr
 299                 else
 300                     tmpa = max( 0.0_r8, adryqvol_box(isize,itype) )*fact_apvolumr
 301                 end if
 302                 vol_distrib(isize,itype,ll)   = tmpa
 303                 sumold(ll) = sumold(ll) + tmpa
 304             end do
 305
 306 !   calculate dry and wet diameters and wet density
 307             if (xxmass .le. 1.01*smallmassbb) then
 308                 dpdry(isize,itype) = dcen_sect(isize,itype)*fact_apdiam
 309                 dpwet(isize,itype) = dpdry(isize,itype)
 310                 denswet(isize,itype) = xxdens
 311             else
 312                 dpdry(isize,itype) = (xxvolu/(xxnumb*piover6))**third
 313                 dpdry(isize,itype) = max( dpdry(isize,itype), dlo_sect(isize,itype)*fact_apdiam )
 314                 l = waterptr_aer(isize,itype)
 315                 if (l .ge. 1) then
 316                     tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr
 317                     xxmasswet = xxmass + tmpa
 318                     xxvoluwet = xxvolu + tmpa/(dens_water_aer*fact_apdens)
 319                     dpwet(isize,itype) = (xxvoluwet/(xxnumb*piover6))**third
 320                     dpwet(isize,itype) = min( dpwet(isize,itype), 30.0_r8*dhi_sect(isize,itype)*fact_apdiam )
 321                     denswet(isize,itype) = (xxmasswet/xxvoluwet)
 322                 else
 323                     dpwet(isize,itype) = dpdry(isize,itype)
 324                     denswet(isize,itype) = xxdens
 325                 end if
 326             end if
 327             densdry(isize,itype) = xxdens
 328
 329             vpdry(isize,itype) = piover6 * (dpdry(isize,itype)**3)
 330
 331         end do isize_loop_aa
 332
 333         end do itype_loop_aa
 334
 335
 336 !
 337 !   make call to coagulation routine
 338 !
 339         if (idiag_coag_onoff > 0) call coag_3d_conserve_check(   &
 340             nsize, maxd_asize, ntype_aer, maxd_atype, ncomp_coag, ncomp_coag_maxd, 1, lunout_coag,   &
 341             dpdry, num_distrib, vol_distrib, sum_num, sum_vol )
 342
 343         write(*,'(/a/)') 'mosaic_coag_3d_1box calling movingcenter_coag_3d'
 344         call movingcenter_coag_3d(   &
 345                 idiagbb, lunout_coag, nsize, maxd_asize, ntype_aer, maxd_atype, iphase, &
 346                 ncomp_coag, ncomp_coag_maxd,   &
 347                 temp_box, press_cgs, dtcoag, nsubstep_coag,   &
 348                 vpcut, vpdry, dpdry, dpwet, densdry, denswet, num_distrib, vol_distrib )
 349         write(*,'(/a/)') 'mosaic_coag_3d_1box backfrm movingcenter_coag_3d'
 350
 351         if (idiag_coag_onoff > 0) call coag_3d_conserve_check(   &
 352             nsize, maxd_asize, ntype_aer, maxd_atype, ncomp_coag, ncomp_coag_maxd, 2, lunout_coag,   &
 353             dpdry, num_distrib, vol_distrib, sum_num, sum_vol )
 354
 355
 356 !
 357 !   unload number/mass/volume from working arrays
 358 !
 359         sumnew(:) = 0.0
 360
 361 itype_loop_yy: &
 362         do itype = 1, ntype_aer
 363
 364 isize_loop_xx: &
 365         do isize = 1, nsize
 366             do ll = 1, ncomp_coag
 367                 sumnew(ll) = sumnew(ll) + max( 0.0_r8, vol_distrib(isize,itype,ll) )
 368             end do
 369
 370             l = numptr_aer(isize,itype,iphase)
 371             rbox(l) = max( 0.0_r8, num_distrib(isize,itype)/(rhoair_box*fact_apnumbmr) )
 372
 373 !   unload mass mixing ratios into rbox; also calculate dry mass and volume
 374             xxmass = 0.0
 375             xxvolu = 0.0
 376             do ll = 1, ncomp_aer(itype)
 377                 l = massptr_aer(ll,isize,itype,iphase)
 378                 if (l .ge. 1) then
 379                     tmpa = max( 0.0_r8, vol_distrib(isize,itype,ll) )
 380                     rbox(l) = tmpa/fact_apmassmr
 381                     xxmass = xxmass + tmpa
 382                     xxvolu = xxvolu + tmpa/(dens_aer(ll,itype)*fact_apdens)
 383                 end if
 384             end do
 385             adryqmas_box(isize,itype) = xxmass/fact_apmassmr
 386             xxvolubb(isize,itype) = xxvolu
 387
 388             ll = (ncomp_coag-3) + 1
 389             l = hyswptr_aer(isize,itype)
 390             if (l .ge. 1) rbox(l) = max( 0.0_r8, vol_distrib(isize,itype,ll) )/fact_apmassmr
 391
 392             ll = (ncomp_coag-3) + 2
 393             l = waterptr_aer(isize,itype)
 394             if (l .ge. 1) rbox(l) = max( 0.0_r8, vol_distrib(isize,itype,ll) )/fact_apmassmr
 395
 396             ll = (ncomp_coag-3) + 3
 397             adryqvol_box(isize,itype) = max( 0.0_r8, vol_distrib(isize,itype,ll) )/fact_apvolumr
 398         end do isize_loop_xx
 399
 400
 401 !
 402 !   calculate new dry density,
 403 !   and check for new mean dry-size within bounds
 404 !
 405 isize_loop_yy: &
 406         do isize = 1, nsize
 407
 408             xxmass = adryqmas_box(isize,itype)*fact_apmassmr
 409             xxvolu = adryqvol_box(isize,itype)*fact_apvolumr
 410             l = numptr_aer(isize,itype,iphase)
 411             xxnumb = rbox(l)*fact_apnumbmr
 412             iconform_numb = 1
 413
 414 !   do a cautious calculation of density, using volume from coag routine
 415             if (xxmass .le. smallmassbb) then
 416                 ncountaa(8) = ncountaa(8) + 1
 417                 xxdens = densdefault
 418                 xxvolu = xxmass/xxdens
 419                 xxnumb = xxmass/(volumcen_sect(isize,itype)*fact_apdia3*xxdens)
 420                 xxdpav = dcen_sect(isize,itype)*fact_apdiam
 421                 wwdens = xxdens
 422                 wwdpav = xxdpav
 423                 l = hyswptr_aer(isize,itype)
 424                 if (l .ge. 1) rbox(l) = 0.0
 425                 l = waterptr_aer(isize,itype)
 426                 if (l .ge. 1) rbox(l) = 0.0
 427                 iconform_numb = 0
 428                 if (idiagbb .ge. 210)   &
 429                     write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-7a',   &
 430                     itstep, itype, isize, xxmass, xxvolu, xxdens
 431             else if (xxmass .gt. 1000.0*xxvolu) then
 432 !   in this case, density is too large.  setting density=1000 forces
 433 !   next IF block while avoiding potential divide by zero or overflow
 434                 xxdens = 1000.0
 435             else
 436                 xxdens = xxmass/xxvolu
 437             end if
 438
 439             if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then
 440 !   (exception) case -- new dry density is unrealistic,
 441 !   so use dry volume from rbox instead of that from coag routine
 442                 ncountaa(7) = ncountaa(7) + 1
 443                 if (idiagbb .ge. 200)   &
 444                     write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-7b',   &
 445                     itstep, itype, isize, xxmass, xxvolu, xxdens
 446                 xxvolu = xxvolubb(isize,itype)
 447                 xxdens = xxmass/xxvolu
 448                 if (idiagbb .ge. 200)   &
 449                     write(lunout_coag,'(a,18x,1p,4e10.2)') 'coagaa-7c',   &
 450                     xxmass, xxvolu, xxdens
 451             end if
 452
 453 !   check for mean size ok, and conform number if not
 454             if (iconform_numb .gt. 0) then
 455                 if (xxnumb .gt. xxvolu/(volumlo_sect(isize,itype)*fact_apdia3)) then
 456                     ncountaa(9) = ncountaa(9) + 1
 457                     tmpa = xxnumb
 458                     xxnumb = xxvolu/(volumlo_sect(isize,itype)*fact_apdia3)
 459                     xxdpav = dlo_sect(isize,itype)*fact_apdiam
 460                     if (idiagbb .ge. 200)   &
 461                         write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-9 ',   &
 462                         itstep, itype, isize, xxvolu, tmpa, xxnumb
 463                 else if (xxnumb .lt. xxvolu/(volumhi_sect(isize,itype)*fact_apdia3)) then
 464                     ncountaa(10) = ncountaa(10) + 1
 465                     tmpa = xxnumb
 466                     xxnumb = xxvolu/(volumhi_sect(isize,itype)*fact_apdia3)
 467                     xxdpav = dhi_sect(isize,itype)*fact_apdiam
 468                     if (idiagbb .ge. 200)   &
 469                         write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-10',   &
 470                         itstep, itype, isize, xxvolu, tmpa, xxnumb
 471                 else
 472                     xxdpav = (xxvolu/(xxnumb*piover6))**third
 473                 end if
 474
 475                 tmpb = 0.0
 476                 l = waterptr_aer(isize,itype)
 477                 if (l .ge. 1) tmpb = max(0.0_r8,rbox(l))*fact_apmassmr
 478                 wwmass = xxmass + tmpb
 479                 wwvolu = xxvolu + tmpb/(dens_water_aer*fact_apdens)
 480                 wwdens = wwmass/wwvolu
 481                 wwdpav = xxdpav*((wwvolu/xxvolu)**third)
 482             end if
 483
 484 !   load number, mass, volume, dry-density back into arrays
 485             l = numptr_aer(isize,itype,iphase)
 486             rbox(l) = xxnumb/fact_apnumbmr
 487             adrydens_box(isize,itype) = xxdens/fact_apdens
 488             adrydpav_box(isize,itype) = xxdpav/fact_apdiam
 489             adryqmas_box(isize,itype) = xxmass/fact_apmassmr
 490             adryqvol_box(isize,itype) = xxvolu/fact_apvolumr
 491             awetdens_box(isize,itype) = wwdens/fact_apdens
 492             awetdpav_box(isize,itype) = wwdpav/fact_apdiam
 493
 494             dpdry(isize,itype) = xxdpav
 495
 496         end do isize_loop_yy
 497
 498         end do itype_loop_yy
 499
 500
 501         if (idiag_coag_onoff > 0) call coag_3d_conserve_check(   &
 502             nsize, maxd_asize, ntype_aer, maxd_atype, ncomp_coag, ncomp_coag_maxd, 3, lunout_coag,   &
 503             dpdry, num_distrib, vol_distrib, sum_num, sum_vol )
 504
 505
 506 !   temporary diagnostics
 507         do ll = 1, ncomp_coag
 508             ! ncountbb > 0 when mass/volume conservation relative-error > 1e-6
 509             tmpa = max( sumold(ll), sumnew(ll), 1.0e-35_r8 )
 510             if (abs(sumold(ll)-sumnew(ll)) .gt. 1.0e-6*tmpa) then
 511                 ncountbb(ll) = ncountbb(ll) + 1
 512                 ncountbb(0) = ncountbb(0) + 1
 513             end if
 514         end do
 515
 516         if (idiagbb .ge. 100) then
 517             write(msg,93030) 'coagbb ncntaa ', ncountaa(1:10)
 518             call peg_message( lunerr, msg )
 519             if (ncountbb(0) .gt. 0) then
 520                 do llx = 1, (ncomp_coag+9)/10
 521                     llb = llx*10
 522                     lla = llb - 9
 523                     llb = min( llb, ncomp_coag)
 524                     write(msg,93032) 'coagbb ncntbb',   &
 525                         mod(llx,10), ncountbb(lla:llb)
 526                     call peg_message( lunerr, msg )
 527                 end do
 528             end if
 529         end if
 530 93020   format( a, 1p, 10e10.2 )
 531 93030   format( a, 1p, 10i10 )
 532 93032   format( a, 1p, i1, 10i10 )
 533
 534
 535         return
 536         end subroutine mosaic_coag_3d_1box
 537
 538
 539 !----------------------------------------------------------------------
 540       subroutine coag_3d_conserve_check(   &
 541         nbin, nbin_maxd, ntype, ntype_maxd, ncomp, ncomp_maxd, iflagaa, lunout,   &
 542         dpdry, num_distrib, vol_distrib, sum_num, sum_vol )
 543
 544       implicit none
 545 !
 546 ! checks on conservation of volume before/after coag routine
 547 !
 548       integer, intent(in) :: nbin            ! actual number of size bins
 549       integer, intent(in) :: nbin_maxd       ! size-bin dimension for input (argument) arrays
 550       integer, intent(in) :: ntype           ! actual number of aerosol types
 551       integer, intent(in) :: ntype_maxd      ! aerosol type dimension for input (argument) arrays
 552       integer, intent(in) :: ncomp           ! actual number of aerosol volume components
 553       integer, intent(in) :: ncomp_maxd      ! volume-component dimension for input (argument) arrays
 554       integer, intent(in) :: lunout          ! logical unit for warning & diagnostic output
 555       integer, intent(in) :: iflagaa         !
 556
 557       real(r8), intent(in) :: dpdry(nbin_maxd,ntype_maxd)
 558       real(r8), intent(in) :: num_distrib(nbin_maxd,ntype_maxd)
 559       real(r8), intent(in) :: vol_distrib(nbin_maxd,ntype_maxd,ncomp_maxd)
 560       real(r8), intent(inout) :: sum_num(2)
 561       real(r8), intent(inout) :: sum_vol(0:ncomp_maxd,2)
 562
 563 ! local variables
 564       integer :: i, itype, l
 565       real(r8), parameter :: pi = 3.14159265358979323846_r8
 566       real(r8) :: voldry(nbin,ntype)
 567
 568
 569       if (iflagaa == 1) then
 570          i = 1
 571       else if (iflagaa == 2) then
 572          i = 2
 573       else if (iflagaa == 3) then
 574          i = 2
 575       else
 576          return
 577       end if
 578
 579
 580 ! sum initial or number and volumes
 581       voldry(1:nbin,1:ntype) = (pi/6.0)*(dpdry(1:nbin,1:ntype)**3)
 582       sum_num(i) = sum( num_distrib(1:nbin,1:ntype) )
 583       sum_vol(0,i) = sum( num_distrib(1:nbin,1:ntype)*voldry(1:nbin,1:ntype) )
 584       do l = 1, ncomp
 585          sum_vol(l,i) = sum( vol_distrib(1:nbin,1:ntype,l) )
 586       end do
 587
 588 ! write diagnostics
 589       if ((iflagaa /= 2) .and. (iflagaa /= 3)) return
 590       write(lunout,'(/a,i3)') 'coag_3d_conserve_check -- iflagaa = ', iflagaa
 591
 592       write(lunout,'(2a,1p,2e14.6)')   &
 593                'number total ', '    initial, final   =',   &
 594                sum_num(1:2)
 595       write(lunout,'(2a,1p,2e14.6)')   &
 596                'number total ', '   (final/initial)   =',   &
 597                sum_num(2)/sum_num(1)
 598
 599 ! with jacobson coag algorithm, dry-size of each bin is unchanged
 600 !    by the coagulation
 601 ! with moving-center coag algorithm, dry-size does change,
 602 !    so sum_vol(0,i) is currently not correct for iflagaa=2, because
 603 !    the incoming voldry() is currently the initial (i=1) value
 604 ! so do not output this when iflagaa=2
 605       if (iflagaa == 3) then
 606       write(lunout,'(2a,1p,2e14.6)')   &
 607                'volume total ', '    initial, final   =',   &
 608                sum_vol(0,1:2)
 609       write(lunout,'(2a,1p,2e14.6)')   &
 610                'volume total ', '   (initial/final)-1 =',   &
 611                (sum_vol(0,1)/sum_vol(0,2))-1.0
 612       end if
 613
 614       do l = 1, ncomp
 615          if (abs(sum_vol(l,2)) > 0.0) then
 616             write(lunout,'(a,i4,a,1p,2e14.6)')   &
 617                'volume l=',  l, '   (initial/final)-1 =',   &
 618                (sum_vol(l,1)/sum_vol(l,2))-1.0
 619          else
 620             write(lunout,'(a,i4,a,1p,2e14.6)')   &
 621                'volume l=',  l, '    initial, final   =',   &
 622                sum_vol(l,1), sum_vol(l,2)
 623          end if
 624       end do
 625
 626       return
 627       end subroutine coag_3d_conserve_check
 628
 629
 630
 631 !----------------------------------------------------------------------
 632       subroutine movingcenter_coag_3d(   &
 633         idiagbb, lunout, nbin, nbin_maxd, ntype, ntype_maxd, iphase, &
 634         ncomp, ncomp_maxd,   &
 635         tempk, press, deltat, nsubstep_in,   &
 636         vpcut, vpdry_in, dpdry, dpwet, densdry, denswet,   &
 637         num_distrib, vol_distrib )
 638
 639       use module_data_mosaic_aero, only:  &
 640          method_bcfrac, method_kappa, msectional_flag2
 641       use module_data_mosaic_asecthp, only:  &
 642          dens_aer, hygro_aer, &
 643          itype_md1_of_itype, itype_md2_of_itype, itype_of_itype_md1md2, &
 644          lptr_bc_aer, massptr_aer, &
 645          ncomp_aer, ntype_aer, ntype_md1_aer, ntype_md2_aer, &
 646          smallmassbb, xcut_atype_md1, xcut_atype_md2
 647
 648       implicit none
 649 !
 650 ! calculates aerosol coagulation for sectional representation
 651 !    and a single "particle type" using the single-type algorithm of
 652 !    jacobson (2002)
 653 !
 654 ! reference
 655 !    jacobson, m. z., 2002, analysis of aerosol interactions with numerical
 656 !       techniques for solving coagulation, nucleation, condensation,
 657 !       dissolution, and reversible chemistry among multiple size
 658 !       distributions, j. geophys. res., 107, d19, 4366
 659 !
 660
 661 !   IN arguments
 662       integer, intent(in) :: idiagbb
 663       integer, intent(in) :: lunout
 664       integer, intent(in) :: nbin            ! actual number of size bins
 665       integer, intent(in) :: nbin_maxd       ! size-bin dimension for input (argument) arrays
 666       integer, intent(in) :: ntype           ! actual number of aerosol types
 667       integer, intent(in) :: ntype_maxd      ! aerosol type dimension for input (argument) arrays
 668       integer, intent(in) :: iphase
 669       integer, intent(in) :: ncomp           ! actual number of aerosol volume components
 670       integer, intent(in) :: ncomp_maxd      ! volume-component dimension for input (argument) arrays
 671       integer, intent(in) :: nsubstep_in     ! number of time sub-steps for the integration
 672
 673       real(r8), intent(in) :: tempk               ! air temperature (K)
 674       real(r8), intent(in) :: press               ! air pressure (dyne/cm2)
 675       real(r8), intent(in) :: deltat              ! integration time (s)
 676       real(r8), intent(in) :: dpdry(nbin_maxd,ntype_maxd)    ! bin current volume-mean dry diameter (cm)
 677       real(r8), intent(in) :: dpwet(nbin_maxd,ntype_maxd)    ! bin wet (ambient) diameter (cm)
 678       real(r8), intent(in) :: densdry(nbin_maxd,ntype_maxd)  ! bin dry density (g/cm3)
 679       real(r8), intent(in) :: denswet(nbin_maxd,ntype_maxd)  ! bin wet (ambient) density (g/cm3)
 680       real(r8), intent(in) :: vpdry_in(nbin_maxd,ntype_maxd) ! bin current mean 1-particle dry volume (cm^3)
 681       real(r8), intent(in) :: vpcut(0:nbin_maxd)  ! 1-particle dry volumes at bin boundaries (cm^3)
 682
 683 !   INOUT arguments
 684       real(r8), intent(inout) :: num_distrib(nbin_maxd,ntype_maxd)
 685 !          num_distrib(i)   = number concentration for bin i (#/cm3)
 686       real(r8), intent(inout) :: vol_distrib(nbin_maxd,ntype_maxd,ncomp_maxd)
 687 !          vol_distrib(i,l) = volume concentration for bin i, component l (cm3/cm3)
 688
 689 !
 690 ! NOTE -- The sectional representation is based on dry particle size.
 691 ! Dry sizes are used to calculate the receiving bin indices and product fractions
 692 !   for each (bin-i1, bin-i2) coagulation combination.
 693 ! Wet sizes and densities are used to calculate the brownian coagulation kernel.
 694 !
 695 ! NOTE -- Units for num_distrib and vol_distrib
 696 !    num_distrib units MUST BE (#/cm3)
 697 !    vol_distrib
 698 !        (cm3/cm3) units are most consistent with the num_distrib units
 699 !        However, the algorithm works with almost any units!  So vol_distrib can
 700 !        be either masses or volumes, and either mixing ratios or concentrations
 701 !        (e.g., g/cm3-air, ug/kg-air, um3/cm3-air, cm3/kg-air, ...).
 702 !        Use whatever is convenient.
 703 !
 704
 705 !
 706 !   local variables
 707 !
 708       character(len=256) :: errmsg
 709       integer :: i, isubstep, itype, itype1, itype2, iwrite193
 710       integer :: j, jtype, jtype1, jtype2
 711       integer :: k, ktype, ktype1, ktype2, ktmp1, ktmp2
 712       integer :: l, method_ktype, nsubstep
 713
 714       real(r8), parameter :: frac_loss_limit = 0.8_r8
 715       real(r8), parameter :: pi = 3.14159265358979323846_r8
 716
 717       real(r8) :: del_num, del_voli, del_volj, deltat_sub
 718       real(r8) :: tmp_bcfrac, tmp_beta
 719       real(r8) :: tmp_densbc
 720       real(r8) :: tmp_kappa, tmp_kappaxvolu
 721       real(r8) :: tmp_massbc, tmp_massdry
 722       real(r8) :: tmp_volu
 723       real(r8) :: tmpa, tmpb, tmpc, tmpi, tmpj
 724       real(r8) :: vpdry_ipj
 725
 726       real(r8) ::   &
 727              betaxh(nbin,nbin),   &
 728              bcfrac_bin(nbin,ntype),   &
 729              cnum(nbin,ntype), cnum_old(nbin,ntype),   &
 730              cvol(nbin,ntype,ncomp), cvol_old(nbin,ntype,ncomp),   &
 731              kappa_bin(nbin,ntype),   &
 732              tmpvecb(nbin),   &
 733              vpdry(nbin,ntype)
 734
 735
 736 ! set method for determining ktype
 737       method_ktype = 0
 738       if ((ntype > 1) .and. (msectional_flag2 > 0)) then
 739          if (ncomp /= ncomp_aer(1) + 3) then
 740             write(errmsg,'(/2a,2(1x,i5)/)') '*** movingcenter_coag_3d - ', &
 741                'mismatch of ncomp & ncomp_aer(1)', ncomp, ncomp_aer(1)
 742             call wrf_error_fatal(trim(adjustl(errmsg)))
 743          end if
 744          if ((method_bcfrac /= 1) .and. (method_bcfrac /= 11)) then
 745             write(errmsg,'(/2a,2(1x,i5)/)') '*** movingcenter_coag_3d - ', &
 746                'method_bcfrac must be 1 or 11 but is', method_bcfrac
 747             call wrf_error_fatal(trim(adjustl(errmsg)))
 748          end if
 749          if ((method_kappa /= 11) .and. (method_kappa /= 12)) then
 750             write(errmsg,'(/2a,2(1x,i5)/)') '*** movingcenter_coag_3d - ', &
 751                'method_kappa must be 11 or 12 but is', method_kappa
 752             call wrf_error_fatal(trim(adjustl(errmsg)))
 753          end if
 754          method_ktype = 1
 755       end if
 756 !     method_ktype = 0
 757       iwrite193 = 0
 758
 759 !     tmp_densbc = dens_aer(ibc_a,1) ! this only works in mosaic box model
 760       tmp_densbc = -2.0
 761       do itype = 1, ntype
 762       do l = 1, ncomp_aer(itype)
 763          if (massptr_aer(l,1,itype,iphase) == lptr_bc_aer(1,itype,iphase)) then
 764             tmp_densbc = dens_aer(l,itype)
 765          end if
 766       end do
 767       end do
 768       if (tmp_densbc < -1.0) then
 769          call wrf_error_fatal('*** movingcenter_coag_3d - cannot find bc')
 770       end if
 771
 772 !
 773 ! check that fractional number loss from any bin over deltat_sub does
 774 ! not exceed frac_loss_limit
 775 !
 776       nsubstep = nsubstep_in
 777       deltat_sub = deltat/nsubstep   ! time substep (s)
 778       tmpa = 0.0
 779       do jtype = 1, ntype
 780          tmpvecb(:) = 0.0
 781
 782          do itype = 1, ntype
 783 ! calculate brownian coagulation kernel for itype,jtype
 784          call brownian_kernel_3dsub(   &
 785             nbin, nbin_maxd, ntype, ntype_maxd, itype, jtype, &
 786             tempk, press, dpwet, denswet, betaxh )
 787
 788          do j = 1, nbin
 789          do i = 1, nbin
 790             if ((i /= j) .or. (itype /= jtype)) then
 791                tmpvecb(j) = tmpvecb(j) + betaxh(i,j)*num_distrib(i,itype)
 792             else
 793                tmpvecb(j) = tmpvecb(j) + betaxh(i,j)*num_distrib(i,itype)*0.5
 794             end if
 795          end do   ! i
 796          end do   ! j
 797          end do   ! itype
 798
 799          do j = 1, nbin
 800             tmpa = max( tmpa, tmpvecb(j) )
 801          end do   ! j
 802       end do   ! jtype
 803
 804       if (idiagbb > 0) write(lunout,'(/a,i10,1p,4e12.4)') 'coag aa nsub, dtsub, fracloss', &
 805             nsubstep, deltat_sub, tmpa*deltat_sub, frac_loss_limit
 806       if (tmpa*deltat_sub > frac_loss_limit) then
 807          tmpb = tmpa*deltat/frac_loss_limit
 808          isubstep = int(tmpb) + 1
 809          tmpb = deltat/isubstep
 810          if (idiagbb > 0) write(lunout,'( a,i10,1p,4e12.4)') 'coag bb nsub, dtsub, fracloss', &
 811             isubstep, tmpb, tmpa*tmpb
 812          nsubstep = isubstep
 813          deltat_sub = deltat/nsubstep
 814       end if
 815
 816
 817 ! initialize working arrays
 818 !   cnum(i)   = current number conc (#/cm3)   for bin i
 819 !   cvol(i,l) = current volume conc (cm3/cm3) for bin i, species l
 820       cnum(1:nbin,1:ntype) = num_distrib(1:nbin,1:ntype)
 821       cvol(1:nbin,1:ntype,1:ncomp) = vol_distrib(1:nbin,1:ntype,1:ncomp)
 822
 823
 824 !
 825 ! solve coagulation over multiple time substeps
 826 !
 827 solve_isubstep_loop:   &
 828       do isubstep = 1, nsubstep
 829 !        write(91,*) 'maa =', 2
 830
 831
 832 ! copy current num, vol values to "old" arrays
 833          cnum_old(1:nbin,1:ntype) = cnum(1:nbin,1:ntype)
 834          cvol_old(1:nbin,1:ntype,1:ncomp) = cvol(1:nbin,1:ntype,1:ncomp)
 835
 836 ! calc 1-particle dry volume
 837          if (isubstep == 1) then
 838              vpdry(1:nbin,1:ntype) = vpdry_in(1:nbin,1:ntype)
 839 !        else
 840 ! ***
 841 ! *** should be updating vpdry(:) after first substep
 842 ! ***
 843 ! *** also densdry(:,:), but only if it is needed in calc of ktype
 844 ! ***    (which is when method_bcfrac=1)
 845 ! ***
 846          end if
 847
 848          if (method_ktype == 1) then
 849 ! calc information needed for getting ktype
 850             do itype = 1, ntype
 851             do i = 1, nbin
 852                tmp_massbc  = 0.0
 853                tmp_massdry = 0.0
 854                tmp_kappaxvolu  = 0.0
 855                tmp_volu = 0.0
 856                do l = 1, ncomp_aer(itype)
 857                   tmpa = vol_distrib(i,itype,l)
 858                   tmpc = tmpa/dens_aer(l,itype)
 859                   if (method_bcfrac == 11) then
 860                      ! bcfrac is dry-volume fraction
 861                      tmpb = tmpc
 862                   else
 863                      ! bcfrac is dry-mass fraction
 864                      tmpb = tmpa
 865                   end if
 866                   tmp_massdry = tmp_massdry + tmpb
 867 !                 if (l == ibc_a) then  ! this only works in mosaic box model, when the ordering
 868 !                                       ! of species in aer(...) and massptr_aer(...) are the same
 869                   if (massptr_aer(l,i,itype,iphase) == lptr_bc_aer(i,itype,iphase)) then
 870                      tmp_massbc = tmpb
 871                      if (method_kappa == 12) cycle
 872                      ! in this case, kappa includes non-bc species only
 873                   end if
 874                   tmp_volu = tmp_volu + tmpc
 875                   tmp_kappaxvolu = tmp_kappaxvolu + tmpc*hygro_aer(l,itype)
 876                end do   ! l
 877
 878                if (tmp_massdry <= smallmassbb) then
 879                   itype1 = itype_md1_of_itype( itype )
 880                   bcfrac_bin(i,itype) = 0.5*( xcut_atype_md1(itype1-1) + xcut_atype_md1(itype1) )
 881                else
 882                   bcfrac_bin(i,itype) = tmp_massbc/tmp_massdry
 883                end if
 884                bcfrac_bin(i,itype) = max( 0.0_r8, xcut_atype_md1(0), &
 885                                           bcfrac_bin(i,itype) )
 886                bcfrac_bin(i,itype) = min( 1.0_r8, xcut_atype_md1(ntype_md1_aer), &
 887                                           bcfrac_bin(i,itype) )
 888
 889                if (tmp_volu <= smallmassbb) then
 890                   itype2 = itype_md2_of_itype( itype )
 891                   kappa_bin(i,itype)  = sqrt( max( 0.0_r8, &
 892                                               xcut_atype_md2(itype2-1) * xcut_atype_md2(itype2) ) )
 893                else
 894                   kappa_bin(i,itype) = tmp_kappaxvolu/tmp_volu
 895                end if
 896                kappa_bin(i,itype) = max( 1.0e-10_r8, xcut_atype_md2(0), &
 897                                           kappa_bin(i,itype) )
 898                kappa_bin(i,itype) = min( 2.0_r8,     xcut_atype_md2(ntype_md2_aer), &
 899                                           kappa_bin(i,itype) )
 900             end do   ! i
 901             end do   ! itype
 902          end if   ! (method_ktype == 1)
 903
 904
 905 ! main loops over i/jtype and i/j(size)
 906          do jtype = 1, ntype
 907          jtype1 = itype_md1_of_itype( jtype )
 908          jtype2 = itype_md2_of_itype( jtype )
 909
 910          do itype = 1, jtype   ! itype <= jtype avoids double counting
 911          itype1 = itype_md1_of_itype( itype )
 912          itype2 = itype_md2_of_itype( itype )
 913
 914 ! calculate brownian coagulation kernel again for itype,jtype
 915 ! problem is that it can be too big !
 916          call brownian_kernel_3dsub(   &
 917             nbin, nbin_maxd, ntype, ntype_maxd, itype, jtype, &
 918             tempk, press, dpwet, denswet, betaxh )
 919
 920          if (iwrite193 > 0) &
 921             write(193,'(/a,2(4x,2i4))') 'i/jtype1 then ...2', itype1, jtype1, itype2, jtype2
 922          do j = 1, nbin
 923          do i = 1, nbin
 924             if (itype == jtype) then
 925                if (i > j) cycle   ! i <= j avoids double counting when itype = jtype
 926             end if
 927
 928             ! determine the bin that receives the "i+j" particles
 929             vpdry_ipj = vpdry(i,itype) + vpdry(j,jtype)
 930             k = max( i, j )
 931             do while (k > -99)
 932                if (vpdry_ipj > vpcut(k)) then
 933                   k = k+1
 934                   if (k >= nbin) exit   ! do while ...
 935                else if (vpdry_ipj < vpcut(k-1)) then
 936                   k = k-1
 937                   if (k <= 1) exit   ! do while ...
 938                else
 939                   exit   ! do while ...
 940                end if
 941             end do   ! while ...
 942             k = max( 1, min( nbin, k ) )
 943
 944             if (method_ktype == 1) then
 945                if (method_bcfrac == 11) then
 946                   tmpi = vpdry(i,itype)
 947                   tmpj = vpdry(j,jtype)
 948                else
 949                   tmpi = vpdry(i,itype)*densdry(i,itype)
 950                   tmpj = vpdry(j,jtype)*densdry(j,jtype)
 951                end if
 952                tmpa = tmpi/max( (tmpi + tmpj), 1.0e-35_r8 )
 953                tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
 954                tmp_bcfrac = bcfrac_bin(i,itype)*tmpa + bcfrac_bin(j,jtype)*(1.0_r8 - tmpa)
 955
 956                ktype1 = ntype_md1_aer
 957                do ktmp1 = 1, ntype_md1_aer - 1
 958                    if (tmp_bcfrac <= xcut_atype_md1(ktmp1)) then
 959                        ktype1 = ktmp1
 960                        exit
 961                    end if
 962                end do
 963
 964                if (method_kappa == 11) then
 965                   tmpi = vpdry(i,itype)
 966                   tmpj = vpdry(j,jtype)
 967                else if (method_bcfrac == 11) then
 968                   tmpa = 1.0_r8 - bcfrac_bin(i,itype)
 969                   tmpi = vpdry(i,itype)*max(tmpa,1.0e-10_r8)
 970                   tmpa = 1.0_r8 - bcfrac_bin(j,jtype)
 971                   tmpj = vpdry(j,jtype)*max(tmpa,1.0e-10_r8)
 972                else
 973                   tmpa = 1.0_r8 - bcfrac_bin(i,itype)*(densdry(i,itype)/tmp_densbc)
 974                   tmpi = vpdry(i,itype)*max(tmpa,1.0e-10_r8)
 975                   tmpa = 1.0_r8 - bcfrac_bin(j,jtype)*(densdry(j,jtype)/tmp_densbc)
 976                   tmpj = vpdry(j,jtype)*max(tmpa,1.0e-10_r8)
 977                end if
 978                tmpa = tmpi/max( (tmpi + tmpj), 1.0e-35_r8 )
 979                tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) )
 980                tmp_kappa = kappa_bin(i,itype)*tmpa + kappa_bin(j,jtype)*(1.0_r8 - tmpa)
 981
 982                ktype2 = ntype_md2_aer
 983                do ktmp2 = 1, ntype_md2_aer - 1
 984                    if (tmp_kappa <= xcut_atype_md2(ktmp2)) then
 985                        ktype2 = ktmp2
 986                        exit
 987                    end if
 988                end do
 989
 990                ktype = itype_of_itype_md1md2(ktype1,ktype2)
 991
 992                if (iwrite193 > 0) then
 993
 994 !              if ( (num_distrib(i,itype) >= 1.0e30) .and. &
 995 !                   (num_distrib(j,jtype) >= 1.0e-1) ) then
 996 !                 write(193,'(a,3(2x,3i3),2f9.4)') 'i(s,t1,t2); j; k; wbc, kap', &
 997 !                    i, itype_md1_of_itype(itype), itype_md2_of_itype(itype), &
 998 !                    j, itype_md1_of_itype(jtype), itype_md2_of_itype(jtype), &
 999 !                    k, ktype1, ktype2, tmp_bcfrac, tmp_kappa
1000 !              end if
1001
1002 !              if ( (j==6 .or. j==12 .or. j==18 .or. j==24) .and. &
1003 !                   (i==j .or. i==j-1 .or. i==j-2 .or. i==j-3 .or. i==j-4 .or. i==j-5 .or. i==j-10) .and. &
1004 !                   ( ( (jtype1==1 .or. jtype1==5 .or. jtype1==10 .or. jtype1==15) .and. &
1005 !                       (jtype1-itype1 <= 1) ) .or. &
1006 !                     ( (jtype1==2 .or. jtype1==6 .or. jtype1==11 .or. jtype1==16) .and. &
1007 !                       (jtype1-itype1 == 1) ) ) .and. &
1008 !                   ( ( (jtype2==1 .or. jtype2==3 .or. jtype2==6 .or. jtype2== 9) .and. &
1009 !                       (jtype2-itype2 <= 1) ) .or. &
1010 !                     ( (jtype2==2 .or. jtype2==4 .or. jtype2==7 .or. jtype2==10) .and. &
1011 !                       (jtype1-itype1 == 1) ) ) ) then
1012
1013                if ( ( j==6 .or. j==12 .or. j==18 .or. j==24 ) .and. &
1014                     ( abs(i-j)<=4 .or. abs(i-j)==8 ) .and. &
1015                     ( jtype1==1 .or. jtype1==5 .or. jtype1==10 .or. jtype1==15 .or. &
1016                       itype1==1 .or. itype1==5 .or. itype1==10 .or. itype1==15 ) .and. &
1017                     ( abs(jtype1-itype1)<=1 .or. abs(jtype1-itype1)==5 .or. abs(jtype1-itype1)==10 ) .and. &
1018                     ( jtype2==1 .or. jtype2==5 .or. jtype2==9 .or. &
1019                       itype2==1 .or. itype2==5 .or. itype2==9 ) .and. &
1020                     ( abs(jtype2-itype2)<=1 .or. abs(jtype2-itype2)==4 .or. abs(jtype2-itype2)==8 ) ) then
1021                   write(193,'(a,3i3,1p,5e9.2,0p,2f6.3,3i3,3f6.3,3i3,3f7.4)') 'ijk&types', &
1022                      i, j, k, dpdry(i,itype), dpdry(j,jtype), vpdry(i,itype), vpdry(j,jtype), vpdry_ipj, &
1023                      densdry(i,itype), densdry(j,jtype), &
1024                      itype1, jtype1, ktype1, bcfrac_bin(i,itype), bcfrac_bin(j,jtype), tmp_bcfrac, &
1025                      itype2, jtype2, ktype2, kappa_bin(i,itype), kappa_bin(j,jtype), tmp_kappa
1026                end if
1027
1028                end if  ! (iwrite193 > 0)
1029
1030             else
1031                if (itype == jtype) then
1032                   ktype = itype
1033                else
1034                   ktype = ntype
1035                end if
1036             end if
1037 !           write(183,'(3i4,1p,5e10.2)') i, j, k, vpdry(i), vpdry(j), vpdry_ipj
1038 !           write(183,'(12x,1p,5e10.2)') vpcut(i), vpcut(j), vpcut(k-1:k)
1039
1040             tmp_beta = betaxh(i,j) * deltat_sub
1041
1042             ! now calc the coagulation changes to number and volume
1043             if ((i == j) .and. (itype == jtype)) then
1044                del_num = 0.5*tmp_beta*cnum_old(i,itype)*cnum_old(i,itype)
1045                if ((i == k) .and. (itype == ktype)) then
1046                   ! here i,itype + i,itype --> i,itype
1047                   ! each coag event consumes 2 and produces 1 "i" particle
1048                   cnum(i,itype) = cnum(i,itype) - del_num
1049                   cycle
1050                else
1051                   ! here i,itype + i,itype --> k,ktype AND i,k and/or itype,ktype differ
1052                   ! each coag event consumes 2 "i" and produces 1 "k" particle
1053                   ! and there is transfer of mass" and produces 1 "k" particle
1054                   cnum(i,itype) = cnum(i,itype) - del_num*2.0
1055                   cnum(k,ktype) = cnum(k,ktype) + del_num
1056                   do l = 1, ncomp
1057                      del_voli = tmp_beta*cnum_old(i,itype)*cvol_old(i,itype,l)
1058                      cvol(i,itype,l) = cvol(i,itype,l) - del_voli
1059                      cvol(k,ktype,l) = cvol(k,ktype,l) + del_voli
1060                   end do
1061                end if
1062
1063             else   ! here i < j OR itype < jtype
1064                del_num = tmp_beta*cnum_old(i,itype)*cnum_old(j,jtype)
1065                cnum(i,itype) = cnum(i,itype) - del_num
1066                if ((j == k) .and. (jtype == ktype)) then
1067                   ! here i,itype and j,jtype are different but j,jtype = k,ktype
1068                   ! i loses number and transfers mass to k=j
1069                   ! j=k does not lose any of its own number or mass
1070                   do l = 1, ncomp
1071                      del_voli = tmp_beta*cnum_old(j,jtype)*cvol_old(i,itype,l)
1072                      cvol(i,itype,l) = cvol(i,itype,l) - del_voli
1073                      cvol(k,ktype,l) = cvol(k,ktype,l) + del_voli
1074                   end do
1075                else
1076                   ! here i,itype; j,jtype; and k,ktype are all different
1077                   ! both i and j lose number and transfer mass to k
1078                   cnum(j,jtype) = cnum(j,jtype) - del_num
1079                   cnum(k,ktype) = cnum(k,ktype) + del_num
1080                   do l = 1, ncomp
1081                      del_voli = tmp_beta*cnum_old(j,jtype)*cvol_old(i,itype,l)
1082                      del_volj = tmp_beta*cnum_old(i,itype)*cvol_old(j,jtype,l)
1083                      cvol(i,itype,l) = cvol(i,itype,l) - del_voli
1084                      cvol(j,jtype,l) = cvol(j,jtype,l) - del_volj
1085                      cvol(k,ktype,l) = cvol(k,ktype,l) + del_voli + del_volj
1086                   end do
1087                end if
1088
1089             end if
1090
1091          end do   ! i
1092          end do   ! j
1093
1094          end do   ! itype
1095          end do   ! jtype
1096
1097       end do solve_isubstep_loop
1098
1099
1100 ! transfer new concentrations from working arrays
1101       num_distrib(1:nbin,1:ntype) = cnum(1:nbin,1:ntype)
1102       do l = 1, ncomp
1103          vol_distrib(1:nbin,1:ntype,l) = cvol(1:nbin,1:ntype,l)
1104       end do
1105
1106
1107 ! all done
1108       return
1109       end subroutine movingcenter_coag_3d
1110
1111
1112
1113 !-----------------------------------------------------------------------
1114       subroutine brownian_kernel_3dsub(   &
1115          nbin, nbin_maxd, ntype, ntype_maxd, itype, jtype, &
1116          tempk, press, dpwet, denswet, bckernel )
1117 !
1118 ! this routine calculates brownian coagulation kernel
1119 ! using on eqn 16.28 of
1120 !    jacobson,  m. z. (1999) fundamentals of atmospheric modeling.
1121 !       cambridge university press, new york, 656 pp.
1122 !
1123       implicit none
1124
1125 ! arguments
1126       integer, intent(in) :: nbin            ! actual number of size bins
1127       integer, intent(in) :: nbin_maxd       ! size-bin dimension for input (argument) arrays
1128       integer, intent(in) :: ntype           ! actual number of aerosol types
1129       integer, intent(in) :: ntype_maxd      ! aerosol type dimension for input (argument) arrays
1130       integer, intent(in) :: itype, jtype    ! the itype & jtype for which to calc the kernal
1131
1132       real(r8), intent(in)  :: tempk            ! air temperature (K)
1133       real(r8), intent(in)  :: press            ! air pressure (dyne/cm2)
1134       real(r8), intent(in)  :: dpwet(nbin_maxd,ntype_maxd)      ! bin wet (ambient) diameter (cm)
1135       real(r8), intent(in)  :: denswet(nbin_maxd,ntype_maxd)    ! bin wet (ambient) density (g/cm3)
1136       real(r8), intent(out) :: bckernel(nbin,nbin)  ! brownian coag kernel (cm3/#/s)
1137
1138 ! local variables
1139       integer i, j
1140
1141       real(r8), parameter :: pi = 3.14159265358979323846_r8
1142
1143       real(r8)   &
1144              apfreepath, avogad,   &
1145              boltz,   &
1146              cunning,   &
1147              deltasq_i, deltasq_j,   &
1148              diffus_i, diffus_j, diffus_ipj,   &
1149              dpwet_i, dpwet_j, dpwet_ipj,   &
1150              gasfreepathx2, gasspeed,   &
1151              knud,   &
1152              mass_i, mass_j, mwair,   &
1153              pi6,   &
1154              rgas, rhoair,   &
1155              speedsq_i, speedsq_j,   &
1156              tmp1,   &
1157              viscosd, viscosk
1158
1159 ! boltz   = boltzmann's constant (erg/K = g*cm2/s/K)
1160 ! avogad  = avogadro's number (molecules/mol)
1161 ! mwair   = molecular weight of air (g/mol)
1162 ! rgas    = gas constant ((dyne/cm2)/(mol/cm3)/K)
1163 ! rhoair  = air density (g/cm3)
1164 ! viscosd = air dynamic viscosity (g/cm/s)
1165 ! viscosk = air kinematic viscosity (cm2/s)
1166 ! gasspeed      = air molecule mean thermal velocity (cm/s)
1167 ! gasfreepathx2 = air molecule mean free path (cm) x 2
1168       pi6 = pi/6.0
1169
1170       boltz = 1.38065e-16
1171       avogad = 6.02214e+23
1172       mwair = 28.966
1173       rgas = 8.31447e7
1174
1175       rhoair = press*mwair/(rgas*tempk)
1176
1177       viscosd = 1.8325e-04*(416.16/(tempk+120.0)) * (tempk/296.16)**1.5
1178       viscosk = viscosd/rhoair
1179       gasspeed = sqrt(8.0*boltz*tempk*avogad/(pi*mwair))
1180       gasfreepathx2 = 4.0*viscosk/gasspeed
1181
1182 ! coagulation kernel from eqn 16.28 of jacobson (1999) text
1183 !
1184 ! diffus_i/j  = particle brownian diffusion coefficient  (cm2/s)
1185 ! speedsq_i/j = square of particle mean thermal velocity (cm/s)
1186 ! apfreepath  = particle mean free path (cm)
1187 ! cunning     = cunningham slip-flow correction factor
1188 ! deltasq_i/j = square of "delta_i" in eqn 16.29
1189
1190        bckernel(:,:) = -2.0_r8
1191
1192        do i = 1, nbin
1193
1194           dpwet_i = dpwet(i,itype)                    ! particle wet diameter (cm)
1195           mass_i = pi6*denswet(i,itype)*(dpwet_i**3)  ! particle wet mass (g)
1196           knud = gasfreepathx2/dpwet_i
1197           cunning = 1.0 + knud*(1.249 + 0.42*exp(-0.87/knud))
1198           diffus_i = boltz*tempk*cunning/(3.0*pi*dpwet_i*viscosd)
1199           speedsq_i = 8.0*boltz*tempk/(pi*mass_i)
1200           apfreepath = 8.0*diffus_i/(pi*sqrt(speedsq_i))
1201           tmp1 = (dpwet_i + apfreepath)**3   &
1202                - (dpwet_i*dpwet_i + apfreepath*apfreepath)**1.5
1203           deltasq_i = ( tmp1/(3.0*dpwet_i*apfreepath) - dpwet_i )**2
1204
1205           do j = 1, nbin
1206
1207            if ( (itype == jtype) .and.   &
1208                 (bckernel(j,i) > -1.0_r8) ) then
1209              ! value already calculated
1210              bckernel(i,j) = bckernel(j,i)
1211
1212            else
1213
1214              dpwet_j = dpwet(j,jtype)                    ! particle wet diameter (cm)
1215              mass_j = pi6*denswet(j,jtype)*(dpwet_j**3)  ! particle wet mass (g)
1216              knud = gasfreepathx2/dpwet_j
1217              cunning = 1.0 + knud*(1.249 + 0.42*exp(-0.87/knud))
1218              diffus_j = boltz*tempk*cunning/(3.0*pi*dpwet_j*viscosd)
1219              speedsq_j = 8.0*boltz*tempk/(pi*mass_j)
1220              apfreepath = 8.0*diffus_j/(pi*sqrt(speedsq_j))
1221              tmp1 = (dpwet_j + apfreepath)**3   &
1222                   - (dpwet_j*dpwet_j + apfreepath*apfreepath)**1.5
1223              deltasq_j = ( tmp1/(3.0*dpwet_j*apfreepath) - dpwet_j )**2
1224
1225              dpwet_ipj = dpwet_i + dpwet_j
1226              diffus_ipj = diffus_i + diffus_j
1227              tmp1 = (dpwet_ipj/(dpwet_ipj + 2.0*sqrt(deltasq_i + deltasq_j)))   &
1228                   + (8.0*diffus_ipj/(dpwet_ipj*sqrt(speedsq_i + speedsq_j)))
1229              bckernel(i,j) = max( 0.0_r8, 2.0*pi*dpwet_ipj*diffus_ipj/tmp1 )
1230
1231            end if
1232
1233           end do   ! j
1234        end do   ! i
1235
1236        return
1237        end subroutine brownian_kernel_3dsub
1238
1239
1240
1241 !-----------------------------------------------------------------------
1242
1243
1244
1245         end module module_mosaic_coag3d