Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 6189)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 6954)
@@ -19,6 +19,10 @@
 !!                 This routine was originaly developed by Patricia deRosnay.
 !!
-!! RECENT CHANGE(S) : None
-!!
+!! RECENT CHANGE(S) : November 2020: It is possible to define soil hydraulic parameters from maps,
+!!                    as needed for the SP-MIP project (Tafasca Salma and Ducharne Agnes).
+!!                    Here, it leads to change dimensions and indices. 
+!!                    We can also impose kfact_root=1 in all soil layers to cancel the effect of
+!!                    roots on ks profile (keyword KFACT_ROOT_TYPE).
+!!                 
 !! REFERENCE(S) :
 !! - de Rosnay, P., J. Polcher, M. Bruen, and K. Laval, Impact of a physically based soil
@@ -43,4 +47,6 @@
 !! of a new soil freezing scheme for a land-surface model with physically-based hydrology.
 !! The Cryosphere, 6, 407-430, doi: 10.5194/tc-6-407-2012. \n
+!! - Tafasca S. (2020). Evaluation de lâimpact des propriÃ©tÃ©s du sol sur lâhydrologie simulee dans le
+!! modÃšle ORCHIDEE, PhD thesis, Sorbonne Universite. \n
 !!
 !! SVN          :
@@ -91,28 +97,8 @@
   ! one dimension array allocated, computed, saved and got in hydrol module
   ! Values per soil type
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: nvan                !! Van Genuchten coeficients n (unitless)
-                                                                          ! RK: 1/n=1-m
-!$OMP THREADPRIVATE(nvan)                                                 
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: avan                !! Van Genuchten coeficients a
-                                                                         !!  @tex $(mm^{-1})$ @endtex
-!$OMP THREADPRIVATE(avan)                                                
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcr                 !! Residual volumetric water content
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
-!$OMP THREADPRIVATE(mcr)                                                 
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcs                 !! Saturated volumetric water content
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex
-!$OMP THREADPRIVATE(mcs)                                                  
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: ks                  !! Hydraulic conductivity at saturation
-                                                                         !!  @tex $(mm d^{-1})$ @endtex
-!$OMP THREADPRIVATE(ks)                                                  
+                                                              
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: pcent               !! Fraction of saturated volumetric soil moisture above 
                                                                          !! which transpir is max (0-1, unitless)
-!$OMP THREADPRIVATE(pcent)
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcfc                !! Volumetric water content at field capacity
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex 
-!$OMP THREADPRIVATE(mcfc)                                                 
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mcw                 !! Volumetric water content at wilting point
-                                                                         !!  @tex $(m^{3} m^{-3})$ @endtex 
-!$OMP THREADPRIVATE(mcw)                                                 
+!$OMP THREADPRIVATE(pcent)                                                      
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)   :: mc_awet             !! Vol. wat. cont. above which albedo is cst
                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex 
@@ -216,5 +202,5 @@
 !$OMP THREADPRIVATE(kfact_root)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: kfact            !! Factor to reduce Ks with depth (unitless)
-                                                                         !! DIM = nslm * nscm
+                                                                         !! DIM = nslm * kjpindex
 !$OMP THREADPRIVATE(kfact)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: zz               !! Depth of nodes [znh in vertical_soil] transformed into (mm) 
@@ -228,21 +214,21 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: mc_lin   !! 50 Vol. Wat. Contents to linearize K and D, for each texture 
                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex
-                                                                 !! DIM = imin:imax * nscm
+                                                                 !! DIM = imin:imax * kjpindex
 !$OMP THREADPRIVATE(mc_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: k_lin    !! 50 values of unsaturated K, for each soil layer and texture
                                                                  !!  @tex $(mm d^{-1})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(k_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: d_lin    !! 50 values of diffusivity D, for each soil layer and texture
                                                                  !!  @tex $(mm^2 d^{-1})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(d_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: a_lin    !! 50 values of the slope in K=a*mc+b, for each soil layer and texture
                                                                  !!  @tex $(mm d^{-1})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(a_lin)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: b_lin    !! 50 values of y-intercept in K=a*mc+b, for each soil layer and texture
                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex
-                                                                 !! DIM = imin:imax * nslm * nscm
+                                                                 !! DIM = imin:imax * nslm * kjpindex
 !$OMP THREADPRIVATE(b_lin)
 
@@ -259,8 +245,8 @@
 !$OMP THREADPRIVATE(kk)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: avan_mod_tab  !! VG parameter a modified from  exponantial profile
-                                                                      !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,nscm)
+                                                                      !! @tex $(mm^{-1})$ @endtex !! DIMENSION (nslm,kjpindex)
 !$OMP THREADPRIVATE(avan_mod_tab)  
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: nvan_mod_tab  !! VG parameter n  modified from  exponantial profile
-                                                                      !! (unitless) !! DIMENSION (nslm,nscm)  
+                                                                      !! (unitless) !! DIMENSION (nslm,kjpindex)  
 !$OMP THREADPRIVATE(nvan_mod_tab)
   
@@ -450,5 +436,7 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE hydrol_initialize ( kjit,           kjpindex,  index,         rest_id,          &
+  SUBROUTINE hydrol_initialize ( ks,             nvan,      avan,          mcr,              &
+                                 mcs,            mcfc,      mcw,           kjit,             &
+                                 kjpindex,       index,     rest_id,                         &
                                  njsc,           soiltile,  veget,         veget_max,        &
                                  humrel,         vegstress, drysoil_frac,                    &
@@ -461,4 +449,16 @@
     !! 0. Variable and parameter declaration
     !! 0.1 Input variables
+
+
+    !salma: added soil params as input variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+    !end salma: added soil params as input variables
+
     INTEGER(i_std), INTENT(in)                         :: kjit             !! Time step number 
     INTEGER(i_std), INTENT(in)                         :: kjpindex         !! Domain size
@@ -495,12 +495,25 @@
     !! 0.4 Local variables
     INTEGER(i_std)                                       :: jsl
+    
+    !salma: added soil params which depend on soil texture
+    REAL(r_std), DIMENSION (nscm)                        :: nvan_text
+    REAL(r_std), DIMENSION (nscm)                        :: avan_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcr_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcs_text
+    REAL(r_std), DIMENSION (nscm)                        :: ks_text
+    REAL(r_std), DIMENSION (nscm)                        :: pcent_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcfc_text
+    REAL(r_std), DIMENSION (nscm)                        :: mcw_text
+    REAL(r_std), DIMENSION (nscm)                        :: mc_awet_text
+    REAL(r_std), DIMENSION (nscm)                        :: mc_adry_text
+     !end salma: added soil params which depend on soil texture
 !_ ================================================================================================================================
 
-    CALL hydrol_init (kjit, kjpindex, index, rest_id, veget_max, soiltile, &
+    CALL hydrol_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc, kjit, kjpindex, index, rest_id, veget_max, soiltile, &
          humrel, vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
          snowdz, snowgrain, snowrho,    snowtemp,   snowheat, &
          drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
     
-    CALL hydrol_var_init (kjpindex, veget, veget_max, &
+    CALL hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget, veget_max, &
          soiltile, njsc, mx_eau_var, shumdiag_perma, &
          drysoil_frac, qsintveg, mc_layh, mcl_layh) 
@@ -557,5 +570,6 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE hydrol_main (kjit, kjpindex, &
+  SUBROUTINE hydrol_main (ks, nvan, avan, mcr, mcs, mcfc, mcw,  &
+       & kjit, kjpindex, &
        & index, indexveg, indexsoil, indexlayer, indexnslm, &
        & temp_sol_new, floodout, runoff, drainage, frac_nobio, totfrac_nobio, vevapwet, veget, veget_max, njsc, &
@@ -607,4 +621,15 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: transpir         !! Transpiration
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinf_slope      !! Slope coef
+
+    !salma: added input soil params:
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+    !end salma: added input soil params
+ 
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot           !! Soil Potential Evaporation
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot_penm      !! Soil Potential Evaporation Correction
@@ -672,8 +697,10 @@
 
     !! 0.4 Local variables
-
+    !salma: added variable kfact_root_type
     INTEGER(i_std)                                     :: jst              !! Index of soil tiles (unitless, 1-3)
     INTEGER(i_std)                                     :: jsl              !! Index of soil layers (unitless)
     INTEGER(i_std)                                     :: ji, jv
+    CHARACTER(LEN=80)                                  :: kfact_root_type  !! read from run.def: when equal to 'cons', it indicates that ks does not increase
+                                                                           !!   in the rootzone, ie, kfact_root=1; else, kfact_root defined as usual
     REAL(r_std),DIMENSION (kjpindex)                   :: soilwet          !! A temporary diagnostic of soil wetness
     REAL(r_std),DIMENSION (kjpindex)                   :: snowdepth_diag   !! Depth of snow layer containing default values, only for diagnostics
@@ -760,5 +787,5 @@
              IF(soiltile(ji,jst) .GT. min_sechiba) THEN
                 kfact_root(ji,jsl,jst) = kfact_root(ji,jsl,jst) * &
-                     & MAX((MAXVAL(ks_usda)/ks(njsc(ji)))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), &
+                     & MAX((MAXVAL(ks_usda)/ks(ji))**(- vegetmax_soil(ji,jv,jst)/2 * (humcste(jv)*zz(jsl)/mille - un)/deux), &
                      un) 
              ENDIF
@@ -767,4 +794,20 @@
     ENDDO
 
+    !salma: added config key of KFACT_ROOT_TYPE
+    !! KFACT_ROOT_TYPE = cons is used to impose that kfact_root = 1 in every soil layer,
+    !! so that ks does not increase in the rootzone; else, kfact_root defined as usual
+    !Config Key   = KFACT_ROOT_TYPE
+    !Config Desc  = keyword added for spmip exp1 and exp4 to get a constant ks over soil depth and rootzone
+    !Config If    = spmip exp1 or exp4
+    !Config Def   = var
+    !Config Help  = can have two values: 'cons' or 'var'. If var then no changes, if cons then kfact_root=un
+    !Config Units = [mm/d]
+    CALL getin_p("KFACT_ROOT_TYPE",kfact_root_type)
+
+    IF (kfact_root_type=='cons') THEN
+    kfact_root(:,:,:) = un
+    ENDIF
+    !end salma: added config key
+
     !
     !! 3.3 computes canopy  ==>hydrol_canop
@@ -778,5 +821,5 @@
     !! 3.5 computes soil hydrology ==>hydrol_soil
 
-    CALL hydrol_soil(kjpindex, veget_max, soiltile, njsc, reinf_slope,  &
+    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope,  &
          transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
          returnflow, reinfiltration, irrigation, &
@@ -869,8 +912,9 @@
     !     mcs etc are identical in all layers (no normalization by vegtot to be comparable to mc)
     DO jsl=1,nslm
-       land_mcs(:,jsl) = mcs(njsc(:))
-       land_mcfc(:,jsl) = mcfc(njsc(:))
-       land_mcw(:,jsl) = mcw(njsc(:))
-       land_mcr(:,jsl) = mcr(njsc(:))
+       !salma: replaced mcs(njsc(:)) by mcs(:) and same for other variables
+       land_mcs(:,jsl) = mcs(:)
+       land_mcfc(:,jsl) = mcfc(:)
+       land_mcw(:,jsl) = mcw(:)
+       land_mcr(:,jsl) = mcr(:)
     ENDDO
     CALL xios_orchidee_send_field("mcs",land_mcs) ! in m3/m3
@@ -878,5 +922,6 @@
     CALL xios_orchidee_send_field("mcw",land_mcw) ! in m3/m3 
     CALL xios_orchidee_send_field("mcr",land_mcr) ! in m3/m3 
-          
+
+      
     CALL xios_orchidee_send_field("water2infilt",water2infilt)   
     CALL xios_orchidee_send_field("mc",mc)
@@ -915,6 +960,6 @@
 
     ! For the soil wetness above wilting point for CMIP6 (mrsow)
-    mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(njsc(:)) ) &
-         / ( zmaxh*mille*( mcs(njsc(:)) - mcw(njsc(:)) ) )
+    mrsow(:) = MAX( zero,humtot(:) - zmaxh*mille*mcw(:) ) &
+         / ( zmaxh*mille*( mcs(:) - mcw(:) ) )
     CALL xios_orchidee_send_field("mrsow",mrsow)
 
@@ -1288,6 +1333,8 @@
 !_ ================================================================================================================================
 !!_ hydrol_init
-
-  SUBROUTINE hydrol_init(kjit, kjpindex, index, rest_id, veget_max, soiltile, &
+!salma: added variables and njsc in arguments and input variables
+
+  SUBROUTINE hydrol_init(ks, nvan, avan, mcr, mcs, mcfc, mcw, njsc,&
+       kjit, kjpindex, index, rest_id, veget_max, soiltile, &
        humrel,  vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
        snowdz,  snowgrain, snowrho,    snowtemp,   snowheat, &
@@ -1298,5 +1345,6 @@
 
     !! 0.1 Input variables
-
+    !salma introduced njsc variable
+    INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)    :: njsc               !! Index of the dominant soil textural class in the grid cell (1-nscm, unitless)
     INTEGER(i_std), INTENT (in)                         :: kjit               !! Time step number 
     INTEGER(i_std), INTENT (in)                         :: kjpindex           !! Domain size
@@ -1305,6 +1353,15 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)   :: veget_max          !! Carte de vegetation max
     REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in)  :: soiltile           !! Fraction of each soil tile within vegtot (0-1, unitless)
-
+   
     !! 0.2 Output variables
+
+    !salma: added input variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)         :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
 
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)  :: humrel             !! Stress hydrique, relative humidity
@@ -1407,32 +1464,12 @@
     !! 2.1 array allocation for soil texture
 
-    ALLOCATE (nvan(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan','','')
-
-    ALLOCATE (avan(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan','','')
-
-    ALLOCATE (mcr(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcr','','')
-
-    ALLOCATE (mcs(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcs','','')
-
-    ALLOCATE (ks(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ks','','')
-
-    ALLOCATE (pcent(nscm),stat=ier)
+   
+    ALLOCATE (pcent(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable pcent','','')
-
-    ALLOCATE (mcfc(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcfc','','')
-
-    ALLOCATE (mcw(nscm),stat=ier)
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mcw','','')
-    
-    ALLOCATE (mc_awet(nscm),stat=ier)
+    
+    ALLOCATE (mc_awet(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_awet','','')
 
-    ALLOCATE (mc_adry(nscm),stat=ier)
+    ALLOCATE (mc_adry(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_adry','','')
        
@@ -1442,26 +1479,12 @@
     CASE (3)
               
-       nvan(:) = nvan_fao(:)       
-       avan(:) = avan_fao(:)
-       mcr(:) = mcr_fao(:)
-       mcs(:) = mcs_fao(:)
-       ks(:) = ks_fao(:)
-       pcent(:) = pcent_fao(:)
-       mcfc(:) = mcf_fao(:)
-       mcw(:) = mcw_fao(:)
-       mc_awet(:) = mc_awet_fao(:)
-       mc_adry(:) = mc_adry_fao(:)
-    CASE (12)
-       
-       nvan(:) = nvan_usda(:)
-       avan(:) = avan_usda(:)
-       mcr(:) = mcr_usda(:)
-       mcs(:) = mcs_usda(:)
-       ks(:) = ks_usda(:)
-       pcent(:) = pcent_usda(:)
-       mcfc(:) = mcf_usda(:)
-       mcw(:) = mcw_usda(:)
-       mc_awet(:) = mc_awet_usda(:)
-       mc_adry(:) = mc_adry_usda(:)
+       pcent(:) = pcent_fao(njsc(:))
+       mc_awet(:) = mc_awet_fao(njsc(:))
+       mc_adry(:) = mc_adry_fao(njsc(:))
+    CASE (13)
+           
+       pcent(:) = pcent_usda(njsc(:))
+       mc_awet(:) = mc_awet_usda(njsc(:))
+       mc_adry(:) = mc_adry_usda(njsc(:))
        
     CASE DEFAULT
@@ -1474,97 +1497,10 @@
     !! 2.3 Read in the run.def the parameters values defined by the user
 
-    !Config Key   = CWRR_N_VANGENUCHTEN
-    !Config Desc  = Van genuchten coefficient n
-    !Config If    = 
-    !Config Def   = 1.89, 1.56, 1.31
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [-]
-    CALL getin_p("CWRR_N_VANGENUCHTEN",nvan)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(nvan(:) <= zero) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for CWRR_N_VANGENUCHTEN.", &
-            &     "This parameter should be positive. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = CWRR_A_VANGENUCHTEN
-    !Config Desc  = Van genuchten coefficient a
-    !Config If    = 
-    !Config Def   = 0.0075, 0.0036, 0.0019
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [1/mm]  
-    CALL getin_p("CWRR_A_VANGENUCHTEN",avan)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(avan(:) <= zero) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for CWRR_A_VANGENUCHTEN.", &
-            &     "This parameter should be positive. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = VWC_RESIDUAL
-    !Config Desc  = Residual soil water content
-    !Config If    = 
-    !Config Def   = 0.065, 0.078, 0.095
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [m3/m3]  
-    CALL getin_p("VWC_RESIDUAL",mcr)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcr(:) < zero) .OR. ANY(mcr(:) > 1.)  ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_RESIDUAL.", &
-            &     "This parameter is ranged between 0 and 1. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-    
-    !Config Key   = VWC_SAT
-    !Config Desc  = Saturated soil water content
-    !Config If    = 
-    !Config Def   = 0.41, 0.43, 0.41
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [m3/m3]  
-    CALL getin_p("VWC_SAT",mcs)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcs(:) < zero) .OR. ANY(mcs(:) > 1.) .OR. ANY(mcs(:) <= mcr(:)) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_SAT.", &
-            &     "This parameter should be greater than VWC_RESIDUAL and less than 1. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = CWRR_KS 
-    !Config Desc  = Hydraulic conductivity Saturation
-    !Config If    = 
-    !Config Def   = 1060.8, 249.6, 62.4
-    !Config Help  = This parameter will be constant over the entire 
-    !Config         simulated domain, thus independent from soil
-    !Config         texture.   
-    !Config Units = [mm/d]   
-    CALL getin_p("CWRR_KS",ks)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(ks(:) <= zero) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for CWRR_KS.", &
-            &     "This parameter should be positive. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
+
+
+
+
+
+
 
 
@@ -1587,39 +1523,5 @@
 
 
-    !Config Key   = VWC_FC 
-    !Config Desc  = Volumetric water content field capacity
-    !Config If    = 
-    !Config Def   = 0.32, 0.32, 0.32
-    !Config Help  = This parameter is independent from soil texture for
-    !Config         the time being.
-    !Config Units = [m3/m3]   
-    CALL getin_p("VWC_FC",mcfc)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcfc(:) > mcs(:)) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_FC.", &
-            &     "This parameter should be less than VWC_SAT. ", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
-
-    !Config Key   = VWC_WP
-    !Config Desc  = Volumetric water content Wilting pt
-    !Config If    = 
-    !Config Def   = 0.10, 0.10, 0.10 
-    !Config Help  = This parameter is independent from soil texture for
-    !Config         the time being.
-    !Config Units = [m3/m3]   
-    CALL getin_p("VWC_WP",mcw)
-
-    !! Check parameter value (correct range)
-    IF ( ANY(mcw(:) > mcfc(:)) .OR. ANY(mcw(:) < mcr(:)) ) THEN
-       CALL ipslerr_p(error_level, "hydrol_init.", &
-            &     "Wrong parameter value for VWC_WP.", &
-            &     "This parameter should be greater or equal than VWC_RESIDUAL and less or equal than VWC_SAT.", &
-            &     "Please, check parameter value in run.def. ")
-    END IF
-
+   
 
     !Config Key   = VWC_MIN_FOR_WET_ALB
@@ -1745,8 +1647,8 @@
     kk(:,:,:) = 276.48
     
-    ALLOCATE (avan_mod_tab(nslm,nscm),stat=ier) 
+    ALLOCATE (avan_mod_tab(nslm,kjpindex),stat=ier) 
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable avan_mod_tab','','')
     
-    ALLOCATE (nvan_mod_tab(nslm,nscm),stat=ier) 
+    ALLOCATE (nvan_mod_tab(nslm,kjpindex),stat=ier) 
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nvan_mod_tab','','')
 
@@ -1932,5 +1834,5 @@
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact_root','','')
 
-    ALLOCATE (kfact(nslm, nscm),stat=ier)
+    ALLOCATE (kfact(nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable kfact','','')
 
@@ -1944,17 +1846,17 @@
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable dh','','')
 
-    ALLOCATE (mc_lin(imin:imax, nscm),stat=ier)
+    ALLOCATE (mc_lin(imin:imax, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc_lin','','')
 
-    ALLOCATE (k_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (k_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable k_lin','','')
 
-    ALLOCATE (d_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (d_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable d_lin','','')
 
-    ALLOCATE (a_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (a_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable a_lin','','')
 
-    ALLOCATE (b_lin(imin:imax, nslm, nscm),stat=ier)
+    ALLOCATE (b_lin(imin:imax, nslm, kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable b_lin','','')
 
@@ -2451,12 +2353,6 @@
 
     ! Allocation for soiltile related parameters
-    IF ( ALLOCATED (nvan)) DEALLOCATE (nvan)
-    IF ( ALLOCATED (avan)) DEALLOCATE (avan)
-    IF ( ALLOCATED (mcr)) DEALLOCATE (mcr)
-    IF ( ALLOCATED (mcs)) DEALLOCATE (mcs)
-    IF ( ALLOCATED (ks)) DEALLOCATE (ks)
+   
     IF ( ALLOCATED (pcent)) DEALLOCATE (pcent)
-    IF ( ALLOCATED (mcfc)) DEALLOCATE (mcfc)
-    IF ( ALLOCATED (mcw)) DEALLOCATE (mcw)
     IF ( ALLOCATED (mc_awet)) DEALLOCATE (mc_awet)
     IF ( ALLOCATED (mc_adry)) DEALLOCATE (mc_adry)
@@ -2811,5 +2707,6 @@
 !_ hydrol_var_init
 
-  SUBROUTINE hydrol_var_init (kjpindex, veget, veget_max, soiltile, njsc, &
+  SUBROUTINE hydrol_var_init (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
+       kjpindex, veget, veget_max, soiltile, njsc, &
        mx_eau_var, shumdiag_perma, &
        drysoil_frac, qsintveg, mc_layh, mcl_layh) 
@@ -2820,4 +2717,12 @@
 
     !! 0.1 Input variables
+    !salma: added the following soil parameters
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
 
     ! input scalar 
@@ -2855,5 +2760,5 @@
     REAL(r_std)                                         :: nvan_mod      !! VG parameter n  modified from  exponantial profile
                                                                          !! (unitless)
-    REAL(r_std), DIMENSION(nslm,nscm)                   :: afact, nfact  !! Multiplicative factor for decay of a and n with depth
+    REAL(r_std), DIMENSION(nslm,kjpindex)               :: afact, nfact  !! Multiplicative factor for decay of a and n with depth
                                                                          !! (unitless)
     ! parameters for "soil densification" with depth
@@ -2883,5 +2788,8 @@
     REAL(r_std), DIMENSION (kjpindex,nslm)              :: nvg             !! VG param n modified with depth at each node
                                                                            !! (unitless)
-    INTEGER(i_std)                                      :: jiref           !! To identify the mc_lins where k_lin and d_lin 
+                                                                           !! need special treatment
+    !salma: added local variable ii and jiref replaced with iiref
+    INTEGER(i_std)                                      :: ii
+    INTEGER(i_std)                                      :: iiref           !! To identify the mc_lins where k_lin and d_lin
                                                                            !! need special treatment
 
@@ -3050,23 +2958,25 @@
     END IF
 
+ 
+
     !-
     !! 3 Compute the profile for a and n
     !-
-
-    ! For every soil texture
-    DO jsc = 1, nscm 
+    !salma: for every grid cell
+    DO ji = 1, kjpindex
        DO jsl=1,nslm
           ! PhD thesis of d'Orgeval, 2006, p81, Eq. 4.38; d'Orgeval et al. 2008, Eq. 2
           ! Calibrated against Hapex-Sahel measurements
-          kfact(jsl,jsc) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un)
-          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14 
-          
-          nfact(jsl,jsc) = ( kfact(jsl,jsc) )**nk_rel
-          afact(jsl,jsc) = ( kfact(jsl,jsc) )**ak_rel
-       ENDDO
-    ENDDO
-
-    ! For every soil texture
-    DO jsc = 1, nscm
+          kfact(jsl,ji) = MIN(MAX(EXP(- f_ks * (zz(jsl)/mille - dp_comp)), un/kfact_max),un)
+          ! PhD thesis of d'Orgeval, 2006, p81, Eqs. 4.39; 4.42, and Fig 4.14
+
+          nfact(jsl,ji) = ( kfact(jsl,ji) )**nk_rel
+          afact(jsl,ji) = ( kfact(jsl,ji) )**ak_rel
+       ENDDO
+    ENDDO
+    !end salma
+    
+    ! For every grid cell
+     DO ji = 1, kjpindex
        !-
        !! 4 Compute the linearized values of k, a, b and d
@@ -3080,75 +2990,77 @@
 
        ! We define 51 bounds for 50 bins of mc between mcr and mcs
-       mc_lin(imin,jsc)=mcr(jsc)
-       mc_lin(imax,jsc)=mcs(jsc)
-       DO ji= imin+1, imax-1 ! ji=2,50
-          mc_lin(ji,jsc) = mcr(jsc) + (ji-imin)*(mcs(jsc)-mcr(jsc))/(imax-imin)
+       mc_lin(imin,ji)=mcr(ji)
+       mc_lin(imax,ji)=mcs(ji)
+       DO ii= imin+1, imax-1 ! ii=2,50
+          mc_lin(ii,ji) = mcr(ji) + (ii-imin)*(mcs(ji)-mcr(ji))/(imax-imin)
        ENDDO
 
        DO jsl = 1, nslm
           ! From PhD thesis of d'Orgeval, 2006, p81, Eq. 4.42
-          nvan_mod = n0 + (nvan(jsc)-n0) * nfact(jsl,jsc)
-          avan_mod = a0 + (avan(jsc)-a0) * afact(jsl,jsc)
+          nvan_mod = n0 + (nvan(ji)-n0) * nfact(jsl,ji)
+          avan_mod = a0 + (avan(ji)-a0) * afact(jsl,ji)
           m = un - un / nvan_mod
           ! Creation of arrays for SP-MIP output by landpoint
-          nvan_mod_tab(jsl,jsc) = nvan_mod
-          avan_mod_tab(jsl,jsc) = avan_mod
-          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(jsc) * kfact(jsl,jsc)
-          DO ji = imax,imin,-1 
-             frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
-             k_lin(ji,jsl,jsc) = ks(jsc) * kfact(jsl,jsc) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
+          nvan_mod_tab(jsl,ji) = nvan_mod
+          avan_mod_tab(jsl,ji) = avan_mod
+          ! We apply Van Genuchten equation for K(theta) based on Ks(z)=ks(ji) * kfact(jsl,ji)
+          DO ii = imax,imin,-1
+             frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
+             k_lin(ii,jsl,ji) = ks(ji) * kfact(jsl,ji) * (frac**0.5) * ( un - ( un - frac ** (un/m)) ** m )**2
           ENDDO
 
           ! k_lin should not be zero, nor too small
-          ! We track jiref, the bin under which mc is too small and we may get zero k_lin     
-          ji=imax-1
-          DO WHILE ((k_lin(ji,jsl,jsc) > 1.e-32) .and. (ji>0))
-             jiref=ji
-             ji=ji-1
+          ! We track iiref, the bin under which mc is too small and we may get zero k_lin
+          !salma: ji replaced with ii and jiref replaced with iiref and jsc with ji
+          ii=imax-1
+          DO WHILE ((k_lin(ii,jsl,ji) > 1.e-32) .and. (ii>0))
+             iiref=ii
+             ii=ii-1
           ENDDO
-          DO ji=jiref-1,imin,-1
-             k_lin(ji,jsl,jsc)=k_lin(ji+1,jsl,jsc)/10.
+          DO ii=iiref-1,imin,-1
+             k_lin(ii,jsl,ji)=k_lin(ii+1,jsl,ji)/10.
           ENDDO
-         
-          DO ji = imin,imax-1 ! ji=1,50
+
+          DO ii = imin,imax-1 ! ii=1,50
              ! We deduce a_lin and b_lin based on continuity between segments k_lin = a_lin*mc-lin+b_lin
-             a_lin(ji,jsl,jsc) = (k_lin(ji+1,jsl,jsc)-k_lin(ji,jsl,jsc)) / (mc_lin(ji+1,jsc)-mc_lin(ji,jsc))
-             b_lin(ji,jsl,jsc)  = k_lin(ji,jsl,jsc) - a_lin(ji,jsl,jsc)*mc_lin(ji,jsc)
+             a_lin(ii,jsl,ji) = (k_lin(ii+1,jsl,ji)-k_lin(ii,jsl,ji)) / (mc_lin(ii+1,ji)-mc_lin(ii,ji))
+             b_lin(ii,jsl,ji)  = k_lin(ii,jsl,ji) - a_lin(ii,jsl,ji)*mc_lin(ii,ji)
 
              ! We calculate the d_lin for each mc bin, from Van Genuchten equation for D(theta)
-             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin 
-             IF (ji.NE.imin .AND. ji.NE.imax-1) THEN
-                frac=MIN(un,(mc_lin(ji,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
-                d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) *  &
-                     ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) * &
+             ! d_lin is constant and taken as the arithmetic mean between the values at the bounds of each bin
+             IF (ii.NE.imin .AND. ii.NE.imax-1) THEN
+                frac=MIN(un,(mc_lin(ii,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
+                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) *  &
+                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) * &
                      (  frac**(-un/m) -un ) ** (-m)
-                frac=MIN(un,(mc_lin(ji+1,jsc)-mcr(jsc))/(mcs(jsc)-mcr(jsc)))
-                d_lin(ji+1,jsl,jsc) =(k_lin(ji+1,jsl,jsc) / (avan_mod*m*nvan_mod))*&
-                     ( (frac**(-un/m))/(mc_lin(ji+1,jsc)-mcr(jsc)) ) * &
+                frac=MIN(un,(mc_lin(ii+1,ji)-mcr(ji))/(mcs(ji)-mcr(ji)))
+                d_lin(ii+1,jsl,ji) =(k_lin(ii+1,jsl,ji) / (avan_mod*m*nvan_mod))*&
+                     ( (frac**(-un/m))/(mc_lin(ii+1,ji)-mcr(ji)) ) * &
                      (  frac**(-un/m) -un ) ** (-m)
-                d_lin(ji,jsl,jsc) = undemi * (d_lin(ji,jsl,jsc)+d_lin(ji+1,jsl,jsc))
-             ELSE IF(ji.EQ.imax-1) THEN
-                d_lin(ji,jsl,jsc) =(k_lin(ji,jsl,jsc) / (avan_mod*m*nvan_mod)) * &
-                     ( (frac**(-un/m))/(mc_lin(ji,jsc)-mcr(jsc)) ) *  &
+                d_lin(ii,jsl,ji) = undemi * (d_lin(ii,jsl,ji)+d_lin(ii+1,jsl,ji))
+             ELSE IF(ii.EQ.imax-1) THEN
+                d_lin(ii,jsl,ji) =(k_lin(ii,jsl,ji) / (avan_mod*m*nvan_mod)) * &
+                     ( (frac**(-un/m))/(mc_lin(ii,ji)-mcr(ji)) ) *  &
                      (  frac**(-un/m) -un ) ** (-m)
              ENDIF
-          ENDDO
-
-          ! Special case for ji=imin
-          d_lin(imin,jsl,jsc) = d_lin(imin+1,jsl,jsc)/1000.
+          ENDDO !Salma end loop over landpoints
+
+          ! Special case for ii=imin
+          d_lin(imin,jsl,ji) = d_lin(imin+1,jsl,ji)/1000.
 
           ! We adjust d_lin where k_lin was previously adjusted otherwise we might get non-monotonous variations
           ! We don't want d_lin = zero
-          DO ji=jiref-1,imin,-1
-             d_lin(ji,jsl,jsc)=d_lin(ji+1,jsl,jsc)/10.
+          DO ii=iiref-1,imin,-1
+             d_lin(ii,jsl,ji)=d_lin(ii+1,jsl,ji)/10.
           ENDDO
 
        ENDDO
     ENDDO
+
 
     ! Output of alphavg and nvg at each node for SP-MIP
     DO jsl = 1, nslm
-       alphavg(:,jsl) = avan_mod_tab(jsl,njsc(:))*1000. ! from mm-1 to m-1
-       nvg(:,jsl) = nvan_mod_tab(jsl,njsc(:))
+       alphavg(:,jsl) = avan_mod_tab(jsl,:)*1000. ! from mm-1 to m-1
+       nvg(:,jsl) = nvan_mod_tab(jsl,:)
     ENDDO
     CALL xios_orchidee_send_field("alphavg",alphavg) ! in m-1
@@ -3165,7 +3077,7 @@
 
     mx_eau_var(:) = zero
-    mx_eau_var(:) = zmaxh*mille*mcs(njsc(:)) 
-
-    DO ji = 1,kjpindex 
+    mx_eau_var(:) = zmaxh*mille*mcs(:)
+
+    DO ji = 1,kjpindex
        IF (vegtot(ji) .LE. zero) THEN
           mx_eau_var(ji) = mx_eau_nobio*zmaxh
@@ -3181,13 +3093,13 @@
 
     ! Loop on soiltiles to compute the variables (ji,jst)
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji = 1, kjpindex
-          tmcs(ji,jst)=zmaxh* mille*mcs(njsc(ji))
-          tmcr(ji,jst)=zmaxh* mille*mcr(njsc(ji))
-          tmcfc(ji,jst)=zmaxh* mille*mcfc(njsc(ji))
-          tmcw(ji,jst)=zmaxh* mille*mcw(njsc(ji))
-       ENDDO
-    ENDDO
-       
+          tmcs(ji,jst)=zmaxh* mille*mcs(ji)
+          tmcr(ji,jst)=zmaxh* mille*mcr(ji)
+          tmcfc(ji,jst)=zmaxh* mille*mcfc(ji)
+          tmcw(ji,jst)=zmaxh* mille*mcw(ji)
+       ENDDO
+    ENDDO
+
     ! The total soil moisture for each soiltile:
     DO jst=1,nstm
@@ -3197,5 +3109,5 @@
     ENDDO
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO jsl=2,nslm-1
           DO ji=1,kjpindex
@@ -3206,5 +3118,5 @@
     ENDDO
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji=1,kjpindex
           tmc(ji,jst) = tmc(ji,jst) +  dz(nslm) * (trois * mc(ji,nslm,jst) + mc(ji,nslm-1,jst))/huit
@@ -3213,5 +3125,5 @@
     END DO
 
-!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.    
+!JG: hydrol_tmc_update should not be called in the initialization phase. Call of hydrol_tmc_update makes the model restart differenlty.
 !    ! If veget has been updated before restart (with LAND USE or DGVM),
 !    ! tmc and mc must be modified with respect to humtot conservation.
@@ -3220,40 +3132,40 @@
     ! The litter variables:
     ! level 1
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji=1,kjpindex
           tmc_litter(ji,jst) = dz(2) * (trois*mcl(ji,1,jst)+mcl(ji,2,jst))/huit
-          tmc_litter_wilt(ji,jst) = dz(2) * mcw(njsc(ji)) / deux
-          tmc_litter_res(ji,jst) = dz(2) * mcr(njsc(ji)) / deux
-          tmc_litter_field(ji,jst) = dz(2) * mcfc(njsc(ji)) / deux
-          tmc_litter_sat(ji,jst) = dz(2) * mcs(njsc(ji)) / deux
-          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(njsc(ji)) / deux
-          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(njsc(ji)) / deux
+          tmc_litter_wilt(ji,jst) = dz(2) * mcw(ji) / deux
+          tmc_litter_res(ji,jst) = dz(2) * mcr(ji) / deux
+          tmc_litter_field(ji,jst) = dz(2) * mcfc(ji) / deux
+          tmc_litter_sat(ji,jst) = dz(2) * mcs(ji) / deux
+          tmc_litter_awet(ji,jst) = dz(2) * mc_awet(ji) / deux
+          tmc_litter_adry(ji,jst) = dz(2) * mc_adry(ji) / deux
        ENDDO
     END DO
     ! sum from level 2 to 4
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO jsl=2,4
           DO ji=1,kjpindex
-             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * & 
+             tmc_litter(ji,jst) = tmc_litter(ji,jst) + dz(jsl) * &
                   & ( trois*mcl(ji,jsl,jst) + mcl(ji,jsl-1,jst))/huit &
                   & + dz(jsl+1)*(trois*mcl(ji,jsl,jst) + mcl(ji,jsl+1,jst))/huit
              tmc_litter_wilt(ji,jst) = tmc_litter_wilt(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))*& 
-                  & mcw(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))*&
+                  & mcw(ji)/deux
              tmc_litter_res(ji,jst) = tmc_litter_res(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))*& 
-                  & mcr(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))*&
+                  & mcr(ji)/deux
              tmc_litter_sat(ji,jst) = tmc_litter_sat(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))* & 
-                  & mcs(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))* &
+                  & mcs(ji)/deux
              tmc_litter_field(ji,jst) = tmc_litter_field(ji,jst) + &
-                  & (dz(jsl)+ dz(jsl+1))* & 
-                  & mcfc(njsc(ji))/deux
+                  & (dz(jsl)+ dz(jsl+1))* &
+                  & mcfc(ji)/deux
              tmc_litter_awet(ji,jst) = tmc_litter_awet(ji,jst) + &
-                  &(dz(jsl)+ dz(jsl+1))* & 
-                  & mc_awet(njsc(ji))/deux
+                  &(dz(jsl)+ dz(jsl+1))* &
+                  & mc_awet(ji)/deux
              tmc_litter_adry(ji,jst) = tmc_litter_adry(ji,jst) + &
-                  & (dz(jsl)+ dz(jsl+1))* & 
-                  & mc_adry(njsc(ji))/deux
+                  & (dz(jsl)+ dz(jsl+1))* &
+                  & mc_adry(ji)/deux
           END DO
        END DO
@@ -3261,8 +3173,8 @@
 
 
-    DO jst=1,nstm 
+    DO jst=1,nstm
        DO ji=1,kjpindex
           ! here we set that humrelv=0 in PFT1
-          humrelv(ji,1,jst) = zero
+         humrelv(ji,1,jst) = zero
        ENDDO
     END DO
@@ -3270,5 +3182,5 @@
 
     ! Calculate shumdiag_perma for thermosoil
-    ! Use resdist instead of soiltile because we here need to have 
+    ! Use resdist instead of soiltile because we here need to have
     ! shumdiag_perma at the value from previous time step.
     ! Here, soilmoist is only used as a temporary variable to calculate shumdiag_perma
@@ -3290,16 +3202,16 @@
     ENDDO
     DO ji=1,kjpindex
-        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average 
-    ENDDO
-    
+        soilmoist(ji,:) = soilmoist(ji,:) * vegtot_old(ji) ! grid cell average
+    ENDDO
+
     ! -- shumdiag_perma for restart
-    !  For consistency with hydrol_soil, we want to calculate a grid-cell average
+   !  For consistency with hydrol_soil, we want to calculate a grid-cell average
     DO jsl = 1, nslm
-       DO ji=1,kjpindex        
-          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji)))
-          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 
-       ENDDO
-    ENDDO
-               
+       DO ji=1,kjpindex
+          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji))
+          shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero)
+       ENDDO
+    ENDDO
+
     ! Calculate drysoil_frac if it was not found in the restart file
     ! For simplicity, we set drysoil_frac to 0.5 in this case
@@ -3310,9 +3222,11 @@
     END IF
 
-    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in 
-    !! thermosoil for the thermal conductivity. 
+    !! Calculate the volumetric soil moisture content (mc_layh and mcl_layh) needed in
+    !! thermosoil for the thermal conductivity.
     ! These values are only used in thermosoil_init in absence of a restart file
+
     mc_layh(:,:) = zero
     mcl_layh(:,:) = zero
+      
     DO jst=1,nstm
        DO jsl=1,nslm
@@ -3668,5 +3582,4 @@
 
   END SUBROUTINE hydrol_flood
-
 
 !! ================================================================================================================================
@@ -3734,6 +3647,6 @@
 !_ ================================================================================================================================
 !_ hydrol_soil
-
-  SUBROUTINE hydrol_soil (kjpindex, veget_max, soiltile, njsc, reinf_slope, &
+  SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
+       kjpindex, veget_max, soiltile, njsc, reinf_slope, &
        & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &
        & returnflow, reinfiltration, irrigation, &
@@ -3748,4 +3661,13 @@
 
     !! 0.1 Input variables
+    !salma added soil params to arguments and input variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+
 
     INTEGER(i_std), INTENT(in)                               :: kjpindex 
@@ -3807,5 +3729,4 @@
                                                                                  !! averaged over the mesh (for thermosoil)
                                                                                  !!  @tex $(m^{3} m^{-3})$ @endtex 
-
     !! 0.3 Modified variables
 
@@ -3953,5 +3874,5 @@
        ! These values will be kept till the end of the prognostic loop
        DO jst=1,nstm
-          CALL hydrol_soil_froz(kjpindex,jst,njsc)
+          CALL hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,jst,njsc)
        ENDDO
 
@@ -4080,5 +4001,5 @@
        !! K and D are computed for the profile of mc before infiltration
        !! They depend on the fraction of soil ice, given by profil_froz_hydro_ns
-       CALL hydrol_soil_coef(kjpindex,jst,njsc)
+       CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc)
 
        !! Infiltration and surface runoff are computed
@@ -4086,5 +4007,5 @@
        !! The conductivity comes from hydrol_soil_coef and relates to the liquid phase only
        !  This seems consistent with ok_freeze
-       CALL hydrol_soil_infilt(kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, &
+       CALL hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, jst, njsc, flux_infilt, qinfilt_ns, ru_infilt_ns, &
             check_infilt_ns)
        ru_ns(:,jst) = ru_infilt_ns(:,jst) 
@@ -4112,5 +4033,5 @@
        !! 2.1 K and D are recomputed after infiltration
        !! They depend on the fraction of soil ice, still given by profil_froz_hydro_ns 
-       CALL hydrol_soil_coef(kjpindex,jst,njsc)
+       CALL hydrol_soil_coef(mcr, mcs,kjpindex,jst,njsc)
  
        !! 2.2 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme
@@ -4121,6 +4042,6 @@
        DO jsl = 1, nslm
           DO ji =1, kjpindex
-             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + &
-                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + &
+                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) )
              ! we always have mcl<=mc
              ! if mc>mcr, then mcl>mcr; if mc=mcr,mcl=mcr; if mc<mcr, then mcl<mcr
@@ -4143,8 +4064,8 @@
           DO ji = 1, kjpindex
              IF (soiltile(ji,jst) .GT. zero) THEN
-                mvg = un - un / nvan_mod_tab(jsl,njsc(ji))
-                avg = avan_mod_tab(jsl,njsc(ji))*1000. ! to convert in m-1
-                mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(njsc(ji)))/(mcs(njsc(ji)) - mcr(njsc(ji))) )**(-un/mvg)
-                psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,njsc(ji))) / avg ! in m 
+                mvg = un - un / nvan_mod_tab(jsl,ji)
+                avg = avan_mod_tab(jsl,ji)*1000. ! to convert in m-1
+                mc_ratio = MAX( 10.**(-14*mvg), (mcl(ji,jsl,jst) - mcr(ji))/(mcs(ji) - mcr(ji)) )**(-un/mvg)
+                psi = - MAX(zero,(mc_ratio - un))**(un/nvan_mod_tab(jsl,ji)) / avg ! in m 
                 psi_moy(ji,jsl) = psi_moy(ji,jsl) + soiltile(ji,jst) * psi ! average across soil tiles
              ENDIF
@@ -4274,5 +4195,5 @@
           DO ji =1, kjpindex
              mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + &
-                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl
           ENDDO
@@ -4286,5 +4207,5 @@
        dr_corr_ns(:,jst) = zero
        ru_corr_ns(:,jst) = zero
-       call hydrol_soil_smooth_over_mcs2(kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns)
+       call hydrol_soil_smooth_over_mcs2(mcs, kjpindex, jst, njsc, is_over_mcs, ru_corr_ns, check_over_ns)
        
        ! In absence of freezing, if F is large enough, the correction of oversaturation is sent to drainage       
@@ -4330,6 +4251,6 @@
                 dmc(ji,jsl) = zero
                 IF ( ( zz(jsl) >= zwt_force(ji,jst)*mille ) ) THEN
-                   dmc(ji,jsl) = mcs(njsc(ji)) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value
-                   mc(ji,jsl,jst) = mcs(njsc(ji))
+                   dmc(ji,jsl) = mcs(ji) - mc(ji,jsl,jst) ! addition to reach mcs (m3/m3) = positive value
+                   mc(ji,jsl,jst) = mcs(ji)
                 ENDIF
              ENDDO
@@ -4366,5 +4287,5 @@
           wtd_ns(ji,jst) = undef_sechiba ! in meters
           jsl=nslm
-          DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(njsc(ji))) .AND. (jsl > 1) )
+          DO WHILE ( (mc(ji,jsl,jst) .EQ. mcs(ji)) .AND. (jsl > 1) )
              wtd_ns(ji,jst) = zz(jsl)/mille ! in meters
              jsl=jsl-1
@@ -4375,5 +4296,5 @@
        !      This routine does not change tmc but decides where we should turn off ET to prevent further mc decrease
        !      Like above, the tests are made on total mc, compared to mcr
-       CALL hydrol_soil_smooth_under_mcr(kjpindex, jst, njsc, is_under_mcr, check_under_ns)
+       CALL hydrol_soil_smooth_under_mcr(mcr, kjpindex, jst, njsc, is_under_mcr, check_under_ns)
  
        !! 4. At the end of the prognostic calculations, we recompute important moisture variables
@@ -4399,6 +4320,6 @@
        DO jsl = 1, nslm
           DO ji =1, kjpindex
-             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + &
-                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + &
+                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc
           ENDDO
@@ -4451,7 +4372,7 @@
        sm(:,1)  = dz(2) * (trois*mcl(:,1,jst) + mcl(:,2,jst))/huit
        smt(:,1) = dz(2) * (trois*mc(:,1,jst) + mc(:,2,jst))/huit
-       smw(:,1) = dz(2) * (quatre*mcw(njsc(:)))/huit
-       smf(:,1) = dz(2) * (quatre*mcfc(njsc(:)))/huit
-       sms(:,1) = dz(2) * (quatre*mcs(njsc(:)))/huit
+       smw(:,1) = dz(2) * (quatre*mcw(:))/huit
+       smf(:,1) = dz(2) * (quatre*mcfc(:))/huit
+       sms(:,1) = dz(2) * (quatre*mcs(:))/huit
        DO jsl = 2,nslm-1
           sm(:,jsl)  = dz(jsl) * (trois*mcl(:,jsl,jst)+mcl(:,jsl-1,jst))/huit &
@@ -4459,19 +4380,19 @@
           smt(:,jsl) = dz(jsl) * (trois*mc(:,jsl,jst)+mc(:,jsl-1,jst))/huit &
                + dz(jsl+1) * (trois*mc(:,jsl,jst)+mc(:,jsl+1,jst))/huit
-          smw(:,jsl) = dz(jsl) * ( quatre*mcw(njsc(:)) )/huit &
-               + dz(jsl+1) * ( quatre*mcw(njsc(:)) )/huit
-          smf(:,jsl) = dz(jsl) * ( quatre*mcfc(njsc(:)) )/huit &
-               + dz(jsl+1) * ( quatre*mcfc(njsc(:)) )/huit
-          sms(:,jsl) = dz(jsl) * ( quatre*mcs(njsc(:)) )/huit &
-               + dz(jsl+1) * ( quatre*mcs(njsc(:)) )/huit
+          smw(:,jsl) = dz(jsl) * ( quatre*mcw(:) )/huit &
+               + dz(jsl+1) * ( quatre*mcw(:) )/huit
+          smf(:,jsl) = dz(jsl) * ( quatre*mcfc(:) )/huit &
+               + dz(jsl+1) * ( quatre*mcfc(:) )/huit
+          sms(:,jsl) = dz(jsl) * ( quatre*mcs(:) )/huit &
+               + dz(jsl+1) * ( quatre*mcs(:) )/huit
        ENDDO
        sm(:,nslm)  = dz(nslm) * (trois*mcl(:,nslm,jst) + mcl(:,nslm-1,jst))/huit  
        smt(:,nslm) = dz(nslm) * (trois*mc(:,nslm,jst) + mc(:,nslm-1,jst))/huit      
-       smw(:,nslm) = dz(nslm) * (quatre*mcw(njsc(:)))/huit
-       smf(:,nslm) = dz(nslm) * (quatre*mcfc(njsc(:)))/huit
-       sms(:,nslm) = dz(nslm) * (quatre*mcs(njsc(:)))/huit
+       smw(:,nslm) = dz(nslm) * (quatre*mcw(:))/huit
+       smf(:,nslm) = dz(nslm) * (quatre*mcfc(:))/huit
+       sms(:,nslm) = dz(nslm) * (quatre*mcs(:))/huit
        ! sm_nostress = soil moisture of each layer at which us reaches 1, here at the middle of [smw,smf]
        DO jsl = 1,nslm
-          sm_nostress(:,jsl) = smw(:,jsl) + pcent(njsc(:)) * (smf(:,jsl)-smw(:,jsl))
+          sm_nostress(:,jsl) = smw(:,jsl) + pcent(:) * (smf(:,jsl)-smw(:,jsl))
        END DO
 
@@ -4642,4 +4563,5 @@
        !! The multiplication by vegtot creates grid-cell average values
        ! *** To be checked for consistency with the use of nobio properties in thermosoil
+           
        DO jsl=1,nslm
           DO ji=1,kjpindex
@@ -4661,5 +4583,5 @@
        DO jsl = 1, nslm
           ksat(:,jsl) = ksat(:,jsl) + soiltile(:,jst) * &
-               ( ks(njsc(:)) * kfact(jsl,njsc(:)) * kfact_root(:,jsl,jst) ) 
+               ( ks(:) * kfact(jsl,:) * kfact_root(:,jsl,jst) ) 
        ENDDO
               
@@ -4671,5 +4593,5 @@
 
     !! 7. Summing 3d variables into 2d variables 
-    CALL hydrol_diag_soil (kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
+    CALL hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
          & evapot, vevapnu, returnflow, reinfiltration, irrigation, &
          & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,tot_melt)
@@ -4762,9 +4684,9 @@
        !! 8.2 Since we estimate bare soile evap for the next time step, we update profil_froz_hydro and mcl
        !     (effect of mc only, the change in temp_hydro is neglected)
-       IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(kjpindex,jst,njsc)
+       IF ( ok_freeze_cwrr ) CALL hydrol_soil_froz(nvan, avan, mcr, mcs,kjpindex,jst,njsc)
         DO jsl = 1, nslm
           DO ji =1, kjpindex
-             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(njsc(ji)) + &
-                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+             mcl(ji,jsl,jst)= MIN( mc(ji,jsl,jst), mcr(ji) + &
+                  (un-profil_froz_hydro_ns(ji,jsl,jst))*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we keep mcl=mc
           ENDDO
@@ -4772,5 +4694,5 @@
 
        !! 8.3 K and D are recomputed for the updated profile of mc/mcl
-       CALL hydrol_soil_coef(kjpindex,jst,njsc)
+       CALL hydrol_soil_coef(mcr, mcs, kjpindex,jst,njsc)
 
        !! 8.4 Set the tridiagonal matrix coefficients for the diffusion/redistribution scheme
@@ -4811,4 +4733,5 @@
              
              flux_top(ji) = zero
+             r_soil_ns(ji,jst) = zero
              
           ENDIF
@@ -4863,5 +4786,5 @@
        DO ji = 1, kjpindex
           ! by construction, mc and mcl are always on the same side of mcr, so we can use mcl here
-          resolv(ji) = (mcl(ji,1,jst).LT.(mcr(njsc(ji))).AND.flux_top(ji).GT.min_sechiba)
+          resolv(ji) = (mcl(ji,1,jst).LT.(mcr(ji)).AND.flux_top(ji).GT.min_sechiba)
        ENDDO
        !! Reset the coefficient for diffusion (tridiag is only solved if resolv(ji) = .TRUE.)O
@@ -4879,5 +4802,5 @@
           tmat(ji,1,2) = un
           tmat(ji,1,3) = zero
-          rhs(ji,1) = mcr(njsc(ji))
+          rhs(ji,1) = mcr(ji)
        ENDDO
        
@@ -4898,5 +4821,5 @@
           DO ji =1, kjpindex
              mc(ji,jsl,jst) = MAX( mcl(ji,jsl,jst), mcl(ji,jsl,jst) + &
-                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(njsc(ji))) )
+                  profil_froz_hydro_ns(ji,jsl,jst)*(mc(ji,jsl,jst)-mcr(ji)) )
              ! if profil_froz_hydro_ns=0 (including NOT ok_freeze_cwrr) we get mc=mcl
           ENDDO
@@ -5050,5 +4973,5 @@
 !_ hydrol_soil_infilt
 
-  SUBROUTINE hydrol_soil_infilt(kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check)
+  SUBROUTINE hydrol_soil_infilt(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, ins, njsc, flux_infilt, qinfilt_ns, ru_infilt, check)
 
     !! 0. Variable and parameter declaration
@@ -5057,4 +4980,14 @@
 
     ! GLOBAL (in or inout)
+    ! salma added input soil hydraulic parameters
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
+
+
     INTEGER(i_std), INTENT(in)                        :: kjpindex        !! Domain size
     INTEGER(i_std), INTENT(in)                        :: ins
@@ -5108,5 +5041,5 @@
     DO ji = 1, kjpindex
        !! First we fill up the first layer (about 1mm) without any resistance and quasi-immediately
-       wat_inf_pot(ji) = MAX((mcs(njsc(ji))-mc(ji,1,ins)) * dz(2) / deux, zero)
+       wat_inf_pot(ji) = MAX((mcs(ji)-mc(ji,1,ins)) * dz(2) / deux, zero)
        wat_inf(ji) = MIN(wat_inf_pot(ji), flux_infilt(ji))
        mc(ji,1,ins) = mc(ji,1,ins) + wat_inf(ji) * deux / dz(2)
@@ -5127,5 +5060,5 @@
           ! This is computed by an simple arithmetic average because 
           ! the time step (30min) is not appropriate for a geometric average (advised by Haverkamp and Vauclin)
-          k_m = (k(ji,jsl) + ks(njsc(ji))*kfact(jsl-1,njsc(ji))*kfact_root(ji,jsl,ins)) / deux 
+          k_m = (k(ji,jsl) + ks(ji)*kfact(jsl-1,ji)*kfact_root(ji,jsl,ins)) / deux 
 
           IF (ok_freeze_cwrr) THEN
@@ -5141,5 +5074,5 @@
 
           !! From which we deduce the time it takes to fill up the layer or to end the time step...
-          wat_inf_pot(ji) =  MAX((mcs(njsc(ji))-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero)
+          wat_inf_pot(ji) =  MAX((mcs(ji)-mc(ji,jsl,ins)) * (dz(jsl) + dz(jsl+1)) / deux, zero)
           IF ( infilt_tmp(ji) > min_sechiba) THEN
              dt_inf(ji) =  MIN(wat_inf_pot(ji)/infilt_tmp(ji), dt_tmp(ji))
@@ -5223,5 +5156,5 @@
 !_ hydrol_soil_smooth_under_mcr
 
-  SUBROUTINE hydrol_soil_smooth_under_mcr(kjpindex, ins, njsc, is_under_mcr, check)
+  SUBROUTINE hydrol_soil_smooth_under_mcr(mcr, kjpindex, ins, njsc, is_under_mcr, check)
 
     !- arguments
@@ -5231,4 +5164,6 @@
     !! 0.1 Input variables
 
+    !salma: added parameter mcr
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: mcr          !! Residual volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                         :: kjpindex     !! Domain size
     INTEGER(i_std), INTENT(in)                         :: ins          !! Soiltile index (1-nstm, unitless)
@@ -5268,5 +5203,5 @@
     DO jsl = 1,nslm-2
        DO ji=1, kjpindex
-          excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+          excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
           mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * &
@@ -5277,5 +5212,5 @@
     jsl = nslm-1
     DO ji=1, kjpindex
-       excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+       excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
        mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) - excess * &
@@ -5285,5 +5220,5 @@
     jsl = nslm
     DO ji=1, kjpindex
-       excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+       excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
        mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * &
@@ -5294,5 +5229,5 @@
     DO jsl = nslm-1,2,-1
        DO ji=1, kjpindex
-          excess = MAX(mcr(njsc(ji))-mc(ji,jsl,ins),zero)
+          excess = MAX(mcr(ji)-mc(ji,jsl,ins),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) + excess
           mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) - excess * &
@@ -5304,5 +5239,5 @@
     ! excess > 0
     DO ji=1, kjpindex
-       excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(njsc(ji))-mc(ji,1,ins),zero)
+       excessji(ji) = mask_soiltile(ji,ins) * MAX(mcr(ji)-mc(ji,1,ins),zero)
     ENDDO
     DO ji=1, kjpindex
@@ -5380,5 +5315,5 @@
 !_ hydrol_soil_smooth_over_mcs
 
-  SUBROUTINE hydrol_soil_smooth_over_mcs(kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
+  SUBROUTINE hydrol_soil_smooth_over_mcs(mcs ,kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
 
     !- arguments
@@ -5386,6 +5321,7 @@
     !! 0. Variable and parameter declaration
 
-    !! 0.1 Input variables
-
+    !! 0.1 Input variables 
+    !salma: added mcs
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: mcs             !! Saturated volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                           :: kjpindex        !! Domain size
     INTEGER(i_std), INTENT(in)                           :: ins             !! Soiltile index (1-nstm, unitless)
@@ -5425,5 +5361,5 @@
     DO jsl = 1, nslm-2
        DO ji=1, kjpindex
-          excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+          excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
           mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * &
@@ -5434,5 +5370,5 @@
     jsl = nslm-1
     DO ji=1, kjpindex
-       excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+       excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
        mc(ji,jsl+1,ins) = mc(ji,jsl+1,ins) + excess * &
@@ -5442,5 +5378,5 @@
     jsl = nslm
     DO ji=1, kjpindex
-       excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+       excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
        mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
        mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * &
@@ -5451,5 +5387,5 @@
     DO jsl = nslm-1,2,-1
        DO ji=1, kjpindex
-          excess = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero)
+          excess = MAX(mc(ji,jsl,ins)-mcs(ji),zero)
           mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess
           mc(ji,jsl-1,ins) = mc(ji,jsl-1,ins) + excess * &
@@ -5461,5 +5397,5 @@
 
     DO ji=1, kjpindex
-       excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(njsc(ji)),zero)
+       excess = mask_soiltile(ji,ins) * MAX(mc(ji,1,ins)-mcs(ji),zero)
        mc(ji,1,ins) = mc(ji,1,ins) - excess ! then mc(1)=mcs
        rudr_corr(ji,ins) = rudr_corr(ji,ins) + excess * dz(2) / deux 
@@ -5512,5 +5448,5 @@
 !_ hydrol_soil_smooth_over_mcs2
 
-  SUBROUTINE hydrol_soil_smooth_over_mcs2(kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
+  SUBROUTINE hydrol_soil_smooth_over_mcs2(mcs, kjpindex, ins, njsc, is_over_mcs, rudr_corr, check)
 
     !- arguments
@@ -5519,5 +5455,6 @@
 
     !! 0.1 Input variables
-
+    !salma: added mcs
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)        :: mcs             !! Saturated volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                           :: kjpindex        !! Domain size
     INTEGER(i_std), INTENT(in)                           :: ins             !! Soiltile index (1-nstm, unitless)
@@ -5562,5 +5499,5 @@
     DO jsl = 1, nslm
        DO ji=1, kjpindex
-          excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(njsc(ji)),zero) ! >=0
+          excess(ji,jsl) = MAX(mc(ji,jsl,ins)-mcs(ji),zero) ! >=0
           mc(ji,jsl,ins) = mc(ji,jsl,ins) - excess(ji,jsl) ! here mc either does not change or decreases
        ENDDO
@@ -5796,6 +5733,6 @@
 !_ ================================================================================================================================
 !_ hydrol_soil_coef
- 
-  SUBROUTINE hydrol_soil_coef(kjpindex,ins,njsc)
+
+  SUBROUTINE hydrol_soil_coef(mcr, mcs, kjpindex,ins,njsc)
 
     IMPLICIT NONE
@@ -5804,5 +5741,7 @@
 
     !! 0.1 Input variables
-
+    !salma: added mcr and mcs
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
     INTEGER(i_std), INTENT(in)                        :: kjpindex         !! Domain size
     INTEGER(i_std), INTENT(in)                        :: ins              !! Index of soil type
@@ -5836,14 +5775,14 @@
              ! It corresponds to a liquid mc, but the expression is different from mcl in hydrol_soil,
              ! to ensure that we get the a, b, d of the first bin when mcl<mcr
-             mc_used = mcr(njsc(ji))+x*MAX((mc(ji,jsl, ins)-mcr(njsc(ji))),zero) 
+             mc_used = mcr(ji)+x*MAX((mc(ji,jsl, ins)-mcr(ji)),zero) 
              !
              ! calcul de k based on mc_liq
              !
-             i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(njsc(ji)))/(mcs(njsc(ji))-mcr(njsc(ji))))))
-             a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d
-             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), &
-                  a_lin(i,jsl,njsc(ji)) * mc_used + b_lin(i,jsl,njsc(ji))) ! in mm/d
+             i= MAX(imin, MIN(imax-1, INT(imin +(imax-imin)*(mc_used-mcr(ji))/(mcs(ji)-mcr(ji)))))
+             a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d
+             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), &
+                  a_lin(i,jsl,ji) * mc_used + b_lin(i,jsl,ji)) ! in mm/d
           ENDDO ! loop on grid
        ENDDO
