Index: trunk/LMDZ.MARS/libf/phymars/iostart.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/iostart.F90	(revision 2911)
+++ trunk/LMDZ.MARS/libf/phymars/iostart.F90	(revision 2913)
@@ -15,5 +15,5 @@
     INTEGER,SAVE :: idim7 ! "Time" dimension
     INTEGER,SAVE :: idim8 ! "nslope" dimension
-    INTEGER,SAVE :: idim9 ! "inter slope" dimension
+    INTEGER,SAVE :: idim9 ! "nslope_plus_1" dimension
     INTEGER,SAVE :: timeindex ! current time index (for time-dependent fields)
     INTEGER,PARAMETER :: length=100 ! size of tab_cntrl array
@@ -564,5 +564,5 @@
       ENDIF
 
-      ierr=NF90_DEF_DIM(nid_restart,"inter slope",nslope+1,idim9)
+      ierr=NF90_DEF_DIM(nid_restart,"inter_slope",nslope+1,idim9)
       IF (ierr/=NF90_NOERR) THEN
         write(*,*)'phyredem: problem defining inter slope dimension'
@@ -660,4 +660,5 @@
   
   REAL                           :: field_glo(klon_glo,field_size)
+  REAL                           :: field_glo_reshape(klon_glo,nsoilmx,nslope,timeindex)
   INTEGER                        :: ierr
   INTEGER                        :: nvarid
@@ -804,5 +805,5 @@
         endif ! of if (.not.present(timeindex))
 
-      ELSE IF (field_size==nsoilmx) THEN
+      ELSE IF (field_size==nsoilmx .or. field_size==nsoilmx*nslope) THEN
         ! input is a 2D "subsurface field" array
         if (.not.present(time)) then ! for a time-independent field
@@ -828,9 +829,19 @@
             ierr=NF90_REDEF(nid_restart)
 #ifdef NC_DOUBLE
+            if(field_name.eq. "tsoil") then
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
+                              (/idim2,idim3,idim8,idim7/),nvarid)
+            else
             ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_DOUBLE,&
                               (/idim2,idim3,idim7/),nvarid)
-#else
+            endif
+#else
+            if(field_name.eq. "tsoil") then
+            ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
+                              (/idim2,idim3,idim8,idim7/),nvarid)
+            else
             ierr=NF90_DEF_VAR(nid_restart,field_name,NF90_FLOAT,&
                               (/idim2,idim3,idim7/),nvarid)
+            endif
 #endif
            if (ierr.ne.NF90_NOERR) then
@@ -842,7 +853,14 @@
           endif
           ! Write the variable
+          if(field_name.eq. "tsoil") then
+
+          field_glo_reshape=RESHAPE(field_glo, (/klon_glo, nsoilmx, nslope, timeindex/))
+
+          ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo_reshape,&
+                            start=(/1,1,1,timeindex/))
+          else
           ierr=NF90_PUT_VAR(nid_restart,nvarid,field_glo,&
                             start=(/1,1,timeindex/))
-
+          endif
         endif ! of if (.not.present(time))
 
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2911)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 2913)
@@ -60,14 +60,14 @@
 
 !  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) :: tsurf(ngrid,nslope) ! surface temperature
+  real,intent(out) :: tsoil(ngrid,nsoil,nslope) ! soil temperature
+  real,intent(out) :: albedo(ngrid,2,nslope) ! surface albedo
+  real,intent(out) :: emis(ngrid,nslope) ! surface emissivity
   real,intent(out) :: q2(ngrid,nlay+1) ! 
-  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
+  real,intent(out) :: qsurf(ngrid,nq,nslope) ! tracers on surface
   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) :: watercap(ngrid) ! h2o_ice_cover
+  real,intent(out) :: watercap(ngrid,nslope) ! h2o_ice_cover
   real, intent(out) :: def_slope(nslope+1) !boundaries for bining of the slopes
   real, intent(out) :: def_slope_mean(nslope)
@@ -143,7 +143,8 @@
       write(*,*)'We take default def_slope and subslope dist'
       allocate(default_def_slope(nslope+1))
-      default_def_slope(1) = 0.
-      default_def_slope(2) = 0.
+      default_def_slope(1) = -90.
+      default_def_slope(2) = 90.
       subslope_dist(:,:)=1.
+      def_slope(:)=default_def_slope(:)
      else
        write(*,*)'Variable def_slope is not present in the start and'
@@ -168,4 +169,5 @@
       default_def_slope(2) = 0.
       subslope_dist(:,:)=1.
+      def_slope(:)=default_def_slope(:)
     else
      write(*,*)'Without starfi, lets first run with nslope=1'
@@ -516,5 +518,5 @@
    endif
 else
-  tsurf(:)=0. ! will be updated afterwards in physiq !
+  tsurf(:,:)=0. ! will be updated afterwards in physiq !
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
@@ -523,15 +525,19 @@
 ! Surface albedo
 if (startphy_file) then
