Index: trunk/LMDZ.MARS/changelog.txt
===================================================================
--- trunk/LMDZ.MARS/changelog.txt	(revision 3179)
+++ trunk/LMDZ.MARS/changelog.txt	(revision 3183)
@@ -4442,2 +4442,5 @@
 == 18/01/2024 == JBC
 Improvement of the error message for tracers initialization with a 1D start file + update of "start1D.txt" in the deftank + small cleanings.
+
+== 25/01/2024 == JBC
+Some changes to prepare the introduction of slopes in 1D: in particular, the subroutine "getslopes.F90" and "param_slope.F90" are now inside the module "slope_mod.F90" + Few small cleanings throughout the code.
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 3179)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 3183)
@@ -1734,5 +1734,5 @@
 
         else if (trim(modif) .eq. 'nslope') then
-          write(*,*) 'set a new number of subgrid scale slope'
+          write(*,*) 'set a new number of subgrid scale slopes'
           write(*,*) 'Current value=', nslope
           write(*,*) 'Enter value for nslope (ex: 1,5,7)?'
@@ -1758,25 +1758,30 @@
 
       allocate(default_def_slope(nslope_new+1))
-      !Sub-grid scale  subslopes
-      if (nslope_new.eq.7) then
-      default_def_slope(1) = -43.
-      default_def_slope(2) = -19.
-      default_def_slope(3) = -9.
-      default_def_slope(4) = -3.
-      default_def_slope(5) = 3.
-      default_def_slope(6) = 9.
-      default_def_slope(7) = 19.
-      default_def_slope(8) = 43.
-      elseif (nslope_new.eq.5) then
-      default_def_slope(1) = -43.
-      default_def_slope(2) = -9.
-      default_def_slope(3) = -3.
-      default_def_slope(4) = 3.
-      default_def_slope(5) = 9.
-      default_def_slope(6) = 43.
-      elseif (nslope_new.eq.1) then
-      default_def_slope(1) = -50.
-      default_def_slope(2) = 50.
-      endif
+      ! Sub-grid scale slopes parameters (minimum/maximun angles)
+      select case(nslope)
+          case(1)
+              default_def_slope(1) = -50.
+              default_def_slope(2) = 50.
+          case(5)
+              default_def_slope(1) = -43.
+              default_def_slope(2) = -9.
+              default_def_slope(3) = -3.
+              default_def_slope(4) = 3.
+              default_def_slope(5) = 9.
+              default_def_slope(6) = 43.
+          case(7)
+              default_def_slope(1) = -43.
+              default_def_slope(2) = -19.
+              default_def_slope(3) = -9.
+              default_def_slope(4) = -3.
+              default_def_slope(5) = 3.
+              default_def_slope(6) = 9.
+              default_def_slope(7) = 19.
+              default_def_slope(8) = 43.
+          case default
+              write(*,*) 'Number of slopes not possible: nslope should 
+     & be 1, 5 or 7!'
+              call abort
+      end select
 
       do islope=1,nslope_new+1
@@ -1785,5 +1790,5 @@
 
       do islope=1,nslope_new
-       def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
+       def_slope_mean(islope)=(def_slope(islope)+def_slope(islope+1))/2.
       enddo
 
@@ -1798,6 +1803,5 @@
        call subslope_mola(ngridmx,nslope_new,def_slope,subslope_dist)
 
-! Surfdat related stuff
-        
+! Surfdat related stuff 
         allocate(tsurf_old_slope(ngridmx,nslope_old))  
         allocate(qsurf_old_slope(ngridmx,nqtot,nslope_old))  
@@ -1805,5 +1809,4 @@
         allocate(watercap_old_slope(ngridmx,nslope_old)) 
         allocate(perennial_co2_old_slope(ngridmx,nslope_old)) 
-
 
         tsurf_old_slope(:,:)=tsurf(:,:)
@@ -1816,5 +1819,4 @@
 
 ! Comsoil related stuff (tsoil)
-
         allocate(tsoil_old_slope(ngridmx,nsoilmx,nslope_old))
         allocate(inertiesoil_old_slope(ngridmx,nsoilmx,nslope_old))
@@ -1829,5 +1831,4 @@
 
 ! Dimradmars related stuff (albedo)
-
         allocate(albedo_old_slope(ngridmx,2,nslope_old))
         albedo_old_slope(:,:,:)=albedo(:,:,:)
Index: trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/comslope_mod.F90	(revision 3183)
@@ -3,5 +3,5 @@
 ! Subject:
 !---------
-!   Module used for the parametrization of subgrid scale slope 
+!   Module used for the parametrization of subgrid scale slope
 !----------------------------------------------------------------------------------------------------------------------!
 ! Reference:
@@ -11,88 +11,97 @@
 !======================================================================================================================!
 
-IMPLICIT NONE
+implicit none
 
-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 (°)
-REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: def_slope_mean        ! mean slope of each bin (°)
-REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: sky_slope             ! portion of the sky view by each slopes (1)
-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)
+integer, save                              :: nslope         ! Number of slopes for the statistic (1)
+integer, save                              :: iflat          ! Flat slope for islope (1)
+real,    save, allocatable, dimension(:)   :: def_slope      ! Bound of the slope statistics / repartitions (°)
+real,    save, allocatable, dimension(:)   :: def_slope_mean ! Mean slope of each bin (°)
+real,    save, allocatable, dimension(:)   :: sky_slope      ! Portion of the sky view by each slopes (1)
+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(nslope,def_slope,def_slope_mean,sky_slope,subslope_dist,iflat,major_slope)
 
-CONTAINS 
-  SUBROUTINE ini_comslope_h(ngrid,nslope_in)
+!=======================================================================
+contains
+!=======================================================================
 
-  IMPLICIT NONE
-  INTEGER, INTENT(IN) :: ngrid
-  INTEGER, INTENT(IN) :: nslope_in
+SUBROUTINE ini_comslope_h(ngrid,nslope_in)
 
-    allocate(def_slope(nslope_in+1))
-    allocate(def_slope_mean(nslope_in))
-    allocate(sky_slope(nslope_in))
-    allocate(subslope_dist(ngrid,nslope_in))
-    allocate(major_slope(ngrid))
-  END SUBROUTINE ini_comslope_h
+implicit none
 
+integer, intent(in) :: ngrid
+integer, intent(in) :: nslope_in
 
-  SUBROUTINE end_comslope_h
+allocate(def_slope(nslope_in+1))
+allocate(def_slope_mean(nslope_in))
+allocate(sky_slope(nslope_in))
+allocate(subslope_dist(ngrid,nslope_in))
+allocate(major_slope(ngrid))
 
-  IMPLICIT NONE
+END SUBROUTINE ini_comslope_h
 
