Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2899)
+++ /trunk/LMDZ.MARS/README	(revision 2900)
@@ -3880,5 +3880,27 @@
 Minor edits then (+ svn update with RV had some issues, so there are some "artefact changes" ...)
 
+== 14/02/2023 == RV
+First modifications to introduce subslope parametrization in surface processes.
+New variables are added in the new module comslope_mod to describe the subslopes: 
+	_ nslope (number of subslope)
+	_ subslope_dist (distribution  of the slopes)
+	_ def_slope (bound of the slope statistics / repartitions) etc..
+These variables are added to the starfi file.
+Other GCM variables, outputs and subroutines are not modified yet.
+Only nslope=1 is accepted so far (ie. same as no subslope parametrization)
+
 == 16/02/2023 == RV
 MARS PEM : Deep cleaning of variables name and allocate.
 All the "dyn to phys" grid change is done in subroutines and not in the main program.
+
+== 21/02/2023 == RV
+Following r2896, further implementation of subslope parametrisation.
+Carefull ! This is a devolpment revision and it still need improvements and tests. 
+However, this commit should not change anything for nslope=1.
+Only nslope=1 is possible for now!
+Utilitaries are not adapted yet.
+Dimension of some variables (30 listed below) is changed to add the nslope dimension.
+Outputs (diagfi and restartfi) are not changed yet.
+nslope is now read and written in the startfi.nc
+
+
Index: /trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90	(revision 2900)
@@ -13,5 +13,5 @@
 IMPLICIT NONE
 
-INTEGER, parameter :: nslope = 1                            ! Number of slopes for the statistic  (1)
+INTEGER, SAVE :: nslope                                ! Number of slopes for the statistic  (1)
 INTEGER, SAVE :: iflat                                      ! Number of slopes for the statistic  (1)
 REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope             ! bound of the slope statistics / repartitions (°)
@@ -20,5 +20,5 @@
 REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: subslope_dist       ! distribution  of the slopes (1)
 INTEGER,SAVE,ALLOCATABLE,DIMENSION(:) :: major_slope        ! Index of the subslope that occupies most of themesh (1)
-!$OMP THREADPRIVATE(def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)
+!$OMP THREADPRIVATE(nslope,def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)
 
 CONTAINS 
Index: /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 2900)
@@ -17,11 +17,11 @@
 
   ! variables (FC: built in firstcall in soil.F)
-  REAL,SAVE,ALLOCATABLE :: tsoil(:,:)       ! sub-surface temperatures (K)
-  real,save,allocatable :: mthermdiff(:,:)  ! (FC) mid-layer thermal diffusivity
-  real,save,allocatable :: thermdiff(:,:)   ! (FC) inter-layer thermal diffusivity
+  REAL,SAVE,ALLOCATABLE :: tsoil(:,:,:)       ! sub-surface temperatures (K)
+  real,save,allocatable :: mthermdiff(:,:,:)  ! (FC) mid-layer thermal diffusivity
+  real,save,allocatable :: thermdiff(:,:,:)   ! (FC) inter-layer thermal diffusivity
   real,save,allocatable :: coefq(:)         ! (FC) q_{k+1/2} coefficients
-  real,save,allocatable :: coefd(:,:)       ! (FC) d_k coefficients
-  real,save,allocatable :: alph(:,:)        ! (FC) alpha_k coefficients
-  real,save,allocatable :: beta(:,:)        ! beta_k coefficients
+  real,save,allocatable :: coefd(:,:,:)       ! (FC) d_k coefficients
+  real,save,allocatable :: alph(:,:,:)        ! (FC) alpha_k coefficients
+  real,save,allocatable :: beta(:,:,:)        ! beta_k coefficients
   real,save :: mu
   real,save,allocatable :: flux_geo(:)       ! Geothermal Flux (W/m^2)
@@ -31,21 +31,21 @@
 contains
 
-  subroutine ini_comsoil_h(ngrid)
+  subroutine ini_comsoil_h(ngrid,nslope)
   
   implicit none
   integer,intent(in) :: ngrid ! number of atmospheric columns
-  
+  integer,intent(in) :: nslope ! number of sub grid slopes
     allocate(layer(nsoilmx)) !soil layer depths
     allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
     allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia
  
-    allocate(tsoil(ngrid,nsoilmx)) ! soil temperatures
+    allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
 
-    allocate(mthermdiff(ngrid,0:nsoilmx-1))
-    allocate(thermdiff(ngrid,nsoilmx-1))
+    allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
+    allocate(thermdiff(ngrid,nsoilmx-1,nslope))
     allocate(coefq(0:nsoilmx-1))
-    allocate(coefd(ngrid,nsoilmx-1))
-    allocate(alph(ngrid,nsoilmx-1))
-    allocate(beta(ngrid,nsoilmx-1))
+    allocate(coefd(ngrid,nsoilmx-1,nslope))
+    allocate(alph(ngrid,nsoilmx-1,nslope))
+    allocate(beta(ngrid,nsoilmx-1,nslope))
     allocate(flux_geo(ngrid))
  
Index: /trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/dimradmars_mod.F90	(revision 2900)
@@ -117,7 +117,7 @@
 !!
   REAL,SAVE,ALLOCATABLE :: dtrad(:,:) ! Net atm. radiative heating rate (K.s-1)
-  REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:) ! rad. flux from sky absorbed by surface (W.m-2)
-  REAL,SAVE,ALLOCATABLE :: fluxrad(:) ! Net radiative surface flux (W.m-2)
-  REAL,SAVE,ALLOCATABLE :: albedo(:,:) ! Surface albedo in each solar band
+  REAL,SAVE,ALLOCATABLE :: fluxrad_sky(:,:) ! rad. flux from sky absorbed by surface (W.m-2)
+  REAL,SAVE,ALLOCATABLE :: fluxrad(:,:) ! Net radiative surface flux (W.m-2)
+  REAL,SAVE,ALLOCATABLE :: albedo(:,:,:) ! Surface albedo in each solar band
   REAL,SAVE,ALLOCATABLE :: totcloudfrac(:) ! total cloud fraction over the column
 ! aerosol (dust or ice) extinction optical depth  at reference wavelength 
@@ -181,5 +181,5 @@
 contains
 
-  subroutine ini_dimradmars_mod(ngrid,nlayer)
+  subroutine ini_dimradmars_mod(ngrid,nlayer,nslope)
   
   implicit none
@@ -187,5 +187,5 @@
   integer,intent(in) :: ngrid ! number of atmospheric columns
   integer,intent(in) :: nlayer ! number of atmospheric layers
-
+  integer,intent(in) :: nslope ! number of subgrid scale slopes
    nflev=nlayer
 !  ndomainsz=ngrid
@@ -195,8 +195,8 @@
    ndlo2=ndlon
 
-   allocate(albedo(ngrid,2))
+   allocate(albedo(ngrid,2,nslope))
    allocate(dtrad(ngrid,nlayer))
-   allocate(fluxrad_sky(ngrid))
-   allocate(fluxrad(ngrid))
+   allocate(fluxrad_sky(ngrid,nslope))
+   allocate(fluxrad(ngrid,nslope))
    allocate(nueffdust(ngrid,nlayer))
    allocate(totcloudfrac(ngrid))
Index: /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/phys_state_var_init_mod.F90	(revision 2900)
@@ -61,4 +61,8 @@
       use dust_rad_adjust_mod, only: ini_dust_rad_adjust_mod, &
                                      end_dust_rad_adjust_mod
+      use comslope_mod, ONLY: nslope
+      use netcdf
+      USE mod_phys_lmdz_para, ONLY: is_master,bcast
+
       IMPLICIT NONE
       
@@ -71,9 +75,32 @@
       INTEGER,INTENT(in) :: dyn_nqperes
       INTEGER,INTENT(in) :: dyn_nqfils(nq)
+      character(len=10) :: filename  ! name of the startfi.nc
+      integer :: ncid, status, nslope_dim_id
+      integer :: nslope_read
+
+      filename="startfi.nc"
+      if(is_master) then
+        status = nf90_open(filename, nf90_nowrite, ncid)
+        if (status /= nf90_noerr) then
+          nslope=1
+        else
+          status = nf90_inq_dimid(ncid, "nslope", nslope_dim_id)
+          if (status /= nf90_noerr) then
+            nslope=1
+          else
+            status = nf90_inquire_dimension(ncid, nslope_dim_id, len = nslope_read)
+            if (status /= nf90_noerr) then
+              call abort_physic("phys_state_var_init","nslope present but not readable",1)
+            else
+              nslope=nslope_read
+            endif
+          endif
+        endif
+      endif
+      call bcast(nslope)
 
       ! set dimension and allocate arrays in tracer_mod
       call end_tracer_mod
       call ini_tracer_mod(nq,tname,dyn_nqperes,dyn_nqfils)! MVals: variables isotopes
-
 
       ! initialize "print_control" constants/flags ("prt_level","lunout", etc.)
@@ -101,5 +128,5 @@
       ! allocate "surfdat_h" arrays
       call end_surfdat_h
