Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90	(revision 6594)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_soil_var.f90	(revision 6954)
@@ -17,4 +17,7 @@
 !! RECENT CHANGE(S): AD: mcw and mcf depend now on soil texture, based on Van Genuchten equations 
 !!                   and classical matric potential values, and pcent is adapted
+!!                   November 2020 by Salma Tafasca and Agnes Ducharne : we introduce a new texture class
+!!                   for clay oxisols (cf. Tafasca, 2020, PhD thesis; Tafasca et al., in prep for GRL).
+!!                   It makes no change if we read a soil texture map with only 12 USDA classes.
 !!
 !! REFERENCE(S)	:
@@ -30,8 +33,12 @@
 !!- Carsel, R. and Parrish, R.: Developing joint probability distributions of soil water
 !! retention characteristics, Water Resour. Res.,24, 755â769, 1988.
-!!- Mualem Y (1976) A new model for predicting the hydraulic conductivity  
+!!- Mualem Y (1976). A new model for predicting the hydraulic conductivity  
 !! of unsaturated porous media. Water Resources Research 12(3):513-522
-!!- Van Genuchten M (1980) A closed-form equation for predicting the  
+!!- Van Genuchten M (1980). A closed-form equation for predicting the  
 !! hydraulic conductivity of unsaturated soils. Soil Sci Soc Am J, 44(5):892-898
+!!- Tafasca S. (2020). Evaluation de lâimpact des propriÃ©tÃ©s du sol sur lâhydrologie simulee dans le
+!! modÃšle ORCHIDEE, PhD thesis, Sorbonne Universite. 
+!!- Tafasca S., Ducharne A. and Valentin C. Accounting for soil structure in pedo-transfer functions:
+!!  swelling vs non swelling clays. In prep for GRL.
 !!
 !! SVN          :
@@ -61,5 +68,5 @@
 !$OMP THREADPRIVATE(soil_classif)
   INTEGER(i_std), PARAMETER :: nscm_fao=3               !! For FAO Classification (unitless)
-  INTEGER(i_std), PARAMETER :: nscm_usda=12             !! For USDA Classification (unitless)
+  INTEGER(i_std), PARAMETER :: nscm_usda=13             !! For USDA Classification (unitless)
   INTEGER(i_std), SAVE      :: nscm=nscm_fao            !! Default value for nscm
 !$OMP THREADPRIVATE(nscm)
@@ -171,5 +178,5 @@
 
   REAL(r_std),DIMENSION(nscm_usda),SAVE :: soilclass_default_usda = &    !! Default soil texture distribution in the above order :
- & (/ 0.28, 0.52, 0.20, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)   !! Thus different from "FAO"'s COARSE, MEDIUM, FINE
+ & (/ 0.28, 0.52, 0.20, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)   !! Thus different from "FAO"'s COARSE, MEDIUM, FINE
                                                                          !! which have indices 3,6,9 in the 12-texture vector
   !$OMP THREADPRIVATE(soilclass_default_usda)
@@ -178,70 +185,84 @@
  & (/ 2.68_r_std, 2.28_r_std, 1.89_r_std, 1.41_r_std, &                   !  RK: 1/n=1-m
  &    1.37_r_std, 1.56_r_std, 1.48_r_std, 1.23_r_std, &
- &    1.31_r_std, 1.23_r_std, 1.09_r_std, 1.09_r_std /)
+ &    1.31_r_std, 1.23_r_std, 1.09_r_std, 1.09_r_std, &
+ &    1.552_r_std    /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: avan_usda = &            !! Van Genuchten coefficient a 
  & (/ 0.0145_r_std, 0.0124_r_std, 0.0075_r_std, 0.0020_r_std, &          !!  @tex $(mm^{-1})$ @endtex
  &    0.0016_r_std, 0.0036_r_std, 0.0059_r_std, 0.0010_r_std, &
- &    0.0019_r_std, 0.0027_r_std, 0.0005_r_std, 0.0008_r_std /)
+ &    0.0019_r_std, 0.0027_r_std, 0.0005_r_std, 0.0008_r_std, &
+ &    0.0132_r_std /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcr_usda = &             !! Residual volumetric water content 
  & (/ 0.045_r_std, 0.057_r_std, 0.065_r_std, 0.067_r_std, &              !!  @tex $(m^{3} m^{-3})$ @endtex
  &    0.034_r_std, 0.078_r_std, 0.100_r_std, 0.089_r_std, &
- &    0.095_r_std, 0.100_r_std, 0.070_r_std, 0.068_r_std /)
+ &    0.095_r_std, 0.100_r_std, 0.070_r_std, 0.068_r_std, &
+ &    0.068_r_std /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcs_usda = &             !! Saturated volumetric water content 
  & (/ 0.43_r_std, 0.41_r_std, 0.41_r_std, 0.45_r_std, &                  !!  @tex $(m^{3} m^{-3})$ @endtex
  &    0.46_r_std, 0.43_r_std, 0.39_r_std, 0.43_r_std, &
- &    0.41_r_std, 0.38_r_std, 0.36_r_std, 0.38_r_std /)
+ &    0.41_r_std, 0.38_r_std, 0.36_r_std, 0.38_r_std, &
+ &    0.503_r_std  /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: ks_usda = &              !! Hydraulic conductivity at saturation
  & (/ 7128.0_r_std, 3501.6_r_std, 1060.8_r_std, 108.0_r_std, &           !!  @tex $(mm d^{-1})$ @endtex
  &    60.0_r_std, 249.6_r_std, 314.4_r_std, 16.8_r_std, &
- &    62.4_r_std, 28.8_r_std, 4.8_r_std, 48.0_r_std /)
+ &    62.4_r_std, 28.8_r_std, 4.8_r_std, 48.0_r_std, &
+ &    6131.4_r_std  /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: pcent_usda = &           !! Fraction of saturated volumetric soil moisture
  & (/ 0.8_r_std, 0.8_r_std, 0.8_r_std, 0.8_r_std, &                      !! above which transpir is max (0-1, unitless)
  &    0.8_r_std, 0.8_r_std, 0.8_r_std, 0.8_r_std, &
- &    0.8_r_std, 0.8_r_std, 0.8_r_std, 0.8_r_std /)
+ &    0.8_r_std, 0.8_r_std, 0.8_r_std, 0.8_r_std, &
+ &    0.8_r_std /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: free_drain_max_usda = &  !! Max=default value of the permeability coeff 
  & (/ 1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std, &                      !! at the bottom of the soil (0-1, unitless)
  &    1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std, &
- &    1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std /)
+ &    1.0_r_std, 1.0_r_std, 1.0_r_std, 1.0_r_std,  &
+ &    1.0_r_std /)
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcf_usda = &             !! Volumetric water content at field capacity
  & (/ 0.0493_r_std, 0.0710_r_std, 0.1218_r_std, 0.2402_r_std, &          !!  @tex $(m^{3} m^{-3})$ @endtex
       0.2582_r_std, 0.1654_r_std, 0.1695_r_std, 0.3383_r_std, &
-      0.2697_r_std, 0.2672_r_std, 0.3370_r_std, 0.3469_r_std /)
+      0.2697_r_std, 0.2672_r_std, 0.3370_r_std, 0.3469_r_std, &
+      0.172_r_std  /) ! oxisols
   
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mcw_usda = &             !! Volumetric water content at wilting point
  & (/ 0.0450_r_std, 0.0570_r_std, 0.0657_r_std, 0.1039_r_std, &          !!  @tex $(m^{3} m^{-3})$ @endtex
       0.0901_r_std, 0.0884_r_std, 0.1112_r_std, 0.1967_r_std, &
-      0.1496_r_std, 0.1704_r_std, 0.2665_r_std, 0.2707_r_std /)
+      0.1496_r_std, 0.1704_r_std, 0.2665_r_std, 0.2707_r_std, &
+      0.075_r_std  /) ! oxisols
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mc_awet_usda = &         !! Vol. wat. cont. above which albedo is cst
  & (/ 0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std, &                  !!  @tex $(m^{3} m^{-3})$ @endtex
  &    0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std, &
- &    0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std /)
+ &    0.25_r_std, 0.25_r_std, 0.25_r_std, 0.25_r_std, &
+ &    0.25_r_std /)
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: mc_adry_usda = &         !! Vol. wat. cont. below which albedo is cst
  & (/ 0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std, &                      !!  @tex $(m^{3} m^{-3})$ @endtex
  &    0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std, &
- &    0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std /)
-
-  REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: SMCMAX_usda = &          !! porosity
+ &    0.1_r_std, 0.1_r_std, 0.1_r_std, 0.1_r_std, &
+ &    0.1_r_std /) ! oxisols
+
+  REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: SMCMAX_usda = &          !! porosity for thermosoil
  & (/ 0.43_r_std, 0.41_r_std, 0.41_r_std, 0.45_r_std, &
  &    0.46_r_std, 0.43_r_std, 0.39_r_std, 0.43_r_std, &
- &    0.41_r_std, 0.38_r_std, 0.36_r_std, 0.38_r_std /)
+ &    0.41_r_std, 0.38_r_std, 0.36_r_std, 0.38_r_std, &
+ &    0.503_r_std  /) ! oxisols
  
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: QZ_usda = &              !! QUARTZ CONTENT (SOIL TYPE DEPENDENT)
- & (/ 0.92_r_std, 0.82_r_std, 0.60_r_std, 0.25_r_std, &
+ & (/ 0.92_r_std, 0.82_r_std, 0.60_r_std, 0.25_r_std, &                  !! Peters et al [1998]
  &    0.10_r_std, 0.40_r_std, 0.60_r_std, 0.10_r_std, &
- &    0.35_r_std, 0.52_r_std, 0.10_r_std, 0.25_r_std /)                  !! Peters et al [1998]
+ &    0.35_r_std, 0.52_r_std, 0.10_r_std, 0.25_r_std, &
+&     0.25_r_std /)  ! oxisols                 
 
   REAL(r_std),PARAMETER,DIMENSION(nscm_usda) :: so_capa_dry_ns_usda = &  !! Dry soil Heat capacity of soils,J.m^{-3}.K^{-1}
- & (/ 1.47e+6_r_std, 1.41e+6_r_std, 1.34e+6_r_std, 1.27e+6_r_std, &
+ & (/ 1.47e+6_r_std, 1.41e+6_r_std, 1.34e+6_r_std, 1.27e+6_r_std, &      !! Pielke [2002, 2013]
  &    1.21e+6_r_std, 1.21e+6_r_std, 1.18e+6_r_std, 1.32e+6_r_std, &
- &    1.23e+6_r_std, 1.18e+6_r_std, 1.15e+6_r_std, 1.09e+6_r_std /)      !! Pielke [2002, 2013]
+ &    1.23e+6_r_std, 1.18e+6_r_std, 1.15e+6_r_std, 1.09e+6_r_std, &
+ &    1.09e+6_r_std /) ! oxisols     
   
   !! Parameters for the numerical scheme used by CWRR
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 6594)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 6954)
@@ -19,6 +19,10 @@
 !!                 This routine was originaly developed by Patricia deRosnay.
 !!
-!! RECENT CHANGE(S) : None
-!!
+!! RECENT CHANGE(S) : November 2020: It is possible to define soil hydraulic parameters from maps,
+!!                    as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
+!!                    Here, it leads to change dimensions and indices. 
+!!                    We can also impose kfact_root=1 in all soil layers to cancel the effect of
+!!                    roots on ks profile (keyword KFACT_ROOT_TYPE).
+!!                 
 !! REFERENCE(S) :
 !! - de Rosnay, P., J. Polcher, M. Bruen, and K. Laval, Impact of a physically based soil
@@ -43,4 +47,6 @@
 !! of a new soil freezing scheme for a land-surface model with physically-based hydrology.
 !! The Cryosphere, 6, 407-430, doi: 10.5194/tc-6-407-2012. \n
+!! - Tafasca S. (2020). Evaluation de lâimpact des propriÃ©tÃ©s du sol sur lâhydrologie simulee dans le
+!! modÃšle ORCHIDEE, PhD thesis, Sorbonne Universite. \n
 !!
 !! SVN          :
@@ -91,28 +97,8 @@
   ! one dimension array allocated, computed, saved and got in hydrol module
   ! Values per soil type
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: nvan                !! Van Genuchten coeficients n (unitless)
-                                                                          ! RK: 1/n=1-m
-!$OMP THREADPRIVATE(nvan)                                                 
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: avan                !! Van Genuchten coeficients a
-                                                                         !!  @tex $(mm^{-1})$ @endtex
-!$OMP THREADPRIVATE(avan)                                                
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcr                 !! Residual volumetric water content
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
-!$OMP THREADPRIVATE(mcr)                                                 
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcs                 !! Saturated volumetric water content
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
-!$OMP THREADPRIVATE(mcs)                                                  
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: ks                  !! Hydraulic conductivity at saturation
-                                                                         !!  @tex $(mm d^{-1})$ @endtex
-!$OMP THREADPRIVATE(ks)                                                  
+                                                              
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pcent               !! Fraction of saturated volumetric soil moisture above 
                                                                          !! which transpir is max (0-1, unitless)
-!$OMP THREADPRIVATE(pcent)
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcfc                !! Volumetric water content at field capacity
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex 
-!$OMP THREADPRIVATE(mcfc)                                                 
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcw                 !! Volumetric water content at wilting point
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex 
-!$OMP THREADPRIVATE(mcw)                                                 
+!$OMP THREADPRIVATE(pcent)                                                      
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mc_awet             !! Vol. wat. cont. above which albedo is cst
                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex 
@@ -216,5 +202,5 @@
 !$OMP THREADPRIVATE(kfact_root)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: kfact            !! Factor to reduce Ks with depth (unitless)
-                                                                         !! DIM = nslm * nscm
+                                                                         !! DIM = nslm * kjpindex
 !$OMP THREADPRIVATE(kfact)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: zz               !! Depth of nodes [znh in vertical_soil] transformed into (mm) 
@@ -228,21 +214,21 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: mc_lin   !! 50 Vol. Wat. Contents to linearize K and D, for each texture 
                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex
-                                                                 !! DIM = imin:imax * nscm
+                                                                 !! DIM = imin:imax * kjpindex
 !$OMP THREADPRIVATE(mc_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: k_lin    !! 50 values of unsaturated K, for each soil layer and texture
                                                                  !!  @tex $(mm d^{-1})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(k_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: d_lin    !! 50 values of diffusivity D, for each soil layer and texture
                                                                  !!  @tex $(mm^2 d^{-1})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(d_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: a_lin    !! 50 values of the slope in K=a*mc+b, for each soil layer and texture
                                                                  !!  @tex $(mm d^{-1})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(a_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: b_lin    !! 50 values of y-intercept in K=a*mc+b, for each soil layer and texture
                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(b_lin)
 
@@ -259,8 +245,8 @@
 !$OMP THREADPRIVATE(kk)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: avan_mod_tab  !! VG parameter a modified from  exponantial profile
-                                                                      !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,nscm)
+                                                                      !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,kjpindex)
 !$OMP THREADPRIVATE(avan_mod_tab)  
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: nvan_mod_tab  !! VG parameter n  modified from  exponantial profile
-                                                                      !! (unitless) !! DIMENSION (nslm,nscm)  
+                                                                      !! (unitless) !! DIMENSION (nslm,kjpindex)  
 !$OMP THREADPRIVATE(nvan_mod_tab)
   
@@ -450,5 +436,7 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          &
+  SUBROUTINE hydrol_initialize ( ks,             nvan,      avan,          mcr,              &
+                                 mcs,            mcfc,      mcw,           kjit,             &
+                                 kjpindex,       index,     rest_id,                         &
                                  njsc,           soiltile,  veget,         veget_max,        &
                                  humrel,         vegstress, drysoil_frac,                    &
@@ -461,4 +449,16 @@
     !! 0. Variable and parameter declaration
     !! 0.1 Input variables
+
+
+    !salma: added soil params as input variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+    !end salma: added soil params as input variables
+
     INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number 
     INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
@@ -495,12 +495,25 @@
     !! 0.4 Local variables
     INTEGER(i_std)                                       :: jsl
+    
+    !salma: added soil params which depend on soil texture
+    REAL(r_std), DIMENSION (nscm)                        :: nvan_text
+    REAL(r_std), DIMENSION (nscm)                        :: avan_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcr_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcs_text
+    REAL(r_std), DIMENSION (nscm)                        :: ks_text
+    REAL(r_std), DIMENSION (nscm)                        :: pcent_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcfc_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcw_text
+    REAL(r_std), DIMENSION (nscm)                        :: mc_awet_text
+    REAL(r_std), DIMENSION (nscm)                        :: mc_adry_text
+     !end salma: added soil params which depend on soil texture
 !_ ================================================================================================================================
 
-    CALL hydrol_init (kjit, kjpindex, index, rest_id, veget_max, soiltile, &
+    CALL hydrol_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc, kjit, kjpindex, index, rest_id, veget_max, soiltile, &
          humrel, vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
          snowdz, snowgrain, snowrho,    snowtemp,   snowheat, &
          drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
     
-    CALL hydrol_var_init (kjpindex, veget, veget_max, &
+    CALL hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget, veget_max, &
          soiltile, njsc, mx_eau_var, shumdiag_perma, &
          drysoil_frac, qsintveg, mc_layh, mcl_layh) 
@@ -557,5 +570,6 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE hydrol_main (kjit, kjpindex, &
+  SUBROUTINE hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw,  &
+       & kjit, kjpindex, &
        & index, indexveg, indexsoil, indexlayer, indexnslm, &
        & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
@@ -607,4 +621,15 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: transpir         !! Transpiration
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinf_slope      !! Slope coef
+
+    !salma: added input soil params:
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+    !end salma: added input soil params
+ 
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot           !! Soil Potential Evaporation
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot_penm      !! Soil Potential Evaporation Correction
@@ -672,8 +697,10 @@
 
     !! 0.4 Local variables
-
+    !salma: added variable kfact_root_type
     INTEGER(i_std)                                     :: jst              !! Index of soil tiles (unitless, 1-3)
     INTEGER(i_std)                                     :: jsl              !! Index of soil layers (unitless)
     INTEGER(i_std)                                     :: ji, jv
+    CHARACTER(LEN=80)                                  :: kfact_root_type  !! read from run.def: when equal to 'cons', it indicates that ks does not increase
+                                                                           !!   in the rootzone, ie, kfact_root=1; else, kfact_root defined as usual
     REAL(r_std),DIMENSION (kjpindex)                   :: soilwet          !! A temporary diagnostic of soil wetness
     REAL(r_std),DIMENSION (kjpindex)                   :: snowdepth_diag   !! Depth of snow layer containing default values, only for diagnostics
@@ -760,5 +787,5 @@
              IF(soiltile(ji,jst) .GT. min_sechiba) THEN
                 kfact_root(ji,jsl,jst) = kfact_root(ji,jsl,jst) * &
-                     & MAX((MAXVAL(ks_usda)/ks(njsc(ji)))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), &
+                     & MAX((MAXVAL(ks_usda)/ks(ji))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), &
                      un) 
              ENDIF
@@ -767,4 +794,20 @@
     ENDDO
 
+    !salma: added config key of KFACT_ROOT_TYPE
+    !! KFACT_ROOT_TYPE = cons is used to impose that kfact_root = 1 in every soil layer,
+    !! so that ks does not increase in the rootzone; else, kfact_root defined as usual
+    !Config Key   = KFACT_ROOT_TYPE
+    !Config Desc  = keyword added for spmip exp1 and exp4 to get a constant ks over soil depth and rootzone
+    !Config If    = spmip exp1 or exp4
+    !Config Def   = var
+    !Config Help  = can have two values: 'cons' or 'var'. If var then no changes, if cons then kfact_root=un
+    !Config Units = [mm/d]
+    CALL getin_p("KFACT_ROOT_TYPE",kfact_root_type)
+
+    IF (kfact_root_type=='cons') THEN
+    kfact_root(:,:,:) = un
+    ENDIF
+    !end salma: added config key
+
     !
     !! 3.3 computes canopy  ==>hydrol_canop
@@ -778,5 +821,5 @@
     !! 3.5 computes soil hydrology ==>hydrol_soil
 
-    CALL hydrol_soil(kjpindex, veget_max, soiltile, njsc, reinf_slope,  &
+    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope,  &
          transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
          returnflow, reinfiltration, irrigation, &
@@ -869,8 +912,9 @@
     !     mcs etc are identical in all layers (no normalization by vegtot to be comparable to mc)
     DO jsl=1,nslm
-       land_mcs(:,jsl) = mcs(njsc(:))
-       land_mcfc(:,jsl) = mcfc(njsc(:))
-       land_mcw(:,jsl) = mcw(njsc(:))
-       land_mcr(:,jsl) = mcr(njsc(:))
+       !salma: replaced mcs(njsc(:)) by mcs(:) and same for other variables
+       land_mcs(:,jsl) = mcs(:)
+       land_mcfc(:,jsl) = mcfc(:)
+       land_mcw(:,jsl) = mcw(:)
+       land_mcr(:,jsl) = mcr(:)
     ENDDO
     CALL xios_orchidee_send_field("mcs",land_mcs) ! in m3/m3
@@ -878,5 +922,6 @@
     CALL xios_orchidee_send_field("mcw",land_mcw) ! in m3/m3 
     CALL xios_orchidee_send_field("mcr",land_mcr) ! in m3/m3 
-          
+
+      
     CALL xios_orchidee_send_field("water2infilt",water2infilt)   
     CALL xios_orchidee_send_field("mc",mc)
@@ -915,6 +960,6 @@
 
     ! For the soil wetness above wilting point for CMIP6 (mrsow)
-    mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(njsc(:)) ) &
-         / ( zmaxh*mille*( mcs(njsc(:)) - mcw(njsc(:)) ) )
+    mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(:) ) &
+         / ( zmaxh*mille*( mcs(:) - mcw(:) ) )
     CALL xios_orchidee_send_field("mrsow",mrsow)
 
@@ -1288,6 +1333,8 @@
 !_ ================================================================================================================================
 !!_ hydrol_init
-
-  SUBROUTINE hydrol_init(kjit, kjpindex, index, rest_id, veget_max, soiltile, &
+!salma: added variables and njsc in arguments and input variables
+
+  SUBROUTINE hydrol_init(ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc,&
+       kjit, kjpindex, index, rest_id, veget_max, soiltile, &
        humrel,  vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
        snowdz,  snowgrain, snowrho,    snowtemp,   snowheat, &
@@ -1298,5 +1345,6 @@
 
     !! 0.1 Input variables
-
+    !salma introduced njsc variable
+    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: njsc               !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
     INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number 
     INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size
@@ -1305,6 +1353,15 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max          !! Carte de vegetation max
     REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in)  :: soiltile           !! Fraction of each soil tile within vegtot (0-1, unitless)
-
+   
     !! 0.2 Output variables
+
+    !salma: added input variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
 
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: humrel             !! Stress hydrique, relative humidity
@@ -1407,32 +1464,12 @@
     !! 2.1 array allocation for soil texture
 
-    ALLOCATE (nvan(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan','','')
-
-    ALLOCATE (avan(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan','','')
-
-    ALLOCATE (mcr(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcr','','')
-
-    ALLOCATE (mcs(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcs','','')
-
-    ALLOCATE (ks(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ks','','')
-
-    ALLOCATE (pcent(nscm),stat=ier)
+   
+    ALLOCATE (pcent(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable pcent','','')
-
-    ALLOCATE (mcfc(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcfc','','')
-
-    ALLOCATE (mcw(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcw','','')
-    
-    ALLOCATE (mc_awet(nscm),stat=ier)
+    
+    ALLOCATE (mc_awet(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_awet','','')
 
-    ALLOCATE (mc_adry(nscm),stat=ier)
+    ALLOCATE (mc_adry(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','')
        
@@ -1442,26 +1479,12 @@
     CASE (3)
               
-       nvan(:) = nvan_fao(:)       
-       avan(:) = avan_fao(:)
-       mcr(:) = mcr_fao(:)
-       mcs(:) = mcs_fao(:)
-       ks(:) = ks_fao(:)
-       pcent(:) = pcent_fao(:)
-       mcfc(:) = mcf_fao(:)
-       mcw(:) = mcw_fao(:)
-       mc_awet(:) = mc_awet_fao(:)
-       mc_adry(:) = mc_adry_fao(:)
-    CASE (12)
-       
-       nvan(:) = nvan_usda(:)
-       avan(:) = avan_usda(:)
-       mcr(:) = mcr_usda(:)
-       mcs(:) = mcs_usda(:)
-       ks(:) = ks_usda(:)
-       pcent(:) = pcent_usda(:)
-       mcfc(:) = mcf_usda(:)
-       mcw(:) = mcw_usda(:)
-       mc_awet(:) = mc_awet_usda(:)
-       mc_adry(:) = mc_adry_usda(:)
+       pcent(:) = pcent_fao(njsc(:))
+       mc_awet(:) = mc_awet_fao(njsc(:))
+       mc_adry(:) = mc_adry_fao(njsc(:))
+    CASE (13)
+           
+       pcent(:) = pcent_usda(njsc(:))
+       mc_awet(:) = mc_awet_usda(njsc(:))
+       mc_adry(:) = mc_adry_usda(njsc(:))
        
     CASE DEFAULT
@@ -1474,97 +1497,10 @@
     !! 2.3 Read in the run.def the parameters values defined by the user
 
-    !Config Key   = CWRR_N_VANGENUCHTEN
-    !Config Desc  = Van genuchten coefficient n
-    !Config If    = 
-    !Config Def   = 1.89, 1.56, 1.31
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [-]
-    CALL getin_p("CWRR_N_VANGENUCHTEN",nvan)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(nvan(:) <= zero) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for CWRR_N_VANGENUCHTEN.", &
-            &     "This parameter should be positive. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = CWRR_A_VANGENUCHTEN
-    !Config Desc  = Van genuchten coefficient a
-    !Config If    = 
-    !Config Def   = 0.0075, 0.0036, 0.0019
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [1/mm]  
-    CALL getin_p("CWRR_A_VANGENUCHTEN",avan)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(avan(:) <= zero) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for CWRR_A_VANGENUCHTEN.", &
-            &     "This parameter should be positive. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = VWC_RESIDUAL
-    !Config Desc  = Residual soil water content
-    !Config If    = 
-    !Config Def   = 0.065, 0.078, 0.095
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [m3/m3]  
-    CALL getin_p("VWC_RESIDUAL",mcr)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcr(:) < zero) .OR. ANY(mcr(:) > 1.)  ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_RESIDUAL.", &
-            &     "This parameter is ranged between 0 and 1. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-    
-    !Config Key   = VWC_SAT
-    !Config Desc  = Saturated soil water content
-    !Config If    = 
-    !Config Def   = 0.41, 0.43, 0.41
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [m3/m3]  
-    CALL getin_p("VWC_SAT",mcs)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcs(:) < zero) .OR. ANY(mcs(:) > 1.) .OR. ANY(mcs(:) <= mcr(:)) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_SAT.", &
-            &     "This parameter should be greater than VWC_RESIDUAL and less than 1. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = CWRR_KS 
-    !Config Desc  = Hydraulic conductivity Saturation
-    !Config If    = 
-    !Config Def   = 1060.8, 249.6, 62.4
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [mm/d]   
-    CALL getin_p("CWRR_KS",ks)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(ks(:) <= zero) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for CWRR_KS.", &
-            &     "This parameter should be positive. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
+
+
+
+
+
+
 
 
@@ -1587,39 +1523,5 @@
 
 
-    !Config Key   = VWC_FC 
-    !Config Desc  = Volumetric water content field capacity
-    !Config If    = 
-    !Config Def   = 0.32, 0.32, 0.32
-    !Config Help  = This parameter is independent from soil texture for
-    !Config         the time being.
-    !Config Units = [m3/m3]   
-    CALL getin_p("VWC_FC",mcfc)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcfc(:) > mcs(:)) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_FC.", &
-            &     "This parameter should be less than VWC_SAT. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = VWC_WP
-    !Config Desc  = Volumetric water content Wilting pt
-    !Config If    = 
-    !Config Def   = 0.10, 0.10, 0.10 
-    !Config Help  = This parameter is independent from soil texture for
-    !Config         the time being.
-    !Config Units = [m3/m3]   
-    CALL getin_p("VWC_WP",mcw)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcw(:) > mcfc(:)) .OR. ANY(mcw(:) < mcr(:)) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_WP.", &
-            &     "This parameter should be greater or equal than VWC_RESIDUAL and less or equal than VWC_SAT.", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
+   
 
     !Config Key   = VWC_MIN_FOR_WET_ALB
@@ -1745,8 +1647,8 @@
     kk(:,:,:) = 276.48
     
-    ALLOCATE (avan_mod_tab(nslm,nscm),stat=ier) 
+    ALLOCATE (avan_mod_tab(nslm,kjpindex),stat=ier) 
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan_mod_tab','','')
     
-    ALLOCATE (nvan_mod_tab(nslm,nscm),stat=ier) 
+    ALLOCATE (nvan_mod_tab(nslm,kjpindex),stat=ier) 
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan_mod_tab','','')
 
@@ -1932,5 +1834,5 @@
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact_root','','')
 
-    ALLOCATE (kfact(nslm, nscm),stat=ier)
+    ALLOCATE (kfact(nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact','','')
 
@@ -1944,17 +1846,17 @@
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dh','','')
 
-    ALLOCATE (mc_lin(imin:imax, nscm),stat=ier)
+    ALLOCATE (mc_lin(imin:imax, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_lin','','')
 
-    ALLOCATE (k_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (k_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k_lin','','')
 
-    ALLOCATE (d_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (d_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d_lin','','')
 
-    ALLOCATE (a_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (a_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a_lin','','')
 
-    ALLOCATE (b_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (b_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b_lin','','')
 
@@ -2451,12 +2353,6 @@
 
     ! Allocation for soiltile related parameters
-    IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
-    IF ( ALLOCATED (avan)) DEALLOCATE (avan)
-    IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
-    IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
-    IF ( ALLOCATED (ks)) DEALLOCATE (ks)
+   
     IF ( ALLOCATED (pcent)) DEALLOCATE (pcent)
-    IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
-    IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
     IF ( ALLOCATED (mc_awet)) DEALLOCATE (mc_awet)
     IF ( ALLOCATED (mc_adry)) DEALLOCATE (mc_adry)
@@ -2811,5 +2707,6 @@
 !_ hydrol_var_init
 
-  SUBROUTINE hydrol_var_init (kjpindex, veget, veget_max, soiltile, njsc, &
+  SUBROUTINE hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
+       kjpindex, veget, veget_max, soiltile, njsc, &
        mx_eau_var, shumdiag_perma, &
        drysoil_frac, qsintveg, mc_layh, mcl_layh) 
@@ -2820,4 +2717,12 @@
 
     !! 0.1 Input variables
+    !salma: added the following soil parameters
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
 
     ! input scalar 
@@ -2855,5 +2760,5 @@
     REAL(r_std)                                         :: nvan_mod      !! VG parameter n  modified from  exponantial profile
                                                                          !! (unitless)
-    REAL(r_std), DIMENSION(nslm,nscm)                   :: afact, nfact  !! Multiplicative factor for decay of a and n with depth
+    REAL(r_std), DIMENSION(nslm,kjpindex)               :: afact, nfact  !! Multiplicative factor for decay of a and n with depth
                                                                          !! (unitless)
     ! parameters for "soil densification" with depth
@@ -2883,5 +2788,8 @@
     REAL(r_std), DIMENSION (kjpindex,nslm)              :: nvg             !! VG param n modified with depth at each node
                                                                            !! (unitless)
-    INTEGER(i_std)                                      :: jiref           !! To identify the mc_lins where k_lin and d_lin 
+                                                                           !! need special treatment
+    !salma: added local variable ii and jiref replaced with iiref
+    INTEGER(i_std)                                      :: ii
+    INTEGER(i_std)                                      :: iiref           !! To identify the mc_lins where k_lin and d_lin
                                                                            !! need special treatment
 
@@ -3050,23 +2958,25 @@
     END IF
 
+ 
+
     !-
     !! 3 Compute the profile for a and n
     !-
-
-    ! For every soil texture
-    DO jsc = 1, nscm 
+    !salma: for every grid cell
+    DO ji = 1, kjpindex
        DO jsl=1,nslm
           ! PhD thesis of d'Orgeval, 2006, p81, Eq. 4.38; d'Orgeval et al. 2008, Eq. 2
           ! Calibrated against Hapex-Sahel measurements
-          kfact(jsl,jsc) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un)
-          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14 
-          
-          nfact(jsl,jsc) = ( kfact(jsl,jsc) )**nk_rel
-          afact(jsl,jsc) = ( kfact(jsl,jsc) )**ak_rel
-       ENDDO
-    ENDDO
-
-    ! For every soil texture
-    DO jsc = 1, nscm
+          kfact(jsl,ji) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un)
+          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14
+
+          nfact(jsl,ji) = ( kfact(jsl,ji) )**nk_rel
+          afact(jsl,ji) = ( kfact(jsl,ji) )**ak_rel
+       ENDDO
+    ENDDO
+    !end salma
+    
+    ! For every grid cell
+     DO ji = 1, kjpindex
        !-
        !! 4 Compute the linearized values of k, a, b and d
@@ -3080,75 +2990,77 @@
 
        ! We define 51 bounds for 50 bins of mc between mcr and mcs
-       mc_lin(imin,jsc)=mcr(jsc)
-       mc_lin(imax,jsc)=mcs(jsc)
-       DO ji= imin+1, imax-1 ! ji=2,50
-          mc_lin(ji,jsc) = mcr(jsc) + (ji-imin)*(mcs(jsc)-mcr(jsc))/(imax-imin)
+       mc_lin(imin,ji)=mcr(ji)
+       mc_lin(imax,ji)=mcs(ji)
+       DO ii= imin+1, imax-1 ! ii=2,50
+          mc_lin(ii,ji) = mcr(ji) + (ii-imin)*(mcs(ji)-mcr(ji))/(imax-imin)
        ENDDO
 
        DO jsl = 1, nslm
           ! From PhD thesis of d'Orgeval, 2006, p81, Eq. 4.42
-          nvan_mod = n0 + (nvan(jsc)-n0) * nfact(jsl,jsc)
-          avan_mod = a0 + (avan(jsc)-a0) * afact(jsl,jsc)
+          nvan_mod = n0 + (nvan(ji)-n0) * nfact(jsl,ji)
+          avan_mod = a0 + (avan(ji)-a0) * afact(jsl,ji)
           m = un - un / nvan_mod
           ! Creation of arrays for SP-MIP output by landpoint
-          nvan_mod_tab(jsl,jsc) = nvan_mod
-          avan_mod_tab(jsl,jsc) = avan_mod
-          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(jsc) * kfact(jsl,jsc)
-          DO ji = imax,imin,-1 
-             frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
-             k_lin(ji,jsl,jsc) = ks(jsc) * kfact(jsl,jsc) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
+          nvan_mod_tab(jsl,ji) = nvan_mod
+          avan_mod_tab(jsl,ji) = avan_mod
+          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(ji) * kfact(jsl,ji)
+          DO ii = imax,imin,-1
+             frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
+             k_lin(ii,jsl,ji) = ks(ji) * kfact(jsl,ji) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
           ENDDO
 
           ! k_lin should not be zero, nor too small
-          ! We track jiref, the bin under which mc is too small and we may get zero k_lin     
-          ji=imax-1
-          DO WHILE ((k_lin(ji,jsl,jsc) > 1.e-32) .and. (ji>0))
-             jiref=ji
-             ji=ji-1
+          ! We track iiref, the bin under which mc is too small and we may get zero k_lin
+          !salma: ji replaced with ii and jiref replaced with iiref and jsc with ji
+          ii=imax-1
+          DO WHILE ((k_lin(ii,jsl,ji) > 1.e-32) .and. (ii>0))
+             iiref=ii
+             ii=ii-1
           ENDDO
-          DO ji=jiref-1,imin,-1
-             k_lin(ji,jsl,jsc)=k_lin(ji+1,jsl,jsc)/10.
+          DO ii=iiref-1,imin,-1
+             k_lin(ii,jsl,ji)=k_lin(ii+1,jsl,ji)/10.
           ENDDO
-         
-          DO ji = imin,imax-1 ! ji=1,50
+
+          DO ii = imin,imax-1 ! ii=1,50
              ! We deduce a_lin and b_lin based on continuity between segments k_lin = a_lin*mc-lin+b_lin
-             a_lin(ji,jsl,jsc) = (k_lin(ji+1,jsl,jsc)-k_lin(ji,jsl,jsc)) / (mc_lin(ji+1,jsc)-mc_lin(ji,jsc))
-             b_lin(ji,jsl,jsc)  = k_lin(ji,jsl,jsc) - a_lin(ji,jsl,jsc)*mc_lin(ji,jsc)
+             a_lin(ii,jsl,ji) = (k_lin(ii+1,jsl,ji)-k_lin(ii,jsl,ji)) / (mc_lin(ii+1,ji)-mc_lin(ii,ji))
+             b_lin(ii,jsl,ji)  = k_lin(ii,jsl,ji) - a_lin(ii,jsl,ji)*mc_lin(ii,ji)
 
              ! We calculate the d_lin for each mc bin, from Van Genuchten equation for D(theta)
-             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin 
-             IF (ji.NE.imin .AND. ji.NE.imax-1) THEN
-                frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
-                d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) *  &
-                     ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) * &
+             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin
+             IF (ii.NE.imin .AND. ii.NE.imax-1) THEN
+                frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
+                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) *  &
+                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) * &
                      (  frac**(-un/m) -un ) ** (-m)
-                frac=MIN(un,(mc_lin(ji+1,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
-                d_lin(ji+1,jsl,jsc) =(k_lin(ji+1,jsl,jsc) / (avan_mod*m*nvan_mod))*&
-                     ( (frac**(-un/m))/(mc_lin(ji+1,jsc)-mcr(jsc)) ) * &
+                frac=MIN(un,(mc_lin(ii+1,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
+                d_lin(ii+1,jsl,ji) =(k_lin(ii+1,jsl,ji) / (avan_mod*m*nvan_mod))*&
+                     ( (frac**(-un/m))/(mc_lin(ii+1,ji)-mcr(ji)) ) * &
                      (  frac**(-un/m) -un ) ** (-m)
-                d_lin(ji,jsl,jsc) = undemi * (d_lin(ji,jsl,jsc)+d_lin(ji+1,jsl,jsc))
-             ELSE IF(ji.EQ.imax-1) THEN
-                d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * &
-                     ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) *  &
+                d_lin(ii,jsl,ji) = undemi * (d_lin(ii,jsl,ji)+d_lin(ii+1,jsl,ji))
+             ELSE IF(ii.EQ.imax-1) THEN
+                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) * &
+                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) *  &
                      (  frac**(-un/m) -un ) ** (-m)
              ENDIF
-          ENDDO
-
-          ! Special case for ji=imin
-          d_lin(imin,jsl,jsc) = d_lin(imin+1,jsl,jsc)/1000.
+          ENDDO !Salma end loop over landpoints
+
+          ! Special case for ii=imin
+          d_lin(imin,jsl,ji) = d_lin(imin+1,jsl,ji)/1000.
 
           ! We adjust d_lin where k_lin was previously adjusted otherwise we might get non-monotonous variations
           ! We don't want d_lin = zero
-          DO ji=jiref-1,imin,-1
-             d_lin(ji,jsl,jsc)=d_lin(ji+1,jsl,jsc)/10.
+          DO ii=iiref-1,imin,-1
+             d_lin(ii,jsl,ji)=d_lin(ii+1,jsl,ji)/10.
           ENDDO
 
        ENDDO
     ENDDO
+
 
     ! Output of alphavg and nvg at each node for SP-MIP
     DO jsl = 1, nslm
-       alphavg(:,jsl) = avan_mod_tab(jsl,njsc(:))*1000. ! from mm-1 to m-1
-       nvg(:,jsl) = nvan_mod_tab(jsl,njsc(:))
+       alphavg(:,jsl) = avan_mod_tab(jsl,:)*1000. ! from mm-1 to m-1
+       nvg(:,jsl) = nvan_mod_tab(jsl,:)
     ENDDO
     CALL xios_orchidee_send_field("alphavg",alphavg) ! in m-1
@@ -3165,7 +3077,7 @@
 
     mx_eau_var(:) = zero
-    mx_eau_var(:) = zmaxh*mille*mcs(njsc(:)) 
-
-    DO ji = 1,kjpindex 
+    mx_eau_var(:) = zmaxh*mille*mcs(:)
+
+    DO ji = 1,kjpindex
        IF (vegtot(ji) .LE. zero) THEN
           mx_eau_var(ji) = mx_eau_nobio*zmaxh
@@ -3181,13 +3093,13 @@
 
     ! Loop on soiltiles to compute the variables (ji,jst)
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji = 1, kjpindex
-          tmcs(ji,jst)=zmaxh* mille*mcs(njsc(ji))
-          tmcr(ji,jst)=zmaxh* mille*mcr(njsc(ji))
-          tmcfc(ji,jst)=zmaxh* mille*mcfc(njsc(ji))
-          tmcw(ji,jst)=zmaxh* mille*mcw(njsc(ji))
-       ENDDO
-    ENDDO
-       
+          tmcs(ji,jst)=zmaxh* mille*mcs(ji)
+          tmcr(ji,jst)=zmaxh* mille*mcr(ji)
+          tmcfc(ji,jst)=zmaxh* mille*mcfc(ji)
+          tmcw(ji,jst)=zmaxh* mille*mcw(ji)
+       ENDDO
+    ENDDO
+
     ! The total soil moisture for each soiltile:
     DO jst=1,nstm
@@ -3197,5 +3109,5 @@
     ENDDO
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO jsl=2,nslm-1
           DO ji=1,kjpindex
@@ -3206,5 +3118,5 @@
     ENDDO
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji=1,kjpindex
           tmc(ji,jst) = tmc(ji,jst) +  dz(nslm) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
@@ -3213,5 +3125,5 @@
     END DO
 
-!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.    
+!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.
 !    ! If veget has been updated before restart (with LAND USE or DGVM),
 !    ! tmc and mc must be modified with respect to humtot conservation.
@@ -3220,40 +3132,40 @@
     ! The litter variables:
     ! level 1
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji=1,kjpindex
           tmc_litter(ji,jst) = dz(2) * (trois*mcl(ji,1,jst)+mcl(ji,2,jst))/huit
-          tmc_litter_wilt(ji,jst) = dz(2) * mcw(njsc(ji)) / deux
-          tmc_litter_res(ji,jst) = dz(2) * mcr(njsc(ji)) / deux
-          tmc_litter_field(ji,jst) = dz(2) * mcfc(njsc(ji)) / deux
-          tmc_litter_sat(ji,jst) = dz(2) * mcs(njsc(ji)) / deux
-          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(njsc(ji)) / deux
-          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(njsc(ji)) / deux
+          tmc_litter_wilt(ji,jst) = dz(2) * mcw(ji) / deux
+          tmc_litter_res(ji,jst) = dz(2) * mcr(ji) / deux
+          tmc_litter_field(ji,jst) = dz(2) * mcfc(ji) / deux
+          tmc_litter_sat(ji,jst) = dz(2) * mcs(ji) / deux
+          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(ji) / deux
+          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(ji) / deux
        ENDDO
     END DO
     ! sum from level 2 to 4
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO jsl=2,4
           DO ji=1,kjpindex
-             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * & 
+             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * &
                   & ( trois*mcl(ji,jsl,jst) + mcl(ji,jsl-1,jst))/huit &
                   & + dz(jsl+1)*(trois*mcl(ji,jsl,jst) + mcl(ji,jsl+1,jst))/huit
              tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))*& 
-                  & mcw(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))*&
+                  & mcw(ji)/deux
              tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))*& 
-                  & mcr(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))*&
+                  & mcr(ji)/deux
              tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))* & 
-                  & mcs(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))* &
+                  & mcs(ji)/deux
              tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + &
-                  & (dz(jsl)+ dz(jsl+1))* & 
-                  & mcfc(njsc(ji))/deux
+                  & (dz(jsl)+ dz(jsl+1))* &
+                  & mcfc(ji)/deux
              tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))* & 
-                  & mc_awet(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))* &
+                  & mc_awet(ji)/deux
              tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + &
-                  & (dz(jsl)+ dz(jsl+1))* & 
-                  & mc_adry(njsc(ji))/deux
+                  & (dz(jsl)+ dz(jsl+1))* &
+                  & mc_adry(ji)/deux
           END DO
        END DO
@@ -3261,8 +3173,8 @@
 
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji=1,kjpindex
           ! here we set that humrelv=0 in PFT1
-          humrelv(ji,1,jst) = zero
+         humrelv(ji,1,jst) = zero
        ENDDO
     END DO
@@ -3270,5 +3182,5 @@
 
     ! Calculate shumdiag_perma for thermosoil
-    ! Use resdist instead of soiltile because we here need to have 
+    ! Use resdist instead of soiltile because we here need to have
     ! shumdiag_perma at the value from previous time step.
     ! Here, soilmoist is only used as a temporary variable to calculate shumdiag_perma
@@ -3290,16 +3202,16 @@
     ENDDO
     DO ji=1,kjpindex
-        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average 
-    ENDDO
-    
+        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average
+    ENDDO
+
     ! -- shumdiag_perma for restart
-    !  For consistency with hydrol_soil, we want to calculate a grid-cell average
+   !  For consistency with hydrol_soil, we want to calculate a grid-cell average
     DO jsl = 1, nslm
-       DO ji=1,kjpindex        
-          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji)))
-          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 
-       ENDDO
-    ENDDO
-               
+       DO ji=1,kjpindex
+          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji))
+          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero)
+       ENDDO
+    ENDDO
+
     ! Calculate drysoil_frac if it was not found in the restart file
     ! For simplicity, we set drysoil_frac to 0.5 in this case
@@ -3310,9 +3222,11 @@
     END IF
 
-    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in 
-    !! thermosoil for the thermal conductivity. 
+    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in
+    !! thermosoil for the thermal conductivity.
     ! These values are only used in thermosoil_init in absence of a restart file
+
     mc_layh(:,:) = zero
     mcl_layh(:,:) = zero
+      
     DO jst=1,nstm
        DO jsl=1,nslm
@@ -3668,5 +3582,4 @@
 
   END SUBROUTINE hydrol_flood
-
 
 !! ================================================================================================================================
@@ -3734,6 +3647,6 @@
 !_ ================================================================================================================================
 !_ hydrol_soil
-
-  SUBROUTINE hydrol_soil (kjpindex, veget_max, soiltile, njsc, reinf_slope, &
+  SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
+       kjpindex, veget_max, soiltile, njsc, reinf_slope, &
        & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &
        & returnflow, reinfiltration, irrigation, &
@@ -3748,4 +3661,13 @@
 
     !! 0.1 Input variables
+    !salma added soil params to arguments and input variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+
 
     INTEGER(i_std), INTENT(in)                               :: kjpindex 
@@ -3807,5 +3729,4 @@
                                                                                  !! averaged over the mesh (for thermosoil)
                                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex 
-
     !! 0.3 Modified variables
 
@@ -3953,5 +3874,5 @@
        ! These values will be kept till the end of the prognostic loop
        DO jst=1,nstm
-          CALL hydrol_soil_froz(kjpindex,jst,njsc)
+          CALL hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,jst,njsc)
        ENDDO
 
@@ -4080,5 +4001,5 @@
        !! K and D are computed for the profile of mc before infiltration
        !! They depend on the fraction of soil ice, given by profil_froz_hydro_ns
-       CALL hydrol_soil_coef(kjpindex,jst,njsc)
+       CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc)
 
        !! Infiltration and surface runoff are computed
@@ -4086,5 +4007,5 @@
        !! The conductivity comes from hydrol_soil_coef and relates to the liquid phase only
        !  This seems consistent with ok_freeze
-       CALL hydrol_soil_infilt(kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, &
+       CALL hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, &
             check_infilt_ns)
        ru_ns(:,jst) = ru_infilt_ns(:,jst) 
@@ -4112,5 +4033,5 @@
        !! 2.1 K and D are recomputed after infiltration
        !! They depend on the fraction of soil ice, still given by profil_froz_hydro_ns 
-       CALL hydrol_soil_coef(kjpindex,jst,njsc)
+       CALL hydrol_soil_coef(mcr, mcs,kjpindex,jst,njsc)
  
        !! 2.2 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme
@@ -4121,6 +4042,6 @@
        DO jsl = 1, nslm
           DO ji =1, kjpindex
-             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + &
-                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + &
+                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) )
              ! we always have mcl<=mc
              ! if mc>mcr, then mcl>mcr; if mc=mcr,mcl=mcr; if mc<mcr, then mcl<mcr
@@ -4143,8 +4064,8 @@
           DO ji = 1, kjpindex
              IF (soiltile(ji,jst) .GT. zero) THEN
-                mvg = un - un / nvan_mod_tab(jsl,njsc(ji))
-                avg = avan_mod_tab(jsl,njsc(ji))*1000. ! to convert in m-1
-                mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(njsc(ji)))/(mcs(njsc(ji)) - mcr(njsc(ji))) )**(-un/mvg)
-                psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,njsc(ji))) / avg ! in m 
+                mvg = un - un / nvan_mod_tab(jsl,ji)
+                avg = avan_mod_tab(jsl,ji)*1000. ! to convert in m-1
+                mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(ji))/(mcs(ji) - mcr(ji)) )**(-un/mvg)
+                psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,ji)) / avg ! in m 
                 psi_moy(ji,jsl) = psi_moy(ji,jsl) + soiltile(ji,jst) * psi ! average across soil tiles
              ENDIF