-        IF (allocated(def_slope)) deallocate(def_slope)
-        IF (allocated(def_slope_mean)) deallocate(def_slope_mean)
-        IF (allocated(sky_slope)) deallocate(sky_slope)
-        IF (allocated(subslope_dist)) deallocate(subslope_dist)
-        IF (allocated(major_slope)) deallocate(major_slope)
+!=======================================================================
 
-  END SUBROUTINE end_comslope_h
+SUBROUTINE end_comslope_h
 
-  SUBROUTINE compute_meshgridavg(ngrid,nq,albedo_slope,emis_slope,tsurf_slope,qsurf_slope,albedo_meshavg,emis_meshavg,tsurf_meshavg, qsurf_meshavg)
-  USE comcstfi_h, only:  pi
-  IMPLICIT NONE
-  INTEGER, INTENT(IN) :: ngrid,nq                  !  # points on the physical grid, tracers  (1)
-  REAL, INTENT(IN) :: albedo_slope(ngrid,2,nslope) !  albedo on each sub-grid slope (1)
-  REAL, INTENT(IN) :: emis_slope(ngrid,nslope)     !  emissivity on each sub-grid slope (1)
-  REAL, INTENT(IN) :: tsurf_slope(ngrid,nslope)    !  Surface Temperature on each sub-grid slope (K)
-  REAL, INTENT(IN) :: qsurf_slope(ngrid,nq,nslope) !  Surface Tracers on each sub-grid slope (kg/m^2 sloped)
-  REAL, INTENT(OUT) :: albedo_meshavg(ngrid,2)     !  grid box average of the albedo (1)
-  REAL, INTENT(OUT) :: emis_meshavg(ngrid)         !  grid box average of the emissivity (1)
-  REAL, INTENT(OUT) :: tsurf_meshavg(ngrid)        !  grid box average of the surface temperature (K)
-  REAL, INTENT(OUT) :: qsurf_meshavg(ngrid,nq)     !  grid box average of the surface tracers (kg/m^2 flat)
+implicit none
+
+if (allocated(def_slope)) deallocate(def_slope)
+if (allocated(def_slope_mean)) deallocate(def_slope_mean)
+if (allocated(sky_slope)) deallocate(sky_slope)
+if (allocated(subslope_dist)) deallocate(subslope_dist)
+if (allocated(major_slope)) deallocate(major_slope)
+
+END SUBROUTINE end_comslope_h
+
+!=======================================================================
+
+SUBROUTINE compute_meshgridavg(ngrid,nq,albedo_slope,emis_slope,tsurf_slope,qsurf_slope,albedo_meshavg,emis_meshavg,tsurf_meshavg, qsurf_meshavg)
+
+use comcstfi_h, only: pi
+
+implicit none
+
+integer, intent(in) :: ngrid, nq                             ! # points on the physical grid, tracers (1)
+real, dimension(ngrid,2,nslope),  intent(in) :: albedo_slope ! albedo on each sub-grid slope (1)
+real, dimension(ngrid,nslope),    intent(in) :: emis_slope   ! emissivity on each sub-grid slope (1)
+real, dimension(ngrid,nslope),    intent(in) :: tsurf_slope  ! Surface Temperature on each sub-grid slope (K)
+real, dimension(ngrid,nq,nslope), intent(in) :: qsurf_slope  ! Surface Tracers on each sub-grid slope (kg/m^2 sloped)
+real, dimension(ngrid,2),  intent(out) :: albedo_meshavg ! grid box average of the albedo (1)
+real, dimension(ngrid),    intent(out) :: emis_meshavg   ! grid box average of the emissivity (1)
+real, dimension(ngrid),    intent(out) :: tsurf_meshavg  ! grid box average of the surface temperature (K)
+real, dimension(ngrid,nq), intent(out) :: qsurf_meshavg  ! grid box average of the surface tracers (kg/m^2 flat)
 ! Local
-  integer :: islope,ig,l,iq
+integer :: islope, ig, l, iq
 
 !-------------------
 
-      albedo_meshavg(:,:) = 0.
-      emis_meshavg(:) = 0.
-      tsurf_meshavg(:) = 0.
-      qsurf_meshavg(:,:) = 0.
+albedo_meshavg = 0.
+emis_meshavg = 0.
+tsurf_meshavg = 0.
+qsurf_meshavg = 0.
 
-  if(nslope.eq.1) then
-      albedo_meshavg(:,:) = albedo_slope(:,:,1)
-      emis_meshavg(:) = emis_slope(:,1)
-      tsurf_meshavg(:) = tsurf_slope(:,1)
-      qsurf_meshavg(:,:) = qsurf_slope(:,:,1)
-  else
+if (nslope == 1) then
+    albedo_meshavg = albedo_slope(:,:,1)
+    emis_meshavg = emis_slope(:,1)
+    tsurf_meshavg = tsurf_slope(:,1)
+    qsurf_meshavg = qsurf_slope(:,:,1)
+else
+    do ig = 1,ngrid
+        do islope = 1,nslope
+            do l = 1,2
+                albedo_meshavg(ig,l) = albedo_meshavg(ig,l) + albedo_slope(ig,l,islope)*subslope_dist(ig,islope)
+            enddo
+            do iq = 1,nq
+                qsurf_meshavg(ig,iq) = qsurf_meshavg(ig,iq) + qsurf_slope(ig,iq,islope)*subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
+            enddo
+            emis_meshavg(ig) = emis_meshavg(ig) + emis_slope(ig,islope)*subslope_dist(ig,islope)
+            tsurf_meshavg(ig) = tsurf_meshavg(ig) + (emis_slope(ig,islope)*tsurf_slope(ig,islope)**4)*subslope_dist(ig,islope)
+        enddo
+        tsurf_meshavg(ig) = (tsurf_meshavg(ig)/emis_meshavg(ig))**(0.25)
+    enddo
+endif
 
-  DO ig = 1,ngrid
-      DO islope = 1,nslope
-          DO l = 1,2
-            albedo_meshavg(ig,l) = albedo_meshavg(ig,l) + albedo_slope(ig,l,islope)*subslope_dist(ig,islope)
-          ENDDO
-          DO iq = 1,nq
-             qsurf_meshavg(ig,iq) = qsurf_meshavg(ig,iq) + qsurf_slope(ig,iq,islope)*subslope_dist(ig,islope)/cos(def_slope_mean(islope)*pi/180.)
-          ENDDO
-          emis_meshavg(ig) = emis_meshavg(ig) + emis_slope(ig,islope)*subslope_dist(ig,islope)
-          tsurf_meshavg(ig) = tsurf_meshavg(ig) + (emis_slope(ig,islope)*tsurf_slope(ig,islope)**4)*subslope_dist(ig,islope)
-      ENDDO
-      tsurf_meshavg(ig) = (tsurf_meshavg(ig)/emis_meshavg(ig))**(0.25)
-  ENDDO
-
-  endif
-
-  END SUBROUTINE compute_meshgridavg
+END SUBROUTINE compute_meshgridavg
 
 END MODULE comslope_mod
