c****************************************************************************
c THIS IS THE ROUTINE THAT CALCULATES THE CHANGES TO THE AEROSOL
c MASS DISTRIBUTION BECAUSE OF COAGULATION FOR Dt.
c
c A HIGH RESOLUTION DISTRIBUTION IS USED FOR THE COAGULATION CALCULATIONS
c BY DIVIDING EACH SECTION OF THE ORIGINAL DISTRIBUTION TO TEN SECTIONS
c THE NEW DISTRIBUTION IS TREATED AS A POPULATION OF MONODISPERSE
c AEROSOL GROUPS. FOR 10 SECTIONS, THE RESULTING CALCULATION USES
c 100 AEROSOL GROUP SOMETHING THAT HAS PROVEN TO BE ACCURATE FOR
c SUCH A CALCULATION
c
c****************************************************************************
c
c PARAMETERS
c   Q(I) : VECTOR OF AEROSOL AND GAS-PHASE CONCENTRATIONS FOR THE
c          COMPUTATIONAL CELL (UNITS: ug/m3 for aerosol)
c   DIAMETROS(i): MEAN DIAMETER OF HIGH RESOLUTION SECTION (um)
c   KCG(I,J): COAGULATION COEFFICIENT BETWEEN PARTICLES OF SECTIONS I AND J
c             THE UNITS ARE cm3/hr
c   CGNEWP(I,J): THE SECTION TO WHICH THE AGGREGATE PARTICLE RESULTING
c                FROM THE COLLISIONS OF PARTICLES FROM SECTION I AND J
c                BELONGS
c   PN(I): PARTICLE NUMBER CONCENTRATION HIGH RESOLUTION (UNITS: cm-3)
c   M(I) : MASS OF ONE PARTICLE IN HIGH RESOLUTION SECTION I (UNITS: ug)
c   DENS : PARTICLE DENSITY (UNITS: ug/m3)
c   DT   : OPERATOR TIMESTEP IN SECONDS
c



      subroutine coagul(q)
c
      include 'dynamic.inc'
      parameter (ndcoag = nsec*nres)
      parameter (dens = 1.0d12)          ! in ug/m3
c
      real*8 q(ntotal), pmtotal(nsec), pm(ndcoag)
      real*8 diametros(ndcoag), pn(ndcoag), m(ndcoag)
      real*8 kcg(ndcoag, ndcoag)
      real*8 dmdt(ndcoag, nsp), dqdt(ntotal)
      real*8 basi, sint, fraction1, fraction2, change
      integer index, index1, index2, isec1, isec2
      integer isec, jres, i ,j, k
      integer cgnewp(ndcoag,ndcoag)
c
c     RETURN IF NO COAGULATION IS REQUIRED
c
      if (icoag .eq. 0) return
c
c CALCULATION OF THE DIAMETERS OF THE HIGH RESOLUTION DISTRIBUTION
c
      basi = dlog10(dsecf(1))
      sint = (dlog10(dsecf(2)) - dlog10(dsecf(1)))/FLOAT(nres)
      basi = basi + sint/2.0d0            ! mid-point of section
      do isec=1, nsec
       do jres=1, nres
        diametros((isec-1)*nres+jres) = 10.0d0**basi
        basi = basi + sint
       enddo
      enddo
c
c CALCULATION OF THE COAGULATION COEFFICIENT AND COLLISION ARRAYS
c
      call kcoag(diametros, kcg, cgnewp)
c
c CALCULATION OF THE NUMBER DISTRIBUTION
c
c* 1. Zeroing of the total mass of each section
c
      do isec = 1, nsec
       pmtotal(isec) = 0.0d0
      enddo
c
c* 2. Calculation of the total mass of each section
c
      do isec = 1, nsec
       do j =1 , nsp
        index = (isec-1)*nsp + j
        pmtotal(isec) = pmtotal(isec) + q(index)
       enddo
      enddo
c
c* 3. Calculation of total mass of each high resolution section
c
      do isec = 1, nsec
       do j = 1, nres
        index = (isec-1)*nres+j
        pm(index) = pmtotal(isec)/FLOAT(nres)
       enddo
      enddo