@@ -5855,12 +5794,12 @@
              
              ! it is impossible to consider a mc<mcr for the binning
-             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(njsc(ji)), zero)/(mcs(njsc(ji))-mcr(njsc(ji)))
+             mc_ratio = MAX(mc(ji,jsl,ins)-mcr(ji), zero)/(mcs(ji)-mcr(ji))
              
              i= MAX(MIN(INT((imax-imin)*mc_ratio)+imin , imax-1), imin)
-             a(ji,jsl) = a_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             b(ji,jsl) = b_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm/d
-             d(ji,jsl) = d_lin(i,jsl,njsc(ji)) * kfact_root(ji,jsl,ins) ! in mm^2/d
-             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,njsc(ji)), &
-                  a_lin(i,jsl,njsc(ji)) * mc(ji,jsl,ins) + b_lin(i,jsl,njsc(ji)))  ! in mm/d
+             a(ji,jsl) = a_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             b(ji,jsl) = b_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm/d
+             d(ji,jsl) = d_lin(i,jsl,ji) * kfact_root(ji,jsl,ins) ! in mm^2/d
+             k(ji,jsl) = kfact_root(ji,jsl,ins) * MAX(k_lin(imin+1,jsl,ji), &
+                  a_lin(i,jsl,ji) * mc(ji,jsl,ins) + b_lin(i,jsl,ji))  ! in mm/d
           END DO 
        END DO