-      call ini_surfdat_h(ngrid,nq)
+      call ini_surfdat_h(ngrid,nq,nslope)
 
       ! allocate "comgeomfi_h" arrays
@@ -109,9 +136,9 @@
       ! allocate "comsoil_h" arrays
       call end_comsoil_h
-      call ini_comsoil_h(ngrid)
+      call ini_comsoil_h(ngrid,nslope)
 
       ! set some variables in "dimradmars_mod"
       call end_dimradmars_mod
-      call ini_dimradmars_mod(ngrid,nlayer)
+      call ini_dimradmars_mod(ngrid,nlayer,nslope)
 
       ! allocate arrays in "yomlw_h"
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2900)
@@ -271,5 +271,5 @@
       real nuiceco2(ngrid,nlayer) ! Estimated effective variance of the 
                                   !   size distribution
-      REAL inertiesoil(ngrid,nsoilmx) ! Time varying subsurface
+      REAL inertiesoil(ngrid,nsoilmx,nslope) ! Time varying subsurface
                                       ! thermal inertia (J.s-1/2.m-2.K-1)
                                       ! (used only when tifeedback=.true.)
@@ -301,6 +301,6 @@
       INTEGER l,ig,ierr,igout,iq,tapphys,isoil
 
-      REAL fluxsurf_lw(ngrid)      !incident LW (IR) surface flux (W.m-2)
-      REAL fluxsurf_dn_sw(ngrid,2) ! Incident SW (solar) surface flux (W.m-2)
+      REAL fluxsurf_lw(ngrid,nslope)      !incident LW (IR) surface flux (W.m-2)
+      REAL fluxsurf_dn_sw(ngrid,2,nslope) ! Incident SW (solar) surface flux (W.m-2)
       REAL fluxsurf_up_sw(ngrid,2) ! Reflected SW (solar) surface flux (W.m-2)
       REAL fluxtop_lw(ngrid)       !Outgoing LW (IR) flux to space (W.m-2)
@@ -333,5 +333,5 @@
 
 c     Tendancies due to various processes:
-      REAL dqsurf(ngrid,nq)  ! tendency for tracers on surface (Kg/m2/s)
+      REAL dqsurf(ngrid,nq,nslope)  ! tendency for tracers on surface (Kg/m2/s)
       REAL zdtlw(ngrid,nlayer)     ! (K/s)
       REAL zdtsw(ngrid,nlayer)     ! (K/s)
@@ -340,16 +340,16 @@
       REAL zdtnirco2(ngrid,nlayer) ! (K/s)
       REAL zdtnlte(ngrid,nlayer)   ! (K/s)
-      REAL zdtsurf(ngrid)            ! (K/s)
+      REAL zdtsurf(ngrid,nslope)            ! (K/s)
       REAL zdtcloud(ngrid,nlayer),zdtcloudco2(ngrid,nlayer)
       REAL zdvdif(ngrid,nlayer),zdudif(ngrid,nlayer)  ! (m.s-2)
-      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid)         ! (K/s)
+      REAL zdhdif(ngrid,nlayer), zdtsdif(ngrid,nslope)         ! (K/s)
       REAL zdvadj(ngrid,nlayer),zduadj(ngrid,nlayer)  ! (m.s-2)
       REAL zdhadj(ngrid,nlayer)                           ! (K/s)
       REAL zdtgw(ngrid,nlayer)                            ! (K/s)
       REAL zdugw(ngrid,nlayer),zdvgw(ngrid,nlayer)    ! (m.s-2)
-      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid)
+      REAL zdtc(ngrid,nlayer),zdtsurfc(ngrid,nslope)
       REAL zdvc(ngrid,nlayer),zduc(ngrid,nlayer)
 
-      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq)
+      REAL zdqdif(ngrid,nlayer,nq), zdqsdif(ngrid,nq,nslope)
       REAL zdqsed(ngrid,nlayer,nq), zdqssed(ngrid,nq)
       REAL zdqdev(ngrid,nlayer,nq), zdqsdev(ngrid,nq)
@@ -357,5 +357,5 @@
       REAL zdqc(ngrid,nlayer,nq)
       REAL zdqcloudco2(ngrid,nlayer,nq)
-      REAL zdqsc(ngrid,nq)
+      REAL zdqsc(ngrid,nq,nslope)
 
       REAL zdteuv(ngrid,nlayer)    ! (K/s)
@@ -366,8 +366,8 @@
       real*8 PhiEscH,PhiEscH2,PhiEscD
 
-      REAL dwatercap(ngrid), dwatercap_dif(ngrid)     ! (kg/m-2)
+      REAL dwatercap(ngrid,nslope), dwatercap_dif(ngrid,nslope)     ! (kg/m-2)
 
 c     Local variable for local intermediate calcul:
-      REAL zflubid(ngrid)
+      REAL zflubid(ngrid,nslope)
       REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
       REAL zdum1(ngrid,nlayer)
@@ -393,5 +393,5 @@
 
       REAL fluxtop_dn_sw_tot(ngrid), fluxtop_up_sw_tot(ngrid)
-      REAL fluxsurf_dn_sw_tot(ngrid), fluxsurf_up_sw_tot(ngrid)
+      REAL fluxsurf_dn_sw_tot(ngrid,nslope), fluxsurf_up_sw_tot(ngrid)
       character*2 str2
 !      character*5 str5
@@ -534,5 +534,12 @@
 
 c Sub-grid scale slopes
+      real :: tsurf_meshavg(ngrid) ! Surface temperature grid box averaged [K]
+      real :: albedo_meshavg(ngrid,2) ! albedo temperature grid box averaged [1]
+      real :: emis_meshavg(ngrid,2) ! emis temperature grid box averaged [1]
+      real :: qsurf_meshavg(ngrid,nq) ! surface tracer mesh averaged [kg/m^2]
+      real :: qsurf_tmp(ngrid,nq) ! temporary qsurf for chimie
+      real :: inertiedat_tmp(ngrid,nsoilmx,nslope) ! temporary inertiedat for soil.F
       integer :: islope
+      real :: zdqsdif_meshavg_tmp(ngrid,nq) ! temporary for dust lifting
 
       logical :: write_restart
@@ -557,5 +564,5 @@
 
 #ifndef MESOSCALE
-         fluxrad(:)=0
+         fluxrad(:,:)=0
          wstar(:)=0.
 #endif
@@ -600,4 +607,8 @@
        PRINT*,'corresponding criterium = ',def_slope_mean(iflat)
 
+       DO islope = 1,nslope
+         inertiedat_tmp(:,:,islope) = inertiedat(:,:)
+       ENDDO
+
 #else
 ! MESOSCALE. Supposedly everything is already set in modules.
@@ -606,13 +617,13 @@
       print*,"check: rad,cpp,g,r,rcp,daysec"
       print*,rad,cpp,g,r,rcp,daysec
-      PRINT*,'check: tsurf ',tsurf(1),tsurf(ngrid)
-      PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngrid,nsoilmx)
+      PRINT*,'check: tsurf ',tsurf(1,:),tsurf(ngrid,:)
+      PRINT*,'check: tsoil ',tsoil(1,1,:),tsoil(ngrid,nsoilmx,:)
       PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngrid,nsoilmx)
       PRINT*,'check: midlayer,layer ', mlayer(:),layer(:)
       PRINT*,'check: tracernames ', noms
-      PRINT*,'check: emis ',emis(1),emis(ngrid)
+      PRINT*,'check: emis ',emis(1,:),emis(ngrid,:)
       PRINT*,'check: q2 ',q2(1,1),q2(ngrid,nlayer+1)
-      PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngrid,nq)
-      PRINT*,'check: co2ice ',qsurf(1,igcm_co2),qsurf(ngrid,igcm_co2)
+      PRINT*,'check: qsurf ',qsurf(1,1,:),qsurf(ngrid,nq,:)
+      PRINT*,'check: co2ice ',qsurf(1,igcm_co2,:),qsurf(ngrid,igcm_co2,:)
       !!!
       day_ini = pday
@@ -645,9 +656,10 @@
            endif
            ! additionnal "academic" initialization of physics
-           write(*,*) "Physiq: initializing tsurf(:) to pt(:,1) !!"
-           tsurf(:)=pt(:,1)
+           do islope = 1,nslope
+             tsurf(:,islope)=pt(:,1)
+           enddo
            write(*,*) "Physiq: initializing tsoil(:) to pt(:,1) !!"
            do isoil=1,nsoilmx
-             tsoil(1:ngrid,isoil)=tsurf(1:ngrid)
+             tsoil(1:ngrid,isoil,:)=tsurf(1:ngrid,:)
            enddo
            write(*,*) "Physiq: initializing inertiedat !!"
@@ -683,9 +695,12 @@
 c           Thermal inertia feedback:
             IF (tifeedback) THEN