Index: trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/comsoil_h.F90	(revision 3183)
@@ -1,121 +1,127 @@
-module comsoil_h
+MODULE comsoil_h
 
 implicit none
+
 ! nsoilmx : number of subterranean layers
-  integer, parameter :: nsoilmx = 57 
+integer, parameter :: nsoilmx = 57
 
-  real,save,allocatable,dimension(:) :: layer      ! soil layer depths
-  real,save,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
-  real,save,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia for present climate
-  real,save,allocatable,dimension(:,:,:) :: inertiesoil ! soil thermal inertia
-  real,save :: volcapa    ! soil volumetric heat capacity
-       ! NB: volcapa is read fromn control(35) from physicq start file
-       !     in physdem (or set via tabfi, or initialized in
-       !                 soil_settings.F)
+real, save, allocatable, dimension(:)     :: layer       ! soil layer depths
+real, save, allocatable, dimension(:)     :: mlayer      ! soil mid-layer depths
+real, save, allocatable, dimension(:,:)   :: inertiedat  ! soil thermal inertia for present climate
+real, save, allocatable, dimension(:,:,:) :: inertiesoil ! soil thermal inertia
+real, save                                :: volcapa     ! soil volumetric heat capacity
+! NB: volcapa is read from control(35) from physics start file
+!     in physdem (or set via tabfi, or initialized in soil_settings.F)
 
 !$OMP THREADPRIVATE(layer,mlayer,inertiedat,inertiesoil,volcapa)
 
-  ! 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 :: 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 :: mu
-  real,save,allocatable :: flux_geo(:,:)      ! Geothermal Flux (W/m^2)
+! Variables (FC: built in firstcall in soil.F)
+real, save, allocatable, dimension(:,:,:) :: tsoil      ! sub-surface temperatures (K)
+real, save, allocatable, dimension(:,:,:) :: mthermdiff ! (FC) mid-layer thermal diffusivity
+real, save, allocatable, dimension(:,:,:) :: thermdiff  ! (FC) inter-layer thermal diffusivity
+real, save, allocatable, dimension(:)     :: coefq      ! (FC) q_{k+1/2} coefficients
+real, save, allocatable, dimension(:,:,:) :: coefd      ! (FC) d_k coefficients
+real, save, allocatable, dimension(:,:,:) :: alph       ! (FC) alpha_k coefficients
+real, save, allocatable, dimension(:,:,:) :: beta       ! beta_k coefficients
+real, save, allocatable, dimension(:,:)   :: flux_geo   ! Geothermal Flux (W/m^2)
+real, save :: mu
 
 !$OMP THREADPRIVATE(tsoil,mthermdiff,thermdiff,coefq,coefd,alph,beta,mu,flux_geo)
 
+! Subsurface tracers:
+logical, save                                  :: adsorption_soil ! boolean to call adosrption (or not)
+real,    save                                  :: choice_ads      ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
+real,    save, allocatable, dimension(:,:,:,:) :: qsoil           ! subsurface tracers (kg/m^3 of regol)
+integer, parameter :: nqsoil = 3                                  ! number of subsurface tracers, only three when working with water
+integer, parameter :: igcm_h2o_vap_soil = 1
+integer, parameter :: igcm_h2o_ice_soil = 2
+integer, parameter :: igcm_h2o_vap_ads  = 3
 
-! Subsurface tracers:
-  logical,save :: adsorption_soil             ! boolean to call adosrption (or not)
-  real,save :: choice_ads                     ! Choice for adsorption isotherm (3 means no adsorption, see soilwater.F90)
-  integer, parameter :: nqsoil = 3            ! number of subsurface tracers, only three when working with water
-  real,save,allocatable :: qsoil(:,:,:,:)     ! subsurface tracers (kg/m^3 of regol)
-  integer, parameter :: igcm_h2o_vap_soil = 1
-  integer, parameter :: igcm_h2o_ice_soil = 2
-  integer, parameter :: igcm_h2o_vap_ads  = 3
 !$OMP THREADPRIVATE(adsorption_soil,qsoil,choice_ads)
 
+!=======================================================================
+contains
+!=======================================================================
 
-contains
+subroutine ini_comsoil_h(ngrid,nslope)
 
-  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 for present climate
-    allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
-    allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
-    allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
-    allocate(thermdiff(ngrid,nsoilmx-1,nslope))
-    allocate(coefq(0:nsoilmx-1))
-    allocate(coefd(ngrid,nsoilmx-1,nslope))
-    allocate(alph(ngrid,nsoilmx-1,nslope))
-    allocate(beta(ngrid,nsoilmx-1,nslope))
-    allocate(flux_geo(ngrid,nslope))
-    allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
- 
-  end subroutine ini_comsoil_h
+implicit none
 
+integer,intent(in) :: ngrid  ! number of atmospheric columns
+integer,intent(in) :: nslope ! number of sub grid slopes
 
-  subroutine end_comsoil_h
+allocate(layer(nsoilmx)) !soil layer depths
+allocate(mlayer(0:nsoilmx - 1)) ! soil mid-layer depths
+allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia for present climate
+allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
+allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
+allocate(mthermdiff(ngrid,0:nsoilmx - 1,nslope))
+allocate(thermdiff(ngrid,nsoilmx - 1,nslope))
+allocate(coefq(0:nsoilmx - 1))
+allocate(coefd(ngrid,nsoilmx - 1,nslope))
+allocate(alph(ngrid,nsoilmx - 1,nslope))
+allocate(beta(ngrid,nsoilmx - 1,nslope))
+allocate(flux_geo(ngrid,nslope))
+allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
 
-  implicit none
+END SUBROUTINE ini_comsoil_h
 
-    if (allocated(layer)) deallocate(layer)
-    if (allocated(mlayer)) deallocate(mlayer)
-    if (allocated(inertiedat)) deallocate(inertiedat)
-    if (allocated(inertiesoil)) deallocate(inertiesoil)
-    if (allocated(tsoil)) deallocate(tsoil)
-    if (allocated(mthermdiff)) deallocate(mthermdiff)
-    if (allocated(thermdiff)) deallocate(thermdiff)
-    if (allocated(coefq)) deallocate(coefq) 
-    if (allocated(coefd)) deallocate(coefd)
-    if (allocated(alph)) deallocate(alph)
-    if (allocated(beta)) deallocate(beta)
-    if (allocated(flux_geo)) deallocate(flux_geo)
-    if (allocated(qsoil))  deallocate(qsoil)
-  end subroutine end_comsoil_h
+!=======================================================================
+SUBROUTINE end_comsoil_h
 