-   call get_field("albedo",albedo(:,1),found,indextime)
+   call get_field("albedo",albedo(:,1,:),found,indextime)
    if (.not.found) then
      write(*,*) "phyetat0: Failed loading <albedo>"
-     albedo(:,1)=albedodat(:)
-   endif
-else
-   albedo(:,1)=albedodat(:)
+     do islope=1,nslope
+       albedo(:,1,islope)=albedodat(:)
+     enddo
+   endif
+else
+   do islope=1,nslope
+     albedo(:,1,islope)=albedodat(:)
+   enddo
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface albedo <albedo> range:", &
-            minval(albedo(:,1)), maxval(albedo(:,1))
-albedo(:,2)=albedo(:,1)
+            minval(albedo(:,1,:)), maxval(albedo(:,1,:))
+albedo(:,2,:)=albedo(:,1,:)
 
 ! Surface emissivity
@@ -547,5 +553,5 @@
   call getin_p("surfemis_without_startfi",surfemis)
   print*,"surfemis_without_startfi",surfemis
-  emis(:)=surfemis
+  emis(:,:)=surfemis
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface emissivity <emis> range:", &
@@ -633,8 +639,8 @@
         ! We first check if co2ice exist in the startfi file (old way)
         ! CO2 ice cover
-        call get_field("co2ice",qsurf(:,iq),found,indextime)
+        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
         ! If not, we load the corresponding tracer in qsurf
         if (.not.found) then
-          call get_field(txt,qsurf(:,iq),found,indextime)
+          call get_field(txt,qsurf(:,iq,:),found,indextime)
           if (.not.found) then
             call abort_physic(modname, &
@@ -643,16 +649,16 @@
         endif
       else ! (not the co2 tracer)
-        call get_field(txt,qsurf(:,iq),found,indextime)
+        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"
-          qsurf(:,iq)=0.
+          qsurf(:,iq,:)=0.
         endif
       endif !endif co2
     else !(not startphy_file)
-      qsurf(:,iq)=0.
+      qsurf(:,iq,:)=0.
     endif ! of if (startphy_file)
     write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
-                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+                 minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
   enddo ! of do iq=1,nq
 
@@ -661,8 +667,8 @@
       ! CO2 ice cover
       if (startphy_file) then
-        call get_field("co2ice",qsurf(:,iq),found,indextime)
+        call get_field("co2ice",qsurf(:,iq,:),found,indextime)
       ! If not, we load the corresponding tracer in qsurf
         if (.not.found) then
-          call get_field(txt,qsurf(:,iq),found,indextime)
+          call get_field(txt,qsurf(:,iq,:),found,indextime)
           if (.not.found) then
             call abort_physic(modname, &
@@ -672,8 +678,8 @@
       else
        ! If we run without startfile, co2ice is set to 0
-        qsurf(:,iq)=0.
+        qsurf(:,iq,:)=0.
       endif !if (startphy_file)
         write(*,*) "phyetat0: CO2 ice cover <co2ice> range:", &
-            minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+            minval(qsurf(:,iq,:)), maxval(qsurf(:,iq,:))
     endif
 
@@ -685,5 +691,5 @@
      write(*,*) "phyetat0: Failed loading <watercap> : ", &
                           "<watercap> is set to zero"
-     watercap(:)=0
+     watercap(:,:)=0
      write(*,*) 'Now transfer negative surface water ice to', &
                 ' watercap !'
@@ -693,8 +699,10 @@
        if (txt.eq."h2o_ice") then
          do ig=1,ngrid
-          if (qsurf(ig,iq).lt.0.0) then
-             watercap(ig) = qsurf(ig,iq)
-             qsurf(ig,iq) = 0.0
-          end if
+          do islope=1,nslope
+            if (qsurf(ig,iq,islope).lt.0.0) then
+              watercap(ig,islope) = qsurf(ig,iq,islope)
+              qsurf(ig,iq,islope) = 0.0
+            end if
+          enddo
          end do
        endif
@@ -702,7 +710,9 @@
        if (txt.eq."hdo_ice") then
          do ig=1,ngrid
-          if (qsurf(ig,iq).lt.0.0) then
-             qsurf(ig,iq) = 0.0
-          end if
+          do islope=1,nslope
+            if (qsurf(ig,iq,islope).lt.0.0) then
+              qsurf(ig,iq,islope) = 0.0
+            end if
+          enddo
          end do
        endif
@@ -712,5 +722,5 @@
    endif ! of if (.not.found)
 else
-   watercap(:)=0
+   watercap(:,:)=0
 endif ! of if (startphy_file) 
 write(*,*) "phyetat0: Surface water ice <watercap> range:", &
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2911)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 2913)
@@ -171,4 +171,5 @@
   use dust_param_mod, only: dustscaling_mode
   use comsoil_h,only: flux_geo
+  use comslope_mod, ONLY: nslope
   implicit none
   
@@ -182,14 +183,14 @@
   real,intent(in) :: phystep
   real,intent(in) :: time