-                CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
+             DO islope = 1,nslope
+                CALL soil_tifeedback(ngrid,nsoilmx,
+     s               qsurf(:,:,islope),inertiesoil(:,:,islope))
+             ENDDO
                 CALL soil(ngrid,nsoilmx,firstcall,inertiesoil,
      s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
             ELSE
-                CALL soil(ngrid,nsoilmx,firstcall,inertiedat,
+                CALL soil(ngrid,nsoilmx,firstcall,inertiedat_tmp,
      s             ptimestep,tsurf,tsoil,capcal,fluxgrd)
             ENDIF ! of IF (tifeedback)
@@ -694,6 +709,6 @@
      &     'PHYSIQ WARNING! Thermal conduction in the soil turned off'
             DO ig=1,ngrid
-               capcal(ig)=1.e5
-               fluxgrd(ig)=0.
+               capcal(ig,:)=1.e5
+               fluxgrd(ig,:)=0.
             ENDDO
          ENDIF
@@ -793,11 +808,11 @@
       pdq(:,:,:)=0
       pdpsrf(:)=0
-      zflubid(:)=0
-      zdtsurf(:)=0
-      dqsurf(:,:)=0
+      zflubid(:,:)=0
+      zdtsurf(:,:)=0
+      dqsurf(:,:,:)=0
       dsodust(:,:)=0.
       dsords(:,:)=0.
       dsotop(:,:)=0.
-      dwatercap(:)=0
+      dwatercap(:,:)=0
       
       ! Dust scenario conversion coefficient from IRabs to VISext
@@ -899,5 +914,8 @@
      &        +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
           end do
-          co2totA = co2totA + qsurf(ig,igcm_co2)
+          do islope = 1,nslope
+          co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
+     &    subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+          enddo
         end do
       else
@@ -909,5 +927,8 @@
      &        +pdq(ig,l,igcm_co2)*ptimestep)
           end do
-          co2totA = co2totA + qsurf(ig,igcm_co2)
+          do islope = 1,nslope
+          co2totA = co2totA + qsurf(ig,igcm_co2,islope)*
+     &    subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)
+          enddo
         end do
       endif ! of if (igcm_co2_ice.ne.0)
@@ -1007,13 +1028,20 @@
            ! callradite for the part with clouds
            clearsky=.false. ! part with clouds for both cases CLFvarying true/false
-           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
-     &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
-     &     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_dn_sw,fluxsurf_up_sw,
+           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
+     &     albedo(:,:,1),
+     &     emis(:,1),mu0,zplev,zplay,pt,tsurf(:,1),fract,dist_sol,igout,
+     &     zdtlw,zdtsw,fluxsurf_lw(:,iflat),fluxsurf_dn_sw(:,:,iflat),
+     &     fluxsurf_up_sw,
      &     fluxtop_lw,fluxtop_dn_sw,fluxtop_up_sw,
      &     tau_pref_scenario,tau_pref_gcm,
      &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,
-     &     qsurf(:,igcm_co2),rstormdust,rtopdust,totstormfract,
+     &     qsurf(:,igcm_co2,1),rstormdust,rtopdust,totstormfract,
      &     clearatm,dsords,dsotop,nohmons,clearsky,totcloudfrac)
+
+           DO islope=1,nslope
+             fluxsurf_lw(:,islope) =fluxsurf_lw(:,iflat)
+             fluxsurf_dn_sw(:,:,islope) =fluxsurf_dn_sw(:,:,iflat)
+           ENDDO
 
            ! case of sub-grid water ice clouds: callradite for the clear case
@@ -1025,5 +1053,6 @@
                clearsky=.true.  
                CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
-     &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
+     &              albedo(:,:,1),emis(:,1),mu0,zplev,zplay,pt,
+     &              tsurf(:,1),fract,
      &              dist_sol,igout,zdtlwclf,zdtswclf,
      &              fluxsurf_lwclf,fluxsurf_dn_swclf,fluxsurf_up_swclf,
@@ -1032,5 +1061,5 @@
      &              dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &              taucloudtesclf,rdust,
-     &              rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2),
+     &              rice,nuice,riceco2, nuiceco2,qsurf(:,igcm_co2,1),
      &              rstormdust,rtopdust,totstormfract,
      &              clearatm,dsords,dsotop,
@@ -1042,9 +1071,11 @@
                   tf_clf=totcloudfrac(ig)
                   ntf_clf=1.-tf_clf
-                  fluxsurf_lw(ig) = ntf_clf*fluxsurf_lwclf(ig) 
-     &                                      + tf_clf*fluxsurf_lw(ig)
-                  fluxsurf_dn_sw(ig,1:2) = 
+                  DO islope=1,nslope
+                  fluxsurf_lw(ig,islope) = ntf_clf*fluxsurf_lwclf(ig) 
+     &                          + tf_clf*fluxsurf_lw(ig,islope)
+                  fluxsurf_dn_sw(ig,1:2,islope) = 
      &                           ntf_clf*fluxsurf_dn_swclf(ig,1:2) 
-     &                          + tf_clf*fluxsurf_dn_sw(ig,1:2)
+     &                          + tf_clf*fluxsurf_dn_sw(ig,1:2,islope)
+                  ENDDO
                   fluxsurf_up_sw(ig,1:2) = 
      &                           ntf_clf*fluxsurf_up_swclf(ig,1:2) 
@@ -1099,9 +1130,14 @@
 c          ---------------------------------------------------------
            IF(callslope) THEN
+           ! assume that in this case, nslope = 1
+            if(nslope.ne.1) then
+              call abort_physic(
+     &      "physiq","callslope=true but nslope.ne.1",1)
+            endif
             print *, 'Slope scheme is on and computing...'
             DO ig=1,ngrid  
               sl_the = theta_sl(ig)
               IF (sl_the .ne. 0.) THEN
-                ztim1=fluxsurf_dn_sw(ig,1)+fluxsurf_dn_sw(ig,2)
+                ztim1=fluxsurf_dn_sw(ig,1,1)+fluxsurf_dn_sw(ig,2,1)
                 DO l=1,2
                  sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
@@ -1109,7 +1145,7 @@
                  sl_lat = 180.*latitude(ig)/pi
                  sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
-                 sl_alb = albedo(ig,l)
+                 sl_alb = albedo(ig,l,1)
                  sl_psi = psi_sl(ig)
-                 sl_fl0 = fluxsurf_dn_sw(ig,l)
+                 sl_fl0 = fluxsurf_dn_sw(ig,l,1)
                  sl_di0 = 0.
                  if ((mu0(ig) .gt. 0.).and.(ztim1.gt.0.)) then
@@ -1117,5 +1153,5 @@
                   sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
                   sl_di0 = sl_di0/ztim1
-                  sl_di0 = fluxsurf_dn_sw(ig,l)*sl_di0
+                  sl_di0 = fluxsurf_dn_sw(ig,l,1)*sl_di0
                  endif
                  ! you never know (roundup concern...)
@@ -1126,12 +1162,55 @@
      &                             sl_di0, sl_fl0, sl_flu )
                  !!!!!!!!!!!!!!!!!!!!!!!!!!
-                 fluxsurf_dn_sw(ig,l) = sl_flu
+                 fluxsurf_dn_sw(ig,l,1) = sl_flu
                 ENDDO
               !!! compute correction on IR flux as well
               sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
-              fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky
+              fluxsurf_lw(ig,:)= fluxsurf_lw(ig,:)*sky
               ENDIF
             ENDDO
-           ENDIF
+c          ELSE        ! not calling subslope, nslope might be > 1
+c          DO islope = 1,nslope
+c            IF(def_slope_mean(islope).ge.0.) THEN
+c              sl_psi = 0. !Northward slope
+c            ELSE
+c             sl_psi = 180. !Southward slope
+c            ENDIF
+c            sl_the=abs(def_slope_mean(islope)) 	
+c            IF (sl_the .gt. 1e-6) THEN
+c              DO ig=1,ngrid  
+c                ztim1=fluxsurf_dn_sw(ig,1,islope)
+c     s                +fluxsurf_dn_sw(ig,2,islope)
+c                DO l=1,2
+c                 sl_lct = ptime*24. + 180.*longitude(ig)/pi/15.
+c                 sl_ra  = pi*(1.0-sl_lct/12.)
+c                 sl_lat = 180.*latitude(ig)/pi
+c                 sl_tau = tau(ig,1) !il faudrait iaerdust(iaer)
+c                 sl_alb = albedo(ig,l,islope)
+c                 sl_psi = psi_sl(ig)
+c                 sl_fl0 = fluxsurf_dn_sw(ig,l,islope)
+c                 sl_di0 = 0.
+c                 if (mu0(ig) .gt. 0.) then
+c		          sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig)))
+c                  sl_di0 = sl_di0*flux_1AU/dist_sol/dist_sol
+c                  sl_di0 = sl_di0/ztim1
+c                  sl_di0 = fluxsurf_dn_sw(ig,l,islope)*sl_di0
+c                 endif
+c                 ! you never know (roundup concern...)
+c                 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0
+c                 !!!!!!!!!!!!!!!!!!!!!!!!!!
+c                 CALL param_slope( mu0(ig), declin, sl_ra, sl_lat, 
+c     &                             sl_tau, sl_alb, sl_the, sl_psi,
+c     &                             sl_di0, sl_fl0, sl_flu )
+c                 !!!!!!!!!!!!!!!!!!!!!!!!!!
+c                 fluxsurf_dn_sw(ig,l,islope) = sl_flu
+c                ENDDO
+c              !!! compute correction on IR flux as well
+c
+c              fluxsurf_lw(ig,islope)= fluxsurf_lw(ig,islope)
+c     &                                *sky_slope(islope)
+c              ENDDO
+c            ENDIF
+c          ENDDO ! islope = 1,nslope
+          ENDIF ! callslope
 
 c          CO2 near infrared absorption