@@ -5887,5 +5826,5 @@
 !_ hydrol_soil_froz
  
-  SUBROUTINE hydrol_soil_froz(kjpindex,ins,njsc)
+  SUBROUTINE hydrol_soil_froz(nvan, avan, mcr, mcs, kjpindex,ins,njsc)
 
     IMPLICIT NONE
@@ -5894,4 +5833,9 @@
 
     !! 0.1 Input variables
+    !salma: added the following soil variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)     :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
 
     INTEGER(i_std), INTENT(in)                        :: kjpindex         !! Domain size
@@ -5925,7 +5869,7 @@
           DO ji=1,kjpindex
              ! Van Genuchten parameter for thermodynamical calculation
-             m = 1. -1./nvan(njsc(ji))
+             m = 1. -1./nvan(ji)
            
-             IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(njsc(ji))+min_sechiba))) THEN
+             IF ((.NOT. ok_thermodynamical_freezing).OR.(mc(ji,jsl, ins).LT.(mcr(ji)+min_sechiba))) THEN
                 ! Linear soil freezing or soil moisture below residual
                 IF (temp_hydro(ji, jsl).GE.(ZeroCelsius+fr_dT/2.)) THEN
@@ -5944,8 +5888,8 @@
                      (temp_hydro(ji,jsl) .LT. (ZeroCelsius+fr_dT/2.)) ) THEN
                    ! Factor 2.2 from the PhD of Isabelle Gouttevin