c
c* 4. Calculation of the number and single particle mass of each high resolution section
c
      do i=1, ndcoag
       m(i) = pi*dens*diametros(i)**3/6.0d18         ! in ug
       pn(i)= 1.0d-6*pm(i)/m(i)                       ! in part/cm3
      enddo
c
c ZEROING OF THE COAGULATION RATE VECTORS
c
      do i=1, ndcoag
       do k=1, nsp
        dmdt(i, k) = 0.0d0
       enddo
      enddo 
c
c CALCULATION OF THE COAGULATION RATES
c
      do i=1, ndcoag
       do j=1, i
         cgrate = kcg(i,j)*pn(i)*pn(j)        ! UNITS cm-3 hr-1
         if ( i .eq. j) cgrate=cgrate/2.0d0   ! to avoid double-counting of i,i collisions
c
       do k=1, nsp
       index1 = ((i-1)/nres)*nsp+ k              ! species position in Q vector
       index2 = ((j-1)/nres)*nsp+ k
       isec1 = ((i-1)/nres) + 1                  ! section in low resolution
       isec2 = ((j-1)/nres) + 1
c
       if (pmtotal(isec1) .le. 0.0d0) then
        fraction1 = 0.0d0
        else
        fraction1 = q(index1) / pmtotal(isec1)    ! fraction of section mass
       endif
c
       if (pmtotal(isec2) .le. 0.0d0) then
        fraction2=0.0d0
        else
        fraction2 = q(index2) / pmtotal(isec2)
       endif
c
       dmdt(i,k)=dmdt(i,k)- cgrate*m(i)*1.d6*fraction1   ! UNITS ug/m3 hr
       dmdt(j,k)=dmdt(j,k)- cgrate*m(j)*1.d6*fraction2   ! UNITS ug/m3 hr
       dmdt(cgnewp(i,j), k)=dmdt(cgnewp(i,j), k) +
     &   cgrate*(m(i)*fraction1+m(j)*fraction2)*1.d6
       enddo
       enddo
      enddo
c
c MASS BALANCE CHECK 1: ** TO BE REMOVED **
c
csp       change = 0.0
csp       do i = 1, ndcoag
csp       do j=1, nsp
csp       change = change + dmdt(i,j)
csp       enddo
csp       enddo
csp       write(3,*)tcom, change
c
c     CALCULATION OF THE COAGULATION RATES FOR THE DIFFERENT AEROSOL SPECIES
c     TRANSITION FROM THE HIGH RESOLUTION TO THE LOW RESOLUTION DISTRIBUTION
c
       do i=1, ntotal
       dqdt(i)=0.0d0
       enddo
c
       do i = 1, ndcoag
       do k = 1, nsp
       index = ((i-1)/nres)*nsp + k
       dqdt(index) = dqdt(index) + dmdt(i,k)        ! in ug m3 hr
       enddo
       enddo
c
c  MASS BALANCE CHECK 2: ** TO BE REMOVED **
c
csp       change = 0.0
csp       do i=1, ntotal
csp       change = change + dqdt(i)
csp       enddo
csp       write(4,*)tcom, change

c
c UPDATING OF THE AEROSOL SIZE/COMPOSITION DISTRIBUTION AFTER COAGULATION
c
c NOTE: IT IS POSSIBLE THAT A GIVEN SECTION MAY DISSAPPEAR COMPLETELY
c       BECAUSE OF OUR ASSUMPTION OF A CONSTANT COAGULATION RATE FOR Dt.
c       THIS CAN RESULT IN A SMALL GAIN OF MASS AND NEGATIVE CONCENTRATIONS
c
c       TO CORRECT THIS THE CONCENTRATION OF THE SPECIES IS SET TO ZERO
c       AND TO CONSERVE MASS THE CONCENTRATION OF THE SAME SPECIES
c       IN THE LARGER SECTION IS REDUCED BY THE SAME AMOUNT
c
       do i = 1, ntotal
       q(i) = q(i) + dqdt(i)*(dt/3600.0d0)