@@ -4274,5 +4195,5 @@
           DO ji =1, kjpindex
              mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + &
-                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl
           ENDDO
@@ -4286,5 +4207,5 @@
        dr_corr_ns(:,jst) = zero
        ru_corr_ns(:,jst) = zero
-       call hydrol_soil_smooth_over_mcs2(kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns)
+       call hydrol_soil_smooth_over_mcs2(mcs, kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns)
        
        ! In absence of freezing, if F is large enough, the correction of oversaturation is sent to drainage       
@@ -4330,6 +4251,6 @@
                 dmc(ji,jsl) = zero
                 IF ( ( zz(jsl) >= zwt_force(ji,jst)*mille ) ) THEN
-                   dmc(ji,jsl) = mcs(njsc(ji)) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value
-                   mc(ji,jsl,jst) = mcs(njsc(ji))
+                   dmc(ji,jsl) = mcs(ji) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value
+                   mc(ji,jsl,jst) = mcs(ji)
                 ENDIF
              ENDDO
@@ -4366,5 +4287,5 @@
           wtd_ns(ji,jst) = undef_sechiba ! in meters
           jsl=nslm
-          DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(njsc(ji))) .AND. (jsl > 1) )
+          DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(ji)) .AND. (jsl > 1) )
              wtd_ns(ji,jst) = zz(jsl)/mille ! in meters
              jsl=jsl-1