-                   x=MIN(((mcs(njsc(ji))-mcr(njsc(ji))) &
-                        *((2.2*1000.*avan(njsc(ji))*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) &
-                        *lhf/ZeroCelsius/10.)**nvan(njsc(ji))+1.)**(-m)) / &
-                        (mc(ji,jsl, ins)-mcr(njsc(ji))),1._r_std)                
+                   x=MIN(((mcs(ji)-mcr(ji)) &
+                        *((2.2*1000.*avan(ji)*(ZeroCelsius+fr_dT/2.-temp_hydro(ji, jsl)) &
+                        *lhf/ZeroCelsius/10.)**nvan(ji)+1.)**(-m)) / &
+                        (mc(ji,jsl, ins)-mcr(ji)),1._r_std)                
                 ELSE
                    x=0._r_std 
@@ -5955,5 +5899,5 @@
              profil_froz_hydro_ns(ji,jsl,ins) = 1._r_std-x
              
-             mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(njsc(ji))
+             mc_ns(ji,jsl)=mc(ji,jsl,ins)/mcs(ji)
 
           ENDDO ! loop on grid
@@ -6397,5 +6341,5 @@
 !_ hydrol_diag_soil
 
-  SUBROUTINE hydrol_diag_soil (kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
+  SUBROUTINE hydrol_diag_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, runoff, drainage, &
        & evapot, vevapnu, returnflow, reinfiltration, irrigation, &
        & shumdiag,shumdiag_perma, k_litt, litterhumdiag, humrel, vegstress, drysoil_frac, tot_melt)