@@ -1146,7 +1225,10 @@
 c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            DO ig=1,ngrid
-               fluxrad_sky(ig)=emis(ig)*fluxsurf_lw(ig)
-     $         +fluxsurf_dn_sw(ig,1)*(1.-albedo(ig,1))
-     $         +fluxsurf_dn_sw(ig,2)*(1.-albedo(ig,2))
+             DO islope = 1,nslope
+               fluxrad_sky(ig,islope) =
+     $           emis(ig,islope)*fluxsurf_lw(ig,islope)
+     $         +fluxsurf_dn_sw(ig,1,islope)*(1.-albedo(ig,1,islope))         
+     $         +fluxsurf_dn_sw(ig,2,islope)*(1.-albedo(ig,2,islope))
+             ENDDO
            ENDDO
 
@@ -1175,13 +1257,20 @@
 c
            DO ig=1,ngrid
-               zplanck(ig)=tsurf(ig)*tsurf(ig)
-               zplanck(ig)=emis(ig)*
+             DO islope = 1,nslope
+               zplanck(ig)=tsurf(ig,islope)*tsurf(ig,islope)
+               zplanck(ig)=emis(ig,islope)*
      $         stephan*zplanck(ig)*zplanck(ig)
-               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
+               fluxrad(ig,islope)=fluxrad_sky(ig,islope)-zplanck(ig)
                IF(callslope) THEN
                  sky= (1.+cos(pi*theta_sl(ig)/180.))/2.
-                 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig)
+                 fluxrad(ig,nslope)=fluxrad(ig,nslope)+
+     $                             (1.-sky)*zplanck(ig)
+               ELSE 
+                 fluxrad(ig,islope)=fluxrad(ig,islope) +
+     $              (1.-sky_slope(iflat))*emis(ig,iflat)*
+     $              stephan*tsurf(ig,iflat)**4  
                ENDIF
            ENDDO
+         ENDDO
 
          DO l=1,nlayer
@@ -1224,5 +1313,5 @@
 c               for radiative transfer
      &                      clearatm,icount,zday,zls,
-     &                      tsurf,qsurf(:,igcm_co2),igout,
+     &                      tsurf,qsurf(:,igcm_co2,iflat),igout,
      &                      totstormfract,tauscaling,
      &                      dust_rad_adjust,IRtoVIScoef,
@@ -1287,5 +1376,6 @@
      &                pq,pdq,pt,pdt,zplev,zplay,zzlev,
      &                zzlay,zdtsw,zdtlw,
-     &                icount,zday,zls,tsurf,qsurf(:,igcm_co2), 
+     &                icount,zday,zls,tsurf(:,iflat),
+     &                qsurf(:,igcm_co2,iflat), 
      &                igout,aerosol,tauscaling,dust_rad_adjust,
      &                IRtoVIScoef,totstormfract,clearatm,
@@ -1349,7 +1439,9 @@
       IF (calldifv) THEN
          DO ig=1,ngrid
-            zflubid(ig)=fluxrad(ig)+fluxgrd(ig)
+           DO islope = 1,nslope
+            zflubid(ig,islope)=fluxrad(ig,islope)
+     &                        +fluxgrd(ig,islope)
+           ENDDO
          ENDDO
-
          zdum1(:,:)=0
          zdum2(:,:)=0
@@ -1371,5 +1463,5 @@
 
           DO ig=1, ngrid
-             IF (zh(ig,1) .lt. tsurf(ig)) THEN
+             IF (zh(ig,1) .lt. tsurf(ig,1)) THEN
                wstar(ig)=1.
                hfmax_th(ig)=0.2
@@ -1390,5 +1482,5 @@
 
 c        Calling vdif (Martian version WITH CO2 condensation)
-         dwatercap_dif(:) = 0.
+         dwatercap_dif(:,:) = 0.
          zcdh(:) = 0.
          zcdv(:) = 0.
@@ -1405,8 +1497,9 @@
 
           DO ig=1,ngrid
-             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
-             dwatercap(ig)=dwatercap(ig)+dwatercap_dif(ig)
+             zdtsurf(ig,:)=zdtsurf(ig,:)+zdtsdif(ig,:)
+             dwatercap(ig,:)=dwatercap(ig,:)+dwatercap_dif(ig,:)
           ENDDO
-
+c          call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,zdqsdif,
+c     &   albedo_meshavg,emis_meshavg,tsurf_meshavg,zdqsdif_meshavg_tmp)
          IF (.not.turb_resolved) THEN
           DO l=1,nlayer
@@ -1429,5 +1522,5 @@
            DO iq=1, nq
               DO ig=1,ngrid
-                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+                 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
               ENDDO
            ENDDO
@@ -1446,23 +1539,25 @@
              zdqdif(1:ngrid,1,1:nq)=0.
              zdqdif(1:ngrid,1,igcm_dust_number) = 
-     .                  -zdqsdif(1:ngrid,igcm_dust_number)
+     .                  -zdqsdif(1:ngrid,igcm_dust_number,iflat)
              zdqdif(1:ngrid,1,igcm_dust_mass) = 
-     .                  -zdqsdif(1:ngrid,igcm_dust_mass)
+     .                  -zdqsdif(1:ngrid,igcm_dust_mass,iflat)
              zdqdif(1:ngrid,2:nlayer,1:nq) = 0.
              DO iq=1, nq
                IF ((iq .ne. igcm_dust_mass)
      &          .and. (iq .ne. igcm_dust_number)) THEN
-                 zdqsdif(:,iq)=0.
+                 zdqsdif(:,iq,:)=0.
                ENDIF
              ENDDO
            ELSE
              zdqdif(1:ngrid,1:nlayer,1:nq) = 0.
-             zdqsdif(1:ngrid,1:nq) = 0.
+             zdqsdif(1:ngrid,1:nq,1:nslope) = 0.
            ENDIF
          ENDIF
       ELSE    
          DO ig=1,ngrid
-            zdtsurf(ig)=zdtsurf(ig)+
-     s        (fluxrad(ig)+fluxgrd(ig))/capcal(ig)
+           DO islope=1,nslope
+             zdtsurf(ig,islope)=zdtsurf(ig,islope)+
+     s   (fluxrad(ig,islope)+fluxgrd(ig,islope))/capcal(ig,islope)
+           ENDDO
          ENDDO
          IF (turb_resolved) THEN
@@ -1707,7 +1802,10 @@
           DO iq=1, nq
             DO ig=1,ngrid
-              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqssed_ccn(ig,iq)
+             DO islope = 1,nslope
+              dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope)+
+     &         zdqssed_ccn(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+             ENDDO !(islope)
             ENDDO  ! (ig)
-           ENDDO    ! (iq)
+           ENDDO    ! (iq)q)
 c Temperature variation due to latent heat release
                pdt(1:ngrid,1:nlayer) =
@@ -1884,5 +1982,8 @@
               do iq=1,nq
                  DO ig=1,ngrid
-                    dqsurf(ig,iq)= dqsurf(ig,iq) + zdqsdev(ig,iq)
+                   DO islope = 1,nslope
+                    dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
+     &          zdqsdev(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+                   ENDDO
                  ENDDO
               enddo
@@ -1933,5 +2034,8 @@
            DO iq=1, nq
              DO ig=1,ngrid
-                dqsurf(ig,iq)= dqsurf(ig,iq) + zdqssed(ig,iq)
+               DO islope = 1,nslope
+                    dqsurf(ig,iq,islope)= dqsurf(ig,iq,islope) +
+     &          zdqssed(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+               ENDDO
              ENDDO
            ENDDO
@@ -1950,5 +2054,5 @@
            DO iq=1, nq
               DO ig=1,ngrid
-                 dqsurf(ig,iq)=dqsurf(ig,iq) + zdqsdif(ig,iq)
+                 dqsurf(ig,iq,:)=dqsurf(ig,iq,:) + zdqsdif(ig,iq,:)
               ENDDO
            ENDDO
@@ -1974,8 +2078,12 @@
      $                       surfdust, surfice)
 !           call photochemistry
+            DO ig = 1,ngrid
+              qsurf_tmp(ig,:) = qsurf(ig,:,major_slope(ig))
+            ENDDO
             call calchim(ngrid,nlayer,nq,
      &                   ptimestep,zplay,zplev,pt,pdt,dist_sol,mu0,
      $                   zzlev,zzlay,zday,pq,pdq,zdqchim,zdqschim,
-     $                   zdqcloud,zdqscloud,tau(:,1),qsurf(:,igcm_co2),
+     $                   zdqcloud,zdqscloud,tau(:,1),
+     $                   qsurf_tmp(:,igcm_co2),
      $                   pu,pdu,pv,pdv,surfdust,surfice)
 
@@ -2006,5 +2114,8 @@
                       ! tracers is zero anyways
              DO ig=1,ngrid
-               dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq)
+              DO islope = 1,nslope
+               dqsurf(ig,iq,islope)=dqsurf(ig,iq,islope) +
+     &           zdqschim(ig,iq)*cos(pi*def_slope_mean(islope)/180.)
+              ENDDO
              ENDDO
            ENDDO ! of DO iq=1,nq