-  subroutine ini_comsoil_h_slope_var(ngrid,nslope)
-  
-  implicit none
-  integer,intent(in) :: ngrid ! number of atmospheric columns
-  integer,intent(in) :: nslope ! number of sub grid slopes
- 
-    allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
-    allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
-    allocate(mthermdiff(ngrid,0:nsoilmx-1,nslope))
-    allocate(thermdiff(ngrid,nsoilmx-1,nslope))
-    allocate(coefd(ngrid,nsoilmx-1,nslope))
-    allocate(alph(ngrid,nsoilmx-1,nslope))
-    allocate(beta(ngrid,nsoilmx-1,nslope))
-    allocate(flux_geo(ngrid,nslope))
-    allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
- 
-  end subroutine ini_comsoil_h_slope_var
+implicit none
 
+if (allocated(layer)) deallocate(layer)
+if (allocated(mlayer)) deallocate(mlayer)
+if (allocated(inertiedat)) deallocate(inertiedat)
+if (allocated(inertiesoil)) deallocate(inertiesoil)
+if (allocated(tsoil)) deallocate(tsoil)
+if (allocated(mthermdiff)) deallocate(mthermdiff)
+if (allocated(thermdiff)) deallocate(thermdiff)
+if (allocated(coefq)) deallocate(coefq)
+if (allocated(coefd)) deallocate(coefd)
+if (allocated(alph)) deallocate(alph)
+if (allocated(beta)) deallocate(beta)
+if (allocated(flux_geo)) deallocate(flux_geo)
+if (allocated(qsoil))  deallocate(qsoil)
 
-  subroutine end_comsoil_h_slope_var
+END SUBROUTINE end_comsoil_h
 
-  implicit none
+!=======================================================================
+SUBROUTINE ini_comsoil_h_slope_var(ngrid,nslope)
 
-    if (allocated(tsoil)) deallocate(tsoil)
-    if (allocated(inertiesoil)) deallocate(inertiesoil)
-    if (allocated(mthermdiff)) deallocate(mthermdiff)
-    if (allocated(thermdiff)) deallocate(thermdiff)
-    if (allocated(coefd)) deallocate(coefd)
-    if (allocated(alph)) deallocate(alph)
-    if (allocated(beta)) deallocate(beta)
-    if (allocated(flux_geo)) deallocate(flux_geo)
-    if (allocated(qsoil))  deallocate(qsoil)
+implicit none
 
-  end subroutine end_comsoil_h_slope_var
+integer,intent(in) :: ngrid  ! number of atmospheric columns
+integer,intent(in) :: nslope ! number of sub grid slopes
 
-end module comsoil_h
+allocate(tsoil(ngrid,nsoilmx,nslope)) ! soil temperatures
+allocate(inertiesoil(ngrid,nsoilmx,nslope)) ! soil thermal inertia
+allocate(mthermdiff(ngrid,0:nsoilmx - 1,nslope))
+allocate(thermdiff(ngrid,nsoilmx - 1,nslope))
+allocate(coefd(ngrid,nsoilmx - 1,nslope))
+allocate(alph(ngrid,nsoilmx - 1,nslope))
+allocate(beta(ngrid,nsoilmx - 1,nslope))
+allocate(flux_geo(ngrid,nslope))
+allocate(qsoil(ngrid,nsoilmx,nqsoil,nslope))
+
+END SUBROUTINE ini_comsoil_h_slope_var
+
+!=======================================================================
+SUBROUTINE end_comsoil_h_slope_var
+
+implicit none
+
+if (allocated(tsoil)) deallocate(tsoil)
+if (allocated(inertiesoil)) deallocate(inertiesoil)
+if (allocated(mthermdiff)) deallocate(mthermdiff)
+if (allocated(thermdiff)) deallocate(thermdiff)
+if (allocated(coefd)) deallocate(coefd)
+if (allocated(alph)) deallocate(alph)
+if (allocated(beta)) deallocate(beta)
+if (allocated(flux_geo)) deallocate(flux_geo)
+if (allocated(qsoil))  deallocate(qsoil)
+
+END SUBROUTINE end_comsoil_h_slope_var
+
+END MODULE comsoil_h
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/init_testphys1d_mod.F90	(revision 3183)
@@ -41,4 +41,5 @@
 use conf_phys_mod,            only: conf_phys
 use paleoclimate_mod,         only: paleoclimate
+use comslope_mod,             only: nslope, subslope_dist, ini_comslope_h, end_comslope_h
 ! Mostly for XIOS outputs:
 use mod_const_mpi,            only: COMM_LMDZ
@@ -208,5 +209,5 @@
 allocate(nqfils(nqtot))
 nqperes = 0
-nqfils(:) = 0
+nqfils = 0
 do iq = 1,nqtot
     if (tnom_transp(iq) == 'air') then
@@ -373,7 +374,6 @@
 longitude = longitude*pi/180.
 
-! Some initializations (some of which have already been
-! done above!) and loads parameters set in callphys.def
-! and allocates some arrays
+! Some initializations (some of which have already been done above!)
+! and loads parameters set in callphys.def and allocates some arrays
 ! Mars possible matter with dttestphys in input and include!!!
 ! Initializations below should mimick what is done in iniphysiq for 3D GCM
@@ -383,8 +383,7 @@
 call init_geometry_cell_area_for_outputs(1,cell_area)
 ! Ehouarn: init_vertial_layers called later (because disvert not called yet)
-!      call init_vertical_layers(nlayer,preff,scaleheight,
-!     &                      ap,bp,aps,bps,presnivs,pseudoalt)
+! call init_vertical_layers(nlayer,preff,scaleheight,ap,bp,aps,bps,presnivs,pseudoalt)
 call init_dimphy(1,nlayer) ! Initialize dimphy module
-call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils)! MVals: variables isotopes
+call phys_state_var_init(1,llm,nq,tname,day0,dayn,time,daysec,dttestphys,rad,g,r,cpp,nqperes,nqfils) ! MVals: variables isotopes
 call ini_fillgeom(1,latitude,longitude,(/1.0/))
 call conf_phys(1,llm,nq)
@@ -461,6 +460,6 @@
 
 ! For the non-orographic gravity waves scheme
-du_nonoro_gwd(1,:)=0
-dv_nonoro_gwd(1,:)=0
+du_nonoro_gwd(1,:) = 0
+dv_nonoro_gwd(1,:) = 0
 
 ! For the slope wind scheme
@@ -559,5 +558,5 @@
 if (.not. therestart1D) then
     tsurf(:,1) = tmp2(0)