@@ -4375,5 +4296,5 @@
        !      This routine does not change tmc but decides where we should turn off ET to prevent further mc decrease
        !      Like above, the tests are made on total mc, compared to mcr
-       CALL hydrol_soil_smooth_under_mcr(kjpindex, jst, njsc, is_under_mcr, check_under_ns)
+       CALL hydrol_soil_smooth_under_mcr(mcr, kjpindex, jst, njsc, is_under_mcr, check_under_ns)
  
        !! 4. At the end of the prognostic calculations, we recompute important moisture variables
@@ -4399,6 +4320,6 @@
        DO jsl = 1, nslm
           DO ji =1, kjpindex
-             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + &
-                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + &
+                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc
           ENDDO
@@ -4451,7 +4372,7 @@
        sm(:,1)  = dz(2) * (trois*mcl(:,1,jst) + mcl(:,2,jst))/huit
        smt(:,1) = dz(2) * (trois*mc(:,1,jst) + mc(:,2,jst))/huit
-       smw(:,1) = dz(2) * (quatre*mcw(njsc(:)))/huit
-       smf(:,1) = dz(2) * (quatre*mcfc(njsc(:)))/huit
-       sms(:,1) = dz(2) * (quatre*mcs(njsc(:)))/huit
+       smw(:,1) = dz(2) * (quatre*mcw(:))/huit
+       smf(:,1) = dz(2) * (quatre*mcfc(:))/huit
+       sms(:,1) = dz(2) * (quatre*mcs(:))/huit
        DO jsl = 2,nslm-1
           sm(:,jsl)  = dz(jsl) * (trois*mcl(:,jsl,jst)+mcl(:,jsl-1,jst))/huit &
@@ -4459,19 +4380,19 @@
           smt(:,jsl) = dz(jsl) * (trois*mc(:,jsl,jst)+mc(:,jsl-1,jst))/huit &
                + dz(jsl+1) * (trois*mc(:,jsl,jst)+mc(:,jsl+1,jst))/huit
-          smw(:,jsl) = dz(jsl) * ( quatre*mcw(njsc(:)) )/huit &
-               + dz(jsl+1) * ( quatre*mcw(njsc(:)) )/huit
-          smf(:,jsl) = dz(jsl) * ( quatre*mcfc(njsc(:)) )/huit &
-               + dz(jsl+1) * ( quatre*mcfc(njsc(:)) )/huit
-          sms(:,jsl) = dz(jsl) * ( quatre*mcs(njsc(:)) )/huit &
-               + dz(jsl+1) * ( quatre*mcs(njsc(:)) )/huit
+          smw(:,jsl) = dz(jsl) * ( quatre*mcw(:) )/huit &
+               + dz(jsl+1) * ( quatre*mcw(:) )/huit
+          smf(:,jsl) = dz(jsl) * ( quatre*mcfc(:) )/huit &
+               + dz(jsl+1) * ( quatre*mcfc(:) )/huit
+          sms(:,jsl) = dz(jsl) * ( quatre*mcs(:) )/huit &
+               + dz(jsl+1) * ( quatre*mcs(:) )/huit
        ENDDO
        sm(:,nslm)  = dz(nslm) * (trois*mcl(:,nslm,jst) + mcl(:,nslm-1,jst))/huit  
        smt(:,nslm) = dz(nslm) * (trois*mc(:,nslm,jst) + mc(:,nslm-1,jst))/huit      
-       smw(:,nslm) = dz(nslm) * (quatre*mcw(njsc(:)))/huit
-       smf(:,nslm) = dz(nslm) * (quatre*mcfc(njsc(:)))/huit
-       sms(:,nslm) = dz(nslm) * (quatre*mcs(njsc(:)))/huit
+       smw(:,nslm) = dz(nslm) * (quatre*mcw(:))/huit
+       smf(:,nslm) = dz(nslm) * (quatre*mcfc(:))/huit
+       sms(:,nslm) = dz(nslm) * (quatre*mcs(:))/huit
        ! sm_nostress = soil moisture of each layer at which us reaches 1, here at the middle of [smw,smf]
        DO jsl = 1,nslm
-          sm_nostress(:,jsl) = smw(:,jsl) + pcent(njsc(:)) * (smf(:,jsl)-smw(:,jsl))
+          sm_nostress(:,jsl) = smw(:,jsl) + pcent(:) * (smf(:,jsl)-smw(:,jsl))
        END DO
 
@@ -4642,4 +4563,5 @@
        !! The multiplication by vegtot creates grid-cell average values
        ! *** To be checked for consistency with the use of nobio properties in thermosoil
+           
        DO jsl=1,nslm
           DO ji=1,kjpindex
@@ -4661,5 +4583,5 @@
        DO jsl = 1, nslm
           ksat(:,jsl) = ksat(:,jsl) + soiltile(:,jst) * &
-               ( ks(njsc(:)) * kfact(jsl,njsc(:)) * kfact_root(:,jsl,jst) ) 
+               ( ks(:) * kfact(jsl,:) * kfact_root(:,jsl,jst) ) 
        ENDDO
               
@@ -4671,5 +4593,5 @@
 
     !! 7. Summing 3d variables into 2d variables 
-    CALL hydrol_diag_soil (kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
+    CALL hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
          & evapot, vevapnu, returnflow, reinfiltration, irrigation, &
          & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,tot_melt)
@@ -4762,9 +4684,9 @@
        !! 8.2 Since we estimate bare soile evap for the next time step, we update profil_froz_hydro and mcl
        !     (effect of mc only, the change in temp_hydro is neglected)
-       IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(kjpindex,jst,njsc)
+       IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(nvan, avan, mcr, mcs,kjpindex,jst,njsc)
         DO jsl = 1, nslm
           DO ji =1, kjpindex
-             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + &
-                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + &
+                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc
           ENDDO
@@ -4772,5 +4694,5 @@
 
        !! 8.3 K and D are recomputed for the updated profile of mc/mcl
-       CALL hydrol_soil_coef(kjpindex,jst,njsc)
+       CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc)
 
        !! 8.4 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme
@@ -4811,4 +4733,5 @@
              
              flux_top(ji) = zero
+             r_soil_ns(ji,jst) = zero
              
           ENDIF
@@ -4863,5 +4786,5 @@
        DO ji = 1, kjpindex
           ! by construction, mc and mcl are always on the same side of mcr, so we can use mcl here
-          resolv(ji) = (mcl(ji,1,jst).LT.(mcr(njsc(ji))).AND.flux_top(ji).GT.min_sechiba)
+          resolv(ji) = (mcl(ji,1,jst).LT.(mcr(ji)).AND.flux_top(ji).GT.min_sechiba)
        ENDDO
        !! Reset the coefficient for diffusion (tridiag is only solved if resolv(ji) = .TRUE.)O
@@ -4879,5 +4802,5 @@
           tmat(ji,1,2) = un
           tmat(ji,1,3) = zero
-          rhs(ji,1) = mcr(njsc(ji))
+          rhs(ji,1) = mcr(ji)
        ENDDO
        
@@ -4898,5 +4821,5 @@
           DO ji =1, kjpindex
              mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + &
-                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl
           ENDDO
@@ -5050,5 +4973,5 @@
 !_ hydrol_soil_infilt
 
-  SUBROUTINE hydrol_soil_infilt(kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check)
+  SUBROUTINE hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check)
 
     !! 0. Variable and parameter declaration
@@ -5057,4 +4980,14 @@
 
     ! GLOBAL (in or inout)
+    ! salma added input soil hydraulic parameters
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+
+
     INTEGER(i_std), INTENT(in)                        :: kjpindex        !! Domain size
     INTEGER(i_std), INTENT(in)                        :: ins
@@ -5108,5 +5041,5 @@
     DO ji = 1, kjpindex
        !! First we fill up the first layer (about 1mm) without any resistance and quasi-immediately
-       wat_inf_pot(ji) = MAX((mcs(njsc(ji))-mc(ji,1,ins)) * dz(2) / deux, zero)
+       wat_inf_pot(ji) = MAX((mcs(ji)-mc(ji,1,ins)) * dz(2) / deux, zero)
        wat_inf(ji) = MIN(wat_inf_pot(ji), flux_infilt(ji))
        mc(ji,1,ins) = mc(ji,1,ins) + wat_inf(ji) * deux / dz(2)
@@ -5127,5 +5060,5 @@
           ! This is computed by an simple arithmetic average because 
           ! the time step (30min) is not appropriate for a geometric average (advised by Haverkamp and Vauclin)
-          k_m = (k(ji,jsl) + ks(njsc(ji))*kfact(jsl-1,njsc(ji))*kfact_root(ji,jsl,ins)) / deux 
+          k_m = (k(ji,jsl) + ks(ji)*kfact(jsl-1,ji)*kfact_root(ji,jsl,ins)) / deux 
 
           IF (ok_freeze_cwrr) THEN
@@ -5141,5 +5074,5 @@
 
           !! From which we deduce the time it takes to fill up the layer or to end the time step...
-          wat_inf_pot(ji) =  MAX((mcs(njsc(ji))-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero)
+          wat_inf_pot(ji) =  MAX((mcs(ji)-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero)
           IF ( infilt_tmp(ji) > min_sechiba) THEN
              dt_inf(ji) =  MIN(wat_inf_pot(ji)/infilt_tmp(ji), dt_tmp(ji))
@@ -5223,5 +5156,5 @@
 !_ hydrol_soil_smooth_under_mcr
 
-  SUBROUTINE hydrol_soil_smooth_under_mcr(kjpindex, ins, njsc, is_under_mcr, check)
+  SUBROUTINE hydrol_soil_smooth_under_mcr(mcr, kjpindex, ins, njsc, is_under_mcr, check)
 
     !- arguments
@@ -5231,4 +5164,6 @@
     !! 0.1 Input variables
 
+    !salma: added parameter mcr
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr          !! Residual volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                         :: kjpindex     !! Domain size
     INTEGER(i_std), INTENT(in)                         :: ins          !! Soiltile index (1-nstm, unitless)
@@ -5268,5 +5203,5 @@
     DO jsl = 1,nslm-2
        DO ji=1, kjpindex
-          excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+          excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
           mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * &
@@ -5277,5 +5212,5 @@
     jsl = nslm-1
     DO ji=1, kjpindex
-       excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+       excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
        mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * &
@@ -5285,5 +5220,5 @@
     jsl = nslm
     DO ji=1, kjpindex
-       excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+       excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
        mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * &
@@ -5294,5 +5229,5 @@
     DO jsl = nslm-1,2,-1
        DO ji=1, kjpindex
-          excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+          excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
           mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * &
@@ -5304,5 +5239,5 @@
     ! excess > 0
     DO ji=1, kjpindex
-       excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(njsc(ji))-mc(ji,1,ins),zero)
+       excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(ji)-mc(ji,1,ins),zero)
     ENDDO
     DO ji=1, kjpindex
@@ -5380,5 +5315,5 @@
 !_ hydrol_soil_smooth_over_mcs
 
-  SUBROUTINE hydrol_soil_smooth_over_mcs(kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
+  SUBROUTINE hydrol_soil_smooth_over_mcs(mcs ,kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
 
     !- arguments
@@ -5386,6 +5321,7 @@
     !! 0. Variable and parameter declaration
 
-    !! 0.1 Input variables
-
+    !! 0.1 Input variables 
+    !salma: added mcs
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: mcs             !! Saturated volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                           :: kjpindex        !! Domain size
     INTEGER(i_std), INTENT(in)                           :: ins             !! Soiltile index (1-nstm, unitless)
@@ -5425,5 +5361,5 @@
     DO jsl = 1, nslm-2
        DO ji=1, kjpindex
-          excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+          excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
           mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * &
@@ -5434,5 +5370,5 @@
     jsl = nslm-1
     DO ji=1, kjpindex
-       excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+       excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
        mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * &
@@ -5442,5 +5378,5 @@
     jsl = nslm
     DO ji=1, kjpindex
-       excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+       excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
        mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * &
@@ -5451,5 +5387,5 @@
     DO jsl = nslm-1,2,-1
        DO ji=1, kjpindex
-          excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+          excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
           mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * &
@@ -5461,5 +5397,5 @@
 
     DO ji=1, kjpindex
-       excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(njsc(ji)),zero)
+       excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(ji),zero)
        mc(ji,1,ins) = mc(ji,1,ins) - excess ! then mc(1)=mcs
        rudr_corr(ji,ins) = rudr_corr(ji,ins) + excess * dz(2) / deux 
@@ -5512,5 +5448,5 @@
 !_ hydrol_soil_smooth_over_mcs2
 
-  SUBROUTINE hydrol_soil_smooth_over_mcs2(kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
+  SUBROUTINE hydrol_soil_smooth_over_mcs2(mcs, kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
 
     !- arguments
@@ -5519,5 +5455,6 @@
 
     !! 0.1 Input variables
-
+    !salma: added mcs
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: mcs             !! Saturated volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                           :: kjpindex        !! Domain size
     INTEGER(i_std), INTENT(in)                           :: ins             !! Soiltile index (1-nstm, unitless)
@@ -5562,5 +5499,5 @@
     DO jsl = 1, nslm
        DO ji=1, kjpindex
-          excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) ! >=0
+          excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(ji),zero) ! >=0
           mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess(ji,jsl) ! here mc either does not change or decreases
        ENDDO
@@ -5796,6 +5733,6 @@
 !_ ================================================================================================================================
 !_ hydrol_soil_coef
- 
-  SUBROUTINE hydrol_soil_coef(kjpindex,ins,njsc)
+
+  SUBROUTINE hydrol_soil_coef(mcr, mcs, kjpindex,ins,njsc)
 
     IMPLICIT NONE
@@ -5804,5 +5741,7 @@
 
     !! 0.1 Input variables
-
+    !salma: added mcr and mcs
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                        :: kjpindex         !! Domain size
     INTEGER(i_std), INTENT(in)                        :: ins              !! Index of soil type
@@ -5836,14 +5775,14 @@
              ! It corresponds to a liquid mc, but the expression is different from mcl in hydrol_soil,
              ! to ensure that we get the a, b, d of the first bin when mcl<mcr
-             mc_used = mcr(njsc(ji))+x*MAX((mc(ji,jsl, ins)-mcr(njsc(ji))),zero) 
+             mc_used = mcr(ji)+x*MAX((mc(ji,jsl, ins)-mcr(ji)),zero) 
              !
              ! calcul de k based on mc_liq
              !
-             i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(njsc(ji)))/(mcs(njsc(ji))-mcr(njsc(ji))))))
-             a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d
-             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), &
-                  a_lin(i,jsl,njsc(ji)) * mc_used + b_lin(i,jsl,njsc(ji))) ! in mm/d
+             i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(ji))/(mcs(ji)-mcr(ji)))))
+             a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d
+             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), &
+                  a_lin(i,jsl,ji) * mc_used + b_lin(i,jsl,ji)) ! in mm/d
           ENDDO ! loop on grid
        ENDDO
@@ -5855,12 +5794,12 @@
              
              ! it is impossible to consider a mc<mcr for the binning
-             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(njsc(ji)), zero)/(mcs(njsc(ji))-mcr(njsc(ji)))
+             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(ji), zero)/(mcs(ji)-mcr(ji))
              
              i= MAX(MIN(INT((imax-imin)*mc_ratio)+imin , imax-1), imin)
-             a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d
-             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), &
-                  a_lin(i,jsl,njsc(ji)) * mc(ji,jsl,ins) + b_lin(i,jsl,njsc(ji)))  ! in mm/d
+             a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d
+             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), &
+                  a_lin(i,jsl,ji) * mc(ji,jsl,ins) + b_lin(i,jsl,ji))  ! in mm/d
           END DO 
        END DO
