module module_mosaic_coag1d use module_data_mosaic_kind use module_peg_util implicit none contains !----------------------------------------------------------------------- subroutine mosaic_coag_1box( istat_coag, & idiagbb_in, itstep, & dtcoag_in, & temp_box, patm_box, rhoair_box, rbox, & fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam, & adrydens_box, awetdens_box, & adrydpav_box, awetdpav_box, adryqmas_box, mcoag_flag1 ) ! ! calculates aerosol particle coagulation for a single grid box ! over timestep dtcoag_in ! ! uses the single-type coagulation algorithm of jacobson (2002) ! use module_data_mosaic_main, only: ntot_used, pi, piover6, third use module_data_mosaic_aero, only: ifreq_coag use module_data_mosaic_asecthp, only: & ai_phase, hyswptr_aer, & lunerr, lunout, & maxd_acomp, maxd_asize, maxd_atype, & ncomp_aer, ncomp_plustracer_aer, nsize_aer, ntype_aer, & massptr_aer, numptr_aer, waterptr_aer, & dens_aer, dens_water_aer, & dcen_sect, dhi_sect, dlo_sect, & smallmassbb, & volumcen_sect, volumcut_sect, volumhi_sect, volumlo_sect ! subr arguments integer, intent(inout) :: istat_coag ! =0 if no problems integer, intent(in) :: idiagbb_in ! controls diagnostics integer, intent(in) :: itstep ! host code time step index integer, intent(in) :: mcoag_flag1 real(r8), intent(in) :: dtcoag_in ! time step (s) real(r8), intent(in) :: temp_box ! air temp (K) real(r8), intent(in) :: patm_box ! air pressure (atm) real(r8), intent(in) :: rhoair_box ! air density (g/cm^3) real(r8), intent(inout) :: rbox(ntot_used) ! rbox - holds aerosol number and component mass mixing ratios ! units for number mr are such that (rbox*fact_apnumbmr) gives (#/g-air) ! units for mass mr are such that (rbox*fact_apmassmr) gives (g/g-air) real(r8), intent(in) :: fact_apmassmr, fact_apnumbmr, fact_apdens, fact_apdiam ! these are units conversion factors ! the module_data_mosaic_asecthp densities are such that ! dens_aer*fact_apdens gives (g/cm^3) ! the module_data_mosaic_asecthp diameters are such that ! d[lo,cen,hi]_sect*fact_apdiam gives (cm) ! the module_data_mosaic_asecthp volumes are such that ! volum[lo,cen,hi]_sect**fact_apdiam**3) gives (cm^3) real(r8), intent(inout) :: adrydens_box(maxd_asize,maxd_atype) real(r8), intent(inout) :: awetdens_box(maxd_asize,maxd_atype) real(r8), intent(inout) :: adrydpav_box(maxd_asize,maxd_atype) real(r8), intent(inout) :: awetdpav_box(maxd_asize,maxd_atype) real(r8), intent(inout) :: adryqmas_box(maxd_asize,maxd_atype) ! adryqmas_box = aerosol total-dry-mass mixing ratio (g/mol-air) ! adrydens_box = aerosol dry density (units same as dens_aer) ! awetdens_box = aerosol wet density (units same as dens_aer) ! adrydpav_box = aerosol mean dry diameter (units same as dlo_sect) ! awetdpav_box = aerosol mean wet diameter (units same as dlo_sect) ! adryqmas_box = aerosol total-dry-mass mixing ratio (units same as rbox) ! local variables integer, parameter :: coag_method = 1 integer, parameter :: ncomp_coag_maxd = maxd_acomp + 3 integer :: l, ll, lla, llb, llx integer :: isize, itype, iphase integer :: iconform_numb integer :: idiagbb, idiag_coag_onoff integer :: lunout_coag integer :: ncomp_coag, nsize, nsubstep_coag, ntau_coag integer,save :: ncountaa(10), ncountbb(0:ncomp_coag_maxd) real(r8), parameter :: factsafe = 1.00001 ! units for local working varibles: ! size = cm ! density = g/cm^3 ! number mixing ratio = #/g-air ! EXCEPT num_distrib = #/cm^3-air ! mass mixing ratio = g/g-air ! volume mixing ratio = cm^3/g-air ! real(r8) :: densdefault, dtcoag real(r8) :: fact_apdia3, fact_apvolumr real(r8) :: press_cgs ! air pressure (dyne/cm2) real(r8) :: tmpa, tmpb real(r8) :: wwdens, wwdpav, wwmass, wwvolu real(r8) :: xxdens, xxdpav, xxmass, xxnumb, xxvolu real(r8) :: xxmasswet, xxvoluwet real(r8) :: adryqvol_box(maxd_asize,maxd_atype) ! aerosol total-dry-volume mixing ratio (cm3/mol-air) real(r8) :: dpdry(maxd_asize), dpwet(maxd_asize), denswet(maxd_asize) real(r8) :: num_distrib(maxd_asize) real(r8) :: sumnew(ncomp_coag_maxd), sumold(ncomp_coag_maxd) real(r8) :: sum_num(2), sum_vol(0:ncomp_coag_maxd,2) real(r8) :: vol_distrib(maxd_asize,ncomp_coag_maxd) real(r8) :: vpcut(0:maxd_asize), vpdry(maxd_asize) real(r8) :: xxvolubb(maxd_asize) character(len=120) :: msg if (mcoag_flag1 <= 0) return if (ifreq_coag > 1) then if (mod(itstep,ifreq_coag) /= 0) return dtcoag = dtcoag_in*ifreq_coag else dtcoag = dtcoag_in end if istat_coag = 0 lunout_coag = 6 if (lunout .gt. 0) lunout_coag = lunout ! set variables that do not change idiagbb = idiagbb_in idiag_coag_onoff = +1 if (idiagbb <= 0) idiag_coag_onoff = 0 iphase = ai_phase fact_apdia3 = fact_apdiam*fact_apdiam*fact_apdiam fact_apvolumr = fact_apmassmr/fact_apdens nsubstep_coag = 1 press_cgs = patm_box*1.01325e6_r8 ! temporary diagnostics ncountaa(:) = 0 ncountbb(:) = 0 do 2800 itype = 1, ntype_aer ! do 2800 itype = 1, 1 nsize = nsize_aer(itype) ncomp_coag = ncomp_plustracer_aer(itype) + 3 densdefault = dens_aer(1,itype)*fact_apdens vpcut(0:nsize) = volumcut_sect(0:nsize,itype)*fact_apdia3 ncountaa(1) = ncountaa(1) + 1 ! do initial calculation of dry mass, volume, and density, ! and initial number conforming (as needed) vol_distrib(:,:) = 0.0 sumold(:) = 0.0 isize_loop_aa: & do isize = 1, nsize l = numptr_aer(isize,itype,iphase) xxnumb = rbox(l)*fact_apnumbmr xxmass = adryqmas_box(isize,itype)*fact_apmassmr xxdens = adrydens_box(isize,itype)*fact_apdens iconform_numb = 1 if ( (xxdens .lt. 0.1) .or. (xxdens .gt. 20.0) ) then ! (exception) case of drydensity not valid ! so compute dry mass and volume from rbox ncountaa(2) = ncountaa(2) + 1 xxvolu = 0.0 if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-2a', & itstep, itype, isize, xxmass, xxvolu, xxdens xxmass = 0.0 do ll = 1, ncomp_aer(itype) l = massptr_aer(ll,isize,itype,iphase) if (l .ge. 1) then tmpa = max( 0.0_r8, rbox(l) ) xxmass = xxmass + tmpa xxvolu = xxvolu + tmpa/dens_aer(ll,itype) end if end do xxmass = xxmass*fact_apmassmr xxvolu = xxvolu*fact_apvolumr if (xxmass .gt. 0.99*smallmassbb) then xxdens = xxmass/xxvolu xxvolu = xxmass/xxdens if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-2b', & itstep, itype, isize, xxmass, xxvolu, xxdens end if else xxvolu = xxmass/xxdens end if if (xxmass .le. 1.01*smallmassbb) then ! when drymass extremely small, use default density and bin center size, ! and zero out water ncountaa(3) = ncountaa(3) + 1 xxdens = densdefault xxvolu = xxmass/xxdens xxnumb = xxmass/(volumcen_sect(isize,itype)*fact_apdia3*xxdens) l = hyswptr_aer(isize,itype) if (l .ge. 1) rbox(l) = 0.0 l = waterptr_aer(isize,itype) if (l .ge. 1) rbox(l) = 0.0 iconform_numb = 0 if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-2c', & itstep, itype, isize, xxmass, xxvolu, xxdens else xxvolu = xxmass/xxdens end if ! check for mean dry-size 'safely' within bounds, and conform number if not if (iconform_numb .gt. 0) then if (xxnumb .gt. & xxvolu/(factsafe*volumlo_sect(isize,itype)*fact_apdia3)) then ncountaa(4) = ncountaa(4) + 1 tmpa = xxnumb xxnumb = xxvolu/(factsafe*volumlo_sect(isize,itype)*fact_apdia3) if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-4 ', & itstep, itype, isize, xxvolu, tmpa, xxnumb else if (xxnumb .lt. & xxvolu*factsafe/(volumhi_sect(isize,itype)*fact_apdia3)) then ncountaa(5) = ncountaa(5) + 1 tmpa = xxnumb xxnumb = xxvolu*factsafe/(volumhi_sect(isize,itype)*fact_apdia3) if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-5 ', & itstep, itype, isize, xxvolu, tmpa, xxnumb if (idiagbb .ge. 200) & write(lunout_coag,'(a,10x, 8x,1p,4e12.4)') 'coagcc-5 ', & dlo_sect(isize,itype), dlo_sect(isize,itype)*fact_apdiam, & volumlo_sect(isize,itype), volumlo_sect(isize,itype)*fact_apdia3 if (idiagbb .ge. 200) & write(lunout_coag,'(a,10x, 8x,1p,4e12.4)') 'coagcc-5 ', & dhi_sect(isize,itype), dhi_sect(isize,itype)*fact_apdiam, & volumhi_sect(isize,itype), volumhi_sect(isize,itype)*fact_apdia3 end if end if ! load numb, mass, volu, dens back into mosaic_asecthp arrays l = numptr_aer(isize,itype,iphase) rbox(l) = xxnumb/fact_apnumbmr adrydens_box(isize,itype) = xxdens/fact_apdens adryqmas_box(isize,itype) = xxmass/fact_apmassmr adryqvol_box(isize,itype) = xxvolu/fact_apvolumr ! ! load number/mass/volume into working arrays ! ! *** NOTE *** ! num_distrib must be (#/cm3) ! vol_distrib units do not matter - they can be either masses or volumes, ! and either mixing ratios or concentrations ! (e.g., g/cm3-air, ug/kg-air, cm3/cm3-air, cm3/kg-air, ...) ! num_distrib(isize) = xxnumb*rhoair_box do ll = 1, ncomp_plustracer_aer(itype) l = massptr_aer(ll,isize,itype,iphase) tmpa = 0.0 if (l .ge. 1) tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr vol_distrib(isize,ll) = tmpa sumold(ll) = sumold(ll) + tmpa end do do llx = 1, 3 ll = (ncomp_coag-3) + llx tmpa = 0.0 if (llx .eq. 1) then l = hyswptr_aer(isize,itype) if (l .ge. 1) tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr else if (llx .eq. 2) then l = waterptr_aer(isize,itype) if (l .ge. 1) tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr else tmpa = max( 0.0_r8, adryqvol_box(isize,itype) )*fact_apvolumr end if vol_distrib(isize,ll) = tmpa sumold(ll) = sumold(ll) + tmpa end do ! calculate dry and wet diameters and wet density if (xxmass .le. 1.01*smallmassbb) then dpdry(isize) = dcen_sect(isize,itype)*fact_apdiam dpwet(isize) = dpdry(isize) denswet(isize) = xxdens else dpdry(isize) = (xxvolu/(xxnumb*piover6))**third dpdry(isize) = max( dpdry(isize), dlo_sect(isize,itype)*fact_apdiam ) l = waterptr_aer(isize,itype) if (l .ge. 1) then tmpa = max( 0.0_r8, rbox(l) )*fact_apmassmr xxmasswet = xxmass + tmpa xxvoluwet = xxvolu + tmpa/(dens_water_aer*fact_apdens) dpwet(isize) = (xxvoluwet/(xxnumb*piover6))**third dpwet(isize) = min( dpwet(isize), 30.0_r8*dhi_sect(isize,itype)*fact_apdiam ) denswet(isize) = (xxmasswet/xxvoluwet) else dpwet(isize) = dpdry(isize) denswet(isize) = xxdens end if end if vpdry(isize) = piover6 * (dpdry(isize)**3) end do isize_loop_aa ! ! make call to coagulation routine ! if (idiag_coag_onoff > 0) call coag_conserve_check( & nsize, maxd_asize, ncomp_coag, ncomp_coag_maxd, 1, lunout, & dpdry, num_distrib, vol_distrib, sum_num, sum_vol ) if ((mcoag_flag1 >= 10) .and. (mcoag_flag1 <= 19)) then call movingcenter_singletype_coag( & idiagbb, lunout_coag, nsize, maxd_asize, ncomp_coag, ncomp_coag_maxd, & temp_box, press_cgs, dtcoag, nsubstep_coag, & vpcut, vpdry, dpdry, dpwet, denswet, num_distrib, vol_distrib ) else call jacobson2002_singletype_coag( & nsize, maxd_asize, ncomp_coag, ncomp_coag_maxd, & temp_box, press_cgs, dtcoag, nsubstep_coag, & dpdry, dpwet, denswet, num_distrib, vol_distrib ) end if if (idiag_coag_onoff > 0) call coag_conserve_check( & nsize, maxd_asize, ncomp_coag, ncomp_coag_maxd, 2, lunout, & dpdry, num_distrib, vol_distrib, sum_num, sum_vol ) ! ! unload number/mass/volume from working arrays ! sumnew(:) = 0.0 isize_loop_xx: & do isize = 1, nsize do ll = 1, ncomp_coag sumnew(ll) = sumnew(ll) + max( 0.0_r8, vol_distrib(isize,ll) ) end do l = numptr_aer(isize,itype,iphase) rbox(l) = max( 0.0_r8, num_distrib(isize)/(rhoair_box*fact_apnumbmr) ) ! unload mass mixing ratios into rbox; also calculate dry mass and volume xxmass = 0.0 xxvolu = 0.0 do ll = 1, ncomp_aer(itype) l = massptr_aer(ll,isize,itype,iphase) if (l .ge. 1) then tmpa = max( 0.0_r8, vol_distrib(isize,ll) ) rbox(l) = tmpa/fact_apmassmr xxmass = xxmass + tmpa xxvolu = xxvolu + tmpa/(dens_aer(ll,itype)*fact_apdens) end if end do adryqmas_box(isize,itype) = xxmass/fact_apmassmr xxvolubb(isize) = xxvolu ll = (ncomp_coag-3) + 1 l = hyswptr_aer(isize,itype) if (l .ge. 1) rbox(l) = max( 0.0_r8, vol_distrib(isize,ll) )/fact_apmassmr ll = (ncomp_coag-3) + 2 l = waterptr_aer(isize,itype) if (l .ge. 1) rbox(l) = max( 0.0_r8, vol_distrib(isize,ll) )/fact_apmassmr ll = (ncomp_coag-3) + 3 adryqvol_box(isize,itype) = max( 0.0_r8, vol_distrib(isize,ll) )/fact_apvolumr end do isize_loop_xx ! check for mass and volume conservation do ll = 1, ncomp_coag tmpa = max( sumold(ll), sumnew(ll), 1.0e-35_r8 ) if (abs(sumold(ll)-sumnew(ll)) .gt. 1.0e-6*tmpa) then ncountbb(ll) = ncountbb(ll) + 1 ncountbb(0) = ncountbb(0) + 1 end if end do ! ! calculate new dry density, ! and check for new mean dry-size within bounds ! isize_loop_yy: & do isize = 1, nsize xxmass = adryqmas_box(isize,itype)*fact_apmassmr xxvolu = adryqvol_box(isize,itype)*fact_apvolumr l = numptr_aer(isize,itype,iphase) xxnumb = rbox(l)*fact_apnumbmr iconform_numb = 1 ! do a cautious calculation of density, using volume from coag routine if (xxmass .le. smallmassbb) then ncountaa(8) = ncountaa(8) + 1 xxdens = densdefault xxvolu = xxmass/xxdens xxnumb = xxmass/(volumcen_sect(isize,itype)*fact_apdia3*xxdens) xxdpav = dcen_sect(isize,itype)*fact_apdiam wwdens = xxdens wwdpav = xxdpav l = hyswptr_aer(isize,itype) if (l .ge. 1) rbox(l) = 0.0 l = waterptr_aer(isize,itype) if (l .ge. 1) rbox(l) = 0.0 iconform_numb = 0 if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-7a', & itstep, itype, isize, xxmass, xxvolu, xxdens else if (xxmass .gt. 1000.0*xxvolu) then ! in this case, density is too large. setting density=1000 forces ! next IF block while avoiding potential divide by zero or overflow xxdens = 1000.0 else xxdens = xxmass/xxvolu end if if ((xxdens .lt. 0.1) .or. (xxdens .gt. 20.0)) then ! (exception) case -- new dry density is unrealistic, ! so use dry volume from rbox instead of that from coag routine ncountaa(7) = ncountaa(7) + 1 if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,4e10.2)') 'coagaa-7b', & itstep, itype, isize, xxmass, xxvolu, xxdens xxvolu = xxvolubb(isize) xxdens = xxmass/xxvolu if (idiagbb .ge. 200) & write(lunout_coag,'(a,18x,1p,4e10.2)') 'coagaa-7c', & xxmass, xxvolu, xxdens end if ! check for mean size ok, and conform number if not if (iconform_numb .gt. 0) then if (xxnumb .gt. xxvolu/(volumlo_sect(isize,itype)*fact_apdia3)) then ncountaa(9) = ncountaa(9) + 1 tmpa = xxnumb xxnumb = xxvolu/(volumlo_sect(isize,itype)*fact_apdia3) xxdpav = dlo_sect(isize,itype)*fact_apdiam if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-9 ', & itstep, itype, isize, xxvolu, tmpa, xxnumb else if (xxnumb .lt. xxvolu/(volumhi_sect(isize,itype)*fact_apdia3)) then ncountaa(10) = ncountaa(10) + 1 tmpa = xxnumb xxnumb = xxvolu/(volumhi_sect(isize,itype)*fact_apdia3) xxdpav = dhi_sect(isize,itype)*fact_apdiam if (idiagbb .ge. 200) & write(lunout_coag,'(a,i10,2i4,1p,3e12.4)') 'coagcc-10', & itstep, itype, isize, xxvolu, tmpa, xxnumb else xxdpav = (xxvolu/(xxnumb*piover6))**third end if tmpb = 0.0 l = waterptr_aer(isize,itype) if (l .ge. 1) tmpb = max(0.0_r8,rbox(l))*fact_apmassmr wwmass = xxmass + tmpb wwvolu = xxvolu + tmpb/(dens_water_aer*fact_apdens) wwdens = wwmass/wwvolu wwdpav = xxdpav*((wwvolu/xxvolu)**third) end if ! load number, mass, volume, dry-density back into arrays l = numptr_aer(isize,itype,iphase) rbox(l) = xxnumb/fact_apnumbmr adrydens_box(isize,itype) = xxdens/fact_apdens adrydpav_box(isize,itype) = xxdpav/fact_apdiam adryqmas_box(isize,itype) = xxmass/fact_apmassmr adryqvol_box(isize,itype) = xxvolu/fact_apvolumr awetdens_box(isize,itype) = wwdens/fact_apdens awetdpav_box(isize,itype) = wwdpav/fact_apdiam end do isize_loop_yy ! temporary diagnostics if (idiagbb .ge. 100) then write(msg,93030) 'coagbb ncntaa ', ncountaa(1:10) call peg_message( lunerr, msg ) if (ncountbb(0) .gt. 0) then do llx = 1, (ncomp_coag+9)/10 llb = llx*10 lla = llb - 9 llb = min( llb, ncomp_coag) write(msg,93032) 'coagbb ncntbb', & mod(llx,10), ncountbb(lla:llb) call peg_message( lunerr, msg ) end do end if end if 93020 format( a, 1p, 10e10.2 ) 93030 format( a, 1p, 10i10 ) 93032 format( a, 1p, i1, 10i10 ) 2800 continue ! itype = 1, ntype_aer return end subroutine mosaic_coag_1box !---------------------------------------------------------------------- subroutine coag_conserve_check( & nbin, nbin_maxd, ncomp, ncomp_maxd, iflagaa, lunout, & dpdry, num_distrib, vol_distrib, sum_num, sum_vol ) implicit none ! ! checks on conservation of volume before/after coag routine ! integer, intent(in) :: nbin ! actual number of size bins integer, intent(in) :: nbin_maxd ! size-bin dimension for input (argument) arrays integer, intent(in) :: ncomp ! actual number of aerosol volume components integer, intent(in) :: ncomp_maxd ! volume-component dimension for input (argument) arrays integer, intent(in) :: lunout ! logical unit for warning & diagnostic output integer, intent(in) :: iflagaa ! real(r8), intent(in) :: dpdry(nbin_maxd) real(r8), intent(in) :: num_distrib(nbin_maxd) real(r8), intent(in) :: vol_distrib(nbin_maxd,ncomp_maxd) real(r8), intent(inout) :: sum_num(2) real(r8), intent(inout) :: sum_vol(0:ncomp_maxd,2) ! local variables integer :: i, l real(r8), parameter :: pi = 3.14159265358979323846_r8 real(r8) :: voldry(nbin) if (iflagaa == 1) then i = 1 else if (iflagaa == 2) then i = 2 else return end if ! sum initial or number and volumes voldry(1:nbin) = (pi/6.0)*(dpdry(1:nbin)**3) sum_num(i) = sum( num_distrib(1:nbin) ) sum_vol(0,i) = sum( num_distrib(1:nbin)*voldry(1:nbin) ) do l = 1, ncomp sum_vol(l,i) = sum( vol_distrib(1:nbin,l) ) end do ! write diagnostics if (iflagaa /= 2) return write(lunout,'(a)') write(lunout,'(2a,1p,2e14.6)') & 'number total ', ' initial, final =', & sum_num(1:2) write(lunout,'(2a,1p,2e14.6)') & 'number total ', ' (final/initial) =', & sum_num(2)/sum_num(1) write(lunout,'(2a,1p,2e14.6)') & 'volume total ', ' initial, final =', & sum_vol(0,1:2) write(lunout,'(2a,1p,2e14.6)') & 'volume total ', ' (initial/final)-1 =', & (sum_vol(0,1)/sum_vol(0,2))-1.0 do l = 1, ncomp if (abs(sum_vol(l,2)) > 0.0) then write(lunout,'(a,i4,a,1p,2e14.6)') & 'volume l=', l, ' (initial/final)-1 =', & (sum_vol(l,1)/sum_vol(l,2))-1.0 else write(lunout,'(a,i4,a,1p,2e14.6)') & 'volume l=', l, ' initial, final =', & sum_vol(l,1), sum_vol(l,2) end if end do return end subroutine coag_conserve_check !---------------------------------------------------------------------- subroutine movingcenter_singletype_coag( & idiagbb, lunout, nbin, nbin_maxd, ncomp, ncomp_maxd, & tempk, press, deltat, nsubstep_in, & vpcut, vpdry_in, dpdry, dpwet, denswet, & num_distrib, vol_distrib ) implicit none ! ! calculates aerosol coagulation for sectional representation ! and a single "particle type" using the single-type algorithm of ! jacobson (2002) ! ! reference ! jacobson, m. z., 2002, analysis of aerosol interactions with numerical ! techniques for solving coagulation, nucleation, condensation, ! dissolution, and reversible chemistry among multiple size ! distributions, j. geophys. res., 107, d19, 4366 ! ! IN arguments integer, intent(in) :: idiagbb integer, intent(in) :: lunout integer, intent(in) :: nbin ! actual number of size bins integer, intent(in) :: nbin_maxd ! size-bin dimension for input (argument) arrays integer, intent(in) :: ncomp ! actual number of aerosol volume components integer, intent(in) :: ncomp_maxd ! volume-component dimension for input (argument) arrays integer, intent(in) :: nsubstep_in ! number of time sub-steps for the integration real(r8), intent(in) :: tempk ! air temperature (K) real(r8), intent(in) :: press ! air pressure (dyne/cm2) real(r8), intent(in) :: deltat ! integration time (s) real(r8), intent(in) :: dpdry(nbin_maxd) ! bin current volume-mean dry diameter (cm) real(r8), intent(in) :: dpwet(nbin_maxd) ! bin wet (ambient) diameter (cm) real(r8), intent(in) :: denswet(nbin_maxd) ! bin wet (ambient) density (g/cm3) real(r8), intent(in) :: vpdry_in(nbin_maxd) ! bin current mean 1-particle dry volume (cm^3) real(r8), intent(in) :: vpcut(0:nbin_maxd) ! 1-particle dry volumes at bin boundaries (cm^3) ! INOUT arguments real(r8), intent(inout) :: num_distrib(nbin_maxd) ! num_distrib(i) = number concentration for bin i (#/cm3) real(r8), intent(inout) :: vol_distrib(nbin_maxd,ncomp_maxd) ! vol_distrib(i,l) = volume concentration for bin i, component l (cm3/cm3) ! ! NOTE -- The sectional representation is based on dry particle size. ! Dry sizes are used to calculate the receiving bin indices and product fractions ! for each (bin-i1, bin-i2) coagulation combination. ! Wet sizes and densities are used to calculate the brownian coagulation kernel. ! ! NOTE -- Units for num_distrib and vol_distrib ! num_distrib units MUST BE (#/cm3) ! vol_distrib ! (cm3/cm3) units are most consistent with the num_distrib units ! However, the algorithm works with almost any units! So vol_distrib can ! be either masses or volumes, and either mixing ratios or concentrations ! (e.g., g/cm3-air, ug/kg-air, um3/cm3-air, cm3/kg-air, ...). ! Use whatever is convenient. ! ! ! local variables ! integer :: i, isubstep, j, k, l, nsubstep real(r8), parameter :: frac_loss_limit = 0.5_r8 real(r8), parameter :: pi = 3.14159265358979323846_r8 real(r8) :: & del_num, del_voli, del_volj, deltat_sub, & tmpa, tmpb, tmp_beta, & vpdry_ipj real(r8) :: & betaxh(nbin,nbin), & cnum(nbin), cnum_old(nbin), & cvol(nbin,ncomp), cvol_old(nbin,ncomp), & vpdry(nbin) ! ! calculate brownian coagulation kernel ! call brownian_kernel( & nbin, tempk, press, dpwet, denswet, betaxh ) ! ! check that fractional number loss from any bin over deltat_sub does ! not exceed frac_loss_limit ! nsubstep = nsubstep_in deltat_sub = deltat/nsubstep ! time substep (s) tmpa = 0.0 do j = 1, nbin tmpb = 0.0 do i = 1, nbin tmpb = tmpb + betaxh(i,j)*num_distrib(i) end do tmpb = tmpb - 0.5*betaxh(j,j)*num_distrib(j) tmpa = max( tmpa, tmpb ) end do if (idiagbb > 0) write(lunout,'(/a,i10,1p,4e12.4)') 'coag aa nsub, dtsub, fracloss', & nsubstep, deltat_sub, tmpa*deltat_sub, frac_loss_limit if (tmpa*deltat_sub > frac_loss_limit) then tmpb = tmpa*deltat/frac_loss_limit isubstep = int(tmpb) + 1 tmpb = deltat/isubstep if (idiagbb > 0) write(lunout,'( a,i10,1p,4e12.4)') 'coag bb nsub, dtsub, fracloss', & isubstep, tmpb, tmpa*tmpb nsubstep = isubstep deltat_sub = deltat/nsubstep end if ! multiply coag kernel by deltat_sub ! write(92,'(a)') & ! ' dpdry_i (um) dpdry_j (um) coag coef (cm3/#/s)' do j = 1, nbin do i = 1, nbin ! write(92,'(1p,3e16.7)') & ! (1.0e4*dpwet(i)), (1.0e4*dpwet(j)), betaxh(i,j) betaxh(i,j) = betaxh(i,j) * deltat_sub end do end do ! initialize working arrays ! cnum(i) = current number conc (#/cm3) for bin i ! cvol(i,l) = current volume conc (cm3/cm3) for bin i, species l cnum(1:nbin) = num_distrib(1:nbin) cvol(1:nbin,1:ncomp) = vol_distrib(1:nbin,1:ncomp) ! ! solve coagulation over multiple time substeps ! solve_isubstep_loop: & do isubstep = 1, nsubstep ! write(91,*) 'maa =', 2 cnum_old(1:nbin) = cnum(1:nbin) cvol_old(1:nbin,1:ncomp) = cvol(1:nbin,1:ncomp) if (isubstep == 1) then vpdry(1:nbin) = vpdry_in(1:nbin) ! else ! *** ! *** should be updating vpdry(:) after first substep ! *** end if do j = 1, nbin do i = 1, j ! determine the bin that receives the "i+j" particles vpdry_ipj = vpdry(i) + vpdry(j) k = max( i, j ) do while (k > -99) if (vpdry_ipj > vpcut(k)) then k = k+1 if (k >= nbin) exit ! do while ... else if (vpdry_ipj < vpcut(k-1)) then k = k-1 if (k <= 1) exit ! do while ... else exit ! do while ... end if end do ! while ... k = max( 1, min( nbin, k ) ) ! write(183,'(3i4,1p,5e10.2)') i, j, k, vpdry(i), vpdry(j), vpdry_ipj ! write(183,'(12x,1p,5e10.2)') vpcut(i), vpcut(j), vpcut(k-1:k) tmp_beta = betaxh(i,j) ! now calc the coagulation changes to number and volume if (i == j) then del_num = 0.5*tmp_beta*cnum_old(i)*cnum_old(i) if (k == i) then ! here i + i --> i ! each coag event consumes 2 and produces 1 "i" particle cnum(i) = cnum(i) - del_num cycle else ! here i + i --> k /= i ! each coag event consumes 2 "i" and produces 1 "k" particle ! and there is transfer of mass" and produces 1 "k" particle cnum(i) = cnum(i) - del_num*2.0 cnum(k) = cnum(k) + del_num do l = 1, ncomp del_voli = tmp_beta*cnum_old(i)*cvol_old(i,l) cvol(i,l) = cvol(i,l) - del_voli cvol(k,l) = cvol(k,l) + del_voli end do end if else ! here i < j del_num = tmp_beta*cnum_old(i)*cnum_old(j) cnum(i) = cnum(i) - del_num if (k == j) then ! here i and j are different but j = k ! i loses number and transfers mass to k=j ! j=k does not lose any of its own number or mass do l = 1, ncomp del_voli = tmp_beta*cnum_old(j)*cvol_old(i,l) cvol(i,l) = cvol(i,l) - del_voli cvol(k,l) = cvol(k,l) + del_voli end do else ! here i, j, and k are all different ! both i and j lose number and transfer mass to k cnum(j) = cnum(j) - del_num cnum(k) = cnum(k) + del_num do l = 1, ncomp del_voli = tmp_beta*cnum_old(j)*cvol_old(i,l) del_volj = tmp_beta*cnum_old(i)*cvol_old(j,l) cvol(i,l) = cvol(i,l) - del_voli cvol(j,l) = cvol(j,l) - del_volj cvol(k,l) = cvol(k,l) + del_voli + del_volj end do end if end if end do ! i end do ! j end do solve_isubstep_loop ! transfer new concentrations from working arrays num_distrib(1:nbin) = cnum(1:nbin) do l = 1, ncomp vol_distrib(1:nbin,l) = cvol(1:nbin,l) end do ! all done return end subroutine movingcenter_singletype_coag !---------------------------------------------------------------------- subroutine jacobson2002_singletype_coag( & nbin, nbin_maxd, ncomp, ncomp_maxd, & tempk, press, deltat, nsubstep, & dpdry, dpwet, denswet, & num_distrib, vol_distrib ) implicit none ! ! calculates aerosol coagulation for sectional representation ! and a single "particle type" using the single-type algorithm of ! jacobson (2002) ! ! reference ! jacobson, m. z., 2002, analysis of aerosol interactions with numerical ! techniques for solving coagulation, nucleation, condensation, ! dissolution, and reversible chemistry among multiple size ! distributions, j. geophys. res., 107, d19, 4366 ! ! IN arguments integer, intent(in) :: nbin ! actual number of size bins integer, intent(in) :: nbin_maxd ! size-bin dimension for input (argument) arrays integer, intent(in) :: ncomp ! actual number of aerosol volume components integer, intent(in) :: ncomp_maxd ! volume-component dimension for input (argument) arrays integer, intent(in) :: nsubstep ! number of time sub-steps for the integration real(r8), intent(in) :: tempk ! air temperature (K) real(r8), intent(in) :: press ! air pressure (dyne/cm2) real(r8), intent(in) :: deltat ! integration time (s) real(r8), intent(in) :: dpdry(nbin_maxd) ! bin dry diameter (cm) real(r8), intent(in) :: dpwet(nbin_maxd) ! bin wet (ambient) diameter (cm) real(r8), intent(in) :: denswet(nbin_maxd) ! bin wet (ambient) density (g/cm3) ! INOUT arguments real(r8), intent(inout) :: num_distrib(nbin_maxd) ! num_distrib(i) = number concentration for bin i (#/cm3) real(r8), intent(inout) :: vol_distrib(nbin_maxd,ncomp_maxd) ! vol_distrib(i,l) = volume concentration for bin i, component l (cm3/cm3) ! ! NOTE -- The sectional representation is based on dry particle size. ! Dry sizes are used to calculate the receiving bin indices and product fractions ! for each (bin-i1, bin-i2) coagulation combination. ! Wet sizes and densities are used to calculate the brownian coagulation kernel. ! ! NOTE -- Units for num_distrib and vol_distrib ! num_distrib units MUST BE (#/cm3) ! vol_distrib ! (cm3/cm3) units are most consistent with the num_distrib units ! However, the algorithm works with almost any units! So vol_distrib can ! be either masses or volumes, and either mixing ratios or concentrations ! (e.g., g/cm3-air, ug/kg-air, um3/cm3-air, cm3/kg-air, ...). ! Use whatever is convenient. ! ! ! local variables ! integer :: i, isubstep, j, k, kp1, l integer :: & ilo_of_jk(nbin,nbin), & ihi_of_jk(nbin,nbin), & klo_of_ij(nbin,nbin), & khi_of_ij(nbin,nbin) real(r8), parameter :: pi = 3.14159265358979323846_r8 real(r8) :: & deltat_sub, & fijk_tmp, & t1xh_num, t3xh_num, & tmpa, & voldry_ipj real(r8) :: & betaxh(nbin,nbin), & cvol(nbin,ncomp), & cnum(nbin), cnum_old(nbin), & fijk_ieqk(nbin,nbin), & fijk_keqklo(nbin,nbin), & fijk_keqkhi(nbin,nbin), & t1xh_vol(ncomp), t3xh_vol(ncomp), & voldry(nbin) ! ! calculate brownian coagulation kernel ! call brownian_kernel( & nbin, tempk, press, dpwet, denswet, betaxh ) ! write(92,'(a)') & ! ' dpdry_i (um) dpdry_j (um) coag coef (cm3/#/s)' deltat_sub = deltat/nsubstep ! time substep (s) do i = 1, nbin do j = 1, nbin ! write(92,'(1p,3e16.7)') & ! (1.0e4*dpwet(i)), (1.0e4*dpwet(j)), betaxh(i,j) betaxh(i,j) = betaxh(i,j) * deltat_sub end do end do ! ! calculate the fijk of eqn 8 ! ! for each (i,j) pair, there are at most two k values for which fijk > 0 ! rather than calculating (and using) the full (3d) fijk array, ! only calculate the non-zero portions ! ! klo_of_ij(i,j) is the lowest k for which i + j --> k ! khi_of_ij(i,j) is the highest k for which i + j --> k ! for most (i,j), khi_of_ij = klo_of_ij + 1 ! fijk_keqklo(i,j) = fijk with k=klo_of_ij(i,j) ! fijk_keqkhi(i,j) = fijk with k=khi_of_ij(i,j) ! fijk_ieqk(i,j) = fijk with k=i ! fijk_ieqk(:,:) = 0.0 fijk_keqklo(:,:) = 0.0 fijk_keqkhi(:,:) = 0.0 klo_of_ij(:,:) = 0 khi_of_ij(:,:) = 0 do i = 1, nbin voldry(i) = (pi/6.0)*(dpdry(i)**3) end do do i = 1, nbin do j = 1, nbin voldry_ipj = voldry(i) + voldry(j) if (voldry_ipj >= voldry(nbin)) then if (i == nbin) fijk_ieqk(i,j) = 1.0 klo_of_ij(i,j) = nbin fijk_keqklo(i,j) = 1.0 else do k = 1, nbin - 1 kp1 = k+1 if ((voldry_ipj >= voldry(k)) .and. (voldry_ipj < voldry(kp1))) then fijk_tmp = ((voldry(kp1) - voldry_ipj)/(voldry(kp1) - voldry(k))) & * voldry(k) / voldry_ipj if (i == k) fijk_ieqk(i,j) = fijk_tmp klo_of_ij(i,j) = k fijk_keqklo(i,j) = fijk_tmp khi_of_ij(i,j) = kp1 fijk_keqkhi(i,j) = 1.0 - fijk_tmp endif end do endif end do ! j end do ! i ! for each j and k, determine the lowest and highest i ! for which i + j --> k is possible do k = 1, nbin do j = 1, k ilo_of_jk(j,k) = nbin+1 ihi_of_jk(j,k) = 0 do i = 1, k-1 if ((klo_of_ij(i,j) == k) .or. (khi_of_ij(i,j) == k)) then ilo_of_jk(j,k) = min( ilo_of_jk(j,k), i ) ihi_of_jk(j,k) = max( ihi_of_jk(j,k), i ) end if end do end do ! j end do ! k ! initialize working arrays ! cnum(i) = current number conc (#/cm3) for bin i ! cvol(i,l) = current volume conc (cm3/cm3) for bin i, species l cnum(1:nbin) = num_distrib(1:nbin) do l = 1, ncomp cvol(1:nbin,l) = vol_distrib(1:nbin,l) end do ! ! solve coagulation over multiple time substeps ! solve_isubstep_loop: & do isubstep = 1, nsubstep ! write(91,*) 'maa =', 2 cnum_old(1:nbin) = cnum(1:nbin) solve_k_loop: & do k = 1, nbin ! T3*h of eqns 7 and 9 t3xh_num = 0. if (k < nbin) then do j = 1, nbin t3xh_num = t3xh_num + (1.0-fijk_ieqk(k,j))*betaxh(k,j)*cnum_old(j) end do endif t3xh_vol(:) = t3xh_num ! T1*h of eqn 7 (for number) and eqn 9 (for volume) t1xh_num = 0. t1xh_vol(:) = 0. if (k > 1) then do j = 1, k do i = ilo_of_jk(j,k), ihi_of_jk(j,k) if (klo_of_ij(i,j) == k) then fijk_tmp = fijk_keqklo(i,j) else if (khi_of_ij(i,j) == k) then fijk_tmp = fijk_keqkhi(i,j) else cycle end if tmpa = fijk_tmp*betaxh(j,i)*cnum_old(j) t1xh_num = t1xh_num + tmpa*cnum(i)*voldry(i) ! write(91,'(a,3i4,1p,2e12.4)') & ! 'num source', i, j, k, fijk_tmp, t1xh_num do l = 1, ncomp t1xh_vol(l) = t1xh_vol(l) + tmpa*cvol(i,l) ! write(91,'(a,4i4,1p,2e12.4)') & ! 'vol so', l, i, j, k, fijk_tmp, t1xh_vol(l) end do end do end do t1xh_num = t1xh_num/voldry(k) endif ! new concentrations cnum(k) = (cnum(k) + t1xh_num) / (1.0 + t3xh_num) cvol(k,:) = (cvol(k,:) + t1xh_vol(:))/(1.0 + t3xh_vol(:)) end do solve_k_loop end do solve_isubstep_loop ! transfer new concentrations from working arrays num_distrib(1:nbin) = cnum(1:nbin) do l = 1, ncomp vol_distrib(1:nbin,l) = cvol(1:nbin,l) end do ! all done return end subroutine jacobson2002_singletype_coag !----------------------------------------------------------------------- subroutine brownian_kernel( & nbin, tempk, press, dpwet, denswet, bckernel ) ! ! this routine calculates brownian coagulation kernel ! using on eqn 16.28 of ! jacobson, m. z. (1999) fundamentals of atmospheric modeling. ! cambridge university press, new york, 656 pp. ! implicit none ! arguments integer, intent(in) :: nbin ! actual number of size bins real(r8), intent(in) :: tempk ! air temperature (K) real(r8), intent(in) :: press ! air pressure (dyne/cm2) real(r8), intent(in) :: dpwet(nbin) ! bin wet (ambient) diameter (cm) real(r8), intent(in) :: denswet(nbin) ! bin wet (ambient) density (g/cm3) real(r8), intent(out) :: bckernel(nbin,nbin) ! brownian coag kernel (cm3/#/s) ! local variables integer i, j real(r8), parameter :: pi = 3.14159265358979323846_r8 real(r8) & apfreepath, avogad, & boltz, & cunning, & deltasq_i, deltasq_j, & diffus_i, diffus_j, diffus_ipj, & dpwet_i, dpwet_j, dpwet_ipj, & gasfreepathx2, gasspeed, & knud, & mass_i, mass_j, mwair, & pi6, & rgas, rhoair, & speedsq_i, speedsq_j, & tmp1, & viscosd, viscosk ! boltz = boltzmann's constant (erg/K = g*cm2/s/K) ! avogad = avogadro's number (molecules/mol) ! mwair = molecular weight of air (g/mol) ! rgas = gas constant ((dyne/cm2)/(mol/cm3)/K) ! rhoair = air density (g/cm3) ! viscosd = air dynamic viscosity (g/cm/s) ! viscosk = air kinematic viscosity (cm2/s) ! gasspeed = air molecule mean thermal velocity (cm/s) ! gasfreepathx2 = air molecule mean free path (cm) x 2 pi6 = pi/6.0 boltz = 1.38065e-16 avogad = 6.02214e+23 mwair = 28.966 rgas = 8.31447e7 rhoair = press*mwair/(rgas*tempk) viscosd = 1.8325e-04*(416.16/(tempk+120.0)) * (tempk/296.16)**1.5 viscosk = viscosd/rhoair gasspeed = sqrt(8.0*boltz*tempk*avogad/(pi*mwair)) gasfreepathx2 = 4.0*viscosk/gasspeed ! coagulation kernel from eqn 16.28 of jacobson (1999) text ! ! diffus_i/j = particle brownian diffusion coefficient (cm2/s) ! speedsq_i/j = square of particle mean thermal velocity (cm/s) ! apfreepath = particle mean free path (cm) ! cunning = cunningham slip-flow correction factor ! deltasq_i/j = square of "delta_i" in eqn 16.29 do i = 1, nbin dpwet_i = dpwet(i) ! particle wet diameter (cm) mass_i = pi6*denswet(i)*(dpwet_i**3) ! particle wet mass (g) knud = gasfreepathx2/dpwet_i cunning = 1.0 + knud*(1.249 + 0.42*exp(-0.87/knud)) diffus_i = boltz*tempk*cunning/(3.0*pi*dpwet_i*viscosd) speedsq_i = 8.0*boltz*tempk/(pi*mass_i) apfreepath = 8.0*diffus_i/(pi*sqrt(speedsq_i)) tmp1 = (dpwet_i + apfreepath)**3 & - (dpwet_i*dpwet_i + apfreepath*apfreepath)**1.5 deltasq_i = ( tmp1/(3.0*dpwet_i*apfreepath) - dpwet_i )**2 do j = 1, nbin if (j >= i) then dpwet_j = dpwet(j) ! particle wet diameter (cm) mass_j = pi6*denswet(j)*(dpwet_j**3) ! particle wet mass (g) knud = gasfreepathx2/dpwet_j cunning = 1.0 + knud*(1.249 + 0.42*exp(-0.87/knud)) diffus_j = boltz*tempk*cunning/(3.0*pi*dpwet_j*viscosd) speedsq_j = 8.0*boltz*tempk/(pi*mass_j) apfreepath = 8.0*diffus_j/(pi*sqrt(speedsq_j)) tmp1 = (dpwet_j + apfreepath)**3 & - (dpwet_j*dpwet_j + apfreepath*apfreepath)**1.5 deltasq_j = ( tmp1/(3.0*dpwet_j*apfreepath) - dpwet_j )**2 dpwet_ipj = dpwet_i + dpwet_j diffus_ipj = diffus_i + diffus_j tmp1 = (dpwet_ipj/(dpwet_ipj + 2.0*sqrt(deltasq_i + deltasq_j))) & + (8.0*diffus_ipj/(dpwet_ipj*sqrt(speedsq_i + speedsq_j))) bckernel(i,j) = 2.0*pi*dpwet_ipj*diffus_ipj/tmp1 else bckernel(i,j) = bckernel(j,i) end if end do ! j end do ! i return end subroutine brownian_kernel !----------------------------------------------------------------------- end module module_mosaic_coag1d