c=======================================================================
      subroutine massflowrateCO2(P,T,Sat,Radius,Matm,Ic)
c
c     Determination of the mass transfer rate
c
c     newton-raphson method

c     CLASSICAL  (no SF etc.)
      
c     AUTOMATIC SETTING OF RANGES FOR NEWTON-RAPHSON FOR THE PAPER

c     MASS FLUX Ic 

c=======================================================================
      USE comcstfi_h

      implicit none

      include "microphys.h"
c      include "microphysCO2.h"



c   arguments: INPUT
c   ----------

      REAL T,Matm
      REAL*8 SAT
      real P
      DOUBLE PRECISION Radius
c   arguments: OUTPUT
c   ----------

      DOUBLE PRECISION   Ic

c   Local Variables
c   ----------

      DOUBLE PRECISION   Tcm
      DOUBLE PRECISION   T_inf, T_sup, T_dT
      DOUBLE PRECISION   C0,C1,C2
      DOUBLE PRECISION   kmix,Lsub,cond
      DOUBLE PRECISION   rtsafe
      DOUBLE PRECISION   left, fval, dfval

c    function for newton-raphson iterative method
c    --------------------------

      EXTERNAL classical

         
      Tcm = dble(T)    ! initialize pourquoi 0 et pas t(i)

      T_inf = 0d0
      T_sup = 200d0

      T_dT = 0.1  ! precision - mettre petit et limiter nb iteration? 
      
666   CONTINUE

c      print*, 'Radius ', Radius
c      print*, 'SAT = ', Sat
      call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2)

      if (isnan(C0) .eq. .true.)  C0=0d0

   
c     FIND SURFACE TEMPERATURE (Tc) : iteration sur t 

      cond = 4.*pi*Radius*kmix
 
      Tcm = rtsafe(classical,T_inf,T_sup,T_dT,Radius,C0,C1,C2)


      if (Tcm.LE.0d0) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10

            Tcm = 0d0

      endif


c     THEN COMPUTE MASS FLUX Ic from FINAL Tsurface (Tc)

      Ic = (Tcm-T)

      Ic = cond*Ic/(-Lsub)
c regarder de combien varie la solution Ic entre Tcm et Tcm+T_dT
    
      RETURN

      END


c****************************************************************

      FUNCTION rtsafe(funcd,x1,x2,xacc,Radius,C0,C1,C2)
*
*
*     Newton Raphsen routine (Numerical Recipe)
*
c****************************************************************

      implicit none

      INTEGER MAXIT 
       DOUBLE PRECISION x1,x2,xacc 
      DOUBLE PRECISION rtsafe
      DOUBLE PRECISION Radius
      DOUBLE PRECISION C0,C1,C2


      EXTERNAL funcd

      PARAMETER (MAXIT=10000)	!Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection, 
				!find the root of a function bracketed between x1 and x2. The root, returned as the function value rtsafe,
				!will be refined until its accuracy is known within !!±xacc. funcd is a user-supplied subroutine which 
				!returns both the function value and the first derivative of the function.

      INTEGER j 
      DOUBLE PRECISION df,dx,dxold
      DOUBLE PRECISION f,fh,fl,temp,xh,xl 

      call funcd(x1,fl,df,C0,C1,C2)
      call funcd(x2,fh,df,C0,C1,C2)


      if ((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.) ) then
         
         
         x1=0d0
         x2=500d0
         
         call funcd(x1,fl,df,C0,C1,C2)
         call funcd(x2,fh,df,C0,C1,C2)
         
         write(*,*) 'root must be bracketed in rtsafe'
      endif


      if (fl.eq.0.) then

	 rtsafe=x1
	 return 

      else if (fh.eq.0.) then

	 rtsafe=x2
         return 

      else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.
						 						
	 xl=x1
	 xh=x2 

      else

	xh=x1
	xl=x2 

      endif

      rtsafe = .5*(x1+x2) 	!Initialize the guess for root,
      dxold  = abs(x2-x1)       !the stepsize before last,
      dx     = dxold  	        ! and the last step.


      call funcd(rtsafe,f,df,C0,C1,C2)

      DO 11 j=1,MAXIT		!Loop over allowed iterations. 

         
         !print*, 'iteration:', j
         !print*, rtsafe


         if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).gt.0.   ! Bisect if Newton out of range
     *    .or. abs(2.*f).gt.abs(dxold*df) ) then	       ! or not decreasing fst enough

             dxold=dx 
             dx=0.5*(xh-xl) 
             rtsafe=xl+dx 

             if (xl.eq.rtsafe) return			!Change in root is negligible. Newton step acceptable. Take it.

         else 

	    dxold=dx
	    dx=f/df 
	    temp=rtsafe

	    rtsafe=rtsafe-dx

            if(temp.eq.rtsafe) return 

         endif


        if(abs(dx).lt.xacc) return !Convergence criterion. The one new function evaluation per iteration. Maintain the bracket on the root.

         call funcd(rtsafe,f,df,C0,C1,C2)

        if(f.lt.0.) then
         xl=rtsafe 
        else
         xh=rtsafe 
        endif