@@ -5887,5 +5826,5 @@
 !_ hydrol_soil_froz
  
-  SUBROUTINE hydrol_soil_froz(kjpindex,ins,njsc)
+  SUBROUTINE hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,ins,njsc)
 
     IMPLICIT NONE
@@ -5894,4 +5833,9 @@
 
     !! 0.1 Input variables
+    !salma: added the following soil variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
 
     INTEGER(i_std), INTENT(in)                        :: kjpindex         !! Domain size
@@ -5925,7 +5869,7 @@
           DO ji=1,kjpindex
              ! Van Genuchten parameter for thermodynamical calculation
-             m = 1. -1./nvan(njsc(ji))
+             m = 1. -1./nvan(ji)
            
-             IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(njsc(ji))+min_sechiba))) THEN
+             IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(ji)+min_sechiba))) THEN
                 ! Linear soil freezing or soil moisture below residual
                 IF (temp_hydro(ji, jsl).GE.(ZeroCelsius+fr_dT/2.)) THEN
@@ -5944,8 +5888,8 @@
                      (temp_hydro(ji,jsl) .LT. (ZeroCelsius+fr_dT/2.)) ) THEN
                    ! Factor 2.2 from the PhD of Isabelle Gouttevin
-                   x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) &
-                        *((2.2*1000.*avan(njsc(ji))*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) &
-                        *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m)) / &
-                        (mc(ji,jsl, ins)-mcr(njsc(ji))),1._r_std)                
+                   x=MIN(((mcs(ji)-mcr(ji)) &
+                        *((2.2*1000.*avan(ji)*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) &
+                        *lhf/ZeroCelsius/10.)**nvan(ji)+1.)**(-m)) / &
+                        (mc(ji,jsl, ins)-mcr(ji)),1._r_std)                
                 ELSE
                    x=0._r_std 
@@ -5955,5 +5899,5 @@
              profil_froz_hydro_ns(ji,jsl,ins) = 1._r_std-x
              
-             mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(njsc(ji))
+             mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(ji)
 
           ENDDO ! loop on grid
@@ -6397,5 +6341,5 @@
 !_ hydrol_diag_soil
 
-  SUBROUTINE hydrol_diag_soil (kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
+  SUBROUTINE hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
        & evapot, vevapnu, returnflow, reinfiltration, irrigation, &
        & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac, tot_melt)
@@ -6406,4 +6350,12 @@
 
     !! 0.1 Input variables
+    !salma: added the following soil variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
 
     ! input scalar 
@@ -6558,5 +6510,5 @@
              i= MAX(MIN(INT((imax-imin)*tmc_litter_ratio)+imin, imax-1), imin)
           ENDIF       
-          k_tmp = MAX(k_lin(i,1,njsc(ji))*ks(njsc(ji)), zero)
+          k_tmp = MAX(k_lin(i,1,ji)*ks(ji), zero)
           k_litt(ji) = k_litt(ji) + vegtot(ji)*soiltile(ji,jst) * SQRT(k_tmp) ! grid-cell average
        ENDDO      
@@ -6630,5 +6582,5 @@
           DO ji=1,kjpindex
              shumdiag(ji,jsl) = shumdiag(ji,jsl) + soil_wet_ns(ji,jsl,jst) * soiltile(ji,jst) * &
-                               ((mcs(njsc(ji))-mcw(njsc(ji)))/(mcfc(njsc(ji))-mcw(njsc(ji))))
+                               ((mcs(ji)-mcw(ji))/(mcfc(ji)-mcw(ji)))
              shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero) 
           ENDDO
@@ -6641,5 +6593,5 @@
     DO jsl=1,nslm              
        DO ji=1,kjpindex
-          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji)))
+          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji))
           shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 
        ENDDO
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 6594)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 6954)
@@ -18,5 +18,7 @@
 !! which is called before the modules which calculate them.
 !! 
-!! RECENT CHANGE(S): None 
+!! RECENT CHANGE(S): November 2020: It is possible to define soil hydraulic parameters from maps,
+!!                    as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
+!!                    Here, it leads to declare and allocate global variables. 
 !! 
 !! REFERENCE(S) : None
@@ -113,4 +115,24 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: reinf_slope    !! slope coefficient (reinfiltration) 
 !$OMP THREADPRIVATE(reinf_slope)
+
+
+!salma: adding soil hydraulic params
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: ks              !! Saturated soil conductivity (mm d^{-1})
+!$OMP THREADPRIVATE(ks)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: nvan            !! Van Genushten n parameter (unitless)
+!$OMP THREADPRIVATE(nvan)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: avan            !! Van Genushten alpha parameter (mm ^{-1})
+!$OMP THREADPRIVATE(avan)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcr             !! Residual soil moisture (m^{3} m^{-3})
+!$OMP THREADPRIVATE(mcr)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcs             !! Saturated soil moisture (m^{3} m^{-3})
+!$OMP THREADPRIVATE(mcs)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcfc            !! Volumetric water content at field capacity (m^{3} m^{-3})
+!$OMP THREADPRIVATE(mcfc)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)     :: mcw             !! Volumetric water content at wilting point (m^{3} m^{-3})
+!$OMP THREADPRIVATE(mcw)
+!end salma:adding soil hydraulic params
+
+
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag       !! Mean relative soil moisture in the different levels used 
                                                                      !! by thermosoil.f90 (unitless, 0-1)
@@ -431,5 +453,7 @@
                               index,         indexveg,     lalo,           neighbours,        &
                               resolution,    contfrac,     temp_air,                          &
-                              soiltile,      reinf_slope,  deadleaf_cover, assim_param,       &
+                              soiltile,      reinf_slope,  ks,             nvan,              &
+                              avan,          mcr,          mcs,            mcfc,              &
+                              mcw,           deadleaf_cover,               assim_param,       &
                               lai,           frac_age,     height,         veget,             &
                               frac_nobio,    njsc,         veget_max,      fraclut,           &
@@ -451,5 +475,6 @@
     
     !! 1.7 Initialize remaining hydrological variables
-    CALL hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          &
+    CALL hydrol_initialize (ks,  nvan, avan, mcr, mcs, mcfc, mcw,    &
+         kjit,           kjpindex,  index,         rest_id,          &
          njsc,           soiltile,  veget,         veget_max,        &
          humrel,         vegstress, drysoil_frac,                    &
@@ -718,5 +743,5 @@
     !! 4. Compute hydrology
     !! 4.1 Water balance from CWRR module (11 soil layers)
-    CALL hydrol_main (kjit, kjpindex, &
+    CALL hydrol_main (ks,  nvan, avan, mcr, mcs, mcfc, mcw, kjit, kjpindex, &
          & index, indexveg, indexsoil, indexlayer, indexnslm, &
          & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
@@ -1373,4 +1398,6 @@
                             njsc,       lai,       height,   veget,      &
                             frac_nobio, veget_max, reinf_slope,          &
+                            ks,         nvan,      avan,     mcr,        &
+                            mcs,        mcfc,      mcw,                  &
                             co2_to_bm,  assim_param, frac_age)
     
@@ -1530,4 +1557,30 @@
     ALLOCATE (reinf_slope(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for reinf_slope','','')
+
+    !salma: adding soil params
+    ALLOCATE (ks(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for ks','','')
+
+    ALLOCATE (nvan(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for nvan ','','')
+
+    ALLOCATE (avan(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for avan','','')
+
+    ALLOCATE (mcr(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcr','','')
+
+    ALLOCATE (mcs(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcs','','')
+
+    ALLOCATE (mcfc(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcfc','','')
+    
+    ALLOCATE (mcw(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for mcw','','')
+    !end salma: adding soil params
+
+
+
 
     ALLOCATE (vbeta1(kjpindex),stat=ier)
@@ -1877,4 +1930,15 @@
     IF ( ALLOCATED (njsc)) DEALLOCATE (njsc)
     IF ( ALLOCATED (reinf_slope)) DEALLOCATE (reinf_slope)
+
+    !salma: adding soil hydraulic params
+    IF ( ALLOCATED (ks)) DEALLOCATE (ks)
+    IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
+    IF ( ALLOCATED (avan)) DEALLOCATE (avan)
+    IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
+    IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
+    IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
+    IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
+    !end salma: adding soil hydraulic params
+
     IF ( ALLOCATED (vbeta1)) DEALLOCATE (vbeta1)
     IF ( ALLOCATED (vbeta4)) DEALLOCATE (vbeta4)
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90	(revision 6594)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90	(revision 6954)
@@ -17,7 +17,11 @@
 !!
 !! RECENT CHANGE(S): Allowed reading of USDA map, Nov 2014, ADucharne
+!!                   November 2020: It is possible to define soil hydraulic parameters from maps,
+!!                   as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
+!!                   Changes in slowproc_xios_initialize, slowproc_soilt and slowproc_finalize
 !!
 !! REFERENCE(S)	:
-!!
+!!- Tafasca S. (2020). Evaluation de lâimpact des propriÃ©tÃ©s du sol sur lâhydrologie simulee dans le
+!! modÃšle ORCHIDEE, PhD thesis, Sorbonne Universite. \n
 !! SVN          :
 !! $HeadURL$
@@ -92,4 +96,7 @@
 !!
 !! DESCRIPTION:	  Initialize xios dependant defintion before closing context defintion
+!!  
+!! RECENT CHANGE(S): Initialization of XIOS to read soil hydraulic parameters from maps,
+!!                   as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
 !!
 !! \n
@@ -219,4 +226,30 @@
     ENDIF
 
+    !Salma
+    !! 6. Prepare for reading of soil parameter files
+
+    ! Get the file name from run.def file and set file attributes accordingly
+    filename = 'params_sp_mip.nc'
+    CALL getin_p('PARAM_FILE',filename)
+    name = filename(1:LEN_TRIM(FILENAME)-3)
+    CALL xios_orchidee_set_file_attr("soilparam_file",name=name)
+    ! Determine if the file will be read by XIOS. If not, deactivate reading of the file.
+    IF (xios_interpolation .AND. restname_in=='NONE' .AND. .NOT. impsoilt) THEN
+       ! Reading will be done with XIOS later
+       IF (printlev>=2) WRITE(numout,*) 'Reading of soil hydraulic parameters file will be done later using XIOS. The filename is ', filename
+    ELSE
+       ! No reading by XIOS, deactivate soilparam_file and related variables declared in context_input_orchidee.xml.
+       ! If this is not done, the model will crash if the file is not available in the run directory.
+       IF (printlev>=2) WRITE(numout,*) 'Reading of soil parameter file will not be done with XIOS.'
+       CALL xios_orchidee_set_file_attr("soilparam_file",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilks",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilnvan",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilavan",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilmcr",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilmcs",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilmcfc",enabled=.FALSE.)
+       CALL xios_orchidee_set_field_attr("soilmcw",enabled=.FALSE.)
+    ENDIF
+
     IF (printlev_loc>=3) WRITE(numout,*) 'End slowproc_xios_intialize'
    
@@ -244,5 +277,7 @@
                                   IndexLand,     indexveg,     lalo,           neighbours,        &
                                   resolution,    contfrac,     temp_air,                          &
-                                  soiltile,      reinf_slope,  deadleaf_cover, assim_param,       &
+                                  soiltile,      reinf_slope,  ks,             nvan,              &
+                                  avan,          mcr,          mcs,            mcfc,              &
+                                  mcw,           deadleaf_cover,               assim_param,       &
                                   lai,           frac_age,     height,         veget,             &
                                   frac_nobio,    njsc,         veget_max,      fraclut,           &
@@ -282,4 +317,13 @@
     REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out)    :: nwdFraclut     !! Fraction of non-woody vegetation in each landuse tile (0-1, unitless)
     REAL(r_std),DIMENSION (kjpindex), INTENT(out)          :: reinf_slope    !! slope coef for reinfiltration
+    !Salma: adding soil params
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: nvan           !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: avan           !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3})
+
     REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (out):: assim_param    !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: deadleaf_cover !! Fraction of soil covered by dead leaves (unitless)
@@ -296,4 +340,5 @@
     CALL slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
          rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, fraclut, nwdfraclut, reinf_slope, &
+         ks,  nvan, avan, mcr, mcs, mcfc, mcw, &
          veget_max, tot_bare_soil, njsc, &
          height, lcanop, veget_update, veget_year)
@@ -684,5 +729,8 @@
 !>\BRIEF         Write to restart file variables for slowproc module and call finalization of stomate module
 !!
-!! DESCRIPTION : 
+!! DESCRIPTION :
+!!
+!! RECENT CHANGE(S): Add arrays of soil hydraulic parameters to the restart file.
+!!                   Linked to SP-MIP project (Tafasca Salma and Ducharne Agnes).
 !!
 !! MAIN OUTPUT VARIABLE(S) : 
@@ -697,4 +745,5 @@
                                 njsc,       lai,       height,   veget,      &
                                 frac_nobio, veget_max, reinf_slope,          &
+                                ks,  nvan, avan, mcr, mcs, mcfc, mcw,        &
                                 co2_to_bm,  assim_param, frac_age )
 
@@ -711,4 +760,13 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: veget_max      !! Maximum fraction of vegetation type including none biological fraction (unitless)
     REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: reinf_slope    !! slope coef for reinfiltration
+    !salma: added following soil params
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: ks             !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: nvan           !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: avan           !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT(in)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3})
+
     REAL(r_std),DIMENSION (kjpindex,nvm),INTENT(in)      :: co2_to_bm      !! virtual gpp flux between atmosphere and biosphere
     REAL(r_std),DIMENSION (kjpindex,nvm,npco2),INTENT (in):: assim_param   !! min+max+opt temperatures & vmax for photosynthesis (K, \mumol m^{-2} s^{-1})
@@ -759,5 +817,13 @@
     CALL restput_p (rest_id, 'clay_frac', nbp_glo, 1, 1, kjit, clayfraction, 'scatter',  nbp_glo, index_g)
     CALL restput_p (rest_id, 'sand_frac', nbp_glo, 1, 1, kjit, sandfraction, 'scatter',  nbp_glo, index_g)
-    !
+    !salma: added the following lines for restput of the soil parameters
+    CALL restput_p (rest_id, 'ks', nbp_glo, 1, 1, kjit, ks, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'mcs', nbp_glo, 1, 1, kjit, mcs, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'mcr', nbp_glo, 1, 1, kjit, mcr, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'mcw', nbp_glo, 1, 1, kjit, mcw, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'mcfc', nbp_glo, 1, 1, kjit, mcfc, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'nvan', nbp_glo, 1, 1, kjit, nvan, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'avan', nbp_glo, 1, 1, kjit, avan, 'scatter',  nbp_glo, index_g)
+
     ! The height of the vegetation could in principle be recalculated at the beginning of the run.
     ! However, this is very tedious, as many special cases have to be taken into account. This variable
@@ -803,4 +869,5 @@
   SUBROUTINE slowproc_init (kjit, kjpindex, IndexLand, lalo, neighbours, resolution, contfrac, &
        rest_id, lai, frac_age, veget, frac_nobio, totfrac_nobio, soiltile, fraclut, nwdfraclut, reinf_slope, &
+       ks,  nvan, avan, mcr, mcs, mcfc, mcw, &
        veget_max, tot_bare_soil, njsc, &
        height, lcanop, veget_update, veget_year)
@@ -837,4 +904,14 @@
     REAL(r_std), DIMENSION (kjpindex,nlut), INTENT(out)   :: nwdfraclut     !! Fraction of non woody vegetation in each landuse tile
     REAL(r_std), DIMENSION (kjpindex), INTENT(out)        :: reinf_slope    !! slope coef for reinfiltration
+
+    !salma: added soil parameters    
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: nvan           !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: avan           !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3})
+
     INTEGER(i_std), DIMENSION(kjpindex), INTENT(out)      :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
     
@@ -1015,5 +1092,44 @@
        njsc = NINT(tmp_real)
     END IF
-    
+
+    ! Salma: we define restget for the 7 soil variables following clay_frac
+
+    var_name= 'ks'
+    CALL ioconf_setatt_p('UNITS', 'mm/d')
+    CALL ioconf_setatt_p('LONG_NAME','Soil saturated water content')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., ks, "gather", nbp_glo, index_g)
+
+    var_name= 'mcs'
+    CALL ioconf_setatt_p('UNITS', 'm3/m3')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcs, "gather", nbp_glo, index_g)
+
+    var_name= 'mcr'
+    CALL ioconf_setatt_p('UNITS', 'm3/m3')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcr, "gather", nbp_glo, index_g)
+
+    var_name= 'mcfc'
+    CALL ioconf_setatt_p('UNITS', 'm3/m3')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcfc, "gather", nbp_glo, index_g)
+
+    var_name= 'mcw'
+    CALL ioconf_setatt_p('UNITS', 'm3/m3')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., mcw, "gather", nbp_glo, index_g)
+
+    var_name= 'nvan'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., nvan, "gather", nbp_glo, index_g)
+
+    var_name= 'avan'
+    CALL ioconf_setatt_p('UNITS', 'm-1')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., avan, "gather", nbp_glo, index_g)
+
+    !end salma
+
     var_name= 'clay_frac'
     CALL ioconf_setatt_p('UNITS', '-')
@@ -1206,5 +1322,5 @@
                MINVAL(njsc) .EQ. MAXVAL(njsc) .AND. MAXVAL(njsc) .EQ. undef_int ) THEN
              
-             CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, &
+             CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, lalo, neighbours, resolution, contfrac, soilclass, &
                   clayfraction, sandfraction, siltfraction)
              njsc(:) = 0