-    temp(:) = tmp2(1:)
+    temp = tmp2(1:)
 else
     read(3,*) header, (tsurf(1,:), j = 1,size(tsurf,2)), (temp(ilayer), ilayer = 1,nlayer)
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F90	(revision 3183)
@@ -255,7 +255,7 @@
     ! Compute winds and temperature for next time step
     ! ------------------------------------------------
-    u(:) = u(:) + dttestphys*du(:)
-    v(:) = v(:) + dttestphys*dv(:)
-    temp(:) = temp(:) + dttestphys*dtemp(:)
+    u = u + dttestphys*du
+    v = v + dttestphys*dv
+    temp = temp + dttestphys*dtemp
 
     ! Compute pressure for next time step
Index: trunk/LMDZ.MARS/libf/phymars/getslopes.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/getslopes.F90	(revision 3179)
+++ 	(revision )
@@ -1,105 +1,0 @@
-subroutine getslopes(ngrid,geopot)
-    
-use geometry_mod, only: longitude, latitude ! in radians
-use slope_mod, only: theta_sl, psi_sl
-use comcstfi_h, only: g, rad, pi
-use mod_phys_lmdz_para, only: is_parallel
-use mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
-implicit none
-
-
-! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
-! It works fine with a non-regular grid for zoomed simulations.
-! slope inclination angle (deg)  0 == horizontal, 90 == vertical
-! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 == Southward, 270 == Westward
-! TN 04/1013
-
-integer,intent(in) :: ngrid ! nnumber of atmospheric columns
-real,intent(in) :: geopot(ngrid)     ! geopotential on phy grid
-real topogrid(nbp_lon,nbp_lat) ! topography on lat/lon grid with poles and only one -180/180 point
-real latigrid(nbp_lon,nbp_lat),longgrid(nbp_lon,nbp_lat) ! meshgrid of latitude and longitude values (radians)
-real theta_val ! slope inclination
-real psi_val   ! slope orientation
-real gradx(nbp_lon,nbp_lat) ! x: latitude-wise topography gradient,  increasing northward
-real grady(nbp_lon,nbp_lat) ! y: longitude-wise topography gradient, increasing westward
-integer i,j,ig0
-integer id2,idm1 ! a trick to compile testphys1d with debug option
-
-if (is_parallel) then
-  ! This routine only works in serial mode so stop now.
-  write(*,*) "getslopes Error: this routine is not designed to run in parallel"
-  call abort_physic("getslopes",'cannot be run in parallel',1)
-endif
-
-id2  = 2
-idm1 = nbp_lon-1
-
-! rearrange topography on a 2d array
-do j=2,nbp_lat-1
-   ig0= 1+(j-2)*nbp_lon
-   do i=1,nbp_lon
-      topogrid(i,j)=geopot(ig0+i)/g
-      latigrid(i,j)=latitude(ig0+i)
-      longgrid(i,j)=longitude(ig0+i)
-   enddo
-enddo
-!poles :
-topogrid(:,1) = geopot(1)/g
-latigrid(:,1) = latitude(1)
-longgrid(:,1) = longitude(1)
-topogrid(:,nbp_lat) = geopot(ngrid)/g
-latigrid(:,nbp_lat) = latitude(ngrid)
-longgrid(:,nbp_lat) = longitude(ngrid)
-
-
-
-! compute topography gradient
-! topogrid and rad are both in meters
-do j=2,nbp_lat-1
-   do i=1,nbp_lon
-     gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1))
-     gradx(i,j) = gradx(i,j) / rad
-   enddo
-   grady(1,j) = (topogrid(id2,j) - topogrid(nbp_lon,j)) / (2*pi+longgrid(id2,j)-longgrid(nbp_lon,j)) 
-   grady(1,j) = grady(1,j) / rad
-   grady(nbp_lon,j) = (topogrid(1,j) - topogrid(idm1,j)) / (2*pi+longgrid(1,j)-longgrid(idm1,j)) 
-   grady(nbp_lon,j) = grady(nbp_lon,j) / rad
-   do i=2,nbp_lon-1
-     grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j)) 
-     grady(i,j) = grady(i,j) / rad
-   enddo
-enddo
-! poles :
-gradx(:,1) = 0.
-grady(:,1) = 0.
-gradx(:,nbp_lat) = 0.
-grady(:,nbp_lat) = 0.
-
-
-
-! compute slope inclination and orientation :
-theta_sl(:) = 0.
-psi_sl(:)   = 0.
-do j=2,nbp_lat-1
-   do i=1,nbp_lon
-   
-     ig0= 1+(j-2)*nbp_lon
-   
-     theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 ))
-     
-     psi_val=0.
-     if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j))
-     if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi
-     psi_val = 3*pi/2. - psi_val
-     psi_val = psi_val*180./pi
-     psi_val = MODULO(psi_val,360.)
-     
-     theta_sl(ig0+i) = theta_val
-     psi_sl(ig0+i)   = psi_val
-          
-   enddo
-enddo
-
-
-
-end subroutine getslopes
Index: trunk/LMDZ.MARS/libf/phymars/param_slope.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/param_slope.F90	(revision 3179)
+++ 	(revision )
@@ -1,220 +1,0 @@
-SUBROUTINE param_slope( &
-  !
-  ! INPUTS
-  !
-               &  csza, declin, rho, latitude &
-               &  ,taudust, albedo  &
-               &  ,theta_s, psi_s &
-               &  ,fdir_0, ftot_0  &
-  !
-  ! OUTPUTS
-  !
-               &  ,ftot &
- )
-
-
-!!*****************************************************************************************
-!
-! SUBROUTINE:
-! param_slope
-! 
-!
-! PURPOSE: 
-! computes total solar irradiance on a given Martian slope 
-!
-!
-! INPUTS:
-! csza		cosine solar zenith angle
-! declin	sun declination (rad)
-! rho		sun right ascension (rad)
-! latitude	latitude (deg)
-! taudust	dust optical depth at reference wavelength 0.67 mic. 
-! albedo	spectrally integrated surface Lambertian reflection albedo 
-! theta_s	slope inclination angle (deg)
-!			0 is horizontal, 90 is vertical
-! phi_s		slope azimuth (deg) 
-!       		0 >> Northward
-!       		90 >> Eastward
-!       		180 >> Southward
-!       		270 >> Westward
-! ftot_0	spectrally integrated total irradiance on an horizontal surface (W/m2)
-! fdir_0	spectrally integrated direct irradiance on an horizontal surface (W/m2)
-!
-!
-! OUTPUTS:
-! ftot		spectrally integrated total irradiance on the slope (W/m2)
-!
-! REFERENCE: 
-! "Fast and accurate estimation of irradiance on Martian slopes"
-! A. Spiga & F. Forget
-! .....
-!
-! AUTHOR: 
-! A. Spiga (spiga@lmd.jussieu.fr)
-! March 2008
-!
-!!*****************************************************************************************
-
-IMPLICIT NONE
-
-!!
-!! INPUT
-!!
-REAL, INTENT(IN) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0 
-
-!!
-!! LOCAL
-!!
-REAL :: pi, deg2rad
-REAL :: a
-REAL :: mu_s, sigma_s
-REAL :: fdir, fscat, fscat_0, fref
-REAL, DIMENSION(4,2) :: mat_M, mat_N, mat_T
-REAL, DIMENSION(2) :: g_vector
-REAL, DIMENSION(4) :: s_vector
-REAL :: ratio
-
-!!
-!! OUTPUT
-!!
-REAL, INTENT(OUT) :: ftot
-
-!!*****************************************************************************************
-
-!
-! Prerequisite
-!
-  pi = 2.*asin(1.)
-  deg2rad = pi/180.
-  if ((theta_s > 90.) .or. (theta_s < 0.)) then
-        print *, 'please set theta_s between 0 and 90', theta_s
-        call abort_physic("param_slopes","invalid theta_s",1) 
-  endif
-
-!
-! Solar Zenith angle (radian) 
-!
-  if (csza .lt. 0.01) then
-      !print *, 'sun below horizon' 
-      !fdir_0=0.
-      fdir=0.
-      fscat_0=0.
-      fscat=0.   
-      fref=0.
-  else
-
-!!
-!! Low incidence fix
-!! 
-!  if (csza .lt. 0.15) csza = 0.15
-
-!
-! 'Slope vs Sun' azimuth (radian)
-!
-  if ( ( (cos(declin)*sin(rho)) .eq. 0.0 ) &
-               .and. &
-     ( ( sin(deg2rad*latitude)*cos(declin)*cos(rho) &
-              -cos(deg2rad*latitude)*sin(declin) ) .eq. 0.0 ) ) then
-    a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0)  
-  else
-    a = deg2rad*psi_s + atan2(cos(declin)*sin(rho), &
-               sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
-  end if
-  
-!
-! Cosine of slope-sun phase angle 
-!
-  mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
-  if (mu_s .le. 0.) mu_s=0.
-
-!
-! Sky-view factor
-!
-  sigma_s=0.5*(1.+cos(deg2rad*theta_s))
-
-!
-! Direct flux on the slope
-!
-  fdir = fdir_0 * mu_s/csza
-
-!
-! Reflected flux on the slope
-!
-  fref = albedo * (1-sigma_s) * ftot_0
-
-!
-! Scattered flux on a flat surface
-!
-  fscat_0 = ftot_0 - fdir_0
-
-!
-! Scattering vector (slope vs sky)
-!
-  s_vector=(/ 1., exp(-taudust) , sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)
-
-!
-! Geometry vector (slope vs sun)
-!
-  g_vector=(/ mu_s/csza, 1. /)
-
-!
-! Coupling matrix
-!
-  if (csza .ge. 0.5) then
-        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
-        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
-        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
-        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
-
-  else
-        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
-        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
-        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
-        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
-  endif
-  !
-  mat_T = mat_M + csza*mat_N
-
-
-!
-! Scattered flux slope ratio
-!
-  if (deg2rad*theta_s <= 0.0872664626) then
-  !
-  ! low angles
-  !
-        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
-        ratio = DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector )
-        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
-  else
-  !
-  ! general case
-  !
-        ratio= DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector )  
-                  !
-                  ! NB: ratio= DOT_PRODUCT ( s_vector, MATMUL( mat_T, g_vector ) ) is equivalent
-  endif
-
-!
-! Scattered flux on the slope
-!
-  fscat = ratio * fscat_0
-
-
-endif !! if (csza < 0.01)
-
-!
-! Total flux on the slope
-!
-  ftot = fdir + fref + fscat
-
-!!
-!! Display results
-!!
-!  print *, 'sca component 0 ', ftot_0-fdir_0
-!  print *, 'dir component 0 ', fdir_0
-!  print *, 'scattered component ', fscat
-!  print *, 'direct component ', fdir
-!  print *, 'reflected component ', fref
-
-END SUBROUTINE param_slope
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0_mod.F90	(revision 3183)
@@ -146,58 +146,55 @@
 
 if (startphy_file) then