11    ENDDO

      write(*,*) 'rtsafe exceeding maximum iterations'

      return 

      END


c********************************************************************************

      subroutine classical(x,f,df,C0,C1,C2)
             
c     Function to give as input to RTSAFE (NEWTON-RAPHOEN)


c********************************************************************************

      implicit none

      DOUBLE PRECISION   x
      DOUBLE PRECISION  C0,C1,C2

      DOUBLE PRECISION f
      DOUBLE PRECISION  df
      

     

      f   = x + C0*exp(C1*x) - C2   ! start f
      df  = 1. + C0*C1*exp(C1*x)    ! start df 
  
      return

      END  

c********************************************************************************

      subroutine coefffunc(P,T,S,rc,Matm,kmix,Lsub,C0,C1,C2)

c********************************************************************************
c defini la fonction eq 6 papier 2014 
      use tracer_mod, only: rho_ice_co2
      USE comcstfi_h

      implicit none

      include "microphys.h"
c      include "microphysCO2.h"
    

c   arguments: INPUT
c   ----------------
      REAL P
      real T
      REAL*8 S
      double precision rc
      REAL Matm !g.mol-1 ( = mmean(ig,l) )

c   local:
c   ------


      DOUBLE PRECISION Cpatm,Cpn2,Cpco2
      DOUBLE PRECISION psat, xinf, pco2
      DOUBLE PRECISION Dv     
      DOUBLE PRECISION l0,l1,l2,l3,l4            
      DOUBLE PRECISION knudsen, a, lambda      ! F and S correction
      DOUBLE PRECISION Ak                       ! kelvin factor    
      DOUBLE PRECISION vthatm,lpmt,rhoatm, vthco2 ! for Kn,th

c   arguments: OUTPUT
c   ----------
      DOUBLE PRECISION  C0,C1,C2
      DOUBLE PRECISION  kmix,Lsub

c     DEFINE heat cap. J.kg-1.K-1 and To

      data Cpco2/0.7e3/
      data Cpn2/1e3/

      kmix = 0d0
      Lsub = 0d0

      C0 = 0d0
      C1 = 0d0
      C2 = 0d0

c     Equilibirum pressure over a flat surface

      psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)

c     Compute transport coefficient

      pco2 = psat * dble(S)

c     Latent heat of sublimation if CO2  co2 (J.kg-1)
c     version Azreg_Ainou (J/kg) :

      l0=595594.      
      l1=903.111     
      l2=-11.5959    
      l3=0.0528288 
      l4=-0.000103183
  
      Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * 
     &     dble(T)**3 + l4 * dble(T)**4 ! J/kg

c     atmospheric density

      rhoatm = dble(P*Matm)/(rgp*dble(T))   ! g.m-3
      rhoatm = rhoatm * 1.00e-3 !kg.m-3

      call  KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2
      call  Diffcoeff(P, T, Dv)  

      Dv = Dv * 1.00e-4         !!! cm2.s-1  to m2.s-1

c     ----- FS correction for Diff

      vthco2  = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1

      knudsen = 3*Dv / (vthco2 * rc) 

      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
      
      Dv      = Dv / (1. + lambda * knudsen)
        
c     ----- FS correction for Kth 

      vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg 

      Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1

      lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/
     &     (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer

      knudsen = lpmt / rc

      lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)

      kmix    = kmix /  (1. + lambda * knudsen)

c     --------------------- ASSIGN coeff values for FUNCTION

      xinf = dble(S) * psat / dble(P)

      Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) ))
  
      C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/
     &     (rgp*dble(T)))
 
      C1 = Lsub*mco2/(rgp*dble(T)**2)

      C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
     
      RETURN
      END


c======================================================================

      subroutine Diffcoeff(P, T, Diff)

c     Compute diffusion coefficient CO2/N2
c     cited in Ilona's lecture - from Reid et al. 1987
c======================================================================


       IMPLICIT NONE

       include "microphys.h"
c       include "microphysCO2.h"

c      arguments
c     -----------
      
      REAL P
      REAL Pbar                     !!! has to be in bar for the formula
      REAL T
      
c     output
c     -----------

      DOUBLE PRECISION Diff
    