-  real,intent(in) :: tsurf(ngrid)
-  real,intent(in) :: tsoil(ngrid,nsoil)
-  real,intent(in) :: albedo(ngrid,2)
-  real,intent(in) :: emis(ngrid)
+  real,intent(in) :: tsurf(ngrid,nslope)
+  real,intent(in) :: tsoil(ngrid,nsoil,nslope)
+  real,intent(in) :: albedo(ngrid,2,nslope)
+  real,intent(in) :: emis(ngrid,nslope)
   real,intent(in) :: q2(ngrid,nlay+1)
-  real,intent(in) :: qsurf(ngrid,nq)
+  real,intent(in) :: qsurf(ngrid,nq,nslope)
   real,intent(in) :: tauscaling(ngrid)
   real,intent(in) :: totcloudfrac(ngrid)
   real,intent(in) :: wstar(ngrid)
-  real,intent(in) :: watercap(ngrid)
+  real,intent(in) :: watercap(ngrid,nslope)
   
   integer :: iq
@@ -219,5 +220,5 @@
   
   ! Albedo of the surface
-  call put_field("albedo","Surface albedo",albedo(:,1),time)
+  call put_field("albedo","Surface albedo",albedo(:,1,:),time)
   
   ! Emissivity of the surface
@@ -312,5 +313,5 @@
       end if     
 
-      call put_field(trim(txt),"tracer on surface",qsurf(:,iq),time)
+      call put_field(trim(txt),"tracer on surface",qsurf(:,iq,:),time)
     enddo
   endif
Index: trunk/LMDZ.MARS/libf/phymars/soil_settings.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/soil_settings.F	(revision 2911)
+++ trunk/LMDZ.MARS/libf/phymars/soil_settings.F	(revision 2913)
@@ -5,4 +5,5 @@
       use iostart, only: inquire_field_ndims, get_var, get_field,
      &                   inquire_field, inquire_dimension_length
+      USE comslope_mod, ONLY: nslope
       implicit none
 
@@ -41,8 +42,8 @@
       integer,intent(in) :: ngrid	! # of horizontal grid points
       integer,intent(in) :: nsoil	! # of soil layers
-      real,intent(in) :: tsurf(ngrid)   ! surface temperature
+      real,intent(in) :: tsurf(ngrid,nslope)   ! surface temperature
       integer,intent(in) :: indextime	! position on time axis
 !  output:
-      real,intent(out) :: tsoil(ngrid,nsoil)	! soil temperature
+      real,intent(out) :: tsoil(ngrid,nsoil,nslope)	! soil temperature
 
 !======================================================================
@@ -54,4 +55,5 @@
       integer ndims	! # of dimensions of read <inertiedat> data
       integer ig,iloop	! loop counters
+      integer islope
       
       integer edges(3),corner(3) ! to read a specific time
@@ -64,5 +66,5 @@
       real,dimension(:),allocatable :: oldmlayer
       real,dimension(:,:),allocatable :: oldinertiedat
-      real,dimension(:,:),allocatable :: oldtsoil
+      real,dimension(:,:,:),allocatable :: oldtsoil
       
       ! for interpolation
@@ -317,10 +319,12 @@
 	write(*,*)' => Building <tsoil> from surface values <tsurf>'
 	do iloop=1,nsoil
-	  tsoil(:,iloop)=tsurf(:)
+          do islope=1,nslope
+	    tsoil(:,iloop,islope)=tsurf(:,islope)
+          enddo
 	enddo
       else ! <tsoil> found
        if (interpol) then ! put values in oldtsoil
          if (.not.allocated(oldtsoil)) then
-           allocate(oldtsoil(ngrid,dimlen),stat=ierr)
+           allocate(oldtsoil(ngrid,dimlen,nslope),stat=ierr)
            if (ierr.ne.0) then
              write(*,*) 'soil_settings: failed allocation of ',
@@ -376,6 +380,7 @@
         do ig=1,ngrid
 	  ! copy values
-	  oldval(1)=tsurf(ig)
-	  oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
+          do islope=1,nslope
+	  oldval(1)=tsurf(ig,islope)
+	  oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen,islope)
 	  ! build vertical coordinate
 	  oldgrid(1)=0. ! ground
@@ -385,5 +390,6 @@
 	  call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
 	  ! copy result in tsoil
-	  tsoil(ig,:)=newval(:)
+	  tsoil(ig,:,islope)=newval(:)
+          enddo
 	enddo
 
@@ -419,11 +425,13 @@
 	oldgrid(2:dimlen+1)=oldmlayer(1:dimlen)
         do ig=1,ngrid
+          do islope=1,nslope
 	  ! data
-	  oldval(1)=tsurf(ig)
-	  oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen)
+	  oldval(1)=tsurf(ig,islope)
+	  oldval(2:dimlen+1)=oldtsoil(ig,1:dimlen,islope)
 	  ! interpolate
 	  call interp_line(oldgrid,oldval,dimlen+1,mlayer,newval,nsoil)
 	  ! copy result in inertiedat
-	  tsoil(ig,:)=newval(:)
+	  tsoil(ig,:,islope)=newval(:)
+          enddo
 	enddo
 	