@@ -1343,5 +1459,5 @@
        
        IF (printlev_loc>=4) WRITE (numout,*) 'clayfraction or njcs were not in restart file, call slowproc_soilt'
-       CALL slowproc_soilt(kjpindex, lalo, neighbours, resolution, contfrac, soilclass, &
+       CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, lalo, neighbours, resolution, contfrac, soilclass, &
             clayfraction, sandfraction, siltfraction)
        IF (printlev_loc>=4) WRITE (numout,*) 'After slowproc_soilt'
@@ -2614,4 +2730,6 @@
 !!
 !! RECENT CHANGE(S): Nov 2014, ADucharne
+!!                   Nov 2020, Salma Tafasca and Agnes Ducharne: adding a choice for spmipexp/SPMIPEXP,
+!!                             and everything needed to read all maps and assign parameter values.  
 !!
 !! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction
@@ -2624,5 +2742,5 @@
 !! \n
 !_ ================================================================================================================================
-  SUBROUTINE slowproc_soilt(nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction)
+  SUBROUTINE slowproc_soilt(njsc,  ks,  nvan, avan, mcr, mcs, mcfc, mcw, nbpt, lalo, neighbours, resolution, contfrac, soilclass, clayfraction, sandfraction, siltfraction)
 
     USE interpweight
@@ -2649,4 +2767,14 @@
     !  0.2 OUTPUT
     !
+    !salma: added soil params and njsc because needed in the calculation of the soil params
+    INTEGER(i_std),DIMENSION (nbpt), INTENT (out)      :: njsc           !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: ks             !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: nvan           !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: avan           !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcr            !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcs            !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcfc           !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (nbpt), INTENT (out)         :: mcw            !! Volumetric water content at wilting point (m^{3} m^{-3})
+
     REAL(r_std), INTENT(out)      :: soilclass(nbpt, nscm)  !! Soil type map to be created from the Zobler map
                                                             !! or a map defining the 12 USDA classes (e.g. Reynolds)
@@ -2660,4 +2788,7 @@
     !  0.3 LOCAL
     !
+    !salma: added the following local variable to be used for all the soil hydraulic parameters
+    REAL(r_std), DIMENSION(nbpt)        :: param            !! to be introduced in function: interpweight
+
     CHARACTER(LEN=80) :: filename
     INTEGER(i_std) :: ib, ilf, nbexp, i
@@ -2676,4 +2807,9 @@
     REAL(r_std), DIMENSION(:,:), ALLOCATABLE             :: textrefrac       !! text fractions re-dimensioned
     REAL(r_std), DIMENSION(nbpt)                         :: atext            !! Availability of the texture interpolation
+    !salma added the following 3 variables to control the SP-MIP simulations
+    REAL(r_std), DIMENSION(nbpt)                         :: aparam            !! Availability of the parameter interpolation
+    CHARACTER(LEN=80)                                    :: spmipexp          !! designing the number of sp-mip experiment
+    CHARACTER(LEN=80)                                    :: exp               !! designing the model of experiment 4 (sp_mip)
+
     REAL(r_std)                                          :: vmin, vmax       !! min/max values to use for the 
 
@@ -2705,439 +2841,398 @@
 
     IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt'
-    !
-    !  Needs to be a configurable variable
-    !
-    !
-    !Config Key   = SOILCLASS_FILE
-    !Config Desc  = Name of file from which soil types are read
-    !Config Def   = soils_param.nc
-    !Config If    = NOT(IMPOSE_VEG)
-    !Config Help  = The name of the file to be opened to read the soil types. 
-    !Config         The data from this file is then interpolated to the grid of
-    !Config         of the model. The aim is to get fractions for sand loam and
-    !Config         clay in each grid box. This information is used for soil hydrology
-    !Config         and respiration.
-    !Config Units = [FILE]
-    !
-    ! soils_param.nc file is 1deg soil texture file (Zobler)
-    ! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution)
-
-    filename = 'soils_param.nc'
-    CALL getin_p('SOILCLASS_FILE',filename)
-
-    variablename = 'soiltext'
-
-    !! Variables for interpweight
-    ! Type of calculation of cell fractions
-    fractype = 'default'
-
-    IF (xios_interpolation) THEN
-       IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " &
-            // TRIM(filename) // " for variable " // TRIM(variablename)
+
+    ! Salma: sp_mip_experiment: select experiment:exp1, exp2, exp3 or exp4 in run.def
+    !   exp1: Reading the soil parameter maps of SPMIP
+    !   exp2: Reading a usda soil texture map, with 12 USDA classes + Oxisols
+    !         Works for SP-MIP textur emap based on SoilGrids, Reynolds, and any map with USDA classes
+    !   exp3: Reading the Zobler soil texture map
+    !   exp4: Imposing uniform soil texture over the globe (4 texture options, with parameter values imposed by SP-MIP)
+    ! We read a soil texture map in all experiments but exp4.
+    ! Even with exp1, some parameters are defined based on texture.
+    
+    ! IMPORTANT: if no exp is defined in run.def, the model works as before, by deriving the soil parameters 
+    ! from a soil texture map, itself defined by the SOILTYPE_CLASSIF keyword, and soil_classif variable:
+    !   if soil_classif = zobler, spmipexp = exp3
+    !   if soil_classif = usda, spmipexp = exp2
+    ! For the SP-MIP experiments made by Salma Tafasca, exp1 was made by reading the SPMIP soil texture map (for clayfraction, etc.)
+    ! But to get a uniform texture (exp 4), you need to select a soil texture map using soil_classif, even if it's not read
+    ! WARNING: if you choose the zobler map with usda/exp2, the simulation will run, but not as you probably wish.
+
+    SELECTCASE (soil_classif) ! copied from src-parameters/control.f90
+    CASE ('zobler','none')
+       spmipexp='exp3'
+    CASE ('usda')
+       spmipexp='exp2'
+    CASE DEFAULT
+       WRITE(numout,*) "Unsupported soil type classification: soil_classif=",soil_classif
+       WRITE(numout,*) "Choose between zobler, usda and none according to the map"
+       CALL ipslerr_p(3,'control_initialize','Bad choice of soil_classif','Choose between zobler, usda and none','')
+    ENDSELECT
+
+    !Config Key   = spmipexp
+    !Config Desc  = number of sp_mip experiment
+    !Config Help  = possible values: exp1, exp2, exp3 and exp4
+    !Config Units = [-]
+    CALL getin_p("SPMIPEXP",spmipexp)
+
+    IF (spmipexp == 'exp4') THEN
+       ! case where exp4 is selected: uniform soil parameters
+       ! the values of the hydraulic parameters below come from SP-MIP,
+       ! and correspond to the Rosetta PTF (Schaap et al., 2001)
+
+       ! sp_mip_experiment_4: select another level of experiment: a, b, c or d in run.def
+
+       !Config Key   = EXP4
+       !Config Desc  = number of sp_mip experiment 4
+       !Config Help  = possible values: a, b, c and d
+       !Config Units = [-]
+       CALL getin_p("EXP4",exp)
+
+       SELECTCASE (exp)
+
+       CASE ('a') ! loamy sand
+          clayfraction=0.06
+          sandfraction=0.81
+          siltfraction=0.13
+         ! njsc=2
+          DO ib=1 , nbpt
+             njsc=2
+             mcr(ib) = 0.049
+             mcs(ib) = 0.39
+             ks(ib) = (1.41e-5)*1000*24*3600
+             avan(ib) = 3.475*(1e-3)
+             nvan(ib) = 1.746
+             mcfc(ib) = 0.1039
+             mcw(ib) = 0.05221
+          ENDDO
+          ! definir clayfraction Ã  partir du barycentre du domaine loamy sand
+          ! definir njsc pour cette classe => pcent dans hydrol
+
+       CASE ('b') !loam
+          clayfraction=0.2
+          sandfraction=0.4
+          siltfraction=0.4
+          njsc=6
+          DO ib=1, nbpt
+             mcr(ib) = 0.061
+             mcs(ib) = 0.399
+             ks(ib) = (3.38e-6)*1000*24*3600
+             avan(ib) = 1.112*(1e-3)
+             nvan(ib) = 1.472
+             mcfc(ib) = 0.236
+             mcw(ib) = 0.09115
+          ENDDO
+          ! definir clayfraction Ã  partir du barycentre du domaine
+          ! definir njsc pour cette classe
+
+       CASE ('c') !silt
+          clayfraction=0.1
+          sandfraction=0.06
+          siltfraction=0.84
+          njsc=5
+          DO ib=1, nbpt
+             mcr(ib) = 0.05
+             mcs(ib) = 0.489
+             ks(ib) = (2.81e-6)*1000*24*3600
+             avan(ib) = 0.6577*(1e-3)
+             nvan(ib) = 1.679
+             mcfc(ib) = 0.2854
+             mcw(ib) = 0.06944
+          ENDDO
+
+          ! definir clayfraction Ã  partir du barycentre du domaine
+          ! definir njsc pour cette classe
+       CASE ('d')!clay
+          clayfraction=0.55
+          sandfraction=0.15
+          siltfraction=0.3
+          njsc=12
+          DO ib=1, nbpt
+             mcr(ib) = 0.098
+             mcs(ib) = 0.459
+             ks(ib) = (9.74e-7)*1000*24*3600
+             avan(ib) = 1.496*(1e-3)
+             nvan(ib) = 1.253
+             mcfc(ib) = 0.3329
+             mcw(ib) = 0.1897
+          ENDDO
+          ! definir clayfraction Ã  partir du barycentre du domaine
+          ! definir njsc pour cette classe
+
+       CASE DEFAULT
+
+          WRITE (numout,*) 'Unsupported experiment number. Choose between a, b, c or d according to sp_mip_experiment_4 number'
+          CALL ipslerr_p(3,'hydrol_init','Unsupported experiment number. ',&
+               'Choose between a,b,c or d','')
+       ENDSELECT
+
+    ELSE !if we are here then it's either exp1, exp2 or exp3
        
-       SELECT CASE(soil_classif)
-
-       CASE('none')
-          ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
-          DO ib=1, nbpt
-             soilclass(ib,:) = soilclass_default_fao
-             clayfraction(ib) = clayfraction_default
-          ENDDO
-
-
-       CASE('zobler')
-
-         !
-          soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
-          !
-          IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
-          !
-          ALLOCATE(textrefrac(nbpt,nzobler))
-          ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
-          CALL get_soilcorr_zobler (nzobler, textfrac_table)
-       
-          CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
-          CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
-          CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
-          CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
-          CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
-          CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
-          CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
-
-
-       
-          CALL get_soilcorr_zobler (nzobler, textfrac_table)
-        !
-        !
-          DO ib =1, nbpt
-            soilclass(ib,1)=textrefrac(ib,1)
-            soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
-            soilclass(ib,3)=textrefrac(ib,5)
-            
-            ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
-            ! over the zobler pixels composing the ORCHIDEE grid-cell
-            clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + &
-                               textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + &
-                               textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7)
-
-            sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + &
-                               textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + &
-                               textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7)
-
-            siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + &
-                               textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + &
-                               textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7)
-
-            sgn=SUM(soilclass(ib,1:3))
-            
-            IF (sgn < min_sechiba) THEN
-              soilclass(ib,:) = soilclass_default(:)
-              clayfraction(ib) = clayfraction_default
-              sandfraction(ib) = sandfraction_default
-              siltfraction(ib) = siltfraction_default
-              atext(ib)=0.
-            ELSE
-              atext(ib)=sgn
-              clayfraction(ib) = clayfraction(ib) / sgn
-              sandfraction(ib) = sandfraction(ib) / sgn
-              siltfraction(ib) = siltfraction(ib) / sgn
-              soilclass(ib,1:3) = soilclass(ib,1:3) / sgn
-            ENDIF
-            
-          ENDDO
-          
-          
-       
-       CASE('usda')
-
-           IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda'
- 
-           soilclass_default=soilclass_default_usda
-           !
-           WRITE(numout,*) "Using a soilclass map with usda classification"
-           !
-           ALLOCATE(textrefrac(nbpt,nscm))
-           ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
-           IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
-    
-           CALL get_soilcorr_usda (nscm, textfrac_table)
-    
-           IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
-
-          CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
-          CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
-          CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
-          CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
-          CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
-          CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
-          CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
-          CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8))
-          CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9))
-          CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10))
-          CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11))
-          CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12))
-
-
-       
-          CALL get_soilcorr_usda (nscm, textfrac_table)
-          IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'      
-
-          DO ib =1, nbpt
-            clayfraction(ib) = 0.0
-            DO ilf = 1,nscm
-              soilclass(ib,ilf)=textrefrac(ib,ilf)      
-              clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf)
-              sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf)
-              siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf)
-            ENDDO
-
-
-            sgn=SUM(soilclass(ib,:))
-            
-            IF (sgn < min_sechiba) THEN
-              soilclass(ib,:) = soilclass_default(:)
-              clayfraction(ib) = clayfraction_default
-              sandfraction(ib) = sandfraction_default
-              siltfraction(ib) = siltfraction_default
-              atext(ib)=0
-            ELSE
-              soilclass(ib,:) = soilclass(ib,:) / sgn
-              clayfraction(ib) = clayfraction(ib) / sgn
-              sandfraction(ib) = sandfraction(ib) / sgn
-              siltfraction(ib) = siltfraction(ib) / sgn
-              atext(ib)=sgn
-            ENDIF
-          ENDDO
-
-        CASE DEFAULT
-             WRITE(numout,*) 'slowproc_soilt:'
-             WRITE(numout,*) '  A non supported soil type classification has been chosen'
-             CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
-        END SELECT
-
-
-
-    ELSE              !    xios_interpolation 
-       ! Read and interpolate using stardard method with IOIPSL and aggregate
-    
-       IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
-            // TRIM(filename) // " for variable " // TRIM(variablename)
-
-
-       ! Name of the longitude and latitude in the input file
-       lonname = 'nav_lon'
-       latname = 'nav_lat'
-
-       IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " &
-            // TRIM(filename) // " for variable " // TRIM(variablename)
-
-       IF ( TRIM(soil_classif) /= 'none' ) THEN
-
-          ! Define a variable for the number of soil textures in the input file
-          SELECTCASE(soil_classif)
+       !  Needs to be a configurable variable
+       !
+       !Config Key   = SOILCLASS_FILE
+       !Config Desc  = Name of file from which soil types are read
+       !Config Def   = soils_param.nc
+       !Config If    = NOT(IMPOSE_VEG)
+       !Config Help  = The name of the file to be opened to read the soil types.
+       !Config         The data from this file is then interpolated to the grid of
+       !Config         of the model. The aim is to get fractions for sand loam and
+       !Config         clay in each grid box. This information is used for soil hydrology
+       !Config         and respiration.
+       !Config Units = [FILE]
+       !
+       ! soils_param.nc file is 1deg soil texture file (Zobler)
+       ! The USDA map from Reynolds is soils_param_usda.nc (1/12deg resolution)
+
+       filename = 'soils_param.nc'
+       CALL getin_p('SOILCLASS_FILE',filename)
+
+       variablename = 'soiltext'
+
+       !! Variables for interpweight
+       ! Type of calculation of cell fractions
+       fractype = 'default'
+
+       IF (xios_interpolation) THEN
+          IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Use XIOS to read and interpolate " &
+               // TRIM(filename) // " for variable " // TRIM(variablename)
+
+          SELECT CASE(soil_classif)
+
+          CASE('none')
+             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
+             DO ib=1, nbpt
+                soilclass(ib,:) = soilclass_default_fao
+                clayfraction(ib) = clayfraction_default
+             ENDDO
+
           CASE('zobler')
