Index: trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90	(revision 2985)
+++ trunk/LMDZ.COMMON/libf/evolution/co2glaciers_mod.F90	(revision 2993)
@@ -1,5 +1,6 @@
-      module co2glaciers_mod
-        implicit none
-        LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
+module co2glaciers_mod
+       
+ implicit none 
+ LOGICAL  co2glaciersflow ! True by default, to compute co2 ice flow. Read in  pem.def
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -33,6 +34,6 @@
 
 ! Inputs:
-      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat ! # number of time sample, physical points, subslopes, index of the flat subslope
-      REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x SLopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
+      INTEGER,INTENT(IN) :: timelen,ngrid,nslope,iflat !  number of time sample, physical points, subslopes, index of the flat subslope
+      REAL,INTENT(IN) :: subslope_dist(ngrid,nslope), def_slope_mean(ngrid) ! Physical points x Slopes : Distribution of the subgrid slopes; Slopes: values of the sub grid slope angles
       REAL,INTENT(IN) :: vmr_co2_PEM(ngrid,timelen) ! Physical x Time field : VMR of co2 in the first layer [mol/mol]
       REAL,INTENT(IN) :: ps_GCM(ngrid,timelen)      ! Physical x Time field: surface pressure given by the GCM [Pa]
@@ -55,5 +56,5 @@
       call compute_hmaxglaciers_co2(ngrid,nslope,iflat,Tcond,def_slope_mean,hmax)
 
-      call  transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
+      call transfer_co2ice_duringflow(ngrid,nslope,iflat, subslope_dist,def_slope_mean,hmax,Tcond,co2ice_slope,flag_co2flow,flag_co2flow_mesh)
    RETURN
 end subroutine
@@ -70,5 +71,5 @@
 !!!          before initating flow
 !!!          
-!!! Author: LL, based on theoretical work by A.Grau Galofre (LPG)
+!!! Author: LL,based on  work by A.Grau Galofre (LPG) and Isaac Smith (JGR Planets 2022)
 !!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -87,40 +88,16 @@
       REAL,INTENT(OUT) :: hmax(ngrid,nslope) ! Physical grid x Slope field: maximum co2 thickness before flaw [m]
 ! Local
-      INTEGER,PARAMETER :: n = 7 ! flow law exponent Nye et al., 2000
-      REAL,PARAMETER :: Rg = 8.3145 ! gas constant [J/K/mol]
-      REAL,PARAMETER :: Q = 59000. ! Activation Energy [J/mol], Nye et al., 2000
-      DOUBLE PRECISION,PARAMETER :: C = 1.8138e-21 ! Nye et al., 2000 [s/m^n]
-      DOUBLE PRECISION :: Ad = 1.202e11 ! Softness prefactor [MPa^-n] Nye et al., 2000
-      REAL :: Ro,Ho, S,ratio ! gemoetry from Nye et al., 2000
-      DOUBLE PRECISION :: A,Ao ! softness parameter [s/m^n]
-      DOUBLE PRECISION :: C1 ! intermediate variable
-      DOUBLE PRECISION :: t_0 ! relaxation time (assuming radial symetry) [s]
-      DOUBLE PRECISION :: u ! characteristic horizontal deformation rate [m/s]
-      DOUBLE PRECISION :: tau_d ! characteristic basal drag, understood as the strest that an ice CO2 mass flowing under its weight balanced by viscosity
+      DOUBLE PRECISION :: tau_d = 5.e3 ! characteristic basal drag, understood as the stress that an ice CO2 mass flowing under its weight balanced by viscosity. Value obtained from I.Smith
       REAL :: rho_co2(ngrid) ! co2 ice density [kg/m^3]
       INTEGER :: ig,islope ! loop variables
       REAL :: slo_angle
 
-! 0. Geometry parameters
-      Ro = 200e3
-      Ho = 3000.
-      ratio = 2./3.
-      S = Ho/Ro*1/((2+1./n)*(1+1./n))*(1-ratio**(1+1./n))**(1./(2+1./n)-1)*ratio**(1+1./n-1) 
-! 1. Flow parameters 
-     Ao = 3**(1./(2*n+2))*Ad
-     do ig = 1,n
-     Ao = Ao*1e-6
-     enddo
-! 2. Compute rho_co2
+! 1. Compute rho_co2
     DO ig = 1,ngrid
       rho_co2(ig) = (1.72391 - 2.53e-4*Tcond(ig)-2.87*1e-7*Tcond(ig)**2)*1e3  ! Mangan et al. 2017
     ENDDO
+
 ! 3. Compute max thickness
     DO ig = 1,ngrid
-       A = Ao*exp(-Q/(Rg*Tcond(ig)))
-       C1 = A*(rho_co2(ig)*g)**(float(n))/float(n+2)
-       t_0 = Ro/(C1*(5*n+3))*(float(2*n+1)/float(n+1))**n
-       u = Ro/t_0
-       tau_d = (u*Ho*(rho_co2(ig)*g)**2*S**2/(2*A))**(1./(n+2))      
        DO islope = 1,nslope
           if(islope.eq.iflat) then
@@ -256,8 +233,8 @@
         ENDDO
         Tcond(ig) = ave/timelen
-       ENDDO
+      ENDDO
 RETURN
 
 
-      end subroutine
-      end module
+end subroutine
+end module
