module dust_coagulation_mod

! ------------------------
! Dust coagulation scheme
! Based on Bertrand et al., 2022
! Author: T. Bertrand
! ------------------------

use comcstfi_h, only: pi,g,r,mugaz
use tracer_mod, only: igcm_dust_mass,igcm_dust_number, &
                      igcm_stormdust_mass,igcm_stormdust_number, &
                      igcm_topdust_mass,igcm_topdust_number
use microphys_h, only: kbz,m0co2
implicit none

public ::  dust_coagulation_init

logical,save :: full_coag_equations ! coagulation of dust particles using full equations instead of lookup tables
logical,save :: dust_coag_kernel_b ! flag to activate Brownian coagulation 
logical,save :: dust_coag_kernel_g !                  gravitational
logical,save :: dust_coag_kernel_de !                 diffusion enhancement
logical,save :: dust_coag_kernel_ti ! Not implemented yet. 
!$OMP THREADPRIVATE(full_coag_equations,dust_coag_kernel_b,dust_coag_kernel_g,dust_coag_kernel_de,dust_coag_kernel_ti,coal_coeff_mode) 
integer,save :: coal_coeff_mode ! Flag to choose how to calculate coalescence coefficient
!$OMP THREADPRIVATE(coal_coeff_mode) 

! -------------------------------
! Coagulation Input Parameters
! -------------------------------
real    ::  effb=1.  !  Sticking efficiency for B coagulation
real    ::  coal_fac=1.  ! Overall coalescence factor
! Bin discretization of particle sizes
real*8, save :: r1_coag = 0.01e-6 ! min radius 
real*8, save :: rn_coag = 40.e-6  ! max radius
integer,parameter :: nres_coag = 10  ! nb bins
real*8, save :: vrat_coag   ! volume ratio
real*8, save, allocatable :: rads_coag(:),vols_coag(:),deltar_coag(:)

!$OMP THREADPRIVATE(r1_coag,rn_coag,vrat_coag)
!$OMP THREADPRIVATE(rads_coag,vols_coag,deltar_coag)

! parameters of the look up tables of coagulation coefficients
integer,parameter :: table_numt = 15   ! Nb of temperature bins
integer,parameter :: table_nump = 25   ! Nb of pressure bins
integer,parameter :: table_numm = 25   ! Nb of mass bins
real, save :: table_pt(table_numt) 
real, save :: table_pres(table_nump) 
real, save :: table_mfp(table_numm) 
real, save :: table_b(table_numt,table_nump,nres_coag,nres_coag) ! Brownian
real, save :: table_g(table_numt,table_nump,nres_coag,nres_coag) ! Gravitational
real, save :: table_de(table_numt,table_nump,nres_coag,nres_coag) ! Brownian diffusion enhancement
real, save :: table_ti(table_numt,table_nump,nres_coag,nres_coag) ! Turbulent
real*8, public :: dev_dt = 0.63676    ! standard deviation of dust distribution
real, public ::rho_dust=2500.  ! Mars dust density (kg.m-3)

!$OMP THREADPRIVATE(table_pt,table_pres,table_mfp)
!$OMP THREADPRIVATE(table_b,table_g,table_de,table_ti)

contains

!#######################################################################

subroutine dust_coagulation_init()

! Initialisation of the particle bin discretization

integer  :: i 
allocate(rads_coag(nres_coag))
allocate(vols_coag(nres_coag))
allocate(deltar_coag(nres_coag))

vrat_coag=(rn_coag/r1_coag)**(3./(nres_coag-1))

do i = 1, nres_coag
        rads_coag(i)  = r1_coag*vrat_coag**((i-1)/3.)
        vols_coag(i)  = 4./3.*pi*r1_coag**3*vrat_coag**(i-1)
enddo
! diameter width
deltar_coag(:)=2.*rads_coag(:)*2.**(1/3.)*(vrat_coag**(1./3.)-1)/(1+vrat_coag)**(1./3.)

! Making Lookup Tables
if (.not.full_coag_equations) then
 call make_tables()
endif
!deallocate(rads_coag)
!deallocate(vols_coag)
!deallocate(deltar_coag)

end subroutine dust_coagulation_init

!#######################################################################

SUBROUTINE dust_coagulation_main( &
   ngrid, nlayer, nq, ptime, ptimestep, pq, pdqfi, pt, pdtfi, &
   pplay, pplev, pdqcoag)

