Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 1941)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 1944)
@@ -30,8 +30,11 @@
       use comsoil_h, only: inertiedat, layer, mlayer, nsoilmx, tsoil
       use control_mod, only: day_step, iphysiq, anneeref, planet_type
+      use phyetat0_mod, only: phyetat0
       use phyredem, only: physdem0, physdem1
       use iostart, only: open_startphy
-      use dimradmars_mod, only: tauscaling
-      use turb_mod, only: q2
+      use dimradmars_mod, only: tauscaling, albedo
+      use turb_mod, only: q2, wstar
+      use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,
+     &                        mem_Nccn_co2
       use filtreg_mod, only: inifilr
       USE mod_const_mpi, ONLY: COMM_LMDZ
@@ -418,11 +421,12 @@
         fichnom = 'start.nc'
         CALL dynetat0(fichnom,vcov,ucov,teta,q,masse,
-     .       ps,phis,time)
+     &       ps,phis,time)
 
         write(*,*) 'Reading file STARTFI'
         fichnom = 'startfi.nc'
         CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
-     .        day_ini,time,
-     .        tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling,totcloudfrac)
+     &        day_ini,time,tsurf,tsoil,albedo,emis,
+     &        q2,qsurf,co2ice,tauscaling,totcloudfrac,
+     &        wstar,mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
         
         ! copy albedo and soil thermal inertia
@@ -1640,7 +1644,8 @@
      .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
-     .              dtphys,hour_ini,
-     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
-     .              totcloudfrac)
+     &              dtphys,hour_ini,
+     &              tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,tauscaling,
+     &              totcloudfrac,wstar,
+     &              mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
 
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F	(revision 1941)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F	(revision 1944)
@@ -31,4 +31,5 @@
       USE temps_mod, ONLY: day_ini,hour_ini
       USE iniphysiq_mod, ONLY: iniphysiq
+      USE phyetat0_mod, ONLY: phyetat0
       implicit none
 
@@ -67,4 +68,9 @@
       REAL q2(ngridmx,llm+1)
       REAL emis(ngridmx)
+      REAL albedo(ngridmx,2)
+      REAL wstar(ngridmx)
+      DOUBLE PRECISION mem_Nccn_co2(ngridmx,llm)
+      DOUBLE PRECISION mem_Mccn_co2(ngridmx,llm)
+      DOUBLE PRECISION mem_Mh2o_co2(ngridmx,llm)
       INTEGER start,length
       PARAMETER (length = 100)
@@ -157,6 +163,7 @@
 
       CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
-     &      day_ini_fi,timefi,tsurf,tsoil,emis,q2,qsurf,co2ice,
-     &      tauscaling,totcloudfrac)
+     &      day_ini_fi,timefi,tsurf,tsoil,albedo,emis,q2,qsurf,co2ice,
+     &      tauscaling,totcloudfrac,wstar,mem_Mccn_co2,mem_Nccn_co2,
+     &               mem_Mh2o_co2)
 
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1941)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1944)
@@ -21,4 +21,6 @@
      &                            ecritphy, iphysiq
       use dimradmars_mod, only: tauscaling,tauvis,totcloudfrac
+      use co2cloud_mod, only: mem_Mccn_co2, mem_Mh2o_co2,
+     &                        mem_Nccn_co2
       USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff,sig,
      &                       presnivs,pseudoalt,scaleheight
@@ -92,4 +94,6 @@
       REAL co2ice(1)        ! co2ice layer (kg.m-2)
       REAL emis(1)          ! surface layer
+      REAL albedo(1,1)      ! surface albedo
+      REAL :: wstar(1)=0.    ! Thermals vertical velocity
       REAL q2(nlayer+1)   ! Turbulent Kinetic Energy
       REAL zlay(nlayer)   ! altitude estimee dans les couches (km)
@@ -522,4 +526,5 @@
       call getin("albedo",albedodat(1))
       write(*,*) " albedo = ",albedodat(1)
+      albedo(1,1)=albedodat(1)
 
       inertiedat(1,1)=400 ! default value for inertiedat
@@ -593,5 +598,8 @@
       call getin("co2ice",co2ice)
       write(*,*) " co2ice = ",co2ice
-
+! Initialization for CO2 clouds (could be improved to read initial profiles)
+      mem_Mccn_co2(:,:)=0
+      mem_Mh2o_co2(:,:)=0
+      mem_Nccn_co2(:,:)=0
 c
 c  emissivity
@@ -710,7 +718,9 @@
      &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
       call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