@@ -2013,6 +2124,8 @@
            if (igcm_h2o2.ne.0) then
              DO ig=1,ngrid
-               dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2)
-     &                                +zdqscloud(ig,igcm_h2o2)
+              DO islope = 1,nslope
+              dqsurf(ig,igcm_h2o2,islope)=dqsurf(ig,igcm_h2o2,islope)+
+     &      zdqscloud(ig,igcm_h2o2)*cos(pi*def_slope_mean(islope)/180.)
+              ENDDO
              ENDDO
            endif
@@ -2029,5 +2142,5 @@
       if (callthermos) then
         call thermosphere(ngrid,nlayer,nq,zplev,zplay,dist_sol,
-     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+     $     mu0,ptimestep,ptime,zday,tsurf(:,iflat),zzlev,zzlay,
      &     pt,pq,pu,pv,pdt,pdq,
      $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
@@ -2056,6 +2169,6 @@
          CALL geticecover(ngrid, 180.*zls/pi,
      .                  180.*longitude/pi, 180.*latitude/pi,
-     .                  qsurf(:,igcm_co2) )
-         qsurf(:,igcm_co2) = qsurf(:,igcm_co2) * 10000.
+     .                  qsurf_tmp(:,igcm_co2) )
+         qsurf_tmp(:,igcm_co2) = qsurf_tmp(:,igcm_co2) * 10000.
       ENDIF
       
@@ -2063,13 +2176,13 @@
       IF (callcond) THEN
          zdtc(:,:) = 0.
-         zdtsurfc(:) = 0.
+         zdtsurfc(:,:) = 0.
          zduc(:,:) = 0.
          zdvc(:,:) = 0.
          zdqc(:,:,:) = 0.
-         zdqsc(:,:) = 0.
+         zdqsc(:,:,:) = 0.
          CALL co2condens(ngrid,nlayer,nq,ptimestep,
      $              capcal,zplay,zplev,tsurf,pt,
      $              pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,
-     $              qsurf(:,igcm_co2),albedo,emis,rdust,
+     $              qsurf(:,igcm_co2,:),albedo,emis,rdust,
      $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
      $              fluxsurf_dn_sw,zls,
@@ -2079,5 +2192,5 @@
          DO iq=1, nq
            DO ig=1,ngrid
-              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
+              dqsurf(ig,iq,:)=dqsurf(ig,iq,:)+zdqsc(ig,iq,:)
            ENDDO  ! (ig)
          ENDDO    ! (iq)
@@ -2090,5 +2203,5 @@
          ENDDO
          DO ig=1,ngrid
-           zdtsurf(ig) = zdtsurf(ig) + zdtsurfc(ig)
+           zdtsurf(ig,:) = zdtsurf(ig,:) + zdtsurfc(ig,:)
          ENDDO
 
@@ -2130,7 +2243,8 @@
         DO iq=1, nq
           DO ig=1,ngrid
-
-            qsurf(ig,iq)=qsurf(ig,iq)+ptimestep*dqsurf(ig,iq)
-
+           DO islope = 1,nslope
+            qsurf(ig,iq,islope)=qsurf(ig,iq,islope)+
+     &                ptimestep*dqsurf(ig,iq,islope)
+           ENDDO 
           ENDDO  ! (ig)
         ENDDO    ! (iq)
@@ -2144,5 +2258,8 @@
 
       DO ig=1,ngrid
-         tsurf(ig)=tsurf(ig)+ptimestep*zdtsurf(ig) 
+        DO islope = 1,nslope
+            tsurf(ig,islope)=tsurf(ig,islope)+
+     &            ptimestep*zdtsurf(ig,islope) 
+       ENDDO
       ENDDO
 
@@ -2169,12 +2286,14 @@
 c       -------------------------------------------------------------
          do ig=1,ngrid
-           if ((qsurf(ig,igcm_co2).eq.0).and.
-     &        (qsurf(ig,igcm_h2o_ice).gt.frost_albedo_threshold)) then
+          do islope = 1,nslope
+           if ((qsurf(ig,igcm_co2,islope).eq.0).and.
+     &        (qsurf(ig,igcm_h2o_ice,islope)
+     &       .gt.frost_albedo_threshold)) then
              if ((watercaptag(ig)).and.(cst_cap_albedo)) then
-               albedo(ig,1) = albedo_h2o_cap
-               albedo(ig,2) = albedo_h2o_cap
+               albedo(ig,1,islope) = albedo_h2o_cap
+               albedo(ig,2,islope) = albedo_h2o_cap
              else
-               albedo(ig,1) = albedo_h2o_frost
-               albedo(ig,2) = albedo_h2o_frost
+               albedo(ig,1,islope) = albedo_h2o_frost
+               albedo(ig,2,islope) = albedo_h2o_frost
              endif !((watercaptag(ig)).and.(cst_cap_albedo)) then
 c              write(*,*) "frost thickness", qsurf(ig,igcm_h2o_ice)
@@ -2182,4 +2301,5 @@
 c     &        ,latitude(ig)*180./pi, longitude(ig)*180./pi
            endif
+          enddo ! islope
          enddo  ! of do ig=1,ngrid
       ENDIF  ! of IF (water)
@@ -2191,9 +2311,12 @@
 c       Thermal inertia feedback
         IF (tifeedback) THEN
-         CALL soil_tifeedback(ngrid,nsoilmx,qsurf,inertiesoil)
+         DO islope = 1,nslope
+           CALL soil_tifeedback(ngrid,nsoilmx,
+     s               qsurf(:,:,islope),inertiesoil(:,:,islope))
+         ENDDO
          CALL soil(ngrid,nsoilmx,.false.,inertiesoil,
      s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
         ELSE
-         CALL soil(ngrid,nsoilmx,.false.,inertiedat,
+         CALL soil(ngrid,nsoilmx,.false.,inertiedat_tmp,
      s          ptimestep,tsurf,tsoil,capcal,fluxgrd)
         ENDIF
@@ -2242,5 +2365,8 @@
 
       DO ig=1,ngrid
-         watercap(ig)=watercap(ig)+ptimestep*dwatercap(ig)
+        DO islope = 1,nslope
+         watercap(ig,islope)=watercap(ig,islope)+
+     s   ptimestep*dwatercap(ig,islope)
+        ENDDO
       ENDDO
 
@@ -2248,11 +2374,14 @@
 
         DO ig=1,ngrid
-           if (watercaptag(ig).and.
-     &        (qsurf(ig,igcm_h2o_ice).gt.frost_metam_threshold)) then
-
-                watercap(ig)=watercap(ig)+qsurf(ig,igcm_h2o_ice)
-     &                  - frost_metam_threshold
-                qsurf(ig,igcm_h2o_ice) = frost_metam_threshold
+         DO islope = 1,nslope
+           if (watercaptag(ig).and. (qsurf(ig,igcm_h2o_ice,islope)
+     &        .gt.frost_metam_threshold)) then
+
+                watercap(ig,islope)=watercap(ig,islope)
+     &          +qsurf(ig,igcm_h2o_ice,islope)
+     &         - frost_metam_threshold
+                qsurf(ig,igcm_h2o_ice,islope) = frost_metam_threshold
            endif ! (watercaptag(ig).and.
+          ENDDO
         ENDDO
 
@@ -2263,4 +2392,7 @@
 c  13. Write output files
 c  ----------------------
+
+      call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
+     &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
 
 c    -------------------------------
@@ -2312,6 +2444,7 @@
       fluxtop_up_sw_tot(1:ngrid)=fluxtop_up_sw(1:ngrid,1) +
      &                           fluxtop_up_sw(1:ngrid,2)
-      fluxsurf_dn_sw_tot(1:ngrid)=fluxsurf_dn_sw(1:ngrid,1) +
-     &                            fluxsurf_dn_sw(1:ngrid,2)
+      fluxsurf_dn_sw_tot(1:ngrid,1:nslope)=
+     &                           fluxsurf_dn_sw(1:ngrid,1,1:nslope) +
+     &                           fluxsurf_dn_sw(1:ngrid,2,1:nslope)
       fluxsurf_up_sw_tot(1:ngrid)=fluxsurf_up_sw(1:ngrid,1) +
      &                            fluxsurf_up_sw(1:ngrid,2)
@@ -2344,5 +2477,5 @@
          WRITE(*,'(a6,a10,2a15)') 'Ts','dTs','ps','dps'
          WRITE(*,'(2f10.5,2f15.5)')
-     s   tsurf(igout),zdtsurf(igout)*ptimestep,
+     s   tsurf(igout,:),zdtsurf(igout,:)*ptimestep,
      s   zplev(igout,1),pdpsrf(igout)*ptimestep
          WRITE(*,'(a4,a6,5a10)') 'l','u','du','v','dv','T','dT'
@@ -2374,6 +2507,6 @@
 
            call pbl_parameters(ngrid,nlayer,ps,zplay,z0,
-     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out,
-     & T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
+     & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf(:,iflat),
+     & zh,z_out,n_out,T_out,u_out,ustar,tstar,L_mo,vhf,vvv)
                    ! pourquoi ustar recalcule ici? fait dans vdifc.
 
@@ -2472,7 +2605,8 @@
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
-     .                tsurf,tsoil,albedo,emis,
-     .                q2,qsurf,tauscaling,totcloudfrac,wstar,
-     .                watercap)
+     .                tsurf(:,iflat),tsoil(:,:,iflat),albedo(:,:,iflat),
+     .                emis(:,iflat),
+     .                q2,qsurf(:,:,iflat),tauscaling,totcloudfrac,wstar,
+     .                watercap(:,iflat))
           
          ENDIF ! of IF (write_restart)
@@ -2677,9 +2811,10 @@
          
         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
-        call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
+        call wstats(ngrid,"tsurf","Surface temperature","K",2
+     &      ,tsurf(:,iflat))
         call wstats(ngrid,"co2ice","CO2 ice cover",
-     &                "kg.m-2",2,qsurf(:,igcm_co2))
+     &                "kg.m-2",2,qsurf(:,igcm_co2,iflat))
         call wstats(ngrid,"watercap","H2O ice cover",
-     &                "kg.m-2",2,watercap)
+     &                "kg.m-2",2,watercap(:,iflat))
         call wstats(ngrid,"tau_pref_scenario",
      &              "prescribed visible dod at 610 Pa","NU",
@@ -2690,8 +2825,8 @@
         call wstats(ngrid,"fluxsurf_lw",
      &                "Thermal IR radiative flux to surface","W.m-2",2,
-     &                fluxsurf_lw)
+     &                fluxsurf_lw(:,iflat))
         call wstats(ngrid,"fluxsurf_dn_sw",
      &        "Incoming Solar radiative flux to surface","W.m-2",2,
-     &                fluxsurf_dn_sw_tot)
+     &                fluxsurf_dn_sw_tot(:,iflat))
         call wstats(ngrid,"fluxsurf_up_sw",
      &        "Reflected Solar radiative flux from surface","W.m-2",2,
@@ -2718,5 +2853,5 @@
      &                "m2.s-2",3,q2)
           call wstats(ngrid,"emis","Surface emissivity","w.m-1",2,
-     &                emis)
+     &                emis(:,iflat))
 c          call wstats(ngrid,"ssurf","Surface stress","N.m-2",
 c    &                2,zstress)