!--------------------------------------------------------
! Input Variables
!--------------------------------------------------------
INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
INTEGER, INTENT(IN) :: nq ! number of tracer species
REAL, INTENT(IN) :: ptime
REAL, INTENT(IN) :: ptimestep

REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
REAL, INTENT(IN) :: pt(ngrid,nlayer)      ! temperature at mid-layer (K)
REAL, INTENT(IN) :: pdtfi(ngrid,nlayer)   ! tendancy temperature at mid-layer (K)
REAL, INTENT(IN) :: pplay(ngrid,nlayer)     ! pressure at middle of the layers
REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels

!--------------------------------------------------------
! Output Variables
!--------------------------------------------------------
REAL, INTENT(OUT) :: pdqcoag(ngrid,nlayer,nq) ! tendancy field for dust coagulation

!-----------------------------------------------------------------------
!   Local variables
!-----------------------------------------------------------------------
integer  :: ie, je, id, jd, iq, i, j, k, l, ii, jj 
real dens,dev,cst,sig0 
real mfp0,pt0,rad0,rad10,rad20,rho0,pres0 
real term1,term2,kernel,norm 
real, dimension(ngrid, nlayer) :: r0, rho, rn, ntot, mtot, mtotnew
real, dimension(ngrid, nlayer, nq) :: numb_new, numb_ini, numb_ini_dis, mass_ini, mass_new
real, dimension(ngrid, nlayer, nres_coag) :: ndis, ndis_new
real, dimension(ngrid,nlayer,nres_coag,nq) :: ndis_tab,rat_tab,ndis_tab_new
real, dimension(ngrid, nlayer) :: zt
real, dimension(ngrid, nlayer, nq) :: zq

logical, save :: firstcall = .true.
!$OMP THREADPRIVATE(firstcall)
!for boxinterp function:
integer :: t1,t2,p1,p2,m1,m2

!======================================================================
!   Main routine 

!-----------------------------------------------------------------------
!*** 1) Initialisations
!-----------------------------------------------------------------------

!!! Particles properties 
dens    = rho_dust ! density dust 
sig0     = dev_dt    ! Standard deviation of the dust distribution
cst     = .75 / (pi*dens) * exp( -4.5*sig0**2. )   ! fixed distribution

!!! Radius and Volumes defined in coagul_init
! rads_coag, vols_coag, deltar_coag, vrat_coag

!!! Getting tracers and temperature fields
zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep
zt(:,:)=pt(:,:)+pdtfi(:,:)*ptimestep

!!! Atmospheric state
DO l=1,nlayer
   rho(:,l)=pplay(:,l)/(r*pt(:,l))
ENDDO

!!! Initializing tracers and other fields
pdqcoag(:,:,:)=0.d0
ndis_tab(:,:,:,:)=0.d0
ndis_new(:,:,:)=0.d0
mtot(:,:)=0.d0
rat_tab(:,:,:,:) = 0.d0
numb_ini_dis(:,:,:) = 0.d0
numb_ini(:,:,:) = 0.d0
numb_new(:,:,:) = 0.d0
mass_ini(:,:,:) = 0.d0
mass_new(:,:,:) = 0.d0

!!! Preparing the initial log normal distribution
do iq=1,nq
     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
       mass_ini(:,:,iq)=max(zq(:,:,iq),1e-15)
       numb_ini(:,:,iq)=max(zq(:,:,iq-1),1e-15)

       !!! Get total mass to ensure mass conservation
       mtot(:,:)=mtot(:,:)+mass_ini(:,:,iq)

       !!! Initial Lognormal distribution
       ! sig0 remains constant
       r0(:,:)=(3/4.*mass_ini(:,:,iq)/(pi*numb_ini(:,:,iq)*dens))**(1/3.)*exp(-1.5*sig0**2)

       do i=1,nres_coag ! Sum of the ndis for all dust modes
          ndis_tab(:,:,i,iq)=numb_ini(:,:,iq)*rho(:,:)*deltar_coag(i)/(2.*rads_coag(i)*(2*pi)**0.5*sig0)*exp(-0.5*(log(rads_coag(i)/r0(:,:)))**2/sig0**2)
       enddo
       !!! Initial number mixing ratio from bin discretization nb/kg
       numb_ini_dis(:,:,iq)=sum(ndis_tab(:,:,:,iq),3)/rho(:,:)
       numb_new(:,:,iq)=numb_ini_dis(:,:,iq)
     endif
enddo


