module module_mosaic_coag3d use module_data_mosaic_kind use module_peg_util implicit none contains !----------------------------------------------------------------------- subroutine mosaic_coag_3d_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: mcoag_flag1, ifreq_coag !BSINGH use module_data_mosaic_aero, only: ifreq_coag !BSINGH 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 real(r8) :: 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,maxd_atype), dpwet(maxd_asize,maxd_atype) real(r8) :: densdry(maxd_asize,maxd_atype), denswet(maxd_asize,maxd_atype) real(r8) :: num_distrib(maxd_asize,maxd_atype) 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,maxd_atype,ncomp_coag_maxd) real(r8) :: vpcut(0:maxd_asize), vpdry(maxd_asize,maxd_atype) real(r8) :: xxvolubb(maxd_asize,maxd_atype) 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 write(lunout_coag,'(/a,i10)') 'mosaic_coag_3d_1box - itstep ', itstep ! 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 ncountaa(:) = 0 ! used for temporary diagnostics ncountbb(:) = 0 ! all types have same size cuts and same component species ! so nsize, ncomp_coag, densdefault and vpcut are same for all types itype = 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 itype_loop_aa: & do itype = 1, ntype_aer 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. 210) & 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,itype) = 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,itype,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,itype,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,itype) = dcen_sect(isize,itype)*fact_apdiam dpwet(isize,itype) = dpdry(isize,itype) denswet(isize,itype) = xxdens else dpdry(isize,itype) = (xxvolu/(xxnumb*piover6))**third dpdry(isize,itype) = max( dpdry(isize,itype), 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,itype) = (xxvoluwet/(xxnumb*piover6))**third dpwet(isize,itype) = min( dpwet(isize,itype), 30.0_r8*dhi_sect(isize,itype)*fact_apdiam ) denswet(isize,itype) = (xxmasswet/xxvoluwet) else dpwet(isize,itype) = dpdry(isize,itype) denswet(isize,itype) = xxdens end if end if densdry(isize,itype) = xxdens vpdry(isize,itype) = piover6 * (dpdry(isize,itype)**3) end do isize_loop_aa end do itype_loop_aa ! ! make call to coagulation routine ! if (idiag_coag_onoff > 0) call coag_3d_conserve_check( & nsize, maxd_asize, ntype_aer, maxd_atype, ncomp_coag, ncomp_coag_maxd, 1, lunout_coag, & dpdry, num_distrib, vol_distrib, sum_num, sum_vol ) write(*,'(/a/)') 'mosaic_coag_3d_1box calling movingcenter_coag_3d' call movingcenter_coag_3d( & idiagbb, lunout_coag, nsize, maxd_asize, ntype_aer, maxd_atype, iphase, & ncomp_coag, ncomp_coag_maxd, & temp_box, press_cgs, dtcoag, nsubstep_coag, & vpcut, vpdry, dpdry, dpwet, densdry, denswet, num_distrib, vol_distrib ) write(*,'(/a/)') 'mosaic_coag_3d_1box backfrm movingcenter_coag_3d' if (idiag_coag_onoff > 0) call coag_3d_conserve_check( & nsize, maxd_asize, ntype_aer, maxd_atype, ncomp_coag, ncomp_coag_maxd, 2, lunout_coag, & dpdry, num_distrib, vol_distrib, sum_num, sum_vol ) ! ! unload number/mass/volume from working arrays ! sumnew(:) = 0.0 itype_loop_yy: & do itype = 1, ntype_aer isize_loop_xx: & do isize = 1, nsize do ll = 1, ncomp_coag sumnew(ll) = sumnew(ll) + max( 0.0_r8, vol_distrib(isize,itype,ll) ) end do l = numptr_aer(isize,itype,iphase) rbox(l) = max( 0.0_r8, num_distrib(isize,itype)/(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,itype,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,itype) = xxvolu ll = (ncomp_coag-3) + 1 l = hyswptr_aer(isize,itype) if (l .ge. 1) rbox(l) = max( 0.0_r8, vol_distrib(isize,itype,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,itype,ll) )/fact_apmassmr ll = (ncomp_coag-3) + 3 adryqvol_box(isize,itype) = max( 0.0_r8, vol_distrib(isize,itype,ll) )/fact_apvolumr end do isize_loop_xx ! ! 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. 210) & 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,itype) 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 dpdry(isize,itype) = xxdpav end do isize_loop_yy end do itype_loop_yy if (idiag_coag_onoff > 0) call coag_3d_conserve_check( & nsize, maxd_asize, ntype_aer, maxd_atype, ncomp_coag, ncomp_coag_maxd, 3, lunout_coag, & dpdry, num_distrib, vol_distrib, sum_num, sum_vol ) ! temporary diagnostics do ll = 1, ncomp_coag ! ncountbb > 0 when mass/volume conservation relative-error > 1e-6 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 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 ) return end subroutine mosaic_coag_3d_1box !---------------------------------------------------------------------- subroutine coag_3d_conserve_check( & nbin, nbin_maxd, ntype, ntype_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) :: ntype ! actual number of aerosol types integer, intent(in) :: ntype_maxd ! aerosol type 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,ntype_maxd) real(r8), intent(in) :: num_distrib(nbin_maxd,ntype_maxd) real(r8), intent(in) :: vol_distrib(nbin_maxd,ntype_maxd,ncomp_maxd) real(r8), intent(inout) :: sum_num(2) real(r8), intent(inout) :: sum_vol(0:ncomp_maxd,2) ! local variables integer :: i, itype, l real(r8), parameter :: pi = 3.14159265358979323846_r8 real(r8) :: voldry(nbin,ntype) if (iflagaa == 1) then i = 1 else if (iflagaa == 2) then i = 2 else if (iflagaa == 3) then i = 2 else return end if ! sum initial or number and volumes voldry(1:nbin,1:ntype) = (pi/6.0)*(dpdry(1:nbin,1:ntype)**3) sum_num(i) = sum( num_distrib(1:nbin,1:ntype) ) sum_vol(0,i) = sum( num_distrib(1:nbin,1:ntype)*voldry(1:nbin,1:ntype) ) do l = 1, ncomp sum_vol(l,i) = sum( vol_distrib(1:nbin,1:ntype,l) ) end do ! write diagnostics if ((iflagaa /= 2) .and. (iflagaa /= 3)) return write(lunout,'(/a,i3)') 'coag_3d_conserve_check -- iflagaa = ', iflagaa 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) ! with jacobson coag algorithm, dry-size of each bin is unchanged ! by the coagulation ! with moving-center coag algorithm, dry-size does change, ! so sum_vol(0,i) is currently not correct for iflagaa=2, because ! the incoming voldry() is currently the initial (i=1) value ! so do not output this when iflagaa=2 if (iflagaa == 3) then 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 end if 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_3d_conserve_check !---------------------------------------------------------------------- subroutine movingcenter_coag_3d( & idiagbb, lunout, nbin, nbin_maxd, ntype, ntype_maxd, iphase, & ncomp, ncomp_maxd, & tempk, press, deltat, nsubstep_in, & vpcut, vpdry_in, dpdry, dpwet, densdry, denswet, & num_distrib, vol_distrib ) use module_data_mosaic_aero, only: & method_bcfrac, method_kappa, msectional_flag2 use module_data_mosaic_asecthp, only: & dens_aer, hygro_aer, & itype_md1_of_itype, itype_md2_of_itype, itype_of_itype_md1md2, & lptr_bc_aer, massptr_aer, & ncomp_aer, ntype_aer, ntype_md1_aer, ntype_md2_aer, & smallmassbb, xcut_atype_md1, xcut_atype_md2 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) :: ntype ! actual number of aerosol types integer, intent(in) :: ntype_maxd ! aerosol type dimension for input (argument) arrays integer, intent(in) :: iphase 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,ntype_maxd) ! bin current volume-mean dry diameter (cm) real(r8), intent(in) :: dpwet(nbin_maxd,ntype_maxd) ! bin wet (ambient) diameter (cm) real(r8), intent(in) :: densdry(nbin_maxd,ntype_maxd) ! bin dry density (g/cm3) real(r8), intent(in) :: denswet(nbin_maxd,ntype_maxd) ! bin wet (ambient) density (g/cm3) real(r8), intent(in) :: vpdry_in(nbin_maxd,ntype_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,ntype_maxd) ! num_distrib(i) = number concentration for bin i (#/cm3) real(r8), intent(inout) :: vol_distrib(nbin_maxd,ntype_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 ! character(len=256) :: errmsg integer :: i, isubstep, itype, itype1, itype2, iwrite193 integer :: j, jtype, jtype1, jtype2 integer :: k, ktype, ktype1, ktype2, ktmp1, ktmp2 integer :: l, method_ktype, nsubstep real(r8), parameter :: frac_loss_limit = 0.8_r8 real(r8), parameter :: pi = 3.14159265358979323846_r8 real(r8) :: del_num, del_voli, del_volj, deltat_sub real(r8) :: tmp_bcfrac, tmp_beta real(r8) :: tmp_densbc real(r8) :: tmp_kappa, tmp_kappaxvolu real(r8) :: tmp_massbc, tmp_massdry real(r8) :: tmp_volu real(r8) :: tmpa, tmpb, tmpc, tmpi, tmpj real(r8) :: vpdry_ipj real(r8) :: & betaxh(nbin,nbin), & bcfrac_bin(nbin,ntype), & cnum(nbin,ntype), cnum_old(nbin,ntype), & cvol(nbin,ntype,ncomp), cvol_old(nbin,ntype,ncomp), & kappa_bin(nbin,ntype), & tmpvecb(nbin), & vpdry(nbin,ntype) ! set method for determining ktype method_ktype = 0 if ((ntype > 1) .and. (msectional_flag2 > 0)) then if (ncomp /= ncomp_aer(1) + 3) then write(errmsg,'(/2a,2(1x,i5)/)') '*** movingcenter_coag_3d - ', & 'mismatch of ncomp & ncomp_aer(1)', ncomp, ncomp_aer(1) call wrf_error_fatal(trim(adjustl(errmsg))) end if if ((method_bcfrac /= 1) .and. (method_bcfrac /= 11)) then write(errmsg,'(/2a,2(1x,i5)/)') '*** movingcenter_coag_3d - ', & 'method_bcfrac must be 1 or 11 but is', method_bcfrac call wrf_error_fatal(trim(adjustl(errmsg))) end if if ((method_kappa /= 11) .and. (method_kappa /= 12)) then write(errmsg,'(/2a,2(1x,i5)/)') '*** movingcenter_coag_3d - ', & 'method_kappa must be 11 or 12 but is', method_kappa call wrf_error_fatal(trim(adjustl(errmsg))) end if method_ktype = 1 end if ! method_ktype = 0 iwrite193 = 0 ! tmp_densbc = dens_aer(ibc_a,1) ! this only works in mosaic box model tmp_densbc = -2.0 do itype = 1, ntype do l = 1, ncomp_aer(itype) if (massptr_aer(l,1,itype,iphase) == lptr_bc_aer(1,itype,iphase)) then tmp_densbc = dens_aer(l,itype) end if end do end do if (tmp_densbc < -1.0) then call wrf_error_fatal('*** movingcenter_coag_3d - cannot find bc') end if ! ! 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 jtype = 1, ntype tmpvecb(:) = 0.0 do itype = 1, ntype ! calculate brownian coagulation kernel for itype,jtype call brownian_kernel_3dsub( & nbin, nbin_maxd, ntype, ntype_maxd, itype, jtype, & tempk, press, dpwet, denswet, betaxh ) do j = 1, nbin do i = 1, nbin if ((i /= j) .or. (itype /= jtype)) then tmpvecb(j) = tmpvecb(j) + betaxh(i,j)*num_distrib(i,itype) else tmpvecb(j) = tmpvecb(j) + betaxh(i,j)*num_distrib(i,itype)*0.5 end if end do ! i end do ! j end do ! itype do j = 1, nbin tmpa = max( tmpa, tmpvecb(j) ) end do ! j end do ! jtype 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 ! 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,1:ntype) = num_distrib(1:nbin,1:ntype) cvol(1:nbin,1:ntype,1:ncomp) = vol_distrib(1:nbin,1:ntype,1:ncomp) ! ! solve coagulation over multiple time substeps ! solve_isubstep_loop: & do isubstep = 1, nsubstep ! write(91,*) 'maa =', 2 ! copy current num, vol values to "old" arrays cnum_old(1:nbin,1:ntype) = cnum(1:nbin,1:ntype) cvol_old(1:nbin,1:ntype,1:ncomp) = cvol(1:nbin,1:ntype,1:ncomp) ! calc 1-particle dry volume if (isubstep == 1) then vpdry(1:nbin,1:ntype) = vpdry_in(1:nbin,1:ntype) ! else ! *** ! *** should be updating vpdry(:) after first substep ! *** ! *** also densdry(:,:), but only if it is needed in calc of ktype ! *** (which is when method_bcfrac=1) ! *** end if if (method_ktype == 1) then ! calc information needed for getting ktype do itype = 1, ntype do i = 1, nbin tmp_massbc = 0.0 tmp_massdry = 0.0 tmp_kappaxvolu = 0.0 tmp_volu = 0.0 do l = 1, ncomp_aer(itype) tmpa = vol_distrib(i,itype,l) tmpc = tmpa/dens_aer(l,itype) if (method_bcfrac == 11) then ! bcfrac is dry-volume fraction tmpb = tmpc else ! bcfrac is dry-mass fraction tmpb = tmpa end if tmp_massdry = tmp_massdry + tmpb ! if (l == ibc_a) then ! this only works in mosaic box model, when the ordering ! ! of species in aer(...) and massptr_aer(...) are the same if (massptr_aer(l,i,itype,iphase) == lptr_bc_aer(i,itype,iphase)) then tmp_massbc = tmpb if (method_kappa == 12) cycle ! in this case, kappa includes non-bc species only end if tmp_volu = tmp_volu + tmpc tmp_kappaxvolu = tmp_kappaxvolu + tmpc*hygro_aer(l,itype) end do ! l if (tmp_massdry <= smallmassbb) then itype1 = itype_md1_of_itype( itype ) bcfrac_bin(i,itype) = 0.5*( xcut_atype_md1(itype1-1) + xcut_atype_md1(itype1) ) else bcfrac_bin(i,itype) = tmp_massbc/tmp_massdry end if bcfrac_bin(i,itype) = max( 0.0_r8, xcut_atype_md1(0), & bcfrac_bin(i,itype) ) bcfrac_bin(i,itype) = min( 1.0_r8, xcut_atype_md1(ntype_md1_aer), & bcfrac_bin(i,itype) ) if (tmp_volu <= smallmassbb) then itype2 = itype_md2_of_itype( itype ) kappa_bin(i,itype) = sqrt( max( 0.0_r8, & xcut_atype_md2(itype2-1) * xcut_atype_md2(itype2) ) ) else kappa_bin(i,itype) = tmp_kappaxvolu/tmp_volu end if kappa_bin(i,itype) = max( 1.0e-10_r8, xcut_atype_md2(0), & kappa_bin(i,itype) ) kappa_bin(i,itype) = min( 2.0_r8, xcut_atype_md2(ntype_md2_aer), & kappa_bin(i,itype) ) end do ! i end do ! itype end if ! (method_ktype == 1) ! main loops over i/jtype and i/j(size) do jtype = 1, ntype jtype1 = itype_md1_of_itype( jtype ) jtype2 = itype_md2_of_itype( jtype ) do itype = 1, jtype ! itype <= jtype avoids double counting itype1 = itype_md1_of_itype( itype ) itype2 = itype_md2_of_itype( itype ) ! calculate brownian coagulation kernel again for itype,jtype ! problem is that it can be too big ! call brownian_kernel_3dsub( & nbin, nbin_maxd, ntype, ntype_maxd, itype, jtype, & tempk, press, dpwet, denswet, betaxh ) if (iwrite193 > 0) & write(193,'(/a,2(4x,2i4))') 'i/jtype1 then ...2', itype1, jtype1, itype2, jtype2 do j = 1, nbin do i = 1, nbin if (itype == jtype) then if (i > j) cycle ! i <= j avoids double counting when itype = jtype end if ! determine the bin that receives the "i+j" particles vpdry_ipj = vpdry(i,itype) + vpdry(j,jtype) 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 ) ) if (method_ktype == 1) then if (method_bcfrac == 11) then tmpi = vpdry(i,itype) tmpj = vpdry(j,jtype) else tmpi = vpdry(i,itype)*densdry(i,itype) tmpj = vpdry(j,jtype)*densdry(j,jtype) end if tmpa = tmpi/max( (tmpi + tmpj), 1.0e-35_r8 ) tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) ) tmp_bcfrac = bcfrac_bin(i,itype)*tmpa + bcfrac_bin(j,jtype)*(1.0_r8 - tmpa) ktype1 = ntype_md1_aer do ktmp1 = 1, ntype_md1_aer - 1 if (tmp_bcfrac <= xcut_atype_md1(ktmp1)) then ktype1 = ktmp1 exit end if end do if (method_kappa == 11) then tmpi = vpdry(i,itype) tmpj = vpdry(j,jtype) else if (method_bcfrac == 11) then tmpa = 1.0_r8 - bcfrac_bin(i,itype) tmpi = vpdry(i,itype)*max(tmpa,1.0e-10_r8) tmpa = 1.0_r8 - bcfrac_bin(j,jtype) tmpj = vpdry(j,jtype)*max(tmpa,1.0e-10_r8) else tmpa = 1.0_r8 - bcfrac_bin(i,itype)*(densdry(i,itype)/tmp_densbc) tmpi = vpdry(i,itype)*max(tmpa,1.0e-10_r8) tmpa = 1.0_r8 - bcfrac_bin(j,jtype)*(densdry(j,jtype)/tmp_densbc) tmpj = vpdry(j,jtype)*max(tmpa,1.0e-10_r8) end if tmpa = tmpi/max( (tmpi + tmpj), 1.0e-35_r8 ) tmpa = max( 0.0_r8, min( 1.0_r8, tmpa ) ) tmp_kappa = kappa_bin(i,itype)*tmpa + kappa_bin(j,jtype)*(1.0_r8 - tmpa) ktype2 = ntype_md2_aer do ktmp2 = 1, ntype_md2_aer - 1 if (tmp_kappa <= xcut_atype_md2(ktmp2)) then ktype2 = ktmp2 exit end if end do ktype = itype_of_itype_md1md2(ktype1,ktype2) if (iwrite193 > 0) then ! if ( (num_distrib(i,itype) >= 1.0e30) .and. & ! (num_distrib(j,jtype) >= 1.0e-1) ) then ! write(193,'(a,3(2x,3i3),2f9.4)') 'i(s,t1,t2); j; k; wbc, kap', & ! i, itype_md1_of_itype(itype), itype_md2_of_itype(itype), & ! j, itype_md1_of_itype(jtype), itype_md2_of_itype(jtype), & ! k, ktype1, ktype2, tmp_bcfrac, tmp_kappa ! end if ! if ( (j==6 .or. j==12 .or. j==18 .or. j==24) .and. & ! (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. & ! ( ( (jtype1==1 .or. jtype1==5 .or. jtype1==10 .or. jtype1==15) .and. & ! (jtype1-itype1 <= 1) ) .or. & ! ( (jtype1==2 .or. jtype1==6 .or. jtype1==11 .or. jtype1==16) .and. & ! (jtype1-itype1 == 1) ) ) .and. & ! ( ( (jtype2==1 .or. jtype2==3 .or. jtype2==6 .or. jtype2== 9) .and. & ! (jtype2-itype2 <= 1) ) .or. & ! ( (jtype2==2 .or. jtype2==4 .or. jtype2==7 .or. jtype2==10) .and. & ! (jtype1-itype1 == 1) ) ) ) then if ( ( j==6 .or. j==12 .or. j==18 .or. j==24 ) .and. & ( abs(i-j)<=4 .or. abs(i-j)==8 ) .and. & ( jtype1==1 .or. jtype1==5 .or. jtype1==10 .or. jtype1==15 .or. & itype1==1 .or. itype1==5 .or. itype1==10 .or. itype1==15 ) .and. & ( abs(jtype1-itype1)<=1 .or. abs(jtype1-itype1)==5 .or. abs(jtype1-itype1)==10 ) .and. & ( jtype2==1 .or. jtype2==5 .or. jtype2==9 .or. & itype2==1 .or. itype2==5 .or. itype2==9 ) .and. & ( abs(jtype2-itype2)<=1 .or. abs(jtype2-itype2)==4 .or. abs(jtype2-itype2)==8 ) ) then write(193,'(a,3i3,1p,5e9.2,0p,2f6.3,3i3,3f6.3,3i3,3f7.4)') 'ijk&types', & i, j, k, dpdry(i,itype), dpdry(j,jtype), vpdry(i,itype), vpdry(j,jtype), vpdry_ipj, & densdry(i,itype), densdry(j,jtype), & itype1, jtype1, ktype1, bcfrac_bin(i,itype), bcfrac_bin(j,jtype), tmp_bcfrac, & itype2, jtype2, ktype2, kappa_bin(i,itype), kappa_bin(j,jtype), tmp_kappa end if end if ! (iwrite193 > 0) else if (itype == jtype) then ktype = itype else ktype = ntype end if end if ! 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) * deltat_sub ! now calc the coagulation changes to number and volume if ((i == j) .and. (itype == jtype)) then del_num = 0.5*tmp_beta*cnum_old(i,itype)*cnum_old(i,itype) if ((i == k) .and. (itype == ktype)) then ! here i,itype + i,itype --> i,itype ! each coag event consumes 2 and produces 1 "i" particle cnum(i,itype) = cnum(i,itype) - del_num cycle else ! here i,itype + i,itype --> k,ktype AND i,k and/or itype,ktype differ ! each coag event consumes 2 "i" and produces 1 "k" particle ! and there is transfer of mass" and produces 1 "k" particle cnum(i,itype) = cnum(i,itype) - del_num*2.0 cnum(k,ktype) = cnum(k,ktype) + del_num do l = 1, ncomp del_voli = tmp_beta*cnum_old(i,itype)*cvol_old(i,itype,l) cvol(i,itype,l) = cvol(i,itype,l) - del_voli cvol(k,ktype,l) = cvol(k,ktype,l) + del_voli end do end if else ! here i < j OR itype < jtype del_num = tmp_beta*cnum_old(i,itype)*cnum_old(j,jtype) cnum(i,itype) = cnum(i,itype) - del_num if ((j == k) .and. (jtype == ktype)) then ! here i,itype and j,jtype are different but j,jtype = k,ktype ! 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,jtype)*cvol_old(i,itype,l) cvol(i,itype,l) = cvol(i,itype,l) - del_voli cvol(k,ktype,l) = cvol(k,ktype,l) + del_voli end do else ! here i,itype; j,jtype; and k,ktype are all different ! both i and j lose number and transfer mass to k cnum(j,jtype) = cnum(j,jtype) - del_num cnum(k,ktype) = cnum(k,ktype) + del_num do l = 1, ncomp del_voli = tmp_beta*cnum_old(j,jtype)*cvol_old(i,itype,l) del_volj = tmp_beta*cnum_old(i,itype)*cvol_old(j,jtype,l) cvol(i,itype,l) = cvol(i,itype,l) - del_voli cvol(j,jtype,l) = cvol(j,jtype,l) - del_volj cvol(k,ktype,l) = cvol(k,ktype,l) + del_voli + del_volj end do end if end if end do ! i end do ! j end do ! itype end do ! jtype end do solve_isubstep_loop ! transfer new concentrations from working arrays num_distrib(1:nbin,1:ntype) = cnum(1:nbin,1:ntype) do l = 1, ncomp vol_distrib(1:nbin,1:ntype,l) = cvol(1:nbin,1:ntype,l) end do ! all done return end subroutine movingcenter_coag_3d !----------------------------------------------------------------------- subroutine brownian_kernel_3dsub( & nbin, nbin_maxd, ntype, ntype_maxd, itype, jtype, & 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 integer, intent(in) :: nbin_maxd ! size-bin dimension for input (argument) arrays integer, intent(in) :: ntype ! actual number of aerosol types integer, intent(in) :: ntype_maxd ! aerosol type dimension for input (argument) arrays integer, intent(in) :: itype, jtype ! the itype & jtype for which to calc the kernal real(r8), intent(in) :: tempk ! air temperature (K) real(r8), intent(in) :: press ! air pressure (dyne/cm2) real(r8), intent(in) :: dpwet(nbin_maxd,ntype_maxd) ! bin wet (ambient) diameter (cm) real(r8), intent(in) :: denswet(nbin_maxd,ntype_maxd) ! 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 bckernel(:,:) = -2.0_r8 do i = 1, nbin dpwet_i = dpwet(i,itype) ! particle wet diameter (cm) mass_i = pi6*denswet(i,itype)*(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 ( (itype == jtype) .and. & (bckernel(j,i) > -1.0_r8) ) then ! value already calculated bckernel(i,j) = bckernel(j,i) else dpwet_j = dpwet(j,jtype) ! particle wet diameter (cm) mass_j = pi6*denswet(j,jtype)*(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) = max( 0.0_r8, 2.0*pi*dpwet_ipj*diffus_ipj/tmp1 ) end if end do ! j end do ! i return end subroutine brownian_kernel_3dsub !----------------------------------------------------------------------- end module module_mosaic_coag3d