@@ -2764,8 +2899,8 @@
                call wstats(ngrid,"h2o_ice_s",
      &                    "surface h2o_ice","kg/m2",
-     &                    2,qsurf(1,igcm_h2o_ice))
+     &                    2,qsurf(1,igcm_h2o_ice,iflat))
                call wstats(ngrid,'albedo',
      &                         'albedo',
-     &                         '',2,albedo(1,1))
+     &                         '',2,albedo(1,1,iflat))
                call wstats(ngrid,"mtot",
      &                    "total mass of water vapor","kg/m2",
@@ -2852,5 +2987,5 @@
      &                         'kg/kg',3,zq(1,1,iq))
                  call wstats(ngrid,'qsurf'//str2,'qsurf',
-     &                         'kg.m-2',2,qsurf(1,iq))
+     &                         'kg.m-2',2,qsurf(1,iq,iflat))
                end do
              endif ! (doubleq)
@@ -3019,5 +3154,5 @@
 c        for any variables !!
          call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
-     &                  emis)
+     &                  emis(:,iflat))
 c        call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
 c        call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
@@ -3029,10 +3164,10 @@
      &                    "m2s-2",2,phisfi)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
-     &                  tsurf)
+     &                  tsurf(:,iflat))
          call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
          call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
-     &                              ,"kg.m-2",2,qsurf(:,igcm_co2))
+     &                              ,"kg.m-2",2,qsurf(:,igcm_co2,iflat))
          call WRITEDIAGFI(ngrid,"watercap","Water ice thickness"
-     &                                         ,"kg.m-2",2,watercap)
+     &         ,"kg.m-2",2,watercap(:,iflat))
 
          call WRITEDIAGFI(ngrid,"temp_layer1","temperature in layer 1",
@@ -3041,7 +3176,7 @@
      &                  "K",2,zt(1,7))
          call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2,
-     &                  fluxsurf_lw)
+     &                  fluxsurf_lw(:,iflat))
          call WRITEDIAGFI(ngrid,"fluxsurf_dn_sw","fluxsurf_dn_sw",
-     &                  "W.m-2",2,fluxsurf_dn_sw_tot)
+     &                  "W.m-2",2,fluxsurf_dn_sw_tot(:,iflat))
          call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2,
      &                  fluxtop_lw)
@@ -3079,8 +3214,8 @@
       !!! this is to ensure correct initialisation of mesoscale model
       call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
-     &                 tsurf)
+     &                 tsurf(:,iflat))
       call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps)
       call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2,
-     &                 qsurf(:,igcm_co2))
+     &                 qsurf(:,igcm_co2,iflat))
       call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
       call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
@@ -3089,5 +3224,5 @@
      &                 emis)
       call WRITEDIAGFI(ngrid,"tsoil","Soil temperature",
-     &                 "K",3,tsoil)
+     &                 "K",3,tsoil(:,:,iflat))
       call WRITEDIAGFI(ngrid,"inertiedat","Soil inertia",
      &                 "K",3,inertiedat)
@@ -3170,5 +3305,5 @@
           call WRITEDIAGFI(ngrid,'qsurf02','surface tracer',
      &                     'kg.m-2',2,
-     &                     qsurf(1:ngrid,igcm_h2o_ice))
+     &                     qsurf(1:ngrid,igcm_h2o_ice,iflat))
 #endif
           call WRITEDIAGFI(ngrid,'mtot',
@@ -3265,14 +3400,14 @@
             call WRITEDIAGFI(ngrid,'h2o_ice_s',
      &                       'surface h2o_ice',
-     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice))
+     &                       'kg.m-2',2,qsurf(1,igcm_h2o_ice,iflat))
             if (hdo) then
             call WRITEDIAGFI(ngrid,'hdo_ice_s',
      &                       'surface hdo_ice',
-     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice))
+     &                       'kg.m-2',2,qsurf(1,igcm_hdo_ice,iflat))
 
                 do ig=1,ngrid
-                if (qsurf(ig,igcm_h2o_ice).gt.qperemin) then
-                    DoH_surf(ig) = 0.5*( qsurf(ig,igcm_hdo_ice)/
-     &                  qsurf(ig,igcm_h2o_ice) )/155.76e-6
+                if (qsurf_meshavg(ig,igcm_h2o_ice).gt.qperemin) then
+           DoH_surf(ig) = 0.5*( qsurf_meshavg(ig,igcm_hdo_ice)/
+     &          qsurf_meshavg(ig,igcm_h2o_ice) )/155.76e-6
                 else
                     DoH_surf(ig) = 0.
@@ -3287,5 +3422,5 @@
             CALL WRITEDIAGFI(ngrid,'albedo',
      &                         'albedo',
-     &                         '',2,albedo(1,1))
+     &                         '',2,albedo(1,1,iflat))
               if (tifeedback) then
                  call WRITEDIAGSOIL(ngrid,"soiltemp",
@@ -3294,5 +3429,5 @@
                  call WRITEDIAGSOIL(ngrid,'soilti',
      &                       'Soil Thermal Inertia',
-     &                       'J.s-1/2.m-2.K-1',3,inertiesoil)
+     &                       'J.s-1/2.m-2.K-1',3,inertiesoil(:,:,iflat))
               endif
 !A. Pottier
@@ -3392,5 +3527,5 @@
      &                         'kg/kg',3,zq(1,1,iq))
                call WRITEDIAGFI(ngrid,'qsurf'//str2,'qsurf',
-     &                         'kg.m-2',2,qsurf(1,iq))
+     &                         'kg.m-2',2,qsurf(1,iq,iflat))
              end do
            endif ! (doubleq)
@@ -3419,8 +3554,8 @@
              call WRITEDIAGFI(ngrid,'qsurf',
      &                 'stormdust injection',
-     &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass))
+     &                 'kg.m-2',2,qsurf(:,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'pdqsurf',
      &                  'tendancy stormdust mass at surface',
-     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass))
+     &                  'kg.m-2',2,dqsurf(:,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
      &                        'm/s',3,wspeed(:,1:nlayer))