!!! Defining initial number ratio for each mode
do iq=1,nq
     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
       rat_tab(:,:,:,iq)=ndis_tab(:,:,:,iq)/sum(ndis_tab(:,:,:,:),4)
     endif 
enddo

!!! Normalizing the ratio (make sure sum=1)
do iq=1,nq
     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
       rat_tab(:,:,:,iq)=rat_tab(:,:,:,iq)/sum(rat_tab(:,:,:,:),4)
   endif
enddo

!!! Total initial distribution nb/m3 per bin
ndis(:,:,:)=sum(ndis_tab(:,:,:,:),4) 
!!! Total initial distribution nb/kg
ntot(:,:)=sum(ndis(:,:,:),3)/rho(:,:) 

!-----------------------------------------------------------------------
!*** 2) Coagulation
!-----------------------------------------------------------------------
! Spatial loop => 1D

if (full_coag_equations) then ! computing full coagulation equations 
 do i=1,ngrid
  do j=1,nlayer
      if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold
       pt0=pt(i,j)
       rho0=rho(i,j) 
       mfp0=mfp(pt0,rho0)

       ! Bin loop
       do k=1,nres_coag 

         ! Term 1
         term1=0.
         do jj=1,k
         do ii=1,k-1
             kernel=0.
             if (dust_coag_kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(ii),rads_coag(jj))
             if (dust_coag_kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(ii),rads_coag(jj))
             if (dust_coag_kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj))
             if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj)) 
             term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj)
          enddo
         enddo
         term1=term1*ptimestep/vols_coag(k) 

         ! Term 2
         term2=0.
         do jj=1,nres_coag
            kernel=0.
            if (dust_coag_kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(k),rads_coag(jj))
            if (dust_coag_kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(k),rads_coag(jj))
            if (dust_coag_kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj))
            if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj))
            term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj) 
         enddo
         term2=term2*ptimestep

         ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2)

       enddo

       norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:))
       ndis_new(i,j,:)=ndis_new(i,j,:)*norm 

      endif
  enddo
 enddo

else ! full_coag_equations

 ! Using Lookup tables
 do i=1,ngrid
  do j=1,nlayer
      if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold
       pt0=pt(i,j) 
       rho0=rho(i,j)
       mfp0=mfp(pt0,rho0)
       pres0=pplay(i,j)
       !find the corners for interpolation       
       t2=findval(table_pt,pt0) 
       if (t2==0) t2=table_numt
       if (t2==1) t2=2
       t1=t2-1
       m2=findval(table_mfp,mfp0)
       if (m2==0) m2=table_numm
       if (m2==1) m2=2
       m1=m2-1
       p2=findval(table_pres,pres0)
       if (p2==0) p2=table_numt
       if (p2==1) p2=2
       p1=p2-1

       ! Bin loop
       do k=1,nres_coag

         ! Term 1
         term1=0.
         do jj=1,k
          do ii=1,k-1 
             kernel=0.
             if (dust_coag_kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,ii,jj),       &   
                                                   table_b(t1,m2,ii,jj),       &
                                                   table_b(t2,m1,ii,jj),       &
                                                   table_b(t2,m2,ii,jj),       &
                                                   table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
             if (dust_coag_kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,ii,jj),       &
                                                   table_g(t1,m2,ii,jj),       &
                                                   table_g(t2,m1,ii,jj),       &
                                                   table_g(t2,m2,ii,jj),       &
                                                   table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
             if (dust_coag_kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,ii,jj),    &
                                                   table_de(t1,p2,ii,jj),     &
                                                   table_de(t2,p1,ii,jj),     &
                                                   table_de(t2,p2,ii,jj),     &
                                                   table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0))
             ! if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),ii,jj)  !not implemented yet
             term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj) 
          enddo
         enddo
         term1=term1*ptimestep/vols_coag(k)


         ! Term 2
         term2=0.
         do jj=1,nres_coag 
            kernel=0.
            if (dust_coag_kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,k,jj),       &
                                                  table_b(t1,m2,k,jj),       &
                                                  table_b(t2,m1,k,jj),       &
                                                  table_b(t2,m2,k,jj),       &
                                                  table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
             if (dust_coag_kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,k,jj),       &
                                                   table_g(t1,m2,k,jj),       &
                                                   table_g(t2,m1,k,jj),       &
                                                   table_g(t2,m2,k,jj),       &
                                                   table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
             if (dust_coag_kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,k,jj),       &
                                                   table_de(t1,p2,k,jj),       &
                                                   table_de(t2,p1,k,jj),       &
                                                   table_de(t2,p2,k,jj),       &
                                                   table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0))
            ! if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),k,jj)  !not implemented yet
            term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj) 
         enddo
         term2=term2*ptimestep

         ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2) 
       enddo
      
      !! Volume conservation
      norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:))
      ndis_new(i,j,:)=ndis_new(i,j,:)*norm

      endif 

  enddo
 enddo 
