Changeset 6954

Timestamp:

2020-12-03T18:16:43+01:00 (4 years ago)

Author:

agnes.ducharne

Message:

Introducing the changes of branch SP-MIP r6950.
These developments were made by Salma Tafasca durinh her PhD and allow ORCHIDEE to:
(1) read soil hydraulic parameters instead of soil texture maps (with dimension changes) or impose a uniform soil texture all over land, using keywords SPMIPEXP and EXP4;
(2) cancel the effect of roots on Ks with new keyword KFACT_ROOT_TYPE;
(3) introduce a 13th soil texture class for clay oxisols in addition to the 12 USDA classes (with thermal parameters equal to the ones of clay).
These changes to branch 2.2 have been tested in Salma's PhD simulations, and are fully compatible with former ways to control soil texture and hydraulic parameters.

Location:

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba

Files:

• hydrol.f90 (modified) (108 diffs)
• sechiba.f90 (modified) (8 diffs)
• slowproc.f90 (modified) (24 diffs)

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90

@@ -19 +19 @@
 !!                 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 +47 @@
 !! 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 +97 @@
   ! 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 +202 @@
 !$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 +214 @@
   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 +245 @@
 !$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 +436 @@
 !_ ================================================================================================================================
 
-  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 +449 @@
     !! 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 +495 @@
     !! 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 +570 @@
 !_ ================================================================================================================================
 
-  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 +621 @@
     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 +697 @@
 
     !! 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 +787 @@
              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 +794 @@
     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 +821 @@
     !! 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 +912 @@
     !     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 +922 @@
     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 +960 @@
 
     ! 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 +1333 @@
 !_ ================================================================================================================================
 !!_ 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 +1345 @@
 
     !! 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 +1353 @@
     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 +1464 @@
     !! 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 +1479 @@
     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 +1497 @@
     !! 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 +1523 @@
 
 
-    !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 +1647 @@
     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 +1834 @@
     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 +1846 @@
     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 +2353 @@
 
     ! 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 +2707 @@
 !_ 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 +2717 @@
 
     !! 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 +2760 @@
     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 +2788 @@
     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 +2958 @@
     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 +2990 @@
 
        ! 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 +3077 @@
 
     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 +3093 @@
 
     ! 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 +3109 @@
     ENDDO
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO jsl=2,nslm-1
           DO ji=1,kjpindex
@@ -3206 +3118 @@
     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 +3125 @@
     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 +3132 @@
     ! 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 +3173 @@
 
 
-    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 +3182 @@
 
     ! 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 +3202 @@
     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 +3222 @@
     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 +3582 @@
 
   END SUBROUTINE hydrol_flood
-
 
 !! ================================================================================================================================
@@ -3734 +3647 @@
 !_ ================================================================================================================================
 !_ 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 +3661 @@
 
     !! 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 +3729 @@
                                                                                  !! averaged over the mesh (for thermosoil)
                                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex 
-
     !! 0.3 Modified variables
 
@@ -3953 +3874 @@
        ! 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 +4001 @@
        !! 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 +4007 @@
        !! 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 +4033 @@
        !! 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 +4042 @@
        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 +4064 @@
           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 +4195 @@
           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 +4207 @@
        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 +4251 @@
                 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 +4287 @@
           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 +4296 @@
        !      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 +4320 @@
        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 +4372 @@
        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 +4380 @@
           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 +4563 @@
        !! 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 +4583 @@
        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 +4593 @@
 
     !! 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 +4684 @@
        !! 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 +4694 @@
 
        !! 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 +4733 @@
              
              flux_top(ji) = zero
+             r_soil_ns(ji,jst) = zero
              
           ENDIF
@@ -4863 +4786 @@
        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 +4802 @@
           tmat(ji,1,2) = un
           tmat(ji,1,3) = zero
-          rhs(ji,1) = mcr(njsc(ji))
+          rhs(ji,1) = mcr(ji)
        ENDDO
        
@@ -4898 +4821 @@
           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 +4973 @@
 !_ 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 +4980 @@
 
     ! 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 +5041 @@
     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 +5060 @@
           ! 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 +5074 @@
 
           !! 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 +5156 @@
 !_ 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 +5164 @@
     !! 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 +5203 @@
     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 +5212 @@
     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 +5220 @@
     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 +5229 @@
     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 +5239 @@
     ! 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 +5315 @@
 !_ 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 +5321 @@
     !! 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 +5361 @@
     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 +5370 @@
     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 +5378 @@
     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 +5387 @@
     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 +5397 @@
 
     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 +5448 @@
 !_ 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 +5455 @@
 
     !! 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 +5499 @@
     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 +5733 @@
 !_ ================================================================================================================================
 !_ hydrol_soil_coef
- 
-  SUBROUTINE hydrol_soil_coef(kjpindex,ins,njsc)
+
+  SUBROUTINE hydrol_soil_coef(mcr, mcs, kjpindex,ins,njsc)
 
     IMPLICIT NONE
@@ -5804 +5741 @@
 
     !! 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 +5775 @@
              ! 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 +5794 @@
              
              ! 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 +5826 @@
 !_ hydrol_soil_froz
  
-  SUBROUTINE hydrol_soil_froz(kjpindex,ins,njsc)
+  SUBROUTINE hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,ins,njsc)
 
     IMPLICIT NONE
@@ -5894 +5833 @@
 
     !! 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 +5869 @@
           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 +5888 @@
                      (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 +5899 @@
              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 +6341 @@
 !_ 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 +6350 @@
 
     !! 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 +6510 @@
              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 +6582 @@
           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 +6593 @@
     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

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90

@@ -18 +18 @@
 !! 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 +115 @@
   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 +453 @@
                               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 +475 @@
     
     !! 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 +743 @@
     !! 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 +1398 @@
                             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 +1557 @@
     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 +1930 @@
     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)

branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90

@@ -17 +17 @@
 !!
 !! 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 +96 @@
 !!
 !! 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 +226 @@
     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 +277 @@
                                   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 +317 @@
     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 +340 @@
     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 +729 @@
 !>\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 +745 @@
                                 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 +760 @@
     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 +817 @@
     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 +869 @@
   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 +904 @@
     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 +1092 @@
        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 +1322 @@
                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 +1459 @@
        
        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 +2730 @@
 !!
 !! 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 +2742 @@
 !! \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 +2767 @@
     !  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 +2788 @@
     !  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 +2807 @@
     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 +2841 @@
 
     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 +3240 @@
                    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 +3950 @@
     !! 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 +3984 @@
     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