c
        if (q(i) .lt. 0.0d0) then
csp        write(6,*) 'COAGULATION RESULTED IN NEGATIVE CONCENTRATION'
csp        write(6,*) 'MASS GAIN =', -q(i), 'ug/m3'
cgy the syntax stop 'string' doesn't work well on all platforms
cgy        if (i+nsp.gt.naer) stop'COAGUL: beyond aerosol section' ! bkoo (6/8/00)
        if (i+nsp.gt.naer) then
          write(6,*)'COAGUL: beyond aerosol section'
          stop
        endif
        q(i+nsp)=q(i+nsp)+q(i)
        q(i) =0.0d0
        endif
       enddo
c
c END OF CALCULATION
c
      return
      end




      subroutine kcoag(dp, kcg,cgnewp)
c
c     This routine calculates the coagulation coefficient (K) for Brownian
c     diffusion using the formulas of Fuchs (Seinfeld, 1986, Table 10.1) 
c     and puts them in the KCG array
c     It also determines in what section the new particle will exist and 
c     puts the section number in the CGNEWP array
c
c
      include 'dynamic.inc'
c
      integer ndcoag
      real*8 rho, airl, visc
c
      parameter (ndcoag = nsec*nres)
      parameter (rho = 1.0d0)      ! density of particles in g/cm3
      parameter (airl = 0.0651d0)  ! air mean free path in um
      parameter (visc= 1.83d-4)    ! air viscosity in poise
c
      real*8 kcg(ndcoag,ndcoag),difc(ndcoag),cc(ndcoag),cbar
      real*8 gbar,g(ndcoag),l(ndcoag),mass(ndcoag),knudsen(ndcoag)
      real*8 difsum,dpsum,vol(ndcoag),newvol,volbrk
      real*8 dp(ndcoag)
c
      integer cgnewp(ndcoag,ndcoag)
      integer i,j,k
c     
      do i=1,ndcoag      
       vol(i)=dp(i)**3*pi/6.d0                       ! in um3
       mass(i)=vol(i)*rho*1.0d-12                   ! in grams
       cc(i)=DSQRT(8.d0*temp*1.381d-20/pi/mass(i))   ! in m/s
       knudsen(i)=2.0d0*airl/dp(i)
       difc(i)=1.381d-20*temp/3.d0/pi/visc/dp(i)*1.0d4*
     & ((5.d0+4.d0*knudsen(i)+6.d0*knudsen(i)**2+18.d0*knudsen(i)**3)/
     & (5.0d0-knudsen(i)+(8.0d0+pi)*knudsen(i)**2))    !in m2/s
       l(i)=8.0d0*difc(i)/pi/cc(i)*1.d6               !in um
       g(i)=1.0d0/(3.0d0*dp(i)*l(i))*((dp(i)+l(i))**3-
     &      DSQRT((dp(i)**2+l(i)**2)**3))-dp(i)        !in um
      enddo
c
       do i=1,ndcoag
        do j=1,ndcoag
         cbar=DSQRT(cc(i)**2+cc(j)**2)             ! in m/s
         gbar=DSQRT(g(i)**2+g(j)**2)             ! in um
         difsum = difc(i)+difc(j)            ! in m2/s
         dpsum = dp(i)+dp(j)                     ! in um
         kcg(i,j)=2.0d0*pi*difsum*dpsum*3600.0d0/
     &   (dpsum/(dpsum+2.0d0*gbar)+8.0d0*difsum/cbar/dpsum*1.d6) ! in cm3/hr
c
         newvol=vol(i)+vol(j)
         k=max(i,j)
         cgnewp(i,j)=0
              do while (cgnewp(i,j) .eq. 0)
                if (k .lt. ndcoag) then
                   volbrk = 0.5d0*(vol(k)+vol(k+1))
                else
                   volbrk=2.0d0*vol(ndcoag)+1.0d0
                endif
c
                if(newvol .lt. volbrk) then
                   cgnewp(i,j)=k
                else
                   k = k+1
                endif
              enddo
c
       enddo
      enddo
c
      return
      end