endif ! full_coag_equations

!-----------------------------------------------------------------------
!*** 3) new moments mass and number
!-----------------------------------------------------------------------
mtotnew(:,:)=1.d-15 

do iq=1,nq
     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then 
      !! Estimate new distribution in each dust mode using the initial number ratios 
      ndis_tab_new(:,:,:,iq)=ndis_new(:,:,:)*rat_tab(:,:,:,iq)      

      do i=1,ngrid
       do j=1,nlayer
         !! New number and mass mixing ratio in each mode
         numb_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq))/rho(i,j)
         mass_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq)*vols_coag(:))*dens/rho(i,j)
       enddo
      enddo

      !! New total mass
      mtotnew(:,:)=mtotnew(:,:)+mass_new(:,:,iq)
   endif 
enddo

do iq=1,nq
     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
      !! mass conservation
      mass_new(:,:,iq)=mass_new(:,:,iq)/mtotnew(:,:)*mtot(:,:) 
      !! New tendancy
      pdqcoag(:,:,iq)=(mass_new(:,:,iq)-mass_ini(:,:,iq))/ptimestep 
      !pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini(:,:,iq))/ptimestep 
      pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini_dis(:,:,iq))/ptimestep 
      !! check
      where (zq(:,:,iq)+pdqcoag(:,:,iq)*ptimestep.le.0.) 
        pdqcoag(:,:,iq)=-zq(:,:,iq)/ptimestep+1.e-15
      end where
      where (zq(:,:,iq-1)+pdqcoag(:,:,iq-1)*ptimestep.le.0.)
        pdqcoag(:,:,iq-1)=-zq(:,:,iq-1)/ptimestep+1.e-14
      end where
   endif
enddo

end subroutine dust_coagulation_main


!======================================================================
!======================================================================
real function dynvis(pt)
  implicit none
  ! dynamic viscosity following sutherland's formula for CO2
  ! Pa.s = kg m-1 s-1
  real, intent(in) :: pt
  dynvis=1.37e-5*(273.15 + 222.)/(pt + 222.)*(pt/273.15)**(3./2.)
end function dynvis
  
real function kinvis(pt,rho)
  implicit none
  real, intent(in) :: pt,rho
  kinvis=dynvis(pt)/rho
end function kinvis 

real function thervel(pt,massm)
  implicit none
  ! thermal velocity of an air molecule or a dust particle
  ! m s-1
  real, intent(in) :: pt, massm
  thervel=(8*kbz*pt/(pi*massm))**0.5
end function thervel

real function mfp(pt,rho)
  implicit none
  ! thermal velocity of an air molecule 
  ! m s-1
  real, intent(in) :: pt, rho
  mfp=2.*kinvis(pt,rho)/thervel(pt,m0co2) ! mean free path (m)
end function mfp

real function cun(pt,mfp,rad)
  implicit none
  !! Cunningham slip flow correction
  real, intent(in) :: pt, mfp, rad
  cun=1.+mfp/rad*(1.246+0.42*exp(-0.87/(mfp/rad)))
end function cun

real function diff(pt,mfp,rad)
  implicit none
  !! diffusion coefficient
  real, intent(in) :: pt, mfp, rad
  diff=kbz*pt/(6.*pi*rad*dynvis(pt))*cun(pt,mfp,rad)
end function diff

real function lambda(pt, mfp, rad)
  implicit none
  !! particle mean free path 
  real, intent(in) :: pt,mfp,rad
  real massm,dens
  dens    = rho_dust  ! density dust 
  massm=dens*4./3.*pi*rad**3
  lambda=2.*diff(pt, mfp, rad)/(pi*thervel(pt,massm))
end function lambda

real function delta(pt,mfp,rad)
  implicit none
  !! mean distance 
  real, intent(in) :: pt,mfp,rad
  real :: ltmp

  ltmp=lambda(pt,mfp,rad)
  delta=( (2.*rad+ltmp)**3 - (4.*rad**2+ltmp**2)**1.5 ) /(6*rad*ltmp) - 2*rad
