Index: trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F	(revision 1720)
+++ trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F	(revision 1816)
@@ -1,17 +1,16 @@
-
-
 c=======================================================================
       subroutine massflowrateCO2(P,T,Sat,Radius,Matm,Ic)
 c
-c     Determination of the mass transfer rate
+c     Determination of the mass transfer rate for CO2 condensation & 
+c     sublimation
 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 
+c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
+c           particle radius (Radius), molecular mass of the atm (Matm)
+c   output:  MASS FLUX Ic
+c
+c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
 c=======================================================================
       USE comcstfi_h
@@ -20,5 +19,4 @@
 
       include "microphys.h"
-c      include "microphysCO2.h"
 
 
@@ -26,5 +24,4 @@
 c   arguments: INPUT
 c   ----------
-
       REAL T,Matm
       REAL*8 SAT
@@ -33,10 +30,7 @@
 c   arguments: OUTPUT
 c   ----------
-
       DOUBLE PRECISION   Ic
-
 c   Local Variables
 c   ----------
-
       DOUBLE PRECISION   Tcm
       DOUBLE PRECISION   T_inf, T_sup, T_dT
@@ -45,48 +39,29 @@
       DOUBLE PRECISION   rtsafe
       DOUBLE PRECISION   left, fval, dfval
-
 c    function for newton-raphson iterative method
 c    --------------------------
-
       EXTERNAL classical
-
          
-      Tcm = 80!dble(T)    ! initialize pourquoi 0 et pas t(i)
-
+      Tcm = 110!dble(T)    ! initialize pourquoi 0 et pas t(i)
       T_inf = 0d0
       T_sup = 800d0
-
-      T_dT = 0.01  ! precision - mettre petit et limiter nb iteration? 
+      T_dT = 0.001  ! 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) .eqv. .true.)  C0=0d0
-
-   
+      if (isnan(C0) .eqv. .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
-
+      if (Tcm.LE.0d0 .or. isnan(Tcm)) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10
+         Tcm = T
       endif
-
-
 c     THEN COMPUTE MASS FLUX Ic from FINAL Tsurface (Tc)
-
       Ic = (Tcm-T)
-
-      Ic = cond*Ic/(-Lsub)
+      Ic = cond*Ic/Lsub
 c regarder de combien varie la solution Ic entre Tcm et Tcm+T_dT
-    
+  
       RETURN
 
@@ -95,7 +70,5 @@
 
 c****************************************************************
-
       FUNCTION rtsafe(funcd,x1,x2,xacc,Radius,C0,C1,C2)
-*
 *
 *     Newton Raphsen routine (Numerical Recipe)
@@ -114,5 +87,5 @@
       EXTERNAL funcd
 
-      PARAMETER (MAXIT=50000)	!Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection, 
+      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 
@@ -128,51 +101,30 @@
 
       if ((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.) ) then
-         
-         
          x1=0d0
-         x2=500d0
-         
+         x2=1200d0         
          call funcd(x1,fl,df,C0,C1,C2)
          call funcd(x2,fh,df,C0,C1,C2)
-         
-         write(*,*) 'root must be bracketed in rtsafe'
+c         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.
-						 						
+      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
@@ -182,24 +134,15 @@
              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 
@@ -207,10 +150,8 @@
          xh=rtsafe 
         endif
-
-
 11    ENDDO
 
-      write(*,*) 'rtsafe exceeding maximum iterations'
-
+      rtsafe=0d0
+      write(*,*) 'rtsafe exceeding maximum iterations,Tcm=',rtsafe
       return 
 
@@ -219,10 +160,6 @@
 
 c********************************************************************************
-
       subroutine classical(x,f,df,C0,C1,C2)
-             
 c     Function to give as input to RTSAFE (NEWTON-RAPHOEN)
-
-
 c********************************************************************************
 
@@ -231,10 +168,6 @@
       DOUBLE PRECISION   x
       DOUBLE PRECISION  C0,C1,C2
-
       DOUBLE PRECISION f
       DOUBLE PRECISION  df
-      
-
-     
 
       f   = x + C0*exp(C1*x) - C2   ! start f
@@ -242,5 +175,4 @@
   
       return
-
       END  
 
@@ -250,13 +182,10 @@
 
 c********************************************************************************
-c defini la fonction eq 6 papier 2014 
+c defini la fonction eq 6 papier 2014  (Listowski et al., 2014)
       use tracer_mod, only: rho_ice_co2
       USE comcstfi_h
 
       implicit none
-
       include "microphys.h"
-c      include "microphysCO2.h"
-    
 
 c   arguments: INPUT
@@ -270,6 +199,4 @@
 c   local:
 c   ------
-
-
       DOUBLE PRECISION Cpatm,Cpn2,Cpco2
       DOUBLE PRECISION psat, xinf, pco2
@@ -298,14 +225,9 @@
 
 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     
@@ -313,13 +235,9 @@
       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)  
@@ -328,41 +246,23 @@
 
 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
@@ -370,16 +270,11 @@
 
 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
@@ -425,10 +320,7 @@
 c======================================================================
 
-
        implicit none
        
        include "microphys.h"
-c       include "microphysCO2.h"
-
 c      arguments
 c     -----------
@@ -457,13 +349,10 @@
          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)
@@ -471,13 +360,9 @@
       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) )
@@ -487,10 +372,6 @@
      &                                                          /Gamma2
       
-         
 c     Coefficient of Mason and Saxena
-
-
       epsilon = 1.
-
 
       A11 = 1.
@@ -498,5 +379,4 @@
       A22 = 1.
 
-
       A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
      &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))
@@ -505,5 +385,4 @@
      &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))
 
-
 c     INDIVIDUAL COND.
 
@@ -511,19 +390,13 @@
          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======================================================================
-
+c======================================================================
          subroutine KthN2LemJac(kthn2,T,rho)
-
 c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
 cWITH viscosity
@@ -565,5 +438,4 @@
 
         DOUBLE PRECISION k1, k2
-
 
          N1 = 1.511d0
@@ -612,12 +484,7 @@
          gamma9 = 1.
 
-
-
-
-c----------------------------------------------------------------------  
-         
+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
@@ -626,12 +493,6 @@
          delta = rho/rhoc 
 
-
-
-         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1
-      
-
+         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)         
@@ -642,7 +503,5 @@
      &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) 
 
-
          kthn2 = k1 + k2
-
 
          RETURN
@@ -744,6 +603,4 @@
       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) 
@@ -763,5 +620,4 @@
       g10 = 5.5
 
-
       h1 = 1.
       h2 = 5.
@@ -786,10 +642,6 @@
       n10 = 0.00433043347
 
-
-
       Tr   = T/Tc
       rhor = rho/rhoc
-
-
 
       k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) 
@@ -802,9 +654,7 @@
     
       k2  = exp(-5.*rhor**(2.)) * k2
-        
-        
+                
       kthco2 = (k1 + k2) *  Lambdac   ! mW
 
-
       RETURN
 