@@ -3484,7 +3619,7 @@
      &                        'part/kg',3,nccn)
              call WRITEDIAGFI(ngrid,'surfccnq','Surf nuclei mass mr',
-     &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass))
+     &                        'kg.m-2',2,qsurf(1,igcm_ccn_mass,iflat))
              call WRITEDIAGFI(ngrid,'surfccnN','Surf nuclei number',
-     &                        'kg.m-2',2,qsurf(1,igcm_ccn_number))
+     &                        'kg.m-2',2,qsurf(1,igcm_ccn_number,iflat))
            endif ! (scavenging)
 
@@ -3636,5 +3771,5 @@
          ! Write soil temperature
         call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-     &                     3,tsoil)
+     &                     3,tsoil(:,:,iflat))
          ! Write surface temperature
 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",
@@ -3666,5 +3801,5 @@
 c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
         call writediagsoil(ngrid,"soiltemp","Soil temperature","K",
-     &                     3,tsoil)
+     &                     3,tsoil(:,:,iflat))
 
 ! THERMALS STUFF 1D
@@ -3690,5 +3825,5 @@
          call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
          call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
-     &                  tsurf)
+     &                  tsurf(:,iflat))
          call WRITEDIAGFI(ngrid,"u","u wind","m/s",1,zu)
          call WRITEDIAGFI(ngrid,"v","v wind","m/s",1,zv)
@@ -3704,5 +3839,5 @@
      &                   'w.m-2',1,zdtlw)
          call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness"
-     &                                   ,"kg.m-2",0,qsurf(:,igcm_co2))
+     &       ,"kg.m-2",0,qsurf(:,igcm_co2,iflat))
 
          if (igcm_co2.ne.0) then
@@ -3742,7 +3877,7 @@
      &                        ,1,aerosol(:,:,iaer_stormdust_doubleq))
              call WRITEDIAGFI(ngrid,'dqsdifdustq','diffusion',
-     &                       'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass))
+     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_dust_mass,iflat))
              call WRITEDIAGFI(ngrid,'dqsdifrdsq','diffusion',
-     &                  'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass))
+     &          'kg.m-2.s-1',0,zdqsdif(1,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'mstormdtot',
      &                        'total mass of stormdust only',
@@ -3767,8 +3902,8 @@
              call WRITEDIAGFI(ngrid,'rdsqsurf',
      &                 'stormdust at surface',
-     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass))
+     &                 'kg.m-2',0,qsurf(:,igcm_stormdust_mass,iflat))
              call WRITEDIAGFI(ngrid,'qsurf',
      &                  'dust mass at surface',
-     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass))
+     &                  'kg.m-2',0,qsurf(:,igcm_dust_mass,iflat))
              call WRITEDIAGFI(ngrid,'wspeed','vertical speed stormdust',
      &                        'm/s',1,wspeed)
@@ -3824,9 +3959,9 @@
            mtot = 0
            icetot = 0
-           h2otot = qsurf(1,igcm_h2o_ice)
+           h2otot = qsurf(1,igcm_h2o_ice,iflat)
            if (hdo) THEN
            mtotD = 0
            icetotD = 0
-           hdotot = qsurf(1,igcm_hdo_ice)
+           hdotot = qsurf(1,igcm_hdo_ice,iflat)
            ENDIF !hdo
 
@@ -3860,5 +3995,5 @@
              CALL WRITEDIAGFI(ngrid,'albedo',
      &                         'albedo',
-     &                         '',2,albedo(1,1))
+     &                         '',2,albedo(1,1,iflat))
 
              IF (hdo) THEN
@@ -3970,5 +4105,5 @@
            do iq=1,nq
              call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
-     &            trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq))
+     &       trim(noms(iq))//'_s','kg/kg',0,qsurf(1,iq,iflat))
            end do
      
@@ -4030,5 +4165,5 @@
      &        +(pdq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2_ice))*ptimestep)
           enddo
-          co2totB = co2totB + qsurf(ig,igcm_co2)
+          co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
         enddo
       else
@@ -4039,5 +4174,5 @@
      &             (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
           enddo
-          co2totB = co2totB + qsurf(ig,igcm_co2)
+          co2totB = co2totB + qsurf(ig,igcm_co2,iflat)
         enddo
       endif ! of if (igcm_co2_ice.ne.0)
Index: /trunk/LMDZ.MARS/libf/phymars/soil.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/soil.F	(revision 2900)
@@ -8,5 +8,5 @@
      &                     coefd, alph, beta, mu,flux_geo
       use surfdat_h, only: watercaptag, inert_h2o_ice
-
+      use comslope_mod, ONLY: nslope
       implicit none
 
@@ -29,14 +29,14 @@
       integer nsoil	! number of soil layers 
       logical firstcall ! identifier for initialization call 
-      real therm_i(ngrid,nsoil) ! thermal inertia
+      real therm_i(ngrid,nsoil,nslope) ! thermal inertia
       real timestep	    ! time step
-      real tsurf(ngrid)   ! surface temperature
+      real tsurf(ngrid,nslope)   ! surface temperature
 ! outputs:
-      real tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
-      real capcal(ngrid) ! surface specific heat
-      real fluxgrd(ngrid) ! surface diffusive heat flux
+      real tsoil(ngrid,nsoil,nslope) ! soil (mid-layer) temperature
+      real capcal(ngrid,nslope) ! surface specific heat
+      real fluxgrd(ngrid,nslope) ! surface diffusive heat flux
 
 ! local variables:
-      integer ig,ik
+      integer ig,ik,islope
 
 ! 0. Initialisations and preprocessing step
@@ -46,24 +46,31 @@
 ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
       do ig=1,ngrid
+       do islope = 1,nslope
         if (watercaptag(ig)) then
           do ik=0,nsoil-1
 ! If we have permanent ice, we use the water ice thermal inertia from ground to last layer.
-              mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa
+              mthermdiff(ig,ik,islope)=inert_h2o_ice*
+     &            inert_h2o_ice/volcapa
           enddo
         else
           do ik=0,nsoil-1
-           mthermdiff(ig,ik)=therm_i(ig,ik+1)*therm_i(ig,ik+1)/volcapa
+           mthermdiff(ig,ik,islope)=therm_i(ig,ik+1,islope)*
+     &         therm_i(ig,ik+1,islope)/volcapa
           enddo
         endif
       enddo
+      enddo
 
 #ifdef MESOSCALE
       do ig=1,ngrid
-        if ( therm_i(ig,1) .ge. inert_h2o_ice ) then
-          print *, "limit max TI ", therm_i(ig,1), inert_h2o_ice
+       do islope = 1,nslope
+        if ( therm_i(ig,1,islope) .ge. inert_h2o_ice ) then
+          print *, "limit max TI ", therm_i(ig,1,islope), inert_h2o_ice
           do ik=0,nsoil-1
-               mthermdiff(ig,ik)=inert_h2o_ice*inert_h2o_ice/volcapa
+               mthermdiff(ig,ik,islope)=inert_h2o_ice*
+     &    inert_h2o_ice/volcapa
           enddo
         endif
+       enddo
       enddo
 #endif
@@ -71,10 +78,14 @@
 ! 0.2 Build thermdiff(:), the "interlayer" thermal diffusivities
       do ig=1,ngrid
+       do islope = 1,nslope
         do ik=1,nsoil-1
-      thermdiff(ig,ik)=((layer(ik)-mlayer(ik-1))*mthermdiff(ig,ik)
-     &                +(mlayer(ik)-layer(ik))*mthermdiff(ig,ik-1))
+      thermdiff(ig,ik,islope)=((layer(ik)-mlayer(ik-1))
+     &                        *mthermdiff(ig,ik,islope)
+     &                +(mlayer(ik)-layer(ik))
+     &                *mthermdiff(ig,ik-1,islope))
      &                    /(mlayer(ik)-mlayer(ik-1))
 !	write(*,*),'soil: ik: ',ik,' thermdiff:',thermdiff(ig,ik)
 	enddo
+       enddo
       enddo
 
@@ -92,27 +103,32 @@
 
       do ig=1,ngrid
+       do islope = 1,nslope
 	! d_k
 	do ik=1,nsoil-1
-	  coefd(ig,ik)=thermdiff(ig,ik)/(mlayer(ik)-mlayer(ik-1))
+	  coefd(ig,ik,islope)=thermdiff(ig,ik,islope)
+     &        /(mlayer(ik)-mlayer(ik-1))
 	enddo
 	
 	! alph_{N-1}
-	alph(ig,nsoil-1)=coefd(ig,nsoil-1)/
-     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1))
+	alph(ig,nsoil-1,islope)=coefd(ig,nsoil-1,islope)/
+     &                  (coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
         ! alph_k
         do ik=nsoil-2,1,-1
-	  alph(ig,ik)=coefd(ig,ik)/(coefq(ik)+coefd(ig,ik+1)*
-     &                              (1.-alph(ig,ik+1))+coefd(ig,ik))
+	  alph(ig,ik,islope)=coefd(ig,ik,islope)/
+     &         (coefq(ik)+coefd(ig,ik+1,islope)*
+     &          (1.-alph(ig,ik+1,islope))+coefd(ig,ik,islope))
 	enddo
 
         ! capcal
 ! Cstar
-        capcal(ig)=volcapa*layer(1)+
-     &              (thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
-     &              (timestep*(1.-alph(ig,1)))
+        capcal(ig,islope)=volcapa*layer(1)+
+     &         (thermdiff(ig,1,islope)/(mlayer(1)-mlayer(0)))*
+     &         (timestep*(1.-alph(ig,1,islope)))
 ! Cs
-        capcal(ig)=capcal(ig)/(1.+mu*(1.0-alph(ig,1))*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))
+        capcal(ig,islope)=capcal(ig,islope)/
+     &            (1.+mu*(1.0-alph(ig,1,islope))*
+     &            thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
 !      write(*,*)'soil: ig=',ig,' capcal(ig)=',capcal(ig)
+      enddo ! islope
       enddo ! of do ig=1,ngrid
             
@@ -122,25 +138,28 @@
       IF (.not.firstcall) THEN
 ! First layer:
-      do ig=1,ngrid
-        tsoil(ig,1)=(tsurf(ig)+mu*beta(ig,1)*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))/
-     &              (1.+mu*(1.0-alph(ig,1))*
-     &               thermdiff(ig,1)/mthermdiff(ig,0))
+      do islope = 1,nslope
+      do ig=1,ngrid
+        tsoil(ig,1,islope)=(tsurf(ig,islope)+mu*beta(ig,1,islope)*
+     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))/
+     &      (1.+mu*(1.0-alph(ig,1,islope))*
+     &      thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
       enddo
 ! Other layers:
       do ik=1,nsoil-1
         do ig=1,ngrid
-	  tsoil(ig,ik+1)=alph(ig,ik)*tsoil(ig,ik)+beta(ig,ik)
-	enddo
-      enddo
-      
+	  tsoil(ig,ik+1,islope)=alph(ig,ik,islope)*
+     &    tsoil(ig,ik,islope)+beta(ig,ik,islope)
+	enddo
+      enddo
+       enddo ! islope
       ENDIF! of if (.not.firstcall)
 
 !  2. Compute beta coefficients (preprocessing for next time step)
 ! Bottom layer, beta_{N-1}
-      do ig=1,ngrid
-        beta(ig,nsoil-1)=coefq(nsoil-1)*tsoil(ig,nsoil)
-     &            /(coefq(nsoil-1)+coefd(ig,nsoil-1))
-     &            +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1))
+       do islope = 1,nslope
+      do ig=1,ngrid
+        beta(ig,nsoil-1,islope)=coefq(nsoil-1)*tsoil(ig,nsoil,islope)
+     &      /(coefq(nsoil-1)+coefd(ig,nsoil-1,islope))
+     &      +flux_geo(ig)/(coefq(nsoil-1) + coefd(ig,nsoil-1,islope))
       enddo
    