- call get_var("def_slope",def_slope,found)
-   if(.not.found) then 
-     startphy_slope=.false.
-     write(*,*)'slope_settings: Problem while reading <def_slope>'
-     write(*,*)'Checking that nslope=1'
-     if(nslope.eq.1) then
-      write(*,*)'We take default def_slope and subslope dist'
-      allocate(default_def_slope(nslope+1))
-      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'
-       write(*,*)'you are trying to run with nslope!=1'
-       write(*,*)'This is not possible, you have to go through newstart'
-     endif
+   call get_var("def_slope",def_slope,found)
+   if (.not. found) then 
+       startphy_slope=.false.
+       write(*,*)'slope_settings: Problem while reading <def_slope>'
+       write(*,*)'Checking that nslope=1'
+       if (nslope == 1) then
+           write(*,*)'We take default def_slope and subslope dist'
+           allocate(default_def_slope(nslope + 1))
+           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'
+           write(*,*)'you are trying to run with nslope!=1'
+           write(*,*)'This is not possible, you have to go through newstart'
+       endif
    else
-     startphy_slope=.true.
-     call get_field("subslope_dist",subslope_dist,found,indextime)
-     if(.not.found) then 
-       write(*,*)'slope_settings: Problem while reading <subslope_dist>'
-       write(*,*)'We have to abort.'
-       write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
-       call abort_physic(modname, &
-                "phyetat0: Failed loading <subslope_dist>",1)
-     endif
+       startphy_slope=.true.
+       call get_field("subslope_dist",subslope_dist,found,indextime)
+       if (.not. found) then 
+           write(*,*)'slope_settings: Problem while reading <subslope_dist>'
+           write(*,*)'We have to abort.'
+           write(*,*)'Please check that nslope is coherent, you can modify this parameter with newstart'
+           call abort_physic(modname,"phyetat0: Failed loading <subslope_dist>",1)
+       endif
+   endif
+else ! (no startphy_file)
+    if (nslope == 1) then
+        allocate(default_def_slope(2))
+        default_def_slope = 0.
+        subslope_dist = 1.
+        def_slope = default_def_slope
+    else
+        write(*,*)'Without startfi, lets first run with nslope=1'
+        call abort_physic(modname,"phyetat0: No startfi and nslope!=1",1)
     endif