@@ -6406,4 +6350,12 @@
 
     !! 0.1 Input variables
+    !salma: added the following soil variables
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ks               !! Hydraulic conductivity at saturation (mm {-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: nvan             !! Van Genuchten coeficients n (unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: avan             !! Van Genuchten coeficients a (mm-1})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcr              !! Residual volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcs              !! Saturated volumetric water content (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcfc             !! Volumetric water content at field capacity (m^{3} m^{-3})
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: mcw              !! Volumetric water content at wilting point (m^{3} m^{-3})
 
     ! input scalar 
@@ -6558,5 +6510,5 @@
              i= MAX(MIN(INT((imax-imin)*tmc_litter_ratio)+imin, imax-1), imin)
           ENDIF       
-          k_tmp = MAX(k_lin(i,1,njsc(ji))*ks(njsc(ji)), zero)
+          k_tmp = MAX(k_lin(i,1,ji)*ks(ji), zero)
           k_litt(ji) = k_litt(ji) + vegtot(ji)*soiltile(ji,jst) * SQRT(k_tmp) ! grid-cell average
        ENDDO      
@@ -6630,5 +6582,5 @@
           DO ji=1,kjpindex
              shumdiag(ji,jsl) = shumdiag(ji,jsl) + soil_wet_ns(ji,jsl,jst) * soiltile(ji,jst) * &
-                               ((mcs(njsc(ji))-mcw(njsc(ji)))/(mcfc(njsc(ji))-mcw(njsc(ji))))
+                               ((mcs(ji)-mcw(ji))/(mcfc(ji)-mcw(ji)))
              shumdiag(ji,jsl) = MAX(MIN(shumdiag(ji,jsl), un), zero) 
           ENDDO
@@ -6641,5 +6593,5 @@
     DO jsl=1,nslm              
        DO ji=1,kjpindex
-          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(njsc(ji)))
+          shumdiag_perma(ji,jsl) = soilmoist(ji,jsl) / (dh(jsl)*mcs(ji))
           shumdiag_perma(ji,jsl) = MAX(MIN(shumdiag_perma(ji,jsl), un), zero) 
        ENDDO