c      local
c     -----------

      DOUBLE PRECISION dva, dvb, Mab  ! Mab has to be in g.mol-1
        
        Pbar = P * 1d-5
    
        Mab = 2. / ( 1./mn2 + 1./mco2 ) * 1000.

  	dva = 26.9        ! diffusion volume of CO2,  Reid et al. 1987 (cited in Ilona's lecture)
  	dvb = 18.5        ! diffusion volume of N2
    
        Diff  = 0.00143 * dble(T)**(1.75) / (dble(Pbar) * sqrt(Mab) 
     &       * (dble(dva)**(1./3.) + dble(dvb)**(1./3.))**2.) !!! in cm2.s-1
  
       RETURN

       END


c======================================================================

         subroutine KthMixNEW(Kthmix,T,x,rho)

c        Compute thermal conductivity of CO2/N2 mixture
c         (***WITHOUT*** USE OF VISCOSITY)

c          (Mason & Saxena, 1958 - Wassiljeva 1904)
c======================================================================


       implicit none
       
       include "microphys.h"
c       include "microphysCO2.h"

c      arguments
c     -----------
         
         REAL T
         DOUBLE PRECISION x
         DOUBLE PRECISION rho !kg.m-3

c     outputs
c     -----------

         DOUBLE PRECISION Kthmix

c     local
c    ------------

         DOUBLE PRECISION x1,x2

         DOUBLE PRECISION  Tc1, Tc2, Pc1, Pc2

         DOUBLE PRECISION  A12, A11, A22, A21

         DOUBLE PRECISION  Gamma1, Gamma2, M1, M2
         DOUBLE PRECISION  lambda_trans1, lambda_trans2,epsilon

         DOUBLE PRECISION  kco2, kn2


      x1 = x
      x2 = 1d0 - x


      M1 = mco2
      M2 = mn2


      Tc1 =  304.1282 !(Scalabrin et al. 2006)
      Tc2 =  126.192  ! (Lemmon & Jacobsen 2003)

      Pc1 =  73.773   ! (bars)
      Pc2 =  33.958   ! (bars)
  
  
      Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.)
      Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.)


c Translational conductivities



      lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) )
     &                                                          /Gamma1

      lambda_trans2 = ( exp(0.0464 * T/Tc2) - exp(-0.2412 * T/Tc2) )
     &                                                          /Gamma2
      
         
c     Coefficient of Mason and Saxena


      epsilon = 1.


      A11 = 1.
	 
      A22 = 1.


      A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
     &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))

      A21 = epsilon * (1. + sqrt(lambda_trans2/lambda_trans1)*
     &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))


c     INDIVIDUAL COND.

         call KthCO2Scalab(kco2,T,rho)
         call KthN2LemJac(kn2,T,rho)


c     MIXTURE COND.

        Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
        Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1

        RETURN

        END


c======================================================================

         subroutine KthN2LemJac(kthn2,T,rho)

c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
cWITH viscosity
c======================================================================

       implicit none

        include "microphys.h"
c        include "microphysCO2.h"


c      arguments
c     -----------
         
         REAL T
         DOUBLE PRECISION rho !kg.m-3

c     outputs
c     -----------

         DOUBLE PRECISION kthn2

c     local
c    ------------

        DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
        DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
        DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
        DOUBLE PRECISION d4,d5,d6,d7,d8,d9
        DOUBLE PRECISION l4,l5,l6,l7,l8,l9
        DOUBLE PRECISION t2,t3,t4,t5,t6,t7,t8,t9
        DOUBLE PRECISION gamma4,gamma5,gamma6,gamma7,gamma8,gamma9

        DOUBLE PRECISION Tc,rhoc

        DOUBLE PRECISION tau, delta

        DOUBLE PRECISION visco

        DOUBLE PRECISION k1, k2


         N1 = 1.511d0
         N2 = 2.117d0
         N3 = -3.332d0

         N4 = 8.862
         N5 = 31.11
         N6 = -73.13
         N7 = 20.03
         N8 = -0.7096
         N9 = 0.2672

         t2 = -1.0d0
         t3 = -0.7d0
         t4 = 0.0d0
         t5 = 0.03
         t6 = 0.2
         t7 = 0.8
         t8 = 0.6
         t9 = 1.9
   
         d4 =  1.
         d5 =  2.
         d6 =  3.
         d7 =  4.
         d8 =  8.
         d9 = 10.
   
         l4 = 0. 
         gamma4 = 0.
        
         l5 = 0. 
         gamma5 = 0.
        
         l6 = 1. 
         gamma6 = 1.

         l7 = 2. 
         gamma7 = 1.

         l8 = 2. 
         gamma8 = 1.

         l9 = 2. 
         gamma9 = 1.