end function delta

real function fallv(pt,mfp,rad)
  implicit none
  !! fall velocity
  real, intent(in) :: pt,mfp,rad
  real dens
  dens    = rho_dust  ! density dust
  fallv=2./9.*rad**2*dens*g/dynvis(pt)*cun(pt,mfp,rad)
end function fallv

real function reyn(pt,mfp,rad,rho)
  implicit none
  !! Reynold number
  real, intent(in) :: pt,rad,mfp,rho
  reyn=2.*rad*fallv(pt,mfp,rad)/kinvis(pt,rho)
end function reyn

real function stokes(pt,mfp,rad1,rad2)
  implicit none
  !! Stokes number
  real, intent(in) :: pt,mfp,rad1,rad2
  real :: f1

  f1=fallv(pt,mfp,rad1)
  stokes=f1*abs(fallv(pt,mfp,rad2)-f1)/(rad2*g)
end function stokes

real function schmidt(pt,mfp,rho,rad)
  implicit none
  !! Schmidt number
  real, intent(in) :: pt,mfp,rho,rad
  schmidt=kinvis(pt,rho)/diff(pt,mfp,rad)
end function schmidt

real function coal(pt,mfp,rad1,rad2,rho)
  implicit none
  !! Coalescence efficiency
  real, intent(in) :: pt,mfp,rad1,rad2,rho
  real :: coalea,coalev
  real :: stmp,rtmp

  stmp=stokes(pt,mfp,rad1,rad2)
  rtmp=reyn(pt,mfp,rad2,rho)
  coalea=stmp**2/(stmp+0.5)**2
  if (stmp>1.214) then
    coalev=(1+0.75*log(2.*stmp)/(stmp-1.214))**(-2)
  else
    coalev=0.
  endif
  coal=(60.*coalev+coalea*rtmp)/(60.+rtmp)
end function coal

real function betab(pt,mfp,rad1,rad2)
  implicit none
  real, intent(in) :: pt,mfp,rad1,rad2
  real :: m1,m2,num,den1,den2
  real :: dens
  real :: d1,d2

  dens    = rho_dust  ! density dust
  d1=diff(pt,mfp,rad1)
  d2=diff(pt,mfp,rad2)
  !! kernel transition regime
  ! Mass of 1 particle 1 and 2
  m1=dens*4./3.*pi*rad1**3
  m2=dens*4./3.*pi*rad2**3
  num=4.*pi*(rad1+rad2)*(d1+d2)
  den1=(rad1+rad2)/(rad1+rad2+(delta(pt,mfp,rad1)**2+delta(pt,mfp,rad2)**2)**0.5)
  den2=4./effb*(d1+d2)/((thervel(pt,m1)**2+thervel(pt,m2)**2)**0.5*(rad1+rad2))
  betab=num/(den1+den2)
end function betab

real function betag(pt,mfp,rad1,rad2)
  implicit none
  real, intent(in) :: pt,mfp,rad1,rad2
  real w1,w2,ee
  real :: radt
  !! kernel gravitation
  ! fall vel
  w1=fallv(pt,mfp,rad1)
  w2=fallv(pt,mfp,rad2)
  radt=(rad1+rad2)**2
  if (coal_coeff_mode.eq.0) then ! E=1
    ee=1.
  elseif (coal_coeff_mode.eq.1) then
    ee=1.5*(min(rad1,rad2))**2/radt
  elseif (coal_coeff_mode.eq.2) then
    ee=0.25*(min(rad1,rad2))**2/radt
  endif
  !b=coal(pt,r1,r2)*pi*(r1+r2)**2*abs(w1-w2)
  betag=ee*pi*radt*abs(w1-w2)
end function betag

real function betade(pt,mfp,rho,rad1,rad2)
  implicit none
  real, intent(in) :: pt,mfp,rad1,rad2,rho
  real rd,sc
  !! kernel diffusion enhancement
  rd=reyn(pt,mfp,max(rad1,rad2),rho)
  sc=schmidt(pt,mfp,rho,min(rad1,rad2))
  if (rd.le.1.) then
    betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.)
  else
    betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.)
  endif
end function betade