-             ntextinfile=nzobler
+             !
+             soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
+             !
+             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification, to be read using XIOS"
+             !
+             ALLOCATE(textrefrac(nbpt,nzobler))
+             ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
+             CALL get_soilcorr_zobler (nzobler, textfrac_table)
+             CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
+             CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
+             CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
+             CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
+             CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
+             CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
+             CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
+
+             CALL get_soilcorr_zobler (nzobler, textfrac_table)
+             !             !
+             DO ib =1, nbpt
+                soilclass(ib,1)=textrefrac(ib,1)
+                soilclass(ib,2)=textrefrac(ib,2)+textrefrac(ib,3)+textrefrac(ib,4)+textrefrac(ib,7)
+                soilclass(ib,3)=textrefrac(ib,5)
+
+                ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
+                ! over the zobler pixels composing the ORCHIDEE grid-cell
+                clayfraction(ib) = textfrac_table(1,3) * textrefrac(ib,1)+textfrac_table(2,3) * textrefrac(ib,2) + &
+                     textfrac_table(3,3) * textrefrac(ib,3)+textfrac_table(4,3) * textrefrac(ib,4) + &
+                     textfrac_table(5,3) * textrefrac(ib,5)+textfrac_table(7,3) * textrefrac(ib,7)
+
+                sandfraction(ib) = textfrac_table(1,2) * textrefrac(ib,1)+textfrac_table(2,2) * textrefrac(ib,2) + &
+                     textfrac_table(3,2) * textrefrac(ib,3)+textfrac_table(4,2) * textrefrac(ib,4) + &
+                     textfrac_table(5,2) * textrefrac(ib,5)+textfrac_table(7,2) * textrefrac(ib,7)
+
+                siltfraction(ib) = textfrac_table(1,1) * textrefrac(ib,1)+textfrac_table(2,1) * textrefrac(ib,2) + &
+                     textfrac_table(3,1) * textrefrac(ib,3)+textfrac_table(4,1) * textrefrac(ib,4) + &
+                     textfrac_table(5,1) * textrefrac(ib,5)+textfrac_table(7,1) * textrefrac(ib,7)
+
+                sgn=SUM(soilclass(ib,1:3))
+
+                IF (sgn < min_sechiba) THEN
+                   soilclass(ib,:) = soilclass_default(:)
+                   clayfraction(ib) = clayfraction_default
+                   sandfraction(ib) = sandfraction_default
+                   siltfraction(ib) = siltfraction_default
+                   atext(ib)=0.
+                ELSE
+                   atext(ib)=sgn
+                   clayfraction(ib) = clayfraction(ib) / sgn
+                   sandfraction(ib) = sandfraction(ib) / sgn
+                   siltfraction(ib) = siltfraction(ib) / sgn
+                   soilclass(ib,1:3) = soilclass(ib,1:3) / sgn
+                ENDIF
+
+             ENDDO
+
           CASE('usda')
-             ntextinfile=nscm
+
+             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: start case usda'
+
+             soilclass_default=soilclass_default_usda
+             !
+             WRITE(numout,*) "Using a soilclass map with usda classification, to be read using XIOS"
+             !
+             ALLOCATE(textrefrac(nbpt,nscm))
+             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
+
+             CALL get_soilcorr_usda (nscm, textfrac_table)
+
+             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
+
+             CALL xios_orchidee_recv_field('soiltext1',textrefrac(:,1))
+             CALL xios_orchidee_recv_field('soiltext2',textrefrac(:,2))
+             CALL xios_orchidee_recv_field('soiltext3',textrefrac(:,3))
+             CALL xios_orchidee_recv_field('soiltext4',textrefrac(:,4))
+             CALL xios_orchidee_recv_field('soiltext5',textrefrac(:,5))
+             CALL xios_orchidee_recv_field('soiltext6',textrefrac(:,6))
+             CALL xios_orchidee_recv_field('soiltext7',textrefrac(:,7))
+             CALL xios_orchidee_recv_field('soiltext8',textrefrac(:,8))
+             CALL xios_orchidee_recv_field('soiltext9',textrefrac(:,9))
+             CALL xios_orchidee_recv_field('soiltext10',textrefrac(:,10))
+             CALL xios_orchidee_recv_field('soiltext11',textrefrac(:,11))
+             CALL xios_orchidee_recv_field('soiltext12',textrefrac(:,12))        
+             CALL xios_orchidee_recv_field('soiltext13',textrefrac(:,13))
+
+             CALL get_soilcorr_usda (nscm, textfrac_table)
+             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
+
+             DO ib =1, nbpt
+                clayfraction(ib) = 0.0
+                DO ilf = 1,nscm
+                   soilclass(ib,ilf)=textrefrac(ib,ilf)
+                   clayfraction(ib) = clayfraction(ib) + textfrac_table(ilf,3)*textrefrac(ib,ilf)
+                   sandfraction(ib) = sandfraction(ib) + textfrac_table(ilf,2)*textrefrac(ib,ilf)
+                   siltfraction(ib) = siltfraction(ib) + textfrac_table(ilf,1)*textrefrac(ib,ilf)
+                ENDDO
+
+                sgn=SUM(soilclass(ib,:))
+
+                IF (sgn < min_sechiba) THEN
+                   soilclass(ib,:) = soilclass_default(:)
+                   clayfraction(ib) = clayfraction_default
+                   sandfraction(ib) = sandfraction_default
+                   siltfraction(ib) = siltfraction_default
+                   atext(ib)=0
+                ELSE
+                   soilclass(ib,:) = soilclass(ib,:) / sgn
+                   clayfraction(ib) = clayfraction(ib) / sgn
+                   sandfraction(ib) = sandfraction(ib) / sgn
+                   siltfraction(ib) = siltfraction(ib) / sgn
+                   atext(ib)=sgn
+                ENDIF
+             ENDDO
+
           CASE DEFAULT
              WRITE(numout,*) 'slowproc_soilt:'
              WRITE(numout,*) '  A non supported soil type classification has been chosen'
              CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
-          ENDSELECT
-
-          ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',&
-            '','')
-
-          ! Assigning values to vmin, vmax
-          vmin = un
-          vmax = ntextinfile*un
-          
-          ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','')
-          variabletypevals = -un
-          
-          !! Variables for interpweight
-          ! Should negative values be set to zero from input file?
-          nonegative = .FALSE.
-          ! Type of mask to apply to the input data (see header for more details)
-          maskingtype = 'mabove'
-          ! Values to use for the masking
-          maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
-          ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
-          namemaskvar = ''
-          
-          CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours,        &
-             contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,    & 
-             maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext)
-
-          ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','')
-          ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','')
-          
-          IF (printlev_loc >= 5) THEN
-             WRITE(numout,*)'  slowproc_soilt after interpweight_2D'
-             WRITE(numout,*)'  slowproc_soilt before starting loop nbpt:', nbpt
-             WRITE(numout,*)"  slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..."
+          END SELECT
+
+       ELSE              !    xios_interpolation
+          ! Read and interpolate using stardard method with IOIPSL and aggregate
+
+          IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
+               // TRIM(filename) // " for variable " // TRIM(variablename)
+
+          ! Name of the longitude and latitude in the input file
+          lonname = 'nav_lon'
+          latname = 'nav_lat'
+
+          IF (printlev_loc >= 2) WRITE(numout,*) "slowproc_soilt: Start interpolate " &
+               // TRIM(filename) // " for variable " // TRIM(variablename)
+
+          IF ( TRIM(soil_classif) /= 'none' ) THEN
+
+             ! Define a variable for the number of soil textures in the input file
+             SELECTCASE(soil_classif)
+             CASE('zobler')
+                ntextinfile=nzobler
+             CASE('usda')
+                ntextinfile=nscm
+             CASE DEFAULT
+                WRITE(numout,*) 'slowproc_soilt:'
+                WRITE(numout,*) '  A non supported soil type classification has been chosen'
+                CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
+             ENDSELECT
+
+             ALLOCATE(textrefrac(nbpt,ntextinfile), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable textrefrac',&
+                  '','')
+
+             ! Assigning values to vmin, vmax
+             vmin = un
+             vmax = ntextinfile*un
+
+             ALLOCATE(variabletypevals(ntextinfile), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variabletypevals','','')
+             variabletypevals = -un
+
+             !! Variables for interpweight
+             ! Should negative values be set to zero from input file?
+             nonegative = .FALSE.
+             ! Type of mask to apply to the input data (see header for more details)
+             maskingtype = 'mabove'
+             ! Values to use for the masking
+             maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
+             ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
+
+             namemaskvar = ''
+
+             CALL interpweight_2D(nbpt, ntextinfile, variabletypevals, lalo, resolution, neighbours,        &
+                  contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,    &
+                  maskvals, namemaskvar, 0, 0, -1, fractype, -1., -1., textrefrac, atext)
+
+             ALLOCATE(vecpos(ntextinfile), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable vecpos','','')
+             ALLOCATE(solt(ntextinfile), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'slowproc_soilt','Problem in allocation of variable solt','','')
+
+             IF (printlev_loc >= 5) THEN
+                WRITE(numout,*)'  slowproc_soilt after interpweight_2D'
+                WRITE(numout,*)'  slowproc_soilt before starting loop nbpt:', nbpt
+                WRITE(numout,*)"  slowproc_soilt starting classification '" // TRIM(soil_classif) // "'..."
+             END IF
+          ELSE
+             IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt using default values all points are propertly ' // &
+                  'interpolated atext = 1. everywhere!'
+             atext = 1.
           END IF
-       ELSE
-         IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt using default values all points are propertly ' // &
-           'interpolated atext = 1. everywhere!'
-         atext = 1.
-       END IF
-
-    nbexp = 0
-    SELECTCASE(soil_classif)
-    CASE('none')
-       ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
-       IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
-       DO ib=1, nbpt
-          soilclass(ib,:) = soilclass_default_fao
-          clayfraction(ib) = clayfraction_default
-          sandfraction(ib) = sandfraction_default
-          siltfraction(ib) = siltfraction_default
-       ENDDO
-    CASE('zobler')
-       !
-       soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
-       !
-       IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
-       !
-       ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
-       IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
-       CALL get_soilcorr_zobler (nzobler, textfrac_table)
-       !
-       !
-       IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt after getting table of textures'
-       DO ib =1, nbpt
-          soilclass(ib,:) = zero
-          clayfraction(ib) = zero
-          sandfraction(ib) = zero
-          siltfraction(ib) = zero
-          !
-          ! vecpos: List of positions where textures were not zero
-          ! vecpos(1): number of not null textures found
-          vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq')
-          fopt = vecpos(1)
-
-          IF ( fopt .EQ. 0 ) THEN
-             ! No points were found for current grid box, use default values
-             nbexp = nbexp + 1
-             soilclass(ib,:) = soilclass_default(:)
-             clayfraction(ib) = clayfraction_default
-             sandfraction(ib) = sandfraction_default
-             siltfraction(ib) = siltfraction_default
-          ELSE
-             IF (fopt == nzobler) THEN
-                ! All textures are not zero
-                solt=(/(i,i=1,nzobler)/)
-             ELSE
-               DO ilf = 1,fopt
-                 solt(ilf) = vecpos(ilf+1)
-               END DO
-             END IF
-             !
-             !   Compute the fraction of each textural class
-             !
-             sgn = 0.
-             DO ilf = 1,fopt
-                   !
-                   ! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE
-                   ! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine)
-                   ! For type 6 = glacier, default values are set and it is also taken into account during the normalization 
-                   ! of the fractions (done in interpweight_2D)
-                   ! Note that type 0 corresponds to ocean but it is already removed using the mask above.
-                   !
-                IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. & 
-                     (solt(ilf) .NE. 6) ) THEN
-                   SELECT CASE(solt(ilf))
-                     CASE(1)
-                        soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf))
-                     CASE(2)
-                        soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
-                     CASE(3)
-                        soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
-                     CASE(4)
-                        soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
-                     CASE(5)
-                        soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf))
-                     CASE(7)
-                        soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
-                     CASE DEFAULT
-                        WRITE(numout,*) 'We should not be here, an impossible case appeared'
-                        CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','')
-                   END SELECT
-                   ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
-                   ! over the zobler pixels composing the ORCHIDEE grid-cell
-                   clayfraction(ib) = clayfraction(ib) + &
-                        & textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf))
-                   sandfraction(ib) = sandfraction(ib) + &
-                        & textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf))
-                   siltfraction(ib) = siltfraction(ib) + &
-                        & textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf))
-                   ! Sum the fractions which are not glaciers nor ocean
-                   sgn = sgn + textrefrac(ib,solt(ilf))
-                ELSE
-                   IF (solt(ilf) .GT. nzobler) THEN
-                      WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
-                      CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','')
-                   ENDIF
-                END IF
-             ENDDO
-
-             IF ( sgn .LT. min_sechiba) THEN
-                ! Set default values if grid cells were only covered by glaciers or ocean 
-                ! or if now information on the source grid was found.
-                nbexp = nbexp + 1
-                soilclass(ib,:) = soilclass_default(:)
+
+          nbexp = 0
+          SELECTCASE(soil_classif)
+          CASE('none')
+             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
+             DO ib=1, nbpt
+                soilclass(ib,:) = soilclass_default_fao
                 clayfraction(ib) = clayfraction_default
                 sandfraction(ib) = sandfraction_default
                 siltfraction(ib) = siltfraction_default
-             ELSE
-                ! Normalize using the fraction of surface not including glaciers and ocean
-                soilclass(ib,:) = soilclass(ib,:)/sgn
-                clayfraction(ib) = clayfraction(ib)/sgn
-                sandfraction(ib) = sandfraction(ib)/sgn
-                siltfraction(ib) = siltfraction(ib)/sgn
-             ENDIF
-          ENDIF
-       ENDDO
-      
-       ! The "USDA" case reads a map of the 12 USDA texture classes, 
-       ! such as to assign the corresponding soil properties
-       CASE("usda")
-          IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification"
-
-          soilclass_default=soilclass_default_usda
-
-          ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
-          IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
-
-          CALL get_soilcorr_usda (nscm, textfrac_table)
-
-          IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
-          !
-          DO ib =1, nbpt
-          ! GO through the point we have found
-          !
-          !
-          ! Provide which textures were found
-          ! vecpos: List of positions where textures were not zero
-          !   vecpos(1): number of not null textures found
-          vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq')
-          fopt = vecpos(1)
-         
-          !
-          !    Check that we found some points
-          !
-          soilclass(ib,:) = 0.0
-          clayfraction(ib) = 0.0
-          
-          IF ( fopt .EQ. 0) THEN
-             ! No points were found for current grid box, use default values
-             IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib
-             nbexp = nbexp + 1
-             soilclass(ib,:) = soilclass_default
-             clayfraction(ib) = clayfraction_default
-             sandfraction(ib) = sandfraction_default
-             siltfraction(ib) = siltfraction_default
-          ELSE
-             IF (fopt == nscm) THEN
-                ! All textures are not zero
-                solt(:) = (/(i,i=1,nscm)/)
-             ELSE
-               DO ilf = 1,fopt
-                 solt(ilf) = vecpos(ilf+1) 
-               END DO
-             END IF
-            
+             ENDDO
+          CASE('zobler')
+             !
+             soilclass_default=soilclass_default_fao ! FAO means here 3 final texture classes
+             !
+             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with Zobler classification"
+             !
+             ALLOCATE(textfrac_table(nzobler,ntext), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
+             CALL get_soilcorr_zobler (nzobler, textfrac_table)
+                       
+             IF (printlev_loc >= 5) WRITE(numout,*)'  slowproc_soilt after getting table of textures'
+             DO ib =1, nbpt
+                soilclass(ib,:) = zero
+                clayfraction(ib) = zero
+                sandfraction(ib) = zero
+                siltfraction(ib) = zero
                 !
-                !
-                !   Compute the fraction of each textural class  
-                !
-                !
-                DO ilf = 1,fopt
-                   IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
-                      soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf))
-                      clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) *                &
-                           textrefrac(ib,solt(ilf))
-                      sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * &
-                           textrefrac(ib,solt(ilf))
-                      siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * &
-                        textrefrac(ib,solt(ilf))
-                   ELSE
-                      IF (solt(ilf) .GT. nscm) THEN
-                         WRITE(*,*) 'The file contains a soil color class which is incompatible with this program'
-                         CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','')
-                      ENDIF
-                   ENDIF
-                   !
-                ENDDO
-
-                ! Set default values if the surface in source file is too small
-                IF ( atext(ib) .LT. min_sechiba) THEN
+                ! vecpos: List of positions where textures were not zero
+                ! vecpos(1): number of not null textures found
+                vecpos = interpweight_ValVecR(textrefrac(ib,:),nzobler,zero,'neq')
+                fopt = vecpos(1)
+
+                IF ( fopt .EQ. 0 ) THEN
+                   ! No points were found for current grid box, use default values
                    nbexp = nbexp + 1
                    soilclass(ib,:) = soilclass_default(:)
@@ -3145,35 +3240,305 @@
                    sandfraction(ib) = sandfraction_default
                    siltfraction(ib) = siltfraction_default
