chem/module_mosaic_coag1d.F

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