real function betade_tabfunc(pt,pres,rad1,rad2)
  implicit none
  real, intent(in) :: pt,pres,rad1,rad2
  real rd,sc,rho,mf
  !! kernel diffusion enhancement
  rho=pres/(pt*r)
  mf=mfp(pt,rho)
  rd=reyn(pt,mf,max(rad1,rad2),rho)
  sc=schmidt(pt,mf,rho,min(rad1,rad2))
  if (rd.le.1.) then
    betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.)
  else
    betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.)
  endif
end function betade_tabfunc

real function betati(pt,mfp,rho,rad1,rad2)
  implicit none
  real, intent(in) :: pt,mfp,rad1,rad2,rho
  real w1,w2,eps
  ! fall vel
  w1=fallv(pt,mfp,rad1)
  w2=fallv(pt,mfp,rad2)
  betati=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2)
end function betati

real function betati_tabfunc(pt,pres,rad1,rad2)
  implicit none
  real, intent(in) :: pt,rad1,rad2
  real w1,w2,eps,mf,rho,pres
  rho=pres/(pt*r)
  mf=mfp(pt,rho)
  ! fall vel
  w1=fallv(pt,mf,rad1)
  w2=fallv(pt,mf,rad2)
  betati_tabfunc=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2)
end function betati_tabfunc

real function frac(i,j,k)
  implicit none
  integer, intent(in) :: i,j,k
  real vint 
  ! intermediate volume
  vint=vols_coag(i)+vols_coag(j)
  if (k.lt.nres_coag) then
     if ( (vint.ge.vols_coag(k)) .and. (vint.lt.vols_coag(k+1))) then 
         frac=(vols_coag(k+1)-vint)/(vols_coag(k+1)-vols_coag(k))*vols_coag(k)/vint
         return
     endif
  endif
  if ((k.eq.nres_coag) .and. (vint.ge.vols_coag(k))) then 
     frac=1.
     return
  endif
  if (k.gt.1) then
     if ( (vint.lt.vols_coag(k)) .and. (vint.gt.vols_coag(k-1))) then 
        frac=(-vols_coag(k-1)+vint)/(-vols_coag(k-1)+vols_coag(k))*vols_coag(k)/vint
        return
     endif
  endif

  frac=0.
  return
end function frac

! *********************************************************
! Building lookup tables : called in the firstcall
! *********************************************************

subroutine make_tables()

implicit none
integer :: i,j,k,l
real :: t1,t2,p1,p2 

do i = 1,table_numt
    table_pt(i)=50.+25.*real(i) 
end do

do j = 1,table_nump
    table_pres(j)=(1.e-8)*10.**(real(j)/2.) 
end do

do j = 1,table_numm
    table_mfp(j)=(1.e-8)*10.**(real(j)/2.) 
end do

do i = 1,table_numt
    do j = 1,table_numm
        do k = 1,nres_coag
            do l = 1,nres_coag
                table_b(i,j,k,l)=betab(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l)) 
                table_g(i,j,k,l)=betag(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l)) 
            end do
        end do
    end do
end do

do i = 1,table_numt
    do j = 1,table_nump
        do k = 1,nres_coag
            do l = 1,nres_coag
                table_de(i,j,k,l)=betade_tabfunc(table_pt(i),table_pres(j),rads_coag(k),rads_coag(l)) 
            end do
        end do
    end do
end do

end subroutine make_tables

!======================================================================C
!======================================================================C

real function boxinterp(ft1p1,ft2p1,ft1p2,ft2p2,tref1,tref2,pref1,  &
         pref2,tmid,pmid)
!
!   Calculate 2d interpolation from pressure and temperature
!----------------------------------------------------------------------!
!                   T2 .FT2P1                    .FT2P2
!
!
!                   T1 .FT1P1                    .FT1P2
!                      P1                        P2
!----------------------------------------------------------------------!

implicit none

real :: ft1p1, ft2p1, ft1p2, ft2p2, tref1, tref2
real :: pref1, pref2, tmid, pmid, ans1, ans2, ans

!======================================================================C

ans1 = ft1p1 + (ft2p1 - ft1p1)*(tmid - tref1)/(tref2 - tref1)
ans2 = ft1p2 + (ft2p2 - ft1p2)*(tmid - tref1)/(tref2 - tref1)
boxinterp  = ans1 + (ans2 - ans1)*(pmid - pref1)/(pref2 - pref1)

end function boxinterp

!======================================================================C
!======================================================================C

integer function findval(array,value)
implicit none
real :: array(:)
real :: value

findval=minloc(array,dim=1,mask = (array > value))

end function findval

end module dust_coagulation_mod