-else ! (no startphy_file)
-     if(nslope.eq.1) then
-      allocate(default_def_slope(nslope+1))
-      default_def_slope(1) = 0.
-      default_def_slope(2) = 0.
-      subslope_dist(:,:)=1.
-      def_slope(:)=default_def_slope(:)
-    else
-     write(*,*)'Without starfi, lets first run with nslope=1'
-       call abort_physic(modname, &
-                "phyetat0: No startfi and nslope!=1",1)
-   endif
 endif   
 
-do islope=1,nslope
-    def_slope_mean(islope) =(def_slope(islope)+def_slope(islope+1))/2.
+do islope = 1,nslope
+    def_slope_mean(islope) = (def_slope(islope) + def_slope(islope + 1))/2.
 enddo
  
 DO ig = 1,ngrid
-  sum_dist = 0.
-  DO islope = 1,nslope
-     sum_dist = sum_dist + subslope_dist(ig,islope)
-  ENDDO
-  DO islope = 1,nslope
-    subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
-  ENDDO
+    sum_dist = 0.
+    DO islope = 1,nslope
+        sum_dist = sum_dist + subslope_dist(ig,islope)
+    ENDDO
+    DO islope = 1,nslope
+        subslope_dist(ig,islope) = subslope_dist(ig,islope)/sum_dist
+    ENDDO
 ENDDO
 
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 3183)
@@ -65,5 +65,5 @@
       use nirdata_mod, only: NIR_leedat
       use nirco2abs_mod, only: nirco2abs
-      use slope_mod, only: theta_sl, psi_sl
+      use slope_mod, only: theta_sl, psi_sl, getslopes, param_slope
       use conc_mod, only: rnew, cpnew, mmean
       use time_phylmdz_mod, only: iphysiq, day_step, ecritstart, daysec
Index: trunk/LMDZ.MARS/libf/phymars/slope_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/slope_mod.F90	(revision 3179)
+++ trunk/LMDZ.MARS/libf/phymars/slope_mod.F90	(revision 3183)
@@ -1,32 +1,292 @@
-module slope_mod
-
-implicit none
-
-  real,save,allocatable :: theta_sl(:) ! slope angle versus horizontal (deg)
-  real,save,allocatable :: psi_sl(:)   ! slope orientation (deg)
+MODULE slope_mod
+
+implicit none
+
+real, save, dimension(:), allocatable :: theta_sl ! slope angle versus horizontal (deg)
+real, save, dimension(:), allocatable :: psi_sl   ! slope orientation (deg)
 
 !$OMP THREADPRIVATE(theta_sl,psi_sl)
 
+!=======================================================================
 contains
