!********************************************************************************** ! This computer software was prepared by Battelle Memorial Institute, hereinafter ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY, ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. ! ! MOSAIC module: see module_mosaic_driver.F for references and terms of use !********************************************************************************** module module_mosaic_coag use module_peg_util implicit none contains !----------------------------------------------------------------------- subroutine mosaic_coag_1clm( istat_coag, & it, jt, kclm_calcbgn, kclm_calcend, & idiagbb_in, & dtchem, dtcoag_in, & id, ktau, ktauc, its, ite, jts, jte, kts, kte ) ! ! calculates aerosol particle coagulation for grid points in ! the i=it, j=jt vertical column ! over timestep dtcoag_in ! ! uses a version of subr. coagsolv (see additional disclaimer below) ! that was obtained from mark jacobson in june 2003, ! and then modified to work with mosaic aerosol routines ! use module_data_mosaic_asect use module_data_mosaic_other use module_state_description, only: param_first_scalar ! subr arguments integer, intent(inout) :: istat_coag ! =0 if no problems integer, intent(in) :: & it, jt, kclm_calcbgn, kclm_calcend, & idiagbb_in, & id, ktau, ktauc, its, ite, jts, jte, kts, kte real, intent(in) :: dtchem, dtcoag_in ! NOTE - much information is passed via arrays in ! module_data_mosaic_asect and module_data_mosaic_other ! ! rsub (inout) - aerosol mixing ratios ! aqmassdry_sub, aqvoldry_sub (inout) - ! aerosol dry-mass and dry-volume mixing ratios ! adrydens_sub (inout) - aerosol dry density ! rsub(ktemp,:,:), ptotclm, cairclm (in) - ! air temperature, pressure, and molar density ! local variables integer, parameter :: coag_method = 1 integer, parameter :: ncomp_coag_maxd = maxd_acomp + 3 integer :: k, l, ll, lla, llb, llx, m integer :: isize, itype, iphase integer :: iconform_numb integer :: idiagbb, idiag_coag_onoff, iok integer :: lunout_coag integer :: ncomp_coag, nsize, nsubstep_coag, ntau_coag integer,save :: ncountaa(10), ncountbb(0:ncomp_coag_maxd) integer :: p1st real, parameter :: densdefault = 2.0 real, parameter :: factsafe = 1.00001 real, parameter :: onethird = 1.0/3.0 real, parameter :: piover6 = pi/6.0 real, parameter :: smallmassbb = 1.0e-30 real :: cair_box real :: dtcoag real :: duma, dumb, dumc real :: patm_box real :: temp_box real :: xxdens, xxmass, xxnumb, xxvolu real :: xxmasswet, xxvoluwet real :: dpdry(maxd_asize), dpwet(maxd_asize), denswet(maxd_asize) real :: num_distrib(maxd_asize) real :: sumnew(ncomp_coag_maxd), sumold(ncomp_coag_maxd) real :: vol_distrib(maxd_asize,ncomp_coag_maxd) real :: xxvolubb(maxd_asize) character(len=120) :: msg istat_coag = 0 lunout_coag = 6 if (lunout .gt. 0) lunout_coag = lunout ! when dtcoag_in > dtchem, do not perform coagulation calcs ! on every chemistry time step ntau_coag = nint( dtcoag_in/dtchem ) ntau_coag = max( 1, ntau_coag ) if (mod(ktau,ntau_coag) .ne. 0) return dtcoag = dtchem*ntau_coag ! set variables that do not change p1st = PARAM_FIRST_SCALAR idiagbb = idiagbb_in ! NOTE - code currently just does 1 type itype = 1 iphase = ai_phase nsize = nsize_aer(itype) ncomp_coag = ncomp_plustracer_aer(itype) + 3 ! loop over subareas (currently only 1) and vertical levels do 2900 m = 1, nsubareas do 2800 k = kclm_calcbgn, kclm_calcend ! temporary diagnostics if ((it .eq. its) .and. & (jt .eq. jts) .and. (k .eq. kclm_calcbgn)) then ncountaa(:) = 0 ncountbb(:) = 0 end if ncountaa(1) = ncountaa(1) + 1 if (afracsubarea(k,m) .lt. 1.e-4) goto 2700 cair_box = cairclm(k) temp_box = rsub(ktemp,k,m) patm_box = ptotclm(k)/1.01325e6 nsubstep_coag = 1 idiag_coag_onoff = -1 ! do initial calculation of dry mass, volume, and density, ! and initial number conforming (as needed) vol_distrib(:,:) = 0.0 sumold(:) = 0.0 do isize = 1, nsize l = numptr_aer(isize,itype,iphase) xxnumb = rsub(l,k,m) xxmass = aqmassdry_sub(isize,itype,k,m) xxvolu = aqvoldry_sub( isize,itype,k,m) xxdens = adrydens_sub( isize,itype,k,m) iconform_numb = 1 duma = max( abs(xxmass), abs(xxvolu*xxdens), 0.1*smallmassbb ) dumb = abs(xxmass - xxvolu*xxdens)/duma if ( (xxdens .lt. 0.1) .or. (xxdens .gt. 20.0) & .or. (dumb .gt. 1.0e-4) ) then ! (exception) case of drydensity not valid, or mass /= volu*dens ! so compute dry mass and volume from rsub ncountaa(2) = ncountaa(2) + 1 if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-2a', & ktau, it, jt, k, isize, xxmass, xxvolu, xxdens xxmass = 0.0 xxvolu = 0.0 do ll = 1, ncomp_aer(itype) l = massptr_aer(ll,isize,itype,iphase) if (l .ge. p1st) then duma = max( 0.0, rsub(l,k,m) )*mw_aer(ll,itype) xxmass = xxmass + duma xxvolu = xxvolu + duma/dens_aer(ll,itype) end if end do if (xxmass .gt. 0.99*smallmassbb) then xxdens = xxmass/xxvolu xxvolu = xxmass/xxdens if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-2b', & ktau, it, jt, k, isize, xxmass, xxvolu, xxdens end if 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)*xxdens) l = hyswptr_aer(isize,itype) if (l .ge. p1st) rsub(l,k,m) = 0.0 l = waterptr_aer(isize,itype) if (l .ge. p1st) rsub(l,k,m) = 0.0 iconform_numb = 0 if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-2c', & ktau, it, jt, k, 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))) then ncountaa(4) = ncountaa(4) + 1 duma = xxnumb xxnumb = xxvolu/(factsafe*volumlo_sect(isize,itype)) if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-4 ', & ktau, it, jt, k, isize, xxvolu, duma, xxnumb else if (xxnumb .lt. & xxvolu*factsafe/volumhi_sect(isize,itype)) then ncountaa(5) = ncountaa(5) + 1 duma = xxnumb xxnumb = xxvolu*factsafe/volumhi_sect(isize,itype) if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-5 ', & ktau, it, jt, k, isize, xxvolu, duma, xxnumb end if end if ! load numb, mass, volu, dens back into mosaic_asect arrays l = numptr_aer(isize,itype,iphase) rsub(l,k,m) = xxnumb adrydens_sub( isize,itype,k,m) = xxdens aqmassdry_sub(isize,itype,k,m) = xxmass aqvoldry_sub( isize,itype,k,m) = xxvolu ! ! load coagsolv arrays with number, mass, and volume mixing ratios ! ! *** 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*cair_box do ll = 1, ncomp_plustracer_aer(itype) l = massptr_aer(ll,isize,itype,iphase) duma = 0.0 if (l .ge. p1st) duma = max( 0.0, rsub(l,k,m) ) vol_distrib(isize,ll) = duma sumold(ll) = sumold(ll) + duma end do do llx = 1, 3 ll = (ncomp_coag-3) + llx duma = 0.0 if (llx .eq. 1) then l = hyswptr_aer(isize,itype) if (l .ge. p1st) duma = max( 0.0, rsub(l,k,m) ) else if (llx .eq. 2) then l = waterptr_aer(isize,itype) if (l .ge. p1st) duma = max( 0.0, rsub(l,k,m) ) else duma = max( 0.0, aqvoldry_sub( isize,itype,k,m) ) end if vol_distrib(isize,ll) = duma sumold(ll) = sumold(ll) + duma end do ! calculate dry and wet diameters and wet density if (xxmass .le. 1.01*smallmassbb) then dpdry(isize) = dcen_sect(isize,itype) dpwet(isize) = dpdry(isize) denswet(isize) = xxdens else dpdry(isize) = (xxvolu/(xxnumb*piover6))**onethird dpdry(isize) = max( dpdry(isize), dlo_sect(isize,itype) ) l = waterptr_aer(isize,itype) if (l .ge. p1st) then duma = max( 0.0, rsub(l,k,m) )*mw_water_aer xxmasswet = xxmass + duma xxvoluwet = xxvolu + duma/dens_water_aer dpwet(isize) = (xxvoluwet/(xxnumb*piover6))**onethird dpwet(isize) = min( dpwet(isize), 30.0*dhi_sect(isize,itype) ) denswet(isize) = (xxmasswet/xxvoluwet) else dpwet(isize) = dpdry(isize) denswet(isize) = xxdens end if end if end do ! ! make call to coagulation routine ! call coagsolv( & nsize, maxd_asize, ncomp_coag, ncomp_coag_maxd, & temp_box, patm_box, dtcoag, nsubstep_coag, & lunout_coag, idiag_coag_onoff, iok, & dpdry, dpwet, denswet, num_distrib, vol_distrib ) if (iok .lt. 0) then msg = '*** subr mosaic_coag_1clm -- fatal error in coagsolv' call peg_message( lunout, msg ) call peg_error_fatal( lunout, msg ) else if (iok .gt. 0) then ncountaa(6) = ncountaa(6) + 1 goto 2700 end if ! ! unload coagsolv arrays ! sumnew(:) = 0.0 do isize = 1, nsize do ll = 1, ncomp_coag sumnew(ll) = sumnew(ll) + max( 0.0, vol_distrib(isize,ll) ) end do l = numptr_aer(isize,itype,iphase) rsub(l,k,m) = max( 0.0, num_distrib(isize)/cair_box ) ! unload mass mixing ratios into rsub; 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. p1st) then duma = max( 0.0, vol_distrib(isize,ll) ) rsub(l,k,m) = duma duma = duma*mw_aer(ll,itype) xxmass = xxmass + duma xxvolu = xxvolu + duma/dens_aer(ll,itype) end if end do aqmassdry_sub(isize,itype,k,m) = xxmass xxvolubb(isize) = xxvolu ll = (ncomp_coag-3) + 1 l = hyswptr_aer(isize,itype) if (l .ge. p1st) rsub(l,k,m) = max( 0.0, vol_distrib(isize,ll) ) ll = (ncomp_coag-3) + 2 l = waterptr_aer(isize,itype) if (l .ge. p1st) rsub(l,k,m) = max( 0.0, vol_distrib(isize,ll) ) ll = (ncomp_coag-3) + 3 aqvoldry_sub( isize,itype,k,m) = max( 0.0, vol_distrib(isize,ll) ) end do ! check for mass and volume conservation do ll = 1, ncomp_coag duma = max( sumold(ll), sumnew(ll), 1.0e-35 ) if (abs(sumold(ll)-sumnew(ll)) .gt. 1.0e-6*duma) 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 ! do isize = 1, nsize xxmass = aqmassdry_sub(isize,itype,k,m) xxvolu = aqvoldry_sub( isize,itype,k,m) l = numptr_aer(isize,itype,iphase) xxnumb = rsub(l,k,m) iconform_numb = 1 ! do a cautious calculation of density, using volume from coagsolv if (xxmass .le. smallmassbb) then ncountaa(8) = ncountaa(8) + 1 xxdens = densdefault xxvolu = xxmass/xxdens xxnumb = xxmass/(volumcen_sect(isize,itype)*xxdens) l = hyswptr_aer(isize,itype) if (l .ge. p1st) rsub(l,k,m) = 0.0 l = waterptr_aer(isize,itype) if (l .ge. p1st) rsub(l,k,m) = 0.0 iconform_numb = 0 if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-7a', & ktau, it, jt, k, 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 rsub instead of that from coagsolv ncountaa(7) = ncountaa(7) + 1 if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,4e10.2)') 'coagaa-7b', & ktau, it, jt, k, isize, xxmass, xxvolu, xxdens xxvolu = xxvolubb(isize) xxdens = xxmass/xxvolu if (idiagbb .ge. 200) & write(*,'(a,26x,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)) then ncountaa(9) = ncountaa(9) + 1 duma = xxnumb xxnumb = xxvolu/volumlo_sect(isize,itype) if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-9 ', & ktau, it, jt, k, isize, xxvolu, duma, xxnumb else if (xxnumb .lt. xxvolu/volumhi_sect(isize,itype)) then ncountaa(10) = ncountaa(10) + 1 duma = xxnumb xxnumb = xxvolu/volumhi_sect(isize,itype) if (idiagbb .ge. 200) & write(*,'(a,i10,4i4,1p,3e12.4)') 'coagcc-10', & ktau, it, jt, k, isize, xxvolu, duma, xxnumb end if end if ! load number, mass, volume, dry-density back into arrays l = numptr_aer(isize,itype,iphase) rsub(l,k,m) = xxnumb adrydens_sub( isize,itype,k,m) = xxdens aqmassdry_sub(isize,itype,k,m) = xxmass aqvoldry_sub( isize,itype,k,m) = xxvolu end do 2700 continue ! temporary diagnostics if ((idiagbb .ge. 100) .and. & (it .eq. ite) .and. & (jt .eq. jte) .and. (k .eq. kclm_calcend)) 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 ! k levels 2900 continue ! subareas return end subroutine mosaic_coag_1clm !---------------------------------------------------------------------- subroutine coagsolv( & nbin, nbin_maxd, ncomp, ncomp_maxd, & tkelvin_inp, patm_inp, deltat_inp, nsubstep, & lunout, idiag_onoff, iok, & dpdry_inp, dpwet_inp, denswet_inp, & num_distrib_inp, vol_distrib_inp ) implicit none ! ! ********************************************************************* ! ***************** written by mark jacobson (1993) ******************* ! *** (c) copyright, 1993 by mark z. jacobson *** ! *** this version modified december, 2001 *** ! *** (650) 723-6836 *** ! ********************************************************************* ! ! cccccc ooooo a gggggg ssssss ooooo l v v ! c o o a a g s o o l v v ! c o o a a g ggg sssss o o l v v ! c o o aaaaaaa g g s o o l v v ! cccccc ooooo a a gggggg ssssss ooooo lllllll v ! ! ********************************************************************* ! the use of this module implies that the user agrees not to sell the ! module or any portion of it, not to remove the copyright notice from it, ! and not to change the name of the module, "coagsolv". users may ! modify the module as needed or incorporate it into a larger model. ! any special requests with regard to this module may be directed to ! jacobson@stanford.edu. ! ********************************************************************* ! ! ********************************************************************* ! * this is a box-model version of "coagsolv," a semi-implicit * ! * aerosol coagulation solver. the solver is exactly volume * ! * conserving, unconditionally stable (regardless of time step), * ! * and noniterative. * ! * * ! * this program calculates brownian coagulation kernels and solves * ! * self-coagulation among any number of size bins, one * ! * particle type and any number of volume fractions. * ! * * ! * the program is set up as a box-model. * ! * * ! * the volume ratio of adjacent size bins can be any number > 1 * ! * * ! * the scheme can be used to solve coagulation with any size bin * ! * structure. * ! * * ! * the initial size distribution here is monomer or lognormal * ! * (ifsmoluc = 1 or 0) with a fixed size bin structure * ! * * ! ********************************************************************* ! * references * ! ********************************************************************* ! * semi-implicit scheme using movable or variable bins and with * ! * any volume ratio (vrat) > 1 * ! * * ! * jacobson m. z., turco r. p., jensen, e. j. and toon o. b. (1994) * ! * modeling coagulation among particles of different compostion * ! * and size. atmos. environ. 28a, 1327 - 1338. * ! * * ! * jacobson m. z. (1999) fundamentals of atmospheric modeling. * ! * cambridge university press, new york, 656 pp. * ! * * ! ********************************************************************* ! * semi-implicit scheme using fixed bins with volume ratio >,= 2 * ! * * ! * toon o. b., turco r. p., westphal d., malone r. and liu m. s. * ! * (1988) a multidimensional model for aerosols: description of * ! * computational analogs. j. atmos. sci. 45, 2123 - 2143 * ! * * ! ********************************************************************* ! * orig semi-implicit scheme using fixed bins with volume ratio = 2 * ! * * ! * turco r. p., hamill p., toon o. b., whitten r. c. and kiang c. s. * ! * (1979) the nasa-ames research center stratospheric aerosol * ! * model: i physical processes and computational analogs. nasa * ! * tech. publ. (tp) 1362, iii-94. * ! ********************************************************************* ! ! modified by y. zhang for incorporation into pnnl's mirage ! june-july, 2003 ! ! ********************************************************************* ! ! modified by r.c. easter for incorporation into pnnl's wrf-chem ! feb 2005 (a) ! added "_inp" to all subr arguments ! added iradmaxd_inp & lunout ! define iradmaxd, iaertymaxd, iaeromaxd locally ! no include statements ! changed real*16 to real*8 ! feb 2005 (b) ! in coagsolv, distrib_inp is 2-d array that holds both number and ! volume concentrations; distrib is initialized from it; ! the "fracnaer" stuff is gone ! feb 2005 (c) ! change distrib & cc2 to be 2-d arrays (which simplifies indexing!); ! iaeromaxd and iaero are gone; ! deleted many commented-out lines in coagsolv ! feb 2005 (d) ! changed cc2(i,1) to be number instead of all-species volume; ! prod term for number distrib now follows jgr-2002 article ! [multiply by volume(j), divide by volume(k)] ! apr 2006 (a) ! considerable cleanup (mostly cosmetic) ! pass in the bin sizes directly, rather than vrat & dmin_um ! pass in dry and wet sizes separately ! pass in/out number and volume distributions in separate arrays ! ! ********************************************************************* ! IN arguments integer nbin ! actual number of size bins integer nbin_maxd ! size-bin dimension for input (argument) arrays integer ncomp ! actual number of aerosol volume components integer ncomp_maxd ! volume-component dimension for input (argument) arrays integer lunout ! logical unit for warning & diagnostic output integer idiag_onoff ! if positive, write some mass-conservation diagnostics integer nsubstep ! number of time sub-steps for the integration real tkelvin_inp ! air temperature (K) real patm_inp ! air pressure (atm) real deltat_inp ! integration time (s) real dpdry_inp(nbin_maxd) ! bin dry diameter (cm) real dpwet_inp(nbin_maxd) ! bin wet (ambient) diameter (cm) real denswet_inp(nbin_maxd) ! bin wet (ambient) density (g/cm3) ! IN-OUT arguments real num_distrib_inp(nbin_maxd) ! num_distrib_inp(i) = number concentration for bin i (#/cm3) real vol_distrib_inp(nbin_maxd,ncomp_maxd) ! vol_distrib_inp(i,j) = volume concentration for bin i, ! component j (cm3/cm3) ! OUT arguments integer iok ! status flag (0=success, anything else=failure) ! ! 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 coagulation rate coefficients. ! ! NOTE -- Units for num_distrib and vol_distrib ! num_distrib units MUST BE (#/cm3) ! vol_distrib units do not really matter. They can be either masses ! or volumes, and either mixing ratios or concentrations, ! (e.g., g/cm3-air, ug/kg-air, cm3/m3-air, cm3/kg-air, ...). ! Use whatever is convenient. ! ! ! local variables ! integer iradmaxd_wrk, iaertymaxd_wrk parameter (iradmaxd_wrk=16) parameter (iaertymaxd_wrk=150) integer irad integer iaerty ! irad == nbin = actual number of size bins ! iradmaxd_wrk = size-bin dimension for local (working) arrays ! ! (iaerty-1) == ncomp = actual number of aerosol volume components ! iaertymaxd_wrk = volume-component dimension for local (working) arrays ! ! iaerty = 1 --> calc number concentration only ! > 1 --> calc number concentration + (iaerty-1) volume concentrations integer i, isubstep, j, jaer, jb, k, kout integer jbinik(iradmaxd_wrk,iradmaxd_wrk,iradmaxd_wrk) integer jbins(iradmaxd_wrk,iradmaxd_wrk) !wig, 11-Dec-2006: real*8 causing core dumps on IBM's (i.e. bluesky) ! real*8 aknud, aloss, amu, avg, & real aknud, aloss, amu, avg, & boltg, bpm, & cbr, consmu, cpi, & deltat, deltp1, deltr, & divis, dti1, dti2, & finkhi, finklow, fourpi, fourrsq, & ggr, gmfp, & onepi, & patm, pmfp, prod, & radi, radiust, radj, ratior, & rgas2, rho3, rsuma, rsumsq, & sumdc, sumnaft, sumnbef, & term1, term2, third, tk, tkelvin, tworad, & viscosk, vk1, voltota, vthermg, & wtair !wig, 11-Dec-2006: real*8 causing core dumps on IBM's (i.e. bluesky) ! real*8 cc2(iradmaxd_wrk,iaertymaxd_wrk), & real cc2(iradmaxd_wrk,iaertymaxd_wrk), & conc(iradmaxd_wrk), & denav(iradmaxd_wrk) , & difcof(iradmaxd_wrk), & distrib(iradmaxd_wrk,iaertymaxd_wrk), & fink(iradmaxd_wrk,iradmaxd_wrk,iradmaxd_wrk), & floss(iradmaxd_wrk,iradmaxd_wrk), & radius(iradmaxd_wrk), & radwet(iradmaxd_wrk), & rrate(iradmaxd_wrk,iradmaxd_wrk), & sumdp(iradmaxd_wrk), & sumvt(iradmaxd_wrk), & tvolfin(iaertymaxd_wrk), & tvolinit(iaertymaxd_wrk), & volume(iradmaxd_wrk), & volwet(iradmaxd_wrk), & vthermp(iradmaxd_wrk) ! ********************************************************************* ! set some parameters ! ********************************************************************* irad = nbin iaerty = ncomp + 1 kout = lunout if (irad .gt. iradmaxd_wrk) then write(lunout,*) '*** coagsolv: irad > iradmaxd_wrk: ', & irad, iradmaxd_wrk iok = -1 return endif if (iaerty .gt. iaertymaxd_wrk) then write(lunout,*) '*** coagsolv: iaerty > iaertymaxd_wrk: ', & iaerty, iaertymaxd_wrk iok = -2 return endif ! ! tkelvin = temperature (k) ! denav = particle density (g cm-3) ! patm = air pressure (atm) ! tkelvin = tkelvin_inp patm = patm_inp do i = 1, irad denav(i) = denswet_inp(i) end do ! ! deltat_inp = total integration time (s) ! deltat = individual time step (s) ! nsubstep = number of time steps ! deltat = deltat_inp/nsubstep ! ! boltg = boltzmann's 1.38054e-16 (erg k-1 = g cm**2 sec-1 k-1) ! wtair = molecular weight of air (g mol-1) ! avg = avogadro's number (molec. mol-1) ! rgas2 = gas constant (l-atm mol-1 k-1) ! amu = dynamic viscosity of air (g cm-1 sec-1) ! est value at 20 c = 1.815e-04 ! rho3 = air density (g cm-3) ! viscosk = kinematic viscosity = amu / denair = (cm**2 sec-1) ! vthermg = mean thermal velocity of air molecules (cm sec-1) ! gmfp = mean free path of an air molecule = 2 x viscosk / ! thermal velocity of air molecules (gmfp units of cm) ! tk = tkelvin third = 1. / 3. onepi = 3.14159265358979 fourpi = 4. * onepi boltg = 1.38054e-16 wtair = 28.966 avg = 6.02252e+23 rgas2 = 0.08206 consmu = 1.8325e-04 * (296.16 + 120.0) amu = (consmu / (tk + 120.)) * (tk / 296.16)**1.5 rho3 = patm * wtair * 0.001 / (rgas2 * tk) viscosk = amu / rho3 vthermg = sqrt(8. * boltg * tk * avg / (onepi * wtair)) gmfp = 2.0 * viscosk / vthermg ! ! ********************************************************************* ! * set volume ratio size bin grid * ! ********************************************************************* ! cpi = fourpi / 3. do 30 j = 1, irad radius( j) = dpdry_inp(j) * 0.5 volume( j) = cpi * radius(j) * radius(j) * radius(j) radwet( j) = dpwet_inp(j) * 0.5 volwet( j) = cpi * radwet(j) * radwet(j) * radwet(j) 30 continue ! ! ********************************************************************* ! * determine bins where coagulated terms go * ! ********************************************************************* ! finklow = volume fraction of i+j going to lower (k ) bin ! finkhi = volume fraction of i+j going to higher (k+1) bin ! fink(i,j,k) = volume fraction of i+j going to bin k ! floss = simplified fink term used in loss rates ! no self-coagulation loss out of largest bin ! jbins(i,k) = number of production terms into bin k from bin i ! jbinik(i,k,jb) = identifies each bin j that coagulates with bin i ! to produce bin k ! do 35 i = 1, irad do 34 j = 1, irad jbins(i,j) = 0 floss(i,j) = 0. do 33 k = 1, irad jbinik(i,j,k) = 0 fink( i,j,k) = 0. 33 continue 34 continue 35 continue do 40 i = 1, irad do 39 j = 1, irad voltota = volume(i) + volume(j) if (voltota.ge.volume(irad)) then if (i.eq.irad) then floss(i,j) = 0.0 else floss(i,j) = 1.0 endif if (j.lt.irad) then jbins(i,irad) = jbins(i,irad) + 1 jb = jbins(i,irad) jbinik(i,irad,jb) = j fink( i,jb,irad) = 1.0 endif else do 45 k = 1, irad - 1 vk1 = volume(k+1) if (voltota.ge.volume(k).and.voltota.lt.vk1) then finklow = ((vk1 - voltota)/(vk1 - volume(k))) & * volume(k) / voltota finkhi = 1. - finklow if (i.lt.k) then floss(i,j) = 1. elseif (i.eq.k) then floss(i,j) = finkhi endif if (j.lt.k) then jbins(i,k) = jbins(i,k) + 1 jb = jbins(i,k) jbinik(i,k,jb) = j fink( i,jb,k) = finklow endif jbins(i,k+1) = jbins(i,k+1) + 1 jb = jbins(i,k+1) jbinik(i,k+1,jb) = j fink( i,jb,k+1) = finkhi endif 45 continue endif do 38 k = 1, irad if (jbins(i,k).gt.irad) then write(21,*)'coagsolv: jbins > irad: ',jbins(i,k),i,k iok = -3 return endif 38 continue 39 continue 40 continue ! ! *********************************************************************** ! initialize size distribution ! ! copy initial number concentration (#/cm3-air) and volume concentrations ! (cm3/cm3-air) from num/vol_distrib_inp (real*4) to distrib (real*8) ! *********************************************************************** do i = 1, irad distrib(i,1) = num_distrib_inp(i) end do do jaer = 2, iaerty do i = 1, irad distrib(i,jaer) = vol_distrib_inp(i,jaer-1) end do end do ! ! ********************************************************************* ! coagulation kernel from fuchs equations ! ********************************************************************* ! difcof = brownian particle diffus coef (cm**2 sec-1) ! = boltg * t * bpm / (6 * pi * mu * r(i)) ! = 5.25e-04 (diam = 0.01um); 6.23e-7 (diam = 0.5 um) seinfeld p.325. ! vthermp = avg thermal vel of particle (cm sec-1) ! = (8 * boltg * t / (pi * masmolec)**0.5 ! pmfp = mean free path of particle (cm) ! = 8 * di / (pi * vthermp) ! bpm = correction for particle shape and for effects of low mean ! free path of air. (toon et al., 1989, physical processes in ! polar stratospheric ice clouds, jgr 94, 11,359. f1 and f2 are ! included in the expression below (shape effects correction) ! = 1 + kn[a + bexp(-c/kn)] ! deltp1 = if particles with mean free path pmfp leave the surface of ! an absorbing sphere in all directions, each being probably ! equal, then deltp1 is the mean distance from the surface of the ! sphere reached by particles after covering a distance pmfp. (cm). ! cbr = coag kernel (cm3 partic-1 s-1) due to brownian motion (fuchs, 1964) ! rrate = coag kernal (cm3 partic-1 s-1) * time step (s) ! ! *** use the wet (ambient) radius and volume here *** do 145 i = 1, irad radiust = radwet(i) tworad = radiust + radiust fourrsq = 4. * radiust * radiust aknud = gmfp / radiust bpm = 1. + aknud * (1.257 + 0.42*exp(-1.1/aknud)) difcof(i) = boltg * tk * bpm / (6. * onepi * radiust * amu) ! use bin-varied aerosol density from mirage ! vthermp(i) = sqrt(8. * boltg * tk / (onepi*denav * volume(i))) vthermp(i) = sqrt(8. * boltg * tk / (onepi*denav(i) * & volwet(i))) sumvt(i) = vthermp(i) * vthermp(i) pmfp = 8. * difcof(i) / (onepi * vthermp(i)) dti1 = tworad + pmfp dti2 = (fourrsq + pmfp * pmfp)**1.5 divis = 0.166667 / (radiust * pmfp) deltp1 = divis * (dti1 * dti1 * dti1 - dti2) - tworad sumdp(i) = deltp1 * deltp1 145 continue !yy write(kout,9165) do 155 i = 1, irad do 154 j = 1, irad radi = radwet(i) radj = radwet(j) rsuma = radi + radj rsumsq = rsuma * rsuma ratior = radi / radj sumdc = difcof(i) + difcof(j) ggr = sqrt(sumvt(i) + sumvt(j)) deltr = sqrt(sumdp(i) + sumdp(j)) term1 = rsuma / (rsuma + deltr) term2 = 4. * sumdc / (rsuma * ggr) cbr = fourpi * rsuma * sumdc / (term1 + term2) rrate(i,j) = cbr * deltat !yy write(kout,9190) (2.0e4*radius(i)), (2.0e4*radius(j)), cbr 154 continue 155 continue 9165 format(16x,'coagulation coefficients (cm**3 #-1 sec-1)'/ & 'diam1_um diam2_um brownian ') 9190 format(1pe15.7,7(1x,1pe15.7)) ! ! ********************************************************************* ! ****************** solve coagulation equations ********************** ! ********************************************************************* ! ! ********************************************************************* ! inialize new arrays ! ********************************************************************* ! distrib = initial conc (# cm-3 for num conc.; cm3 cm-3 for volume fractions) ! cc2 = initial, final volume conc (cm3 cm-3) for all particle types. ! conc = initial, final number conc (# cm-3) for number distribution ! volume = volume (cm3 #-1) of particles in bin ! ! do i = 1, irad ! cc2( i,1) = distrib(i,1) * volume(i) cc2( i,1) = distrib(i,1) conc(i) = distrib(i,1) end do do jaer = 2, iaerty do i = 1, irad cc2(i,jaer) = distrib(i,jaer) end do end do ! ********************************************************************* ! solve coagulation over several time steps ! ********************************************************************* ! do 700 isubstep = 1, nsubstep ! do 485 jaer = 1, iaerty do 484 k = 1, irad ! ! production terms ! prod = 0. if (k.gt.1) then if (jaer .eq. 1) then do 465 i = 1, k do 464 jb = 1, jbins(i,k) j = jbinik(i,k,jb) prod = prod + fink(i,jb,k)*rrate(i,j)* & cc2(j,jaer)*volume(j)*conc(i) 464 continue 465 continue prod = prod/volume(k) else do 469 i = 1, k do 468 jb = 1, jbins(i,k) j = jbinik(i,k,jb) prod = prod + & fink(i,jb,k)*rrate(i,j)*cc2(j,jaer)*conc(i) 468 continue 469 continue endif endif ! ! loss terms ! aloss = 0. if (k.lt.irad) then do 475 j = 1, irad aloss = aloss + floss(k,j) * rrate(k,j) * conc(j) 475 continue endif ! ! final concentrations ! cc2(k,jaer) = (cc2(k,jaer) + prod) / (1. + aloss) 484 continue 485 continue ! ! put updated number concentration into conc array do i = 1, irad conc(i) = cc2(i,1) end do 700 continue ! ! ********************************************************* ! ** put the advanced concentration back on the grid ** ! (copy them from the real*8 working array to the ! real*4 caller arrays) ! ********************************************************* ! do i = 1, irad num_distrib_inp(i) = cc2(i,1) end do do jaer = 2, iaerty do i = 1, irad vol_distrib_inp(i,jaer-1) = cc2(i,jaer) end do end do ! ! if no diagnostics, then return ! iok = 0 if (idiag_onoff .le. 0) return ! ! ********************************************************************* ! ************* check whether mass is conserved *********** ! ********************************************************************* ! tvolinit = initial volume conc (cm3 cm-3), summed over all bins ! tvolfin = final volume conc (cm3 cm-3), summed over all bins ! sumnbef = summed number conc (# cm-3) before coagulation ! sumnaft = summed number conc (# cm-3) after coagulation ! do jaer = 1, iaerty tvolinit(jaer) = 0. tvolfin( jaer) = 0. end do sumnbef = 0. sumnaft = 0. ! ! distrib(jaer=1, i=1:nrad) = initial total number conc in # cm-3 ! cc2 (jaer=1, i=1:nrad) = final total number conc in # cm-3 ! do i = 1, irad tvolinit(1) = tvolinit(1) + distrib(i,1) * volume(i) tvolfin( 1) = tvolfin( 1) + cc2( i,1) * volume(i) sumnbef = sumnbef + distrib(i,1) sumnaft = sumnaft + cc2( i,1) end do ! ! distrib(jaer=2:iaerty, i=1:nrad) = initial component volume conc in cm3 cm-3 ! cc2 (jaer=2:iaerty, i=1:nrad) = initial component volume conc in cm3 cm-3 ! do jaer = 2, iaerty do i = 1, irad tvolinit(jaer) = tvolinit(jaer) + distrib(i,jaer) tvolfin( jaer) = tvolfin( jaer) + cc2( i,jaer) end do end do write(kout,*) write(kout,9434) sumnbef, sumnaft, & tvolinit(1)*1.0e+12,tvolfin(1)*1.0e+12 write(kout,9435) sumnaft-sumnbef write(kout,9439) tvolinit(1) / tvolfin(1) do jaer = 2, iaerty if (abs(tvolfin(jaer)) .gt. 0.0) then write(kout,9440) jaer-1, tvolinit(jaer) / tvolfin(jaer) else write(kout,9441) jaer-1, tvolinit(jaer), tvolfin(jaer) end if end do 9434 format('# bef, aft; vol bef, aft =',/4(1pe16.10,1x)/) 9435 format('total partic change in # cm-3 = ',1pe17.10) 9439 format('total partic vol bef / vol aft = ',1pe17.10, & ': if 1.0 -> exact vol conserv') 9440 format('vol comp',i4,' vol bef / vol aft = ',1pe17.10, & ': if 1.0 -> exact vol conserv') 9441 format('vol comp',i4,' vol bef, vol aft = ',2(1pe17.10)) ! ! ************************** print results **************************** ! distrib = initial conc in # cm-3 for total particle (jaer=1, i=1,nrad) ! conc = final conc in # cm-3 for total particle (jaer=1, i=1,nrad) ! ! ! ********************************************************************* ! end of program coagsolv.f ! ********************************************************************* ! return end subroutine coagsolv !----------------------------------------------------------------------- end module module_mosaic_coag