c----------------------------------------------------------------------  
         
         call viscoN2(T,visco)  !! v given in microPa.s

       
         Tc   = 126.192d0
         rhoc = 11.1839  * 1000 * mn2   !!!from mol.dm-3 to kg.m-3

         tau  = Tc / T
         delta = rho/rhoc 



         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1
      

c--------- residual thermal conductivity



         k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4)         
     &  +     N5 * tau**t5 * delta**d5 * exp(-gamma5*delta**l5) 
     &  +     N6 * tau**t6 * delta**d6 * exp(-gamma6*delta**l6) 
     &  +     N7 * tau**t7 * delta**d7 * exp(-gamma7*delta**l7) 
     &  +     N8 * tau**t8 * delta**d8 * exp(-gamma8*delta**l8) 
     &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) 


         kthn2 = k1 + k2


         RETURN

         END


c======================================================================

         subroutine viscoN2(T,visco)

c        Compute viscosity of N2 (Lemmon and Jacobsen, 2003)

c======================================================================

         implicit none

       include "microphys.h"
c       include "microphysCO2.h"
c      arguments
c     -----------
         
      REAL T

c     outputs
c     -----------

      DOUBLE PRECISION visco


c     local
c    ------------

      DOUBLE PRECISION a0,a1,a2,a3,a4 
      DOUBLE PRECISION Tstar,factor,sigma,M2
      DOUBLE PRECISION RGCS
      

c----------------------------------------------------------------------  

  
      factor = 98.94   ! (K)   
  
      sigma  = 0.3656  ! (nm)
  
      a0 =  0.431
      a1 = -0.4623
      a2 =  0.08406
      a3 =  0.005341
      a4 = -0.00331
  
      M2 = mn2 * 1.00e3   !!! to g.mol-1
  
      Tstar = T*1./factor

      RGCS = exp( a0 + a1 * log(Tstar) + a2 * (log(Tstar))**2. + 
     &                a3 * (log(Tstar))**3. + a4 * (log(Tstar))**4. )
  
  
      visco = 0.0266958 * sqrt(M2*T) / ( sigma**2. * RGCS )  !!! microPa.s


      RETURN

      END


c======================================================================

         subroutine KthCO2Scalab(kthco2,T,rho)

c        Compute thermal cond of CO2 (Scalabrin et al. 2006)

c======================================================================

         implicit none



c      arguments
c     -----------

      REAL T
      DOUBLE PRECISION rho

c     outputs
c     -----------

      DOUBLE PRECISION kthco2

c     LOCAL
c     -----------

      DOUBLE PRECISION Tc,Pc,rhoc, Lambdac
      
      DOUBLE PRECISION Tr, rhor, k1, k2

      DOUBLE PRECISION g1,g2,g3,g4,g5,g6,g7,g8,g9,g10
      DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
      DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10

      

      Tc   = 304.1282   !(K) 
      Pc   = 7.3773e6   !(MPa)
      rhoc = 467.6      !(kg.m-3)
      Lambdac = 4.81384 !(mW.m-1K-1)
  
      g1 = 0.
      g2 = 0.
      g3 = 1.5
      g4 = 0.0
      g5 = 1.0
      g6 = 1.5
      g7 = 1.5
      g8 = 1.5
      g9 = 3.5
      g10 = 5.5


      h1 = 1.
      h2 = 5.
      h3 = 1.
      h4 = 1.
      h5 = 2.
      h6 = 0.
      h7 = 5.0
      h8 = 9.0
      h9 = 0.
      h10 = 0.

      n1 = 7.69857587
      n2 = 0.159885811
      n3 = 1.56918621
      n4 = -6.73400790
      n5 = 16.3890156
      n6 = 3.69415242
      n7 = 22.3205514
      n8 = 66.1420950
      n9 = -0.171779133
      n10 = 0.00433043347



      Tr   = T/Tc
      rhor = rho/rhoc



      k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) 
     &     + n3*Tr**(g3)*rhor**(h3)

      k2 = n4*Tr**(g4)*rhor**(h4) + n5*Tr**(g5)*rhor**(h5)  
     &    + n6*Tr**(g6)*rhor**(h6) + n7*Tr**(g7)*rhor**(h7) 
     &    + n8*Tr**(g8)*rhor**(h8) + n9*Tr**(g9)*rhor**(h9) 
     &    + n10*Tr**(g10)*rhor**(h10)
    
      k2  = exp(-5.*rhor**(2.)) * k2
        
        
      kthco2 = (k1 + k2) *  Lambdac   ! mW


      RETURN

      END