+
+                ELSE
+                   IF (fopt == nzobler) THEN
+                      ! All textures are not zero
+                      solt=(/(i,i=1,nzobler)/)
+                   ELSE
+                      DO ilf = 1,fopt
+                         solt(ilf) = vecpos(ilf+1)
+                      END DO
+                   END IF
+                   !
+                   !   Compute the fraction of each textural class
+                   !
+                   sgn = 0.
+                   DO ilf = 1,fopt
+                      !
+                      ! Here we make the correspondance between the 7 zobler textures and the 3 textures in ORCHIDEE
+                      ! and soilclass correspond to surfaces covered by the 3 textures of ORCHIDEE (coase,medium,fine)
+                      ! For type 6 = glacier, default values are set and it is also taken into account during the normalization
+                      ! of the fractions (done in interpweight_2D)
+                      ! Note that type 0 corresponds to ocean but it is already removed using the mask above.
+                     !
+                      IF ( (solt(ilf) .LE. nzobler) .AND. (solt(ilf) .GT. 0) .AND. &
+                           (solt(ilf) .NE. 6) ) THEN
+                         SELECT CASE(solt(ilf))
+                         CASE(1)
+                            soilclass(ib,1) = soilclass(ib,1) + textrefrac(ib,solt(ilf))
+                         CASE(2)
+                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                         CASE(3)
+                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                         CASE(4)
+                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                         CASE(5)
+                            soilclass(ib,3) = soilclass(ib,3) + textrefrac(ib,solt(ilf))
+                         CASE(7)
+                            soilclass(ib,2) = soilclass(ib,2) + textrefrac(ib,solt(ilf))
+                         CASE DEFAULT
+                            WRITE(numout,*) 'We should not be here, an impossible case appeared'
+                            CALL ipslerr_p(3,'slowproc_soilt','Bad value for solt','','')
+                         END SELECT
+                         ! clayfraction is the sum of the % of clay (as a mineral of small granulometry, and not as a texture)
+                         ! over the zobler pixels composing the ORCHIDEE grid-cell
+                         clayfraction(ib) = clayfraction(ib) + &
+                              & textfrac_table(solt(ilf),3) * textrefrac(ib,solt(ilf))
+                         sandfraction(ib) = sandfraction(ib) + &
+                              & textfrac_table(solt(ilf),2) * textrefrac(ib,solt(ilf))
+                         siltfraction(ib) = siltfraction(ib) + &
+                              & textfrac_table(solt(ilf),1) * textrefrac(ib,solt(ilf))
+                         ! Sum the fractions which are not glaciers nor ocean
+                         sgn = sgn + textrefrac(ib,solt(ilf))
+                      ELSE
+                         IF (solt(ilf) .GT. nzobler) THEN
+                            WRITE(numout,*) 'The file contains a soil color class which is incompatible with this program'
+                            CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible','','')
+                         ENDIF
+                      END IF
+                   ENDDO
+
+                   IF ( sgn .LT. min_sechiba) THEN
+                      ! Set default values if grid cells were only covered by glaciers or ocean
+                      ! or if now information on the source grid was found.
+                      nbexp = nbexp + 1
+                      soilclass(ib,:) = soilclass_default(:)
+                      clayfraction(ib) = clayfraction_default
+                      sandfraction(ib) = sandfraction_default
+                      siltfraction(ib) = siltfraction_default
+                   ELSE
+                      ! Normalize using the fraction of surface not including glaciers and ocean
+                      soilclass(ib,:) = soilclass(ib,:)/sgn
+                      clayfraction(ib) = clayfraction(ib)/sgn
+                      sandfraction(ib) = sandfraction(ib)/sgn
+                      siltfraction(ib) = siltfraction(ib)/sgn
+                   ENDIF
                 ENDIF
-             ENDIF
-
-          ENDDO
-          
-          IF (printlev_loc>=4) WRITE (numout,*) '  slowproc_soilt: End case usda'
-          
-       CASE DEFAULT
-          WRITE(numout,*) 'slowproc_soilt _______'
-          WRITE(numout,*) '  A non supported soil type classification has been chosen'
-          CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
-       ENDSELECT
-       IF (printlev_loc >= 5 ) WRITE(numout,*)'  slowproc_soilt end of type classification'
-
-       IF ( nbexp .GT. 0 ) THEN
-          WRITE(numout,*) 'slowproc_soilt:'
-          WRITE(numout,*) '  The interpolation of the bare soil albedo had ', nbexp
-          WRITE(numout,*) '  points without data. This are either coastal points or ice covered land.'
-          WRITE(numout,*) '  The problem was solved by using the default soil types.'
-       ENDIF
-
-       IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals)
-       IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac)
-       IF (ALLOCATED(solt)) DEALLOCATE (solt)
-       IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table)
-    
-    ENDIF        !      xios_interpolation 
-    
+             ENDDO
+
+             ! The "USDA" case reads a map of the 12 USDA texture classes,
+             ! such as to assign the corresponding soil properties
+          CASE("usda")
+             IF (printlev_loc>=2) WRITE(numout,*) "Using a soilclass map with usda classification"
+
+             soilclass_default=soilclass_default_usda
+
+
+             ALLOCATE(textfrac_table(nscm,ntext), STAT=ALLOC_ERR)
+             IF (ALLOC_ERR/=0) CALL ipslerr_p(3,'slowproc_soilt','Error in allocation for textfrac_table','','')
+
+             CALL get_soilcorr_usda (nscm, textfrac_table)
+
+             IF (printlev_loc>=4) WRITE (numout,*) 'slowproc_soilt: After get_soilcorr_usda'
+             !
+             DO ib =1, nbpt
+                ! GO through the point we have found
+                !
+                ! Provide which textures were found
+                ! vecpos: List of positions where textures were not zero
+                !   vecpos(1): number of not null textures found
+                vecpos = interpweight_ValVecR(textrefrac(ib,:),ntextinfile,zero,'neq')
+                fopt = vecpos(1)
+                !
+                !    Check that we found some points
+                !
+                soilclass(ib,:) = 0.0
+                clayfraction(ib) = 0.0
+                sandfraction(ib) = 0.0
+                siltfraction(ib) = 0.0
+
+                IF ( fopt .EQ. 0) THEN
+                   ! No points were found for current grid box, use default values
+                   IF (printlev_loc>=3) WRITE(numout,*)'slowproc_soilt: no soil class in input file found for point=', ib
+                   nbexp = nbexp + 1
+                   soilclass(ib,:) = soilclass_default
+                   clayfraction(ib) = clayfraction_default
+                   sandfraction(ib) = sandfraction_default
+                   siltfraction(ib) = siltfraction_default
+                ELSE
+                   IF (fopt == nscm) THEN
+                      ! All textures are not zero
+                      solt(:) = (/(i,i=1,nscm)/)
+                   ELSE
+                      DO ilf = 1,fopt
+                         solt(ilf) = vecpos(ilf+1)
+                      END DO
+                   END IF
+
+                   !   Compute the fraction of each textural class
+                   DO ilf = 1,fopt
+                      IF ( (solt(ilf) .LE. nscm) .AND. (solt(ilf) .GT. 0) ) THEN
+                         soilclass(ib,solt(ilf)) = textrefrac(ib,solt(ilf))
+                         clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) *                &
+                              textrefrac(ib,solt(ilf))
+                         sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * &
+                              textrefrac(ib,solt(ilf))
+                         siltfraction(ib) = siltfraction(ib) + textfrac_table(solt(ilf),1) * &
+                              textrefrac(ib,solt(ilf))
+                      ELSE
+                         IF (solt(ilf) .GT. nscm) THEN
+                            WRITE(*,*) 'The file contains a soil color class which is incompatible with this program'
+                            CALL ipslerr_p(3,'slowproc_soilt','Problem soil color class incompatible 2','','')
+                         ENDIF
+                      ENDIF
+                      !
+                   ENDDO
+
+                   ! Set default values if the surface in source file is too small
+                   IF ( atext(ib) .LT. min_sechiba) THEN
+                      nbexp = nbexp + 1
+                      soilclass(ib,:) = soilclass_default(:)
+                      clayfraction(ib) = clayfraction_default
+                      sandfraction(ib) = sandfraction_default
+                      siltfraction(ib) = siltfraction_default
+                   ENDIF
+                ENDIF
+
+             ENDDO
+
+             IF (printlev_loc>=4) WRITE (numout,*) '  slowproc_soilt: End case usda'
+
+          CASE DEFAULT
+             WRITE(numout,*) 'slowproc_soilt _______'
+             WRITE(numout,*) '  A non supported soil type classification has been chosen'
+             CALL ipslerr_p(3,'slowproc_soilt','non supported soil type classification','','')
+          ENDSELECT
+          IF (printlev_loc >= 5 ) WRITE(numout,*)'  slowproc_soilt end of type classification'
+
+          IF ( nbexp .GT. 0 ) THEN
+             WRITE(numout,*) 'slowproc_soilt:'
+             WRITE(numout,*) '  The interpolation of variable soiltext had ', nbexp
+             WRITE(numout,*) '  points without data. This are either coastal points or ice covered land.'
+             WRITE(numout,*) '  The problem was solved by using the default soil types.'
+          ENDIF
+
+          IF (ALLOCATED(variabletypevals)) DEALLOCATE (variabletypevals)
+          IF (ALLOCATED(textrefrac)) DEALLOCATE (textrefrac)
+          IF (ALLOCATED(solt)) DEALLOCATE (solt)
+          IF (ALLOCATED(textfrac_table)) DEALLOCATE (textfrac_table)
+
+       ENDIF        !      xios_interpolation
+
+       !salma: calculate njsc
+       njsc(:) = 0
+       DO ib = 1, nbpt
+          njsc(ib) = MAXLOC(soilclass(ib,:),1)
+       ENDDO
+
+       IF (spmipexp == 'exp1') THEN
+              IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt: Read soil hydraulic parameters with IOIPSL'
+
+              ! Read using IOIPSL and interpolate using aggregate tool in ORCHIDEE
+
+              !Config Key   = PARAM_FILE
+              !Config Desc  = Name of file from which soil parameter  values are read
+              !Config Def   = params_sp_mip.nc
+              !Config Help  = The name of the file to be opened to read values of parameters.
+              !Config         The data from this file is then interpolated to the grid of
+              !Config         of the model.
+              !Config Units = [FILE]
+              !
+              ! params_sp_mip.nc file is 0.5 deg soil hydraulic parameters file provided by sp_mip
+
+              filename = 'params_sp_mip.nc'
+              CALL getin_p('PARAM_FILE',filename)
+
+              !! Variables for interpweight
+              ! Type of calculation of cell fractions
+              fractype = 'default'
+              ! Name of the longitude and latitude in the input file
+              lonname = 'nav_lon'
+              latname = 'nav_lat'
+              ! Assigning values to vmin, vmax (there are not types/categories
+              vmin =0.
+              vmax = 99999.
+              !! Variables for interpweight
+              ! Should negative values be set to zero from input file?
+              nonegative = .FALSE.
+              ! Type of mask to apply to the input data (see header for more details)
+              maskingtype = 'mabove'
+              ! Values to use for the masking
+              maskvals = (/ min_sechiba, undef_sechiba, undef_sechiba /)
+              ! Name of the variable with the values for the mask in the input file (only if maskkingtype='var') ( not used)
+              namemaskvar = ''
+
+              variablename = 'ks'
+              IF (printlev_loc >= 1) WRITE(numout,*) "slowproc_soilt: Read and interpolate " &
+                   // TRIM(filename) // " for variable " // TRIM(variablename)
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   ks, aparam)
+              WRITE(numout,*) 'ks map is read _______'
+
+              variablename = 'alpha'
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   avan, aparam)
+              WRITE(numout,*) 'avan map read _______'
+
+              variablename = 'thetar'
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   mcr, aparam)
+              WRITE(numout,*) 'thetar map read _______'
+
+              variablename = 'thetas'
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   mcs, aparam)
+              WRITE(numout,*) 'thetas map read _______'
+
+              variablename = 'thetapwpvg'
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   mcw, aparam)
+              WRITE(numout,*) 'thetapwpvg map read _______'
+
+              variablename = 'thetafcvg'
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   mcfc, aparam)
+              WRITE(numout,*) 'thetafcvg map read _______'
+
+              variablename = 'nvg'
+              CALL interpweight_2Dcont(nbpt, 0, 0, lalo, resolution, neighbours,                                &
+                   contfrac, filename, variablename, lonname, latname, vmin, vmax, nonegative, maskingtype,     &
+                   maskvals, namemaskvar, -1, fractype, 0., 0.,                              &
+                   nvan, aparam)
+              WRITE(numout,*) 'nvan map read _______'
+
+      ELSE
+          IF (spmipexp == 'exp2') THEN
+            ! Texture map from SP-MIP, thus Soilgrids modified
+            nvan(:) = nvan_usda(njsc(:))
+            avan(:) = avan_usda(njsc(:))
+            mcr(:) = mcr_usda(njsc(:))
+            mcs(:) = mcs_usda(njsc(:))
+            ks(:) = ks_usda(njsc(:))
+            mcfc(:) = mcf_usda(njsc(:))
+            mcw(:) = mcw_usda(njsc(:))
+            ! on aura pcent(:) = pcent(njsc(:)) dans hydrol
+         ELSE
+            ! salma: here we are in exp3 -- Zobler map
+            nvan(:) = nvan_fao(njsc(:))
+            avan(:) = avan_fao(njsc(:))
+            mcr(:) = mcr_fao(njsc(:))
+            mcs(:) = mcs_fao(njsc(:))
+            ks(:) = ks_fao(njsc(:))
+            mcfc(:) = mcf_fao(njsc(:))
+            mcw(:) = mcw_fao(njsc(:))
+
+          ENDIF !if EXP2
+       ENDIF !if EXP1
+     ENDIF !if EXP4
+             
     ! Write diagnostics
     CALL xios_orchidee_send_field("atext",atext)
-
     CALL xios_orchidee_send_field("interp_diag_atext",atext)
     CALL xios_orchidee_send_field("interp_diag_soilclass",soilclass)
@@ -3585,5 +3950,5 @@
     !! 0.4 Local variables
 
-    INTEGER(i_std),PARAMETER :: nbtypes_usda = 12                    !! Number of USDA texture classes (unitless)
+    INTEGER(i_std),PARAMETER :: nbtypes_usda = 13                    !! Number of USDA texture classes (unitless)
     INTEGER(i_std) :: n                                              !! Index (unitless)
     
@@ -3619,4 +3984,5 @@
     textfrac_table(11,2:3) = (/ 0.06, 0.46 /) ! Silty Clay
     textfrac_table(12,2:3) = (/ 0.15, 0.55 /) ! Clay
+    textfrac_table(13,2:3) = (/ 0.15, 0.55 /) ! Clay
 
     ! Fraction of silt
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 6594)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/thermosoil.f90	(revision 6954)
@@ -389,5 +389,5 @@
        so_capa_dry_ns(:) = so_capa_dry_ns_fao(:)
        mcs(:) = mcs_fao(:)
-    CASE (12)
+    CASE (13) !Salma changed from 12 to 13 for the new class Oxisols
        SMCMAX(:) = SMCMAX_usda(:)
       QZ(:) = QZ_usda(:)
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/context_input_orchidee.xml
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/context_input_orchidee.xml	(revision 6594)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/context_input_orchidee.xml	(revision 6954)
@@ -5,6 +5,7 @@
 <!-- All specifications specific for the reading and interpolating of input files are regrouped in this file.    -->
 <!-- Find below sections for file_definition, field_definition, domain_definition, grid_definition and           -->
-<!-- variable_definition.                                                                                        -->
-<!--                                                                                                             -->
+<!-- variable_definition.                                                                                        -->   
+<!-- Recent change:  added block to read params_sp_mip.nc (Salma Tafasca and Agnes Ducharne, November 2020)      -->
+<!--                                                                                                             --> 
 <!-- $HeadURL$                -->
 <!-- $Date$                                                      -->
@@ -30,4 +31,16 @@
       <field id="soiltext" name="soiltext" domain_ref="soil" operation="instant" freq_offset="1ts"/>
     </file>
+
+    <!-- Default file name: params_sp_mip.nc for both all soil hydraulic parameters. -->
+    <file id="soilparam_file" mode="read" output_freq="1y" enabled="true">
+      <field id="soilks" name="soilks" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+      <field id="soilnvan" name="soilnvan" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+      <field id="soilavan" name="soilavan" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+      <field id="soilmcr" name="soilmcr" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+      <field id="soilmcs" name="soilmcs" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+      <field id="soilmcfc" name="soilmcfc" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+      <field id="soilmcw" name="soilmcw" domain_ref="soil" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1e10 : 0 </field>
+    </file>
+
     
     <file  id="PFTmap_file"  mode="read" output_freq="1y">