-
-  subroutine ini_slope_mod(ngrid)
-  
-  implicit none
-  integer,intent(in) :: ngrid ! number of atmospheric columns
-  
-  allocate(theta_sl(ngrid))
-  allocate(psi_sl(ngrid))
-  
-  end subroutine ini_slope_mod
-
-
-  subroutine end_slope_mod
-
-  implicit none
-
-  if (allocated(theta_sl)) deallocate(theta_sl)
-  if (allocated(psi_sl)) deallocate(psi_sl)
-
-  end subroutine end_slope_mod
-  
-end module slope_mod
+!=======================================================================
+
+SUBROUTINE getslopes(ngrid,geopot)
+
+use geometry_mod,       only: longitude, latitude ! in radians
+use comcstfi_h,         only: g, rad, pi
+use mod_phys_lmdz_para, only: is_parallel
+use mod_grid_phy_lmdz,  only: nbp_lon, nbp_lat
+
+implicit none
+
+! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
+! It works fine with a non-regular grid for zoomed simulations.
+! slope inclination angle (deg)  0 == horizontal, 90 == vertical
+! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 == Southward, 270 == Westward
+! TN 04/1013
+
+! Input arguments
+integer,                intent(in) :: ngrid  ! nnumber of atmospheric columns
+real, dimension(ngrid), intent(in) :: geopot ! geopotential on phy grid
+
+! Local variables
+real, dimension(nbp_lon,nbp_lat) :: topogrid           ! topography on lat/lon grid with poles and only one -180/180 point
+real, dimension(nbp_lon,nbp_lat) :: latigrid, longgrid ! meshgrid of latitude and longitude values (radians)
+real, dimension(nbp_lon,nbp_lat) :: gradx              ! x: latitude-wise topography gradient,  increasing northward
+real, dimension(nbp_lon,nbp_lat) :: grady              ! y: longitude-wise topography gradient, increasing westward
+real                             :: theta_val          ! slope inclination
+real                             :: psi_val            ! slope orientation
+integer                          :: i, j, ig0
+integer                          :: id2, idm1          ! a trick to compile testphys1d with debug option
+
+if (is_parallel) then
+    ! This routine only works in serial mode so stop now.
+    write(*,*) "getslopes Error: this routine is not designed to run in parallel"
+    call abort_physic("getslopes",'cannot be run in parallel',1)
+endif
+
+id2  = 2
+idm1 = nbp_lon-1
+
+! rearrange topography on a 2d array
+do j = 2,nbp_lat-1
+    ig0 = 1 + (j - 2)*nbp_lon
+    do i = 1,nbp_lon
+        topogrid(i,j) = geopot(ig0 + i)/g
+        latigrid(i,j) = latitude(ig0 + i)
+        longgrid(i,j) = longitude(ig0 + i)
+    enddo
+enddo
+
+! poles:
+topogrid(:,1) = geopot(1)/g
+latigrid(:,1) = latitude(1)
+longgrid(:,1) = longitude(1)
+topogrid(:,nbp_lat) = geopot(ngrid)/g
+latigrid(:,nbp_lat) = latitude(ngrid)
+longgrid(:,nbp_lat) = longitude(ngrid)
+
+! compute topography gradient
+! topogrid and rad are both in meters
+do j = 2,nbp_lat - 1
+    do i=1,nbp_lon
+        gradx(i,j) = (topogrid(i,j + 1) - topogrid(i,j - 1))/(latigrid(i,j + 1)-latigrid(i,j - 1))
+        gradx(i,j) = gradx(i,j)/rad
+    enddo
+    grady(1,j) = (topogrid(id2,j) - topogrid(nbp_lon,j))/(2*pi + longgrid(id2,j) - longgrid(nbp_lon,j))
+    grady(1,j) = grady(1,j) / rad
+    grady(nbp_lon,j) = (topogrid(1,j) - topogrid(idm1,j))/(2*pi + longgrid(1,j) - longgrid(idm1,j))
+    grady(nbp_lon,j) = grady(nbp_lon,j)/rad
+    do i = 2,nbp_lon - 1
+        grady(i,j) = (topogrid(i + 1,j) - topogrid(i-1,j))/(longgrid(i + 1,j) - longgrid(i - 1,j))
+        grady(i,j) = grady(i,j)/rad
+    enddo
+enddo
+
+! poles:
+gradx(:,1) = 0.
+grady(:,1) = 0.
+gradx(:,nbp_lat) = 0.
+grady(:,nbp_lat) = 0.
+
+! compute slope inclination and orientation:
+theta_sl = 0.
+psi_sl  = 0.
+do j = 2,nbp_lat - 1
+    do i = 1,nbp_lon
+        ig0 = 1 + (j - 2)*nbp_lon
+
+        theta_val = atan(sqrt((gradx(i,j))**2 + (grady(i,j))**2))
+
+        psi_val = 0.
+        if (gradx(i,j) /= 0.) psi_val = -pi/2. - atan(grady(i,j)/gradx(i,j))
+        if (gradx(i,j) >= 0.) psi_val = psi_val - pi
+        psi_val = 3*pi/2. - psi_val
+        psi_val = psi_val*180./pi
+        psi_val = modulo(psi_val,360.)
+
+        theta_sl(ig0 + i) = theta_val
+        psi_sl(ig0 + i)   = psi_val
+    enddo
+enddo
+
+end subroutine getslopes
+
+!=======================================================================
+
+SUBROUTINE param_slope(csza,declin,rho,latitude,taudust,albedo,theta_s,psi_s,fdir_0,ftot_0,ftot)
+!***********************************************************************
+!
+! SUBROUTINE:
+! param_slope
+! 
+!
+! PURPOSE: 
+! computes total solar irradiance on a given Martian slope 
+!
+!
+! INPUTS:
+! csza		cosine solar zenith angle
+! declin	sun declination (rad)
+! rho		sun right ascension (rad)
+! latitude	latitude (deg)
+! taudust	dust optical depth at reference wavelength 0.67 mic. 
+! albedo	spectrally integrated surface Lambertian reflection albedo 
+! theta_s	slope inclination angle (deg)
+!			0 is horizontal, 90 is vertical
+! phi_s		slope azimuth (deg) 
+!       		0 >> Northward
+!       		90 >> Eastward
+!       		180 >> Southward
+!       		270 >> Westward
+! ftot_0	spectrally integrated total irradiance on an horizontal surface (W/m2)
+! fdir_0	spectrally integrated direct irradiance on an horizontal surface (W/m2)
+!
+!
+! OUTPUTS:
+! ftot		spectrally integrated total irradiance on the slope (W/m2)
+!
+! REFERENCE: 
+! "Fast and accurate estimation of irradiance on Martian slopes"
+! A. Spiga & F. Forget
+! .....
+!
+! AUTHOR: 
+! A. Spiga (spiga@lmd.jussieu.fr)
+! March 2008
+!
+!***********************************************************************
+
+use comcstfi_h, only: pi
+
+implicit none
+
+! Input arguments
+real, intent(in) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0
+
+! Output arguments
+real, intent(out) :: ftot
+
+! Local variables
+real                 :: deg2rad, a, mu_s, sigma_s, fdir, fscat, fscat_0, fref, ratio
+real, dimension(4,2) :: mat_M, mat_N, mat_T
+real, dimension(2)   :: g_vector
+real, dimension(4)   :: s_vector
+!***********************************************************************
+! Prerequisite
+deg2rad = pi/180.
+if ((theta_s > 90.) .or. (theta_s < 0.)) then
+    write(*,*) 'please set theta_s between 0 and 90', theta_s
+    call abort_physic("param_slope","invalid theta_s",1) 
+endif
+
+! Solar Zenith angle (radian)
+if (csza < 0.01) then
+    !print *, 'sun below horizon' 
+    !fdir_0=0.
+    fdir    = 0.
+    fscat_0 = 0.
+    fscat   = 0.   
+    fref    = 0.
+else
+! Low incidence fix 
+!    if (csza < 0.15) csza = 0.15
+
+! 'Slope vs Sun' azimuth (radian)
+    if (cos(declin)*sin(rho) == 0. .and. &
+        sin(deg2rad*latitude)*cos(declin)*cos(rho) - cos(deg2rad*latitude)*sin(declin) == 0.) then
+        a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0)  
+    else
+        a = deg2rad*psi_s + atan2(cos(declin)*sin(rho),sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
+    endif
+
+! Cosine of slope-sun phase angle 
+    mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
+    if (mu_s <= 0.) mu_s = 0.
+
+! Sky-view factor
+    sigma_s=0.5*(1. + cos(deg2rad*theta_s))
+
+! Direct flux on the slope
+    fdir = fdir_0*mu_s/csza
+
+! Reflected flux on the slope
+    fref = albedo*(1 - sigma_s)*ftot_0
+
+! Scattered flux on a flat surface
+    fscat_0 = ftot_0 - fdir_0
+
+! Scattering vector (slope vs sky)
+    s_vector = (/ 1., exp(-taudust), sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)
+
+! Geometry vector (slope vs sun)
+    g_vector = (/ mu_s/csza, 1. /)
+
+! Coupling matrix
+    if (csza >= 0.5) then
+        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
+        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
+        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
+        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
+    else
+        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
+        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
+        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
+        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
+    endif
+    mat_T = mat_M + csza*mat_N
+
+! Scattered flux slope ratio
+    if (deg2rad*theta_s <= 0.0872664626) then
+    ! low angles
+        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
+        ratio = dot_product(matmul(s_vector, mat_T),g_vector)
+        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
+    else
+    ! general case
+    ratio = dot_product(matmul(s_vector,mat_T),g_vector)  
+    ! NB: ratio = dot_product(s_vector,matmul(mat_T,g_vector)) is equivalent
+    endif
+
+! Scattered flux on the slope
+    fscat = ratio*fscat_0
+endif ! if (csza < 0.01)
+
+! Total flux on the slope
+ftot = fdir + fref + fscat
+
+! Display results
+!  print *, 'sca component 0 ', ftot_0-fdir_0
+!  print *, 'dir component 0 ', fdir_0
+!  print *, 'scattered component ', fscat
+!  print *, 'direct component ', fdir
+!  print *, 'reflected component ', fref
+
+END SUBROUTINE param_slope
+
+!=======================================================================
+
+SUBROUTINE ini_slope_mod(ngrid)
+
+implicit none
+
+integer, intent(in) :: ngrid ! number of atmospheric columns
+
+allocate(theta_sl(ngrid))
+allocate(psi_sl(ngrid))
+
+END SUBROUTINE ini_slope_mod
+
+!=======================================================================
+
+SUBROUTINE end_slope_mod
+
+implicit none
+
+if (allocated(theta_sl)) deallocate(theta_sl)
+if (allocated(psi_sl)) deallocate(psi_sl)
+
+END SUBROUTINE end_slope_mod
+
+END MODULE slope_mod
