!**********************************************************************************  
! 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