module nucleaco2_mod
  implicit none
  contains
subroutine nucleaco2(pco2,temp,sat,n_ccn,nucrate,vo2co2, teta)
USE comcstfi_h

implicit none
!*                                                     *
!*   This subroutine computes the nucleation rate      *
!*   as given in Pruppacher & Klett (1978) in the      *
!*   case of water ice forming on a solid substrate.   *
!*     Definition refined by Keese (jgr,1989)          *
!*   Authors: F. Montmessin                            *
!*     Adapted for the LMD/GCM by J.-B. Madeleine      *
!*     (October 2011)                                  *
!*     Optimisation by A. Spiga (February 2012)        *
!*   CO2 nucleation routine dev. by Constantino        *
!*     Listowski and Joachim Audouard (2016-2017),     *
!*     adapted from the water ice nucleation
!* It computes two different nucleation rates : one
!* on the dust CCN distribution and the other one on
!* the water ice particles distribution
!*******************************************************
 ! nucrate = output
 ! nucrate_h2o en sortie aussi :
!nucleation sur dust et h2o separement ici

include "microphys.h"
include "callkeys.h"

!     Inputs
DOUBLE PRECISION, intent(in) :: pco2,sat,vo2co2
DOUBLE PRECISION, intent(in) :: n_ccn(nbinco2_cld)
REAL, intent(in) :: temp !temperature
REAL, intent(in) :: teta
!     Output
DOUBLE PRECISION, intent(out) :: nucrate(nbinco2_cld)

!     Local variables
DOUBLE PRECISION nco2
DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
!DOUBLE PRECISION fshapeco2   ! function defined at the end of the file
DOUBLE PRECISION deltaf
double precision mtetalocal ! local teta in double precision
double precision fshapeco2simple,zefshapeco2
integer i
!*************************************************

mtetalocal = dble(teta)

nco2   = pco2 / kbz / temp
rstar  = 2. * sigco2 * vo2co2 / (kbz*temp*log(sat))
gstar  = 4. * pi * (rstar * rstar * rstar) / (3.*vo2co2)

fshapeco2simple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal) / 4.

!c       Loop over size bins
do i=1,nbinco2_cld

  if ( n_ccn(i) .lt. 1e-10 ) then
    !c no dust, no need to compute nucleation!
    nucrate(i)=0.
  else
    if (rad_cldco2(i).gt.3000.*rstar) then
      zefshapeco2 = fshapeco2simple
    else
      zefshapeco2 = fshapeco2(mtetalocal, rad_cldco2(i)/rstar)
    endif

    fistar = (4./3.*pi) * sigco2 * (rstar * rstar) * zefshapeco2
    deltaf = (2.*desorpco2-surfdifco2-fistar) / (kbz*temp)
    deltaf = min( max(deltaf, -100.d0), 100.d0)

    if (deltaf.eq.-100.) then
        nucrate(i) = 0.
    else
        nucrate(i) = dble(sqrt ( fistar / (3.*pi*kbz*temp*(gstar*gstar)) ) * kbz * temp * rstar * rstar * 4. * pi &
            * ( nco2*rad_cldco2(i) ) * ( nco2*rad_cldco2(i) ) / ( zefshapeco2 * nusco2 * m0co2 )  * exp (deltaf))
    endif
  endif ! if n_ccn(i) .lt. 1e-10

enddo

end subroutine nucleaco2


!*********************************************************
double precision function fshapeco2(cost,rap)
implicit none
!*        function computing the f(m,x) factor           *
!* related to energy required to form a critical embryo  *
!*********************************************************

double precision cost,rap
double precision yeah

!! PHI
yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
!! FSHAPECO2 = TERM A
fshapeco2 = (1.-cost*rap) / yeah
fshapeco2 = fshapeco2 * fshapeco2 * fshapeco2
fshapeco2 = 1. + fshapeco2
!! ... + TERM B
yeah = (rap-cost)/yeah
fshapeco2 = fshapeco2 + rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
!! ... + TERM C
fshapeco2 = fshapeco2 + 3. * cost * rap * rap * (yeah-1.)
!! FACTOR 1/2
fshapeco2 = 0.5*fshapeco2

end function fshapeco2

end module nucleaco2_mod