@@ -148,8 +167,9 @@
       do ik=nsoil-2,1,-1
         do ig=1,ngrid
-	  beta(ig,ik)=(coefq(ik)*tsoil(ig,ik+1)+
-     &                 coefd(ig,ik+1)*beta(ig,ik+1))/
-     &                 (coefq(ik)+coefd(ig,ik+1)*(1.0-alph(ig,ik+1))
-     &                  +coefd(ig,ik))
+	  beta(ig,ik,islope)=(coefq(ik)*tsoil(ig,ik+1,islope)+
+     &                 coefd(ig,ik+1,islope)*beta(ig,ik+1,islope))/
+     &                 (coefq(ik)+coefd(ig,ik+1,islope)*
+     &                 (1.0-alph(ig,ik+1,islope))
+     &                  +coefd(ig,ik,islope))
 	enddo
       enddo
@@ -162,6 +182,8 @@
 !     &              (timestep*(1.-alph(ig,1)))
 ! Fstar
-        fluxgrd(ig)=(thermdiff(ig,1)/(mlayer(1)-mlayer(0)))*
-     &              (beta(ig,1)+(alph(ig,1)-1.0)*tsoil(ig,1))
+        fluxgrd(ig,islope)=(thermdiff(ig,1,islope)/
+     &              (mlayer(1)-mlayer(0)))*
+     &              (beta(ig,1,islope)+(alph(ig,1,islope)-1.0)
+     &              *tsoil(ig,1,islope))
 
 !        mu=mlayer(ig,0)/(mlayer(ig,1)-mlayer(ig,0))
@@ -169,11 +191,13 @@
 !     &                         thermdiff(ig,1)/mthermdiff(ig,0))
 ! Fs
-        fluxgrd(ig)=fluxgrd(ig)+(capcal(ig)/timestep)*
-     &              (tsoil(ig,1)*(1.+mu*(1.0-alph(ig,1))*
-     &                         thermdiff(ig,1)/mthermdiff(ig,0))
-     &               -tsurf(ig)-mu*beta(ig,1)*
-     &                          thermdiff(ig,1)/mthermdiff(ig,0))
-      enddo
-
+        fluxgrd(ig,islope)=fluxgrd(ig,islope)+(capcal(ig,islope)
+     &              /timestep)*
+     &              (tsoil(ig,1,islope)*
+     &              (1.+mu*(1.0-alph(ig,1,islope))*
+     &               thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
+     &               -tsurf(ig,islope)-mu*beta(ig,1,islope)*
+     &                thermdiff(ig,1,islope)/mthermdiff(ig,0,islope))
+      enddo
+      enddo ! islope
       end
 
Index: /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/surfdat_h.F90	(revision 2900)
@@ -44,10 +44,10 @@
 
   !! variables
-  REAL,SAVE,ALLOCATABLE :: tsurf(:)   ! Surface temperature (K)
-  REAL,SAVE,ALLOCATABLE :: emis(:)    ! Thermal IR surface emissivity
-  REAL,SAVE,ALLOCATABLE :: capcal(:) ! surface heat capacity (J m-2 K-1)
-  REAL,SAVE,ALLOCATABLE :: fluxgrd(:) ! surface conduction flux (W.m-2)
-  REAL,ALLOCATABLE,SAVE :: qsurf(:,:) ! tracer on surface (e.g. kg.m-2)
-  REAL,SAVE,ALLOCATABLE :: watercap(:) ! Surface water ice (kg.m-2)
+  REAL,SAVE,ALLOCATABLE :: tsurf(:,:)   ! Surface temperature (K)
+  REAL,SAVE,ALLOCATABLE :: emis(:,:)    ! Thermal IR surface emissivity
+  REAL,SAVE,ALLOCATABLE :: capcal(:,:) ! surface heat capacity (J m-2 K-1)
+  REAL,SAVE,ALLOCATABLE :: fluxgrd(:,:) ! surface conduction flux (W.m-2)
+  REAL,ALLOCATABLE,SAVE :: qsurf(:,:,:) ! tracer on surface (e.g. kg.m-2)
+  REAL,SAVE,ALLOCATABLE :: watercap(:,:) ! Surface water ice (kg.m-2)
 
 !$OMP THREADPRIVATE(tsurf,emis,capcal,fluxgrd,qsurf,watercap)
@@ -55,10 +55,10 @@
 contains
 
-  subroutine ini_surfdat_h(ngrid,nq)
+  subroutine ini_surfdat_h(ngrid,nq,nslope)
   
   implicit none
   integer,intent(in) :: ngrid ! number of atmospheric columns
   integer,intent(in) :: nq ! number of tracers  
-
+  integer,intent(in) :: nslope ! number of sub-grid scale slope  
     allocate(albedodat(ngrid))
     allocate(phisfi(ngrid))
@@ -71,10 +71,10 @@
     allocate(zthe(ngrid))
     allocate(z0(ngrid))
-    allocate(qsurf(ngrid,nq))
-    allocate(tsurf(ngrid))
-    allocate(watercap(ngrid))
-    allocate(emis(ngrid))
-    allocate(capcal(ngrid))
-    allocate(fluxgrd(ngrid))
+    allocate(qsurf(ngrid,nq,nslope))
+    allocate(tsurf(ngrid,nslope))
+    allocate(watercap(ngrid,nslope))
+    allocate(emis(ngrid,nslope))
+    allocate(capcal(ngrid,nslope))
+    allocate(fluxgrd(ngrid,nslope))
     allocate(hmons(ngrid))
     allocate(summit(ngrid))
Index: /trunk/LMDZ.MARS/libf/phymars/writediagfi.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 2899)
+++ /trunk/LMDZ.MARS/libf/phymars/writediagfi.F	(revision 2900)
@@ -78,5 +78,5 @@
 
       integer,save :: zitau=0
-      character(len=20),save :: firstnom='1234567890'
+      character(len=27),save :: firstnom='1234567890'
 !$OMP THREADPRIVATE(zitau,firstnom)
 