-     .              dtphys,time,
-     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
-     .              totcloudfrac)
+     &              dtphys,time,
+     &              tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,tauscaling,
+     &              totcloudfrac,wstar,
+     &              mem_Mccn_co2,mem_Nccn_co2,
+     &              mem_Mh2o_co2)
 
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0.F90	(revision 1941)
+++ 	(revision )
@@ -1,363 +1,0 @@
-subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
-                     day_ini,time0,tsurf,tsoil,emis,q2,qsurf,co2ice, &
-                     tauscaling,totcloudfrac,mem_Mccn_co2,mem_Nccn_co2,&
-                     mem_Mh2o_co2)
-
-  use tracer_mod, only: noms ! tracer names
-  use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
-                       zmea, zstd, zsig, zgam, zthe
-  use iostart, only: nid_start, open_startphy, close_startphy, &
-                     get_field, get_var, inquire_field, &
-                     inquire_dimension, inquire_dimension_length
-  use ioipsl_getincom, only : getin
-
-  implicit none
-!======================================================================
-! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
-!  Adaptation à Mars : Yann Wanherdrick 
-! Objet: Lecture de l etat initial pour la physique
-! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
-!                      r4 or r8 restarts independently of having compiled
-!                      the GCM in r4 or r8)
-!         June 2013 TN : Possibility to read files with a time axis
-!         November 2013 EM : Enabeling parallel, using iostart module
-!======================================================================
-  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
-  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
-!======================================================================
-!  Arguments:
-!  ---------
-!  inputs:
-  character*(*),intent(in) :: fichnom ! "startfi.nc" file
-  integer,intent(in) :: tab0
-  integer,intent(in) :: Lmodif
-  integer,intent(in) :: nsoil ! # of soil layers
-  integer,intent(in) :: ngrid ! # of atmospheric columns
-  integer,intent(in) :: nlay ! # of atmospheric layers
-  integer,intent(in) :: nq
-  integer :: day_ini
-  real :: time0
-
-!  outputs:
-  real,intent(out) :: tsurf(ngrid) ! surface temperature
-  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
-  real,intent(out) :: emis(ngrid) ! surface emissivity
-  real,intent(out) :: q2(ngrid,nlay+1) ! 
-  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
-  real,intent(out) :: co2ice(ngrid) ! co2 ice cover
-  real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
-  real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
-  real,intent(out) :: mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
-  real,intent(out) :: mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
-  real,intent(out) :: mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
-!======================================================================
-!  Local variables:
-
-      real surffield(ngrid) ! to temporarily store a surface field
-      real xmin,xmax ! to display min and max of a field
-!
-      INTEGER ig,iq,lmax
-      INTEGER nid, nvarid
-      INTEGER ierr, i, nsrf
-!      integer isoil 
-!      INTEGER length
-!      PARAMETER (length=100)
-      CHARACTER*7 str7
-      CHARACTER*2 str2
-      CHARACTER*1 yes
-!
-      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
-      INTEGER nqold
-
-! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
-      logical :: oldtracernames=.false.
-      integer :: count
-      character(len=30) :: txt ! to store some text
-
-! specific for time 
-      REAL,ALLOCATABLE :: time(:) ! times stored in start
-      INTEGER timelen ! number of times stored in the file
-      INTEGER indextime ! index of selected time
-
-      INTEGER :: edges(3),corner(3)
-      LOGICAL :: found
-
-      REAL :: timestart ! to pick which initial state to start from
-
-! open physics initial state file:
-call open_startphy(fichnom)
-
-
-! possibility to modify tab_cntrl in tabfi
-write(*,*)
-write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
-call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
-            p_omeg,p_g,p_mugaz,p_daysec,time0)
-
-
-! Load surface geopotential:
-call get_field("phisfi",phisfi,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <phisfi>"
-  call abort
-else
-  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
-             minval(phisfi), maxval(phisfi)
-endif
-
-
-! Load bare ground albedo:
-call get_field("albedodat",albedodat,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <albedodat>"
-  call abort
-else
-  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
-             minval(albedodat), maxval(albedodat)
-endif
-
-! ZMEA
-call get_field("ZMEA",zmea,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZMEA>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZMEA> range:", &
-             minval(zmea), maxval(zmea)
-endif
-
-
-! ZSTD
-call get_field("ZSTD",zstd,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZSTD>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZSTD> range:", &
-             minval(zstd), maxval(zstd)
-endif
-
-
-! ZSIG
-call get_field("ZSIG",zsig,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZSIG>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZSIG> range:", &
-             minval(zsig), maxval(zsig)
-endif
-
-
-! ZGAM
-call get_field("ZGAM",zgam,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZGAM>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZGAM> range:", &
-             minval(zgam), maxval(zgam)
-endif
-
-
-! ZTHE
-call get_field("ZTHE",zthe,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <ZTHE>"
-  call abort
-else
-  write(*,*) "phyetat0: <ZTHE> range:", &
-             minval(zthe), maxval(zthe)
-endif
-
-      
-! Time axis
-! obtain timestart from run.def
-timestart=-9999 ! default value
-call getin("timestart",timestart)
-
-found=inquire_dimension("Time")
-if (.not.found) then
-  indextime = 1
-  write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
-else
-  write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
-  timelen=inquire_dimension_length("Time")
-  allocate(time(timelen))
-  ! load "Time" array:
-  call get_var("Time",time,found)
-  if (.not.found) then
-    write(*,*) "phyetat0: Failed loading <Time>"
-    call abort
-  endif
-  ! seclect the desired time index
-  IF (timestart .lt. 0) THEN  ! default: we use the last time value
-    indextime = timelen
-  ELSE  ! else we look for the desired value in the time axis
-    indextime = 0
-    DO i=1,timelen
-      IF (abs(time(i) - timestart) .lt. 0.01) THEN
-        indextime = i
-        EXIT
-      ENDIF
-    ENDDO
-    IF (indextime .eq. 0) THEN
-      PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
-      PRINT*, "Stored times are:"
-      DO i=1,timelen
-         PRINT*, time(i)
-      ENDDO
-      CALL abort
-    ENDIF
-  ENDIF ! of IF (timestart .lt. 0)
-  ! In startfi the absolute date is day_ini + time0 + time
-  ! For now on, in the GCM physics, it is day_ini + time0
-  time0 = time(indextime) + time0
-  day_ini = day_ini + INT(time0)
-  time0 = time0 - INT(time0)
-        
-  PRINT*, "phyetat0: Selected time ",time(indextime), &
-          " at index ",indextime
-     
-  DEALLOCATE(time)
-endif ! of if Time not found in file
-
-
-! CO2 ice cover
-call get_field("co2ice",co2ice,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <co2ice>"
-  call abort
-else
-  write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
-             minval(co2ice), maxval(co2ice)
-endif
-
-! Memory of the origin of the co2 particles
-call get_field("mem_Mccn_co2",mem_Mccn_co2,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: <mem_Mccn_co2> not in file"
-  mem_Mccn_co2(:,:)=0
-else
-  write(*,*) "phyetat0: Memory of CCN mass of H2O and dust used by CO2"
-  write(*,*) " <mem_Mccn_co2> range:", &
-             minval(mem_Mccn_co2), maxval(mem_Mccn_co2)
-endif
-
-call get_field("mem_Nccn_co2",mem_Nccn_co2,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: <mem_Nccn_co2> not in file"
-  mem_Nccn_co2(:,:)=0
-else
-  write(*,*) "phyetat0: Memory of CCN number of H2O and dust used by CO2"
-  write(*,*) " <mem_Nccn_co2> range:", &
-             minval(mem_Nccn_co2), maxval(mem_Nccn_co2)
-endif
-
-call get_field("mem_Mh2o_co2",mem_Mh2o_co2,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: <mem_Mh2o_co2> not in file"
-  mem_Mh2o_co2(:,:)=0
-else
-  write(*,*) "phyetat0: Memory of H2O mass integred into CO2 crystal"
-  write(*,*) " <mem_Mh2o_co2> range:", &
-             minval(mem_Mh2o_co2), maxval(mem_Mh2o_co2)
-endif
-
-! Dust conversion factor
-call get_field("tauscaling",tauscaling,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: <tauscaling> not in file"
-  tauscaling(:) = -1
-else
-  write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
-             minval(tauscaling), maxval(tauscaling)
-endif
-
-! Sub-grid cloud fraction
-call get_field("totcloudfrac",totcloudfrac,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
-  totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
-else
-  write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
-             minval(totcloudfrac), maxval(totcloudfrac)
-endif
-
-! Surface temperature :
-call get_field("tsurf",tsurf,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <tsurf>"
-  call abort
-else
-  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
-             minval(tsurf), maxval(tsurf)
-endif
-
-! Surface emissivity
-call get_field("emis",emis,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <emis>"
-  call abort
-else
-  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
-             minval(emis), maxval(emis)
-endif
-
-
-! surface roughness length (NB: z0 is a common in surfdat_h)
-call get_field("z0",z0,found)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <z0>"
-  write(*,*) 'will use constant value of z0_default:',z0_default
-  z0(:)=z0_default
-else
-  write(*,*) "phyetat0: Surface roughness <z0> range:", &
-             minval(z0), maxval(z0)
-endif
-
-
-! pbl wind variance
-call get_field("q2",q2,found,indextime)
-if (.not.found) then
-  write(*,*) "phyetat0: Failed loading <q2>"
-  call abort
-else
-  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
-             minval(q2), maxval(q2)
-endif
-
-
-! tracer on surface
-if (nq.ge.1) then
-  do iq=1,nq
-    txt=noms(iq)
-    if (txt.eq."h2o_vap") then
-      ! There is no surface tracer for h2o_vap;
-      ! "h2o_ice" should be loaded instead
-      txt="h2o_ice"
-      write(*,*) 'phyetat0: loading surface tracer', &
-                           ' h2o_ice instead of h2o_vap'
-    endif
-    call get_field(txt,qsurf(:,iq),found,indextime)
-    if (.not.found) then
-      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
-      write(*,*) "         ",trim(txt)," is set to zero"
-    else
-      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
-                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
-    endif
-  enddo
-endif ! of if (nq.ge.1)
-
-! Call to soil_settings, in order to read soil temperatures,
-! as well as thermal inertia and volumetric heat capacity
-
-      call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
-
-!
-! close file:
-!
-      call close_startphy
-
-      END
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 1944)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 1944)
@@ -0,0 +1,393 @@
+module phyetat0_mod
+
+implicit none
+
+contains
+
+subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
+                     day_ini,time0,tsurf,tsoil,albedo,emis,q2,qsurf,co2ice, &
+                     tauscaling,totcloudfrac,wstar,mem_Mccn_co2,mem_Nccn_co2,&
+                     mem_Mh2o_co2)
+
+  use tracer_mod, only: noms ! tracer names
+  use surfdat_h, only: phisfi, albedodat, z0, z0_default,&
+                       zmea, zstd, zsig, zgam, zthe
+  use iostart, only: nid_start, open_startphy, close_startphy, &
+                     get_field, get_var, inquire_field, &
+                     inquire_dimension, inquire_dimension_length
+  use ioipsl_getincom, only : getin
+
+  implicit none
+!======================================================================
+! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+!  Adaptation à Mars : Yann Wanherdrick 
+! Objet: Lecture de l etat initial pour la physique
+! Modifs: Aug.2010 EM : use NetCDF90 to load variables (enables using
+!                      r4 or r8 restarts independently of having compiled
+!                      the GCM in r4 or r8)
+!         June 2013 TN : Possibility to read files with a time axis
+!         November 2013 EM : Enabeling parallel, using iostart module
+!======================================================================
+  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
+  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
+!======================================================================
+!  Arguments:
+!  ---------
+!  inputs:
+  character*(*),intent(in) :: fichnom ! "startfi.nc" file
+  integer,intent(in) :: tab0
+  integer,intent(in) :: Lmodif
+  integer,intent(in) :: nsoil ! # of soil layers
+  integer,intent(in) :: ngrid ! # of atmospheric columns
+  integer,intent(in) :: nlay ! # of atmospheric layers
+  integer,intent(in) :: nq
+  integer :: day_ini
+  real :: time0
+
+!  outputs:
+  real,intent(out) :: tsurf(ngrid) ! surface temperature
+  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
+  real,intent(out) :: albedo(ngrid,2) ! surface albedo
+  real,intent(out) :: emis(ngrid) ! surface emissivity
+  real,intent(out) :: q2(ngrid,nlay+1) ! 
+  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
+  real,intent(out) :: co2ice(ngrid) ! co2 ice cover
+  real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
+  real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
+  real,intent(out) :: wstar(ngrid) ! Max vertical velocity in thermals (m/s)
+  real,intent(out) :: mem_Mccn_co2(ngrid,nlay) ! Memory of CCN mass of H2O and dust used by CO2
+  real,intent(out) :: mem_Nccn_co2(ngrid,nlay) ! Memory of CCN number of H2O and dust used by CO2
+  real,intent(out) :: mem_Mh2o_co2(ngrid,nlay) ! Memory of H2O mass integred into CO2 crystal
+!======================================================================
+!  Local variables:
+
+      real surffield(ngrid) ! to temporarily store a surface field
+      real xmin,xmax ! to display min and max of a field
+!
+      INTEGER ig,iq,lmax
+      INTEGER nid, nvarid
+      INTEGER ierr, i, nsrf
+!      integer isoil 
+!      INTEGER length
+!      PARAMETER (length=100)
+      CHARACTER*7 str7
+      CHARACTER*2 str2
+      CHARACTER*1 yes
+!
+      REAL p_rad,p_omeg,p_g,p_mugaz,p_daysec
+      INTEGER nqold
+
+! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
+      logical :: oldtracernames=.false.
+      integer :: count
+      character(len=30) :: txt ! to store some text
+
+! specific for time 
+      REAL,ALLOCATABLE :: time(:) ! times stored in start
+      INTEGER timelen ! number of times stored in the file
+      INTEGER indextime ! index of selected time
+
+      INTEGER :: edges(3),corner(3)
+      LOGICAL :: found
+
+      REAL :: timestart ! to pick which initial state to start from
+
+! open physics initial state file:
+call open_startphy(fichnom)
+
+
+! possibility to modify tab_cntrl in tabfi
+write(*,*)
+write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
+call tabfi (nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
+            p_omeg,p_g,p_mugaz,p_daysec,time0)
+
+
+! Load surface geopotential:
+call get_field("phisfi",phisfi,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <phisfi>"
+  call abort
+else
+  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
+             minval(phisfi), maxval(phisfi)
+endif
+
+
+! Load bare ground albedo:
+call get_field("albedodat",albedodat,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <albedodat>"
+  call abort
+else
+  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
+             minval(albedodat), maxval(albedodat)
+endif
+
+! ZMEA
+call get_field("ZMEA",zmea,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZMEA>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZMEA> range:", &
+             minval(zmea), maxval(zmea)
+endif
+
+
+! ZSTD
+call get_field("ZSTD",zstd,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZSTD>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZSTD> range:", &
+             minval(zstd), maxval(zstd)
+endif
+
+
+! ZSIG
+call get_field("ZSIG",zsig,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZSIG>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZSIG> range:", &
+             minval(zsig), maxval(zsig)
+endif
+
+
+! ZGAM
+call get_field("ZGAM",zgam,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZGAM>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZGAM> range:", &
+             minval(zgam), maxval(zgam)
+endif
+
+
+! ZTHE
+call get_field("ZTHE",zthe,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZTHE>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZTHE> range:", &
+             minval(zthe), maxval(zthe)
+endif
+
+      
+! Time axis
+! obtain timestart from run.def
+timestart=-9999 ! default value
+call getin("timestart",timestart)
+
+found=inquire_dimension("Time")
+if (.not.found) then
+  indextime = 1
+  write(*,*) "phyetat0: No time axis found in "//trim(fichnom)
+else
+  write(*,*) "phyetat0: Time axis found in "//trim(fichnom)
+  timelen=inquire_dimension_length("Time")
+  allocate(time(timelen))
+  ! load "Time" array:
+  call get_var("Time",time,found)
+  if (.not.found) then
+    write(*,*) "phyetat0: Failed loading <Time>"
+    call abort
+  endif
+  ! seclect the desired time index
+  IF (timestart .lt. 0) THEN  ! default: we use the last time value
+    indextime = timelen
+  ELSE  ! else we look for the desired value in the time axis
+    indextime = 0
+    DO i=1,timelen
+      IF (abs(time(i) - timestart) .lt. 0.01) THEN
+        indextime = i
+        EXIT
+      ENDIF
+    ENDDO
+    IF (indextime .eq. 0) THEN
+      PRINT*, "Time", timestart," is not in "//trim(fichnom)//"!!"
+      PRINT*, "Stored times are:"
+      DO i=1,timelen
+         PRINT*, time(i)
+      ENDDO
+      CALL abort
+    ENDIF
+  ENDIF ! of IF (timestart .lt. 0)
+  ! In startfi the absolute date is day_ini + time0 + time
+  ! For now on, in the GCM physics, it is day_ini + time0
+  time0 = time(indextime) + time0
+  day_ini = day_ini + INT(time0)
+  time0 = time0 - INT(time0)
+        
+  PRINT*, "phyetat0: Selected time ",time(indextime), &
+          " at index ",indextime
+     
+  DEALLOCATE(time)
+endif ! of if Time not found in file
+
+
+! CO2 ice cover
+call get_field("co2ice",co2ice,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <co2ice>"
+  call abort
+else
+  write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
+             minval(co2ice), maxval(co2ice)
+endif
+
+! Memory of the origin of the co2 particles
+call get_field("mem_Mccn_co2",mem_Mccn_co2,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <mem_Mccn_co2> not in file"
+  mem_Mccn_co2(:,:)=0
+else
+  write(*,*) "phyetat0: Memory of CCN mass of H2O and dust used by CO2"
+  write(*,*) " <mem_Mccn_co2> range:", &
+             minval(mem_Mccn_co2), maxval(mem_Mccn_co2)
+endif
+
+call get_field("mem_Nccn_co2",mem_Nccn_co2,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <mem_Nccn_co2> not in file"
+  mem_Nccn_co2(:,:)=0
+else
+  write(*,*) "phyetat0: Memory of CCN number of H2O and dust used by CO2"
+  write(*,*) " <mem_Nccn_co2> range:", &
+             minval(mem_Nccn_co2), maxval(mem_Nccn_co2)
+endif
+
+call get_field("mem_Mh2o_co2",mem_Mh2o_co2,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <mem_Mh2o_co2> not in file"
+  mem_Mh2o_co2(:,:)=0
+else
+  write(*,*) "phyetat0: Memory of H2O mass integred into CO2 crystal"
+  write(*,*) " <mem_Mh2o_co2> range:", &
+             minval(mem_Mh2o_co2), maxval(mem_Mh2o_co2)
+endif
+
+! Dust conversion factor
+call get_field("tauscaling",tauscaling,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <tauscaling> not in file"
+  tauscaling(:) = -1
+else
+  write(*,*) "phyetat0: dust conversion factor <tauscaling> range:", &
+             minval(tauscaling), maxval(tauscaling)
+endif
+
+! Sub-grid cloud fraction
+call get_field("totcloudfrac",totcloudfrac,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
+  totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
+else
+  write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
+             minval(totcloudfrac), maxval(totcloudfrac)
+endif
+
+! Max vertical velocity in thermals
+call get_field("wstar",wstar,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <totcloudfrac> not in file! Set to zero"
+  wstar(:)=0
+else
+  write(*,*) "phyetat0: Max vertical velocity in thermals <wstar> range:", &
+             minval(wstar),maxval(wstar)
+endif
+
+! Surface temperature :
+call get_field("tsurf",tsurf,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <tsurf>"
+  call abort
+else
+  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
+             minval(tsurf), maxval(tsurf)
+endif
+
+! Surface albedo
+call get_field("albedo",albedo(:,1),found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <albedo>"
+  albedo(:,1)=albedodat(:)
+else
+  write(*,*) "phyetat0: Surface albedo <albedo> range:", &
+             minval(albedo(:,1)), maxval(albedo(:,1))
+endif
+albedo(:,2)=albedo(:,1)
+
+! Surface emissivity
+call get_field("emis",emis,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <emis>"
+  call abort
+else
+  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
+             minval(emis), maxval(emis)
+endif
+
+
+! surface roughness length (NB: z0 is a common in surfdat_h)
+call get_field("z0",z0,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <z0>"
+  write(*,*) 'will use constant value of z0_default:',z0_default
+  z0(:)=z0_default
+else
+  write(*,*) "phyetat0: Surface roughness <z0> range:", &
+             minval(z0), maxval(z0)
+endif
+
+
+! pbl wind variance
+call get_field("q2",q2,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <q2>"
+  call abort
+else
+  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
+             minval(q2), maxval(q2)
+endif
+
+
+! tracer on surface
+if (nq.ge.1) then
+  do iq=1,nq
+    txt=noms(iq)
+    if (txt.eq."h2o_vap") then
+      ! There is no surface tracer for h2o_vap;
+      ! "h2o_ice" should be loaded instead
+      txt="h2o_ice"
+      write(*,*) 'phyetat0: loading surface tracer', &
+                           ' h2o_ice instead of h2o_vap'
+    endif
+    call get_field(txt,qsurf(:,iq),found,indextime)
+    if (.not.found) then
+      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+      write(*,*) "         ",trim(txt)," is set to zero"
+    else
+      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
+                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+    endif
+  enddo
+endif ! of if (nq.ge.1)
+
+! Call to soil_settings, in order to read soil temperatures,
+! as well as thermal inertia and volumetric heat capacity
+call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
+
+!
+! close file:
+!
+call close_startphy
+
+end subroutine phyetat0
+
+end module phyetat0_mod
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1941)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1944)
@@ -144,7 +144,7 @@
 
 subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
-                    phystep,time,tsurf,tsoil,co2ice,emis,q2,qsurf,&
-                    tauscaling,totcloudfrac,mem_Mccn_co2,mem_Nccn_co2,&
-                    mem_Mh2o_co2)
+                    phystep,time,tsurf,tsoil,co2ice,albedo,emis,q2,qsurf,&
+                    tauscaling,totcloudfrac,wstar, &
+                    mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
   ! write time-dependent variable to restart file
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -166,4 +166,5 @@
   real,intent(in) :: tsoil(ngrid,nsoil)
   real,intent(in) :: co2ice(ngrid)
+  real,intent(in) :: albedo(ngrid,2)
   real,intent(in) :: emis(ngrid)
   real,intent(in) :: q2(ngrid,nlay+1)
@@ -171,4 +172,5 @@
   real,intent(in) :: tauscaling(ngrid)
   real,intent(in) :: totcloudfrac(ngrid)
+  real,intent(in) :: wstar(ngrid)
   real,intent(in) :: mem_Mccn_co2(ngrid,nlay) ! CCN mass of H2O and dust used by CO2
   real,intent(in) :: mem_Nccn_co2(ngrid,nlay) ! CCN number of H2O and dust used by CO2
@@ -197,4 +199,7 @@
   ! Soil temperature
   call put_field("tsoil","Soil temperature",tsoil,time)
+  
+  ! Albedo of the surface
+  call put_field("albedo","Surface albedo",albedo(:,1),time)
   
   ! Emissivity of the surface
@@ -211,4 +216,8 @@
   ! Or by any user who wants to reconstruct dust, opacity from the start files.
   call put_field("tauscaling","dust conversion factor",tauscaling,time)
+
+  if (calltherm) then
+    call put_field("wstar","Max vertical velocity in thermals",wstar,time)
+  endif
 
   ! Tracers on the surface
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1941)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1944)
@@ -59,4 +59,5 @@
       use planetwide_mod
       use phyredem, only: physdem0, physdem1
+      use phyetat0_mod, only: phyetat0
       use eofdump_mod, only: eofdump
       USE vertical_layers_mod, ONLY: ap,bp,aps,bps
@@ -435,6 +436,7 @@
      &         nsoilmx,ngrid,nlayer,nq,
      &         day_ini,time_phys,
-     &         tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling,
-     &         totcloudfrac,mem_Mccn_co2,mem_Nccn_co2,
+     &         tsurf,tsoil,albedo,emis,
+     &         q2,qsurf,co2ice,tauscaling,totcloudfrac,wstar,
+     &         mem_Mccn_co2,mem_Nccn_co2,
      &         mem_Mh2o_co2)
 
@@ -476,5 +478,5 @@
 c        Initialize albedo and orbital calculation
 c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         CALL surfini(ngrid,co2ice,qsurf,albedo)
+         CALL surfini(ngrid,co2ice,qsurf)
          CALL iniorbit(aphelie,periheli,year_day,peri_day,obliquit)
 
@@ -1844,7 +1846,7 @@
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
-     .                tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
-     .                totcloudfrac,mem_Mccn_co2,mem_Nccn_co2,
-     .                mem_Mh2o_co2)
+     .                tsurf,tsoil,co2ice,albedo,emis,
+     .                q2,qsurf,tauscaling,totcloudfrac,wstar,
+     .                mem_Mccn_co2,mem_Nccn_co2,mem_Mh2o_co2)
           
          ENDIF
Index: trunk/LMDZ.MARS/libf/phymars/surfini.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 1941)
+++ trunk/LMDZ.MARS/libf/phymars/surfini.F	(revision 1944)
@@ -1,3 +1,3 @@
-      SUBROUTINE surfini(ngrid,piceco2,qsurf,psolaralb)
+      SUBROUTINE surfini(ngrid,piceco2,qsurf)
 
       USE ioipsl_getincom, only: getin
@@ -12,5 +12,5 @@
       use mod_phys_lmdz_para, only : is_master, gather, scatter
 #endif
-      USE comcstfi_h
+      USE comcstfi_h, ONLY: pi
       use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
       use datafile_mod, only: datadir
@@ -29,5 +29,4 @@
       real,intent(in) :: piceco2(ngrid) ! CO2 ice thickness
       real,intent(inout) :: qsurf(ngrid,nqmx) ! tracer on surface (kg/m2)
-      real,intent(out) :: psolaralb(ngrid,2) ! albedo
 
       INTEGER ig,icap,iq,alternate
@@ -458,106 +457,9 @@
        call scatter(watercaptag_glo,watercaptag)
        
-! end of #else of #ifdef MESOSCALE
-       ENDIF ! (caps & water)
-#endif       
-
-c ===============================================================
-c      INITIAL ALBEDO
-c ===============================================================
-
-         write(*,*)"surfini: water frost thickness",
-     s     frost_albedo_threshold
-         write(*,*)"surfini: water ice albedo:", albedo_h2o_ice
-         write(*,*)"surfini: water ice TI:", inert_h2o_ice
-
-c        To start with : Initial albedo = observed dataset 
-c        -------------------------------------------------
-         DO ig=1,ngrid
-              psolaralb(ig,1)=albedodat(ig)
-              psolaralb(ig,2)=albedodat(ig)
-         END DO
-         PRINT*,'surfini: minimum albedo without water caps',
-     &          minval(albedodat)
-         PRINT*,'surfini: maximum albedo without water caps',
-     &          maxval(albedodat)
-
-c        initial albedo if permanent H2O ice is present 
-c        ------------------------------------------------
-         IF ((water) .and. (caps)) THEN
-           DO ig=1,ngrid
-            IF (watercaptag(ig)) THEN
-              psolaralb(ig,1) = albedo_h2o_ice
-              psolaralb(ig,2) = albedo_h2o_ice
-            ENDIF
-           END DO
-           PRINT*,'surfini: minimum albedo with water caps',
-     &            minval(albedodat)
-           PRINT*,'surfini: maximum albedo with water caps',
-     &            maxval(albedodat)
-
-         ENDIF
-
-c      changing initial albedo if CO2 ice is present 
-c      -------------------------------------------
-
-       DO ig=1,ngrid
-         IF (piceco2(ig) .GT. 0.) THEN
-             IF(latitude(ig).LT. 0.) THEN
-                icap=2 ! Southern hemisphere
-             ELSE
-                icap=1 ! Northern hemisphere
-             ENDIF
-             psolaralb(ig,1) = albedice(icap)
-             psolaralb(ig,2) = albedice(icap)    
-         END IF
-       END DO
-
-c      changing initial albedo if water ice frost is present
-c      -------------------------------------------
-       IF (water) THEN
-          do iq=1,nqmx
-c          if there is frost and surface albedo is set to albedo_h2o_ice
-           if(noms(iq).eq."h2o_ice") then
-             do ig=1,ngrid
-             
-              if ((watercaptag(ig).eqv..false.) 
-     &     .and. (qsurf(ig,iq).lt.-frost_albedo_threshold)) then
-              print*, ''
-              print*, '!!! PROBLEM in SURFINI !!!!'
-              print*, 'FOUND NEGATIVE SURFACE ICE VALUE WHERE
-     & WATERCAPTAG IS FALSE'
-              print*, ''
-              print*, 'ig,qsurf,threshold' , 
-     &         ig, qsurf(ig,iq), -frost_albedo_threshold
-              print*, ''
-              print*, '1) Check h2o_ice in startfi and ice
-     & distribution in surfini'
-              print*, '2) Use ini_h2osurf option in newstart'
-              print*, ''
-#ifndef MESOINI
-              CALL ABORT
-#else
-              qsurf(ig,iq) = 0.
-#endif
-            endif
-             
-              if ((piceco2(ig) .eq. 0.).and.
-     &          (qsurf(ig,iq).gt.frost_albedo_threshold)) then
-                     psolaralb(ig,1) = albedo_h2o_ice
-                     psolaralb(ig,2) = albedo_h2o_ice
-c                     PRINT*,'surfini.F frost',
-c     &                  latitude(ig)*180./pi, longitude(ig)*180./pi
-               endif 
-              enddo
-           endif 
-          end do
-          PRINT*,'surfini: minimum albedo with frost and co2',
-     &            minval(albedodat)
-          PRINT*,'surfini: maximum albedo with frost and co2',
-     &            maxval(albedodat)
        ELSE
          watercaptag(:) = .false.
-       END IF ! OF IF(water)
-         
-      RETURN
+       ENDIF ! (caps & water)
+! end of #else of #ifndef MESOSCALE
+#endif       
+
       END
Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 1941)
+++ trunk/LMDZ.MARS/libf/phymars/vdifc.F	(revision 1944)
@@ -36,6 +36,6 @@
 c   -------------
 
-#include "callkeys.h"
-#include "microphys.h"
+      include "callkeys.h"
+      include "microphys.h"
 
 c
@@ -66,5 +66,5 @@
 c    Argument added to account for subgrid gustiness :
 
-      REAL wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
+      REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid)
 
 c    Traceurs :
