Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes.f90	(revision 8418)
@@ -143,5 +143,24 @@
        !Config Units = BOOLEAN    
        CALL getin_p('SPINUP_ANALYTIC',spinup_analytic)
-
+       !
+       !Config Key   = OK_MOYANO_SOILHUMSAT
+       !Config Desc  = Activation of Moyano equation for control_moist in stomate_litter.
+       !Config If    = OK_STOMATE
+       !Config Def   = n
+       !Config Help  = Activate this option if you want to use Moyano equation
+       !to define control_moist in stomate_litter.
+       !Config Units = [FLAG]    
+       CALL getin_p('OK_MOYANO_SOILHUMSAT',ok_moyano_soilhumsat)
+       !
+       !Config Key   = OK_ORGA
+       !Config Desc  = Activation of Moyano equation for control_moist in
+       !stomate_litter.
+       !Config If    = OK_STOMATE and OK_MOYANO_SOILHUMSAT
+       !Config Def   = n
+       !Config Help  = Activate this option if you want to use the organic part
+       !of the Moyano equation for soil rich in C.
+       !to define control_moist in stomate_litter.
+       !Config Units = [FLAG]    
+       CALL getin_p('OK_ORGA',ok_orga)
     ENDIF
 
@@ -1862,5 +1881,172 @@
     !Config Units = [-]
     CALL getin_p('MOISTCONT_MIN',moistcont_min)
-
+    !
+    !Config Key   = BETA1
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = -0.26
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('BETA1',beta1)
+    !
+    !Config Key   = BETA1_ORGA
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(3)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT
+    !Config Def   = -0.67
+    !Config Help  =
+    !Config Units = [-]
+    CALL getin_p('BETA1_ORGA',beta1_orga)
+    !
+    !Config Key   = BETA2
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.32
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('BETA2',beta2)
+    !
+    !Config Key   = BETA2_ORGA
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(3)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT
+    !Config Def   = 1.08
+    !Config Help  =
+    !Config Units = [-]
+    CALL getin_p('BETA2_ORGA',beta2_orga)    
+    !
+    !Config Key   = BETA3
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = -0.15
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('BETA3',beta3)
+    !
+    !Config Key   = BETA3_ORGA
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(3)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT
+    !Config Def   = -0.57
+    !Config Help  =
+    !Config Units = [-]
+    CALL getin_p('BETA3_ORGA',beta3_orga)    
+    !
+    !Config Key   = BETA4
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.08
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('BETA4',beta4)
+    !
+    !Config Key   = BETA5
+    !Config Desc  = Parameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = -0.09
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('BETA5',beta5)
+    !
+    !Config Key   = BETA6
+    !Config Desc  = Prameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.57
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('BETA6',beta6)
+    !
+    !Config Key   = INTERCEPT
+    !Config Desc  = Prameter for Moyano et al 2012 PRSR model(2)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 1.059
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('INTERCEPT',intercept)
+    !
+    !Config Key   = INTERCEPT_ORGA
+    !Config Desc  = Prameter for Moyano et al 2012 PRSR model(3)
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 1.134
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('INTERCEPT_ORGA',intercept_orga)    
+    !
+    !Config Key   = SRO
+    !Config Desc  = Initial respiration value (SR0) arbitrary defined at 1.0 
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 1.0
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('SRO',SRo)
+    !
+    !Config Key   = MOISTCONSTSAT_MIN
+    !Config Desc  = minimum soil wetness to limit the heterotrophic respiration 
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.0
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('MOISTCONSTSAT_MIN',moistcontSAT_min)
+    !    
+    !Config Key   = SOILHEIGHT
+    !Config Desc  = soilheight to converte gC/m2soil to gC/gsoil 
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.1
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('SOILHEIGHT',soilheight)
+    !    
+    !Config Key   = CINI_MOYANO
+    !Config Desc  = Carbon content to initialize Moyano equation 
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 1000.0
+    !Config Help  = 
+    !Config Units = [gC/m3]
+    CALL getin_p('CINI_MOYANO',Cini_Moyano)
+    !    
+    !Config Key   = LITTERINI_MOYANO
+    !Config Desc  = Carbon litter content to initialize Moyano equation 
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 1000.0
+    !Config Help  = 
+    !Config Units = [gC/m3]
+    CALL getin_p('LITTERINI_MOYANO',Litterini_Moyano)
+    !    
+    !Config Key   = MAX_CARBON_MOYANO
+    !Config Desc  = Maximum carbon concentration in the database used by Moyano et al
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.35
+    !Config Help  = 
+    !Config Units = [gC/g soil]
+    CALL getin_p('MAX_CARBON_MOYANO',max_carbon_moyano)
+    !    
+    !Config Key   = MIN_CARBON_MOYANO
+    !Config Desc  = Minimum carbon concentration in the database used by Moyano et al
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.01
+    !Config Help  = 
+    !Config Units = [gC/g soil]
+    CALL getin_p('MIN_CARBON_MOYANO',min_carbon_moyano)
+    !    
+    !Config Key   = LIMIT_CARBON_ORGA
+    !Config Desc  = Minimum carbon concentration in the database for organic soil used by Moyano et al
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.06
+    !Config Help  = 
+    !Config Units = [gC/g soil]
+    CALL getin_p('LIMIT_CARBON_ORGA',limit_carbon_orga)    
+    !    
+    !Config Key   = MAX_CLAY_MOYANO
+    !Config Desc  = Maximum clay fraction in the database used by Moyano et al
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.58
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('MAX_CLAY_MOYANO',max_clay_moyano)
+    !    
+    !Config Key   = MIN_CLAY_MOYANO
+    !Config Desc  = Minimum clay fraction in the database used by Moyano et al
+    !Config If    = OK_STOMATE, OK_MOYANO_SOILHUMSAT 
+    !Config Def   = 0.03
+    !Config Help  = 
+    !Config Units = [-]
+    CALL getin_p('MIN_CLAY_MOYANO',min_clay_moyano)
     !-
     ! lpj parameters
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_parameters/constantes_var.f90	(revision 8418)
@@ -702,4 +702,6 @@
 !$OMP THREADPRIVATE(veget_year_orig)
   
+ REAL(r_std), SAVE :: bulk_default = 1000.0           !! Default value for bulk density of soil (kg/m3)
+!$OMP THREADPRIVATE(bulk_default)
   REAL(r_std), SAVE :: min_vegfrac = 0.001           !! Minimal fraction of mesh a vegetation type can occupy (0-1, unitless)
 !$OMP THREADPRIVATE(min_vegfrac)
@@ -1022,5 +1024,68 @@
 !$OMP THREADPRIVATE(moistcont_min)
 
-
+  ! 4. Moyano et la. 2012 Biogeosciences
+
+  LOGICAL, SAVE :: ok_moyano_soilhumsat = .FALSE.  !! Use soilhumSAT and Moyano formulation for control_moist
+!$OMP THREADPRIVATE(ok_moyano_soilhumsat)
+  LOGICAL, SAVE :: ok_orga = .TRUE.                !! Do we use also the equation designed for organic soils by Moyano et al.
+                                                   !! for control_moist
+!$OMP THREADPRIVATE(ok_orga)
+  REAL(r_std), SAVE :: beta1 = -0.26               !! 1st(Beta1*M) term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation Beta1= [-0.28,-0.24]
+!$OMP THREADPRIVATE(beta1)
+  REAL(r_std), SAVE :: beta2 = 0.32                !! 2nd(Beta2*M^2) term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation Beta2= [0.28,0.36]
+!$OMP THREADPRIVATE(beta2)
+  REAL(r_std), SAVE :: beta3 = -0.15               !! 3rd(Beta3*M^3) term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation Beta3= [-0.18,-0.12]
+!$OMP THREADPRIVATE(beta3)
+  REAL(r_std), SAVE :: beta1_orga = -0.67          !! 1st(Beta1*M) term of Moyano' PRSR model 3
+                                                   !! for fraction of saturation Beta1= [-0.65,-0.69]
+!$OMP THREADPRIVATE(beta1_orga)
+  REAL(r_std), SAVE :: beta2_orga = 1.08           !! 2nd(Beta2*M^2) term of Moyano' PRSR model 3
+                                                   !! for fraction of saturation Beta2= [1.03,1.13]
+!$OMP THREADPRIVATE(beta2_orga)
+  REAL(r_std), SAVE :: beta3_orga = -0.57          !! 3rd(Beta3*M^3) term of Moyano' PRSR model 3
+                                                   !! for fraction of saturation Beta3= [-0.54,-0.60]
+!$OMP THREADPRIVATE(beta3_orga)
+  REAL(r_std), SAVE :: beta4 = 0.08                !! 4th(Beta4*clay) term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation Beta4= [0.07,0.09]
+!$OMP THREADPRIVATE(beta4)
+  REAL(r_std), SAVE :: beta5 = -0.09               !! 5th(Beta5*M*clay) term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation Beta5= [-0.10,-0.08]
+!$OMP THREADPRIVATE(beta5)
+  REAL(r_std), SAVE :: beta6 = 0.57                !! 6th(Beta6*SOC) term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation Beta6= [0.53,0.61]
+!$OMP THREADPRIVATE(beta6)
+  REAL(r_std), SAVE :: intercept = 1.059           !! intercept term of Moyano' PRSR model 2
+                                                   !! for fraction of saturation intercept= [1.056,1.062]
+!$OMP THREADPRIVATE(intercept)
+  REAL(r_std), SAVE :: intercept_orga = 1.134      !! intercept term of Moyano' PRSR model 3
+                                                   !! for fraction of saturation intercept= [1.131,1.137]
+!$OMP THREADPRIVATE(intercept_orga)
+  REAL(r_std), SAVE :: SRo = 1.0                   !! Initial respiration value (SR0) arbitrary defined at 1.0 
+!$OMP THREADPRIVATE(SRo)
+  REAL(r_std), SAVE :: moistcontSAT_min = 0.0      !! minimum soil wetness to limit the heterotrophic respiration
+!$OMP THREADPRIVATE(moistcontSAT_min)
+  REAL(r_std), SAVE :: soilheight = 0.10           !! soil height in meter to converte totalSOC 
+                                                   !! from gC/m2soil to gC/gsoil [0.10, 0.20, 0.30, 0.50, 1m] 
+!$OMP THREADPRIVATE(soilheight)
+  REAL(r_std), SAVE :: Cini_Moyano = 1000_r_std !!Initial carbon content for Moyano equation  (kgC/m3)
+!$OMP THREADPRIVATE(Cini_Moyano)
+  REAL(r_std), SAVE :: Litterini_Moyano = 1000_r_std !!Initial litter carbon content for Moyano equation  (kgC/m3)
+!$OMP THREADPRIVATE(Litterini_Moyano)
+  REAL(r_std), SAVE :: min_carbon_moyano = 0.01   !! Minimum carbon concentration in the database 
+                                                  !! used by Moyano et al (gC/gSoil)
+!$OMP THREADPRIVATE(min_carbon_moyano)
+  REAL(r_std), SAVE :: max_carbon_moyano = 0.35   !! Maximum carbon concentration in the database 
+                                                  !! used by Moyano et al (gC/gSoil)
+!$OMP THREADPRIVATE(max_carbon_moyano)
+  REAL(r_std), SAVE :: limit_carbon_orga = 0.06   !! Minimum carbon concentration in the database 
+                                                  !! for organic soil used by Moyano et al (gC/gSoil)
+!$OMP THREADPRIVATE(limit_carbon_orga)
+  REAL(r_std), SAVE :: min_clay_moyano = 0.03     !! Minimum clay fraction in the database used by Moyano et al (unitless)
+!$OMP THREADPRIVATE(min_clay_moyano)
+  REAL(r_std), SAVE :: max_clay_moyano = 0.58     !! Maximum clay fraction in the database used by Moyano et al (unitless)
+!$OMP THREADPRIVATE(max_clay_moyano)
   !
   ! stomate_lpj.f90
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 8418)
@@ -415,4 +415,7 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: soil_wet_litter  !! Soil wetness aove mvw in the litter (0-1, unitless)
 !$OMP THREADPRIVATE(soil_wet_litter)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: soil_wet_litterSAT!!Soil wetness above mvw in the litter (0-1, unitless)
+                                                                          !! with respect to(tmc_litter_sat-tmc_litter_wilt)
+!$OMP THREADPRIVATE(soil_wet_litterSAT)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: qflux_ns         !! Diffusive water fluxes between soil layers
                                                                          !! (at lower interface)
@@ -572,5 +575,5 @@
        & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
        & humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, evap_bare_lim_ns, &
-       & flood_frac, flood_res, &
+       & flood_frac, flood_res,shumdiagSAT,litterhumdiagSAT, &
        & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil, rest_id, hist_id, hist2_id,&
        & contfrac, stempdiag, &
@@ -650,7 +653,13 @@
                                                                            !! with respect to (mcfc-mcw)
                                                                            !! (unitless; can be out of 0-1)
+    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out):: shumdiagSAT      !! Relative soil moisture in each soil layer 
+                                                                           !! with respect to (mcs-mcw)
+                                                                           !! (unitless; can be out of 0-1)
+                                                                           !! used in Stomate, for Moyano et al 2012 equation
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out):: shumdiag_perma   !! Percent of porosity filled with water (mc/mcs) used for the thermal computations 
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: k_litt           !! litter approximate conductivity
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: litterhumdiag    !! litter humidity
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: litterhumdiagSAT !!litter humidity 
+                                                                           !! with respect to(tmc_litter_sat-tmc_litter_wilt) 
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: tot_melt         !! Total melt    
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: floodout         !! Flux out of floodplains
@@ -806,6 +815,6 @@
          transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
          returnflow, reinfiltration, irrigation, &
-         tot_melt,evap_bare_lim,evap_bare_lim_ns, shumdiag, shumdiag_perma, &
-         k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,&
+         tot_melt,evap_bare_lim,evap_bare_lim_ns,shumdiagSAT, shumdiag, shumdiag_perma, &
+         k_litt,litterhumdiagSAT, litterhumdiag, humrel, vegstress, drysoil_frac,&
          stempdiag,snow,snowdz, tot_bare_soil,  u, v, tq_cdrag, &
          mc_layh, mcl_layh, root_deficit, veget)
@@ -1772,4 +1781,7 @@
     ALLOCATE (soil_wet_litter(kjpindex,nstm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soil_wet_litter','','')
+
+    ALLOCATE (soil_wet_litterSAT(kjpindex,nstm),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable soil_wet_litterSAT','','') 
 
     ALLOCATE (qflux_ns(kjpindex,nslm,nstm),stat=ier) 
@@ -2391,4 +2403,5 @@
     IF (ALLOCATED  (soil_wet_ns)) DEALLOCATE (soil_wet_ns)
     IF (ALLOCATED  (soil_wet_litter)) DEALLOCATE (soil_wet_litter)
+    IF (ALLOCATED  (soil_wet_litterSAT)) DEALLOCATE (soil_wet_litterSAT)
     IF (ALLOCATED  (qflux_ns)) DEALLOCATE (qflux_ns)
     IF (ALLOCATED  (tmat)) DEALLOCATE (tmat)
@@ -3661,6 +3674,6 @@
        & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &
        & returnflow, reinfiltration, irrigation, &
-       & tot_melt, evap_bare_lim, evap_bare_lim_ns, shumdiag, shumdiag_perma,&
-       & k_litt, litterhumdiag, humrel,vegstress, drysoil_frac, &
+       & tot_melt, evap_bare_lim, evap_bare_lim_ns,shumdiagSAT, shumdiag, shumdiag_perma,&
+       & k_litt,litterhumdiagSAT, litterhumdiag, humrel,vegstress, drysoil_frac, &
        & stempdiag,snow, &
        & snowdz, tot_bare_soil, u, v, tq_cdrag, mc_layh, mcl_layh, root_deficit, veget)
@@ -3724,4 +3737,7 @@
     REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: shumdiag         !! Relative soil moisture in each diag soil layer 
                                                                                  !! with respect to (mcfc-mcw) (unitless, [0-1])
+    REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: shumdiagSAT      !! Relative soil moisture in each diag soil layer 
+                                                                                 !! with respect to(mcs-mcw)(unitless,[0-1]) 
+                                                                                 !! used in stomate for Moyano et al 2012 equation
     REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: shumdiag_perma   !! Percent of porosity filled with water (mc/mcs)
                                                                                  !! in each diag soil layer (for the thermal computations)
@@ -3731,4 +3747,7 @@
     REAL(r_std), DIMENSION (kjpindex), INTENT (out)          :: litterhumdiag    !! Mean of soil_wet_litter across soil tiles
                                                                                  !! (unitless, [0-1])
+    REAL(r_std), DIMENSION (kjpindex), INTENT (out)          ::litterhumdiagSAT  !! Mean of soil_wet_litterSAT across soil tiles
+                                                                                 !!(unitless, [0-1])
+                                                                                 !! with respect to(tmc_litter_sat-tmc_litter_wilt) 
     REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out)      :: vegstress        !! Veg. moisture stress (only for vegetation 
                                                                                  !! growth) (unitless, [0-1])
@@ -4404,4 +4423,13 @@
        END DO
 
+       ! Subsequent calculation of soil_wet_litterSAT (tmc-tmcw)/(tmcs-tmcw)
+       ! Based on liquid water content
+       DO ji=1,kjpindex
+          soil_wet_litterSAT(ji,jst) = MIN(un, MAX(zero,&
+               & (tmc_litter(ji,jst)-tmc_litter_wilt(ji,jst)) / &
+               & (tmc_litter_sat(ji,jst)-tmc_litter_wilt(ji,jst)) ))
+       END DO
+
+
        ! Preliminary calculation of various soil moistures (for each layer, in kg/m2)
        sm(:,1)  = dz(2) * (trois*mcl(:,1,jst) + mcl(:,2,jst))/huit
@@ -4624,5 +4652,5 @@
        
        ! For consistency in stomate, we also set moderwilt and soil_wet_ns to zero in this case. 
-       ! They are used later for shumdiag and shumdiag_perma
+       ! They are used later for shumdiag, shumdiagSAT and shumdiag_perma
        DO jsl = 1,nslm
           WHERE (is_under_mcr(:,jst))
@@ -4671,6 +4699,6 @@
     !! 7. Summing 3d variables into 2d variables 
     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)
+         & evapot, vevapnu, returnflow, reinfiltration, irrigation, shumdiagSAT,&
+         & shumdiag,shumdiag_perma, k_litt,litterhumdiagSAT, litterhumdiag, humrel, vegstress, drysoil_frac,tot_melt)
 
     ! Means of wtd, runoff and drainage corrections, across soiltiles    
@@ -6415,6 +6443,6 @@
 
   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)
+       & evapot, vevapnu, returnflow, reinfiltration, irrigation,shumdiagSAT, &
+       & shumdiag,shumdiag_perma, k_litt,litterhumdiagSAT, litterhumdiag, humrel, vegstress, drysoil_frac, tot_melt)
     ! 
     ! interface description
@@ -6448,7 +6476,10 @@
     REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: drainage        !! Drainage
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)      :: shumdiag        !! relative soil moisture
+    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)      :: shumdiagSAT     !! relative soil moisture with respect to (mcs-mcw)
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (out)      :: shumdiag_perma  !! Percent of porosity filled with water (mc/mcs) used for the thermal computations
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: k_litt          !! litter cond.
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: litterhumdiag   !! litter humidity
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)           :: litterhumdiagSAT!! litter humidity
+                                                                                !! with respect to(tmc_litter_sat-tmc_litter_wilt)
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: humrel          !! Relative humidity
     REAL(r_std), DIMENSION (kjpindex, nvm), INTENT(out)      :: vegstress       !! Veg. moisture stress (only for vegetation growth)
@@ -6496,7 +6527,9 @@
     humtot(:) = zero
     shumdiag(:,:)= zero
+    shumdiagSAT(:,:)= zero
     shumdiag_perma(:,:)=zero
     k_litt(:) = zero
     litterhumdiag(:) = zero
+    litterhumdiagSAT(:) = zero
     tmc_litt_dry_mea(:) = zero
     tmc_litt_wet_mea(:) = zero
@@ -6571,5 +6604,5 @@
     ! BUT THIS IS NOT USED ANYMORE WITH THE NEW BACKGROUNG ALBEDO
     !! k_litt is calculated here as a grid-cell average (for consistency with drainage)   
-    !! litterhumdiag, like shumdiag, is averaged over the soiltiles for transmission to stomate
+    !! litterhumdiag, litterhumdiagSAT, like shumdiag,shumdiagSAT is averaged over the soiltiles for transmission to stomate
     DO jst=1,nstm       
        DO ji=1,kjpindex
@@ -6588,4 +6621,7 @@
           litterhumdiag(ji) = litterhumdiag(ji) + &
                & soil_wet_litter(ji,jst) * soiltile(ji,jst)
+
+          litterhumdiagSAT(ji) = litterhumdiagSAT(ji) + &
+               & soil_wet_litterSAT(ji,jst) * soiltile(ji,jst)
 
           tmc_litt_wet_mea(ji) =  tmc_litt_wet_mea(ji) + & 
@@ -6659,5 +6695,16 @@
        ENDDO
     ENDDO
-    
+
+    ! ShumdiagSAT: we start from soil_wet_ns,
+    !  do a spatial average, excluding the nobio fraction on which stomate doesn't act
+    DO jst=1,nstm
+       DO jsl=1,nslm
+          DO ji=1,kjpindex
+             shumdiagSAT(ji,jsl) = shumdiagSAT(ji,jsl) + soil_wet_ns(ji,jsl,jst) * soiltile(ji,jst) 
+             shumdiagSAT(ji,jsl) = MAX(MIN(shumdiagSAT(ji,jsl), un), zero)
+          ENDDO
+       ENDDO
+    ENDDO
+
     ! Shumdiag_perma is based on soilmoist / moisture at saturation in the layer
     ! Her we start from grid averages by hydrol soil layer and transform it to the diag levels
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 8418)
@@ -140,4 +140,7 @@
                                                                      !! by thermosoil.f90 (unitless, 0-1)
 !$OMP THREADPRIVATE(shumdiag)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiagSAT    !! Mean relative soil moisture (0-1, unitless) 
+                                                                     !! with respect to(mcs-mcw)
+!$OMP THREADPRIVATE(shumdiagSAT) 
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: shumdiag_perma !! Saturation degree of the soil 
 !$OMP THREADPRIVATE(shumdiag_perma)
@@ -146,4 +149,6 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: litterhumdiag  !! Litter dryness factor (unitless, 0-1)
 !$OMP THREADPRIVATE(litterhumdiag)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: litterhumdiagSAT!! Litter dryness factor (unitless, 0-1)
+!$OMP THREADPRIVATE(litterhumdiagSAT)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: stempdiag      !! Temperature which controls canopy evolution (K)
 !$OMP THREADPRIVATE(stempdiag)
@@ -765,5 +770,5 @@
          & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
          & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, evap_bare_lim_ns, flood_frac, flood_res, &
-         & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil,&
+         & shumdiagSAT,litterhumdiagSAT,shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope,&
          & rest_id, hist_id, hist2_id,&
          & contfrac, stempdiag, &
@@ -820,5 +825,5 @@
     CALL slowproc_main (kjit, kjpij, kjpindex, &
          index, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
-         temp_air, temp_sol, stempdiag, &
+         temp_air, temp_sol, stempdiag,shumdiagSAT,litterhumdiagSAT, &
          vegstress, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
          deadleaf_cover, &
@@ -1765,4 +1770,7 @@
     ALLOCATE (shumdiag(kjpindex,nslm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiag','','')
+
+    ALLOCATE (shumdiagSAT(kjpindex,nslm),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for shumdiagSAT','','')
     
     ALLOCATE (shumdiag_perma(kjpindex,nslm),stat=ier)
@@ -1771,4 +1779,7 @@
     ALLOCATE (litterhumdiag(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiag','','')
+
+    ALLOCATE (litterhumdiagSAT(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for litterhumdiagSAT','','')
 
     ALLOCATE (ptnlev1(kjpindex),stat=ier)
@@ -2021,6 +2032,8 @@
     IF ( ALLOCATED (co2_flux)) DEALLOCATE (co2_flux)
     IF ( ALLOCATED (shumdiag)) DEALLOCATE (shumdiag)
+    IF ( ALLOCATED (shumdiagSAT)) DEALLOCATE (shumdiagSAT)
     IF ( ALLOCATED (shumdiag_perma)) DEALLOCATE (shumdiag_perma)
     IF ( ALLOCATED (litterhumdiag)) DEALLOCATE (litterhumdiag)
+    IF ( ALLOCATED (litterhumdiagSAT)) DEALLOCATE (litterhumdiagSAT)
     IF ( ALLOCATED (ptnlev1)) DEALLOCATE (ptnlev1)
     IF ( ALLOCATED (k_litt)) DEALLOCATE (k_litt)
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90	(revision 8418)
@@ -72,4 +72,6 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: siltfraction        !! Siltfraction (0-1, unitless)
 !$OMP THREADPRIVATE(siltfraction)  
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)      :: bulk                !!Bulk density (kg/m**3) 
+!$OMP THREADPRIVATE(bulk)   
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:)  :: laimap              !! LAI map when the LAI is prescribed and not calculated by STOMATE
 !$OMP THREADPRIVATE(laimap)
@@ -137,4 +139,44 @@
     
     !! 2. Prepare for reading of PFTmap file
+         !! 
+         !! 2.1 Prepare for reading of bulk variable 
+         !! 
+     IF (ok_moyano_soilhumsat) THEN 
+         ! Get the file name from run.def file and set file attributes accordingly 
+         filename = 'soil_bulk_and_ph.nc'  
+         CALL getin_p('SOIL_BULK_FILE',filename)  
+         name = filename(1:LEN_TRIM(FILENAME)-3) 
+         CALL xios_orchidee_set_file_attr("soilbulk_file",name=name) 
+      
+         ! Set variables that can be used in the xml files 
+         lerr=xios_orchidee_setvar('bulk_default',bulk_default) 
+          
+         ! 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 soilbulk file will be done later using XIOS. The filename is ', filename 
+        ELSE 
+            ! No reading by XIOS, deactivate soilbulk 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_bulk file will not be done with XIOS.' 
+            CALL xios_orchidee_set_file_attr("soilbulk_file",enabled=.FALSE.) 
+            CALL xios_orchidee_set_field_attr("soilbulk",enabled=.FALSE.) 
+            CALL xios_orchidee_set_field_attr("soilbulk_mask",enabled=.FALSE.) 
+         END IF 
+     ELSE
+         ! Not needed if the flag is not activated, deactivate soilbulk 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_bulk file will not be done with XIOS.'
+         CALL xios_orchidee_set_file_attr("soilbulk_file",enabled=.FALSE.)
+         CALL xios_orchidee_set_field_attr("soilbulk",enabled=.FALSE.)
+         CALL xios_orchidee_set_field_attr("soilbulk_mask",enabled=.FALSE.)
+      END IF
+          
+    !! 
+    !! 2.2 Prepare for reading of PFTmap file 
+    !! 
+
     filename = 'PFTmap.nc'
     CALL getin_p('VEGETATION_FILE',filename)
@@ -377,6 +419,6 @@
              rest_id_stom,   hist_id_stom,           hist_id_stom_IPCC,               &
              indexLand,      lalo,                   neighbours,   resolution,        &
-             contfrac,       totfrac_nobio,          clayfraction, temp_air,          &
-             lai,            veget,                  veget_max,                       &
+             contfrac,       totfrac_nobio,          clayfraction, bulk,              &
+             temp_air,       lai,                    veget,        veget_max,         &
              deadleaf_cover,         assim_param,  temp_growth )
     ENDIF
@@ -478,5 +520,5 @@
   SUBROUTINE slowproc_main (kjit, kjpij, kjpindex, &
        IndexLand, indexveg, lalo, neighbours, resolution, contfrac, soiltile, fraclut, nwdFraclut, &
-       temp_air, temp_sol, stempdiag, &
+       temp_air, temp_sol, stempdiag, shumdiagSAT,litterhumdiagSAT, &
        humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
        deadleaf_cover, &
@@ -511,5 +553,9 @@
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)  :: stempdiag           !! Soil temperature (K)
     REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)  :: shumdiag            !! Relative soil moisture (0-1, unitless)
+    REAL(r_std),DIMENSION (kjpindex,nslm), INTENT (in)  :: shumdiagSAT         !! Relative soil moisture (0-1, unitless)
+                                                                               !! with respect to(mcs-mcw)
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: litterhumdiag       !! Litter humidity  (0-1, unitless)
+    REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: litterhumdiagSAT    !! Litter humidity  (0-1, unitless) 
+                                                                               !!with respect to(tmc_litter_sat-tmc_litter_wilt)
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: precip_rain         !! Rain precipitation (mm dt_stomate^{-1})
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)       :: precip_snow         !! Snow precipitation (mm dt_stomate^{-1})
@@ -627,8 +673,9 @@
 
        !! 4.1 Call stomate main routine that will call all c-cycle routines       !
+
        CALL stomate_main (kjit, kjpij, kjpindex, &
             IndexLand, lalo, neighbours, resolution, contfrac, totfrac_nobio, clayfraction, &
-            temp_air, temp_sol, stempdiag, &
-            humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, gpp, &
+            bulk, temp_air, temp_sol, stempdiag, humrel, shumdiag, shumdiagSAT, litterhumdiag, &
+            litterhumdiagSAT, precip_rain, precip_snow, gpp, &
             deadleaf_cover, &
             assim_param, &
@@ -638,5 +685,4 @@
             co2_flux, fco2_lu, fco2_wh, fco2_ha, &
             resp_maint, resp_hetero, resp_growth, temp_growth)
-
 
        !! 4.2 Output the respiration terms and the net primary
@@ -865,5 +911,5 @@
     ! counter, vegetation fraction, max vegetation fraction, LAI
     ! variable from stomate, fraction of bare soil, soiltype
-    ! fraction, clay fraction, height of vegetation, map of LAI
+    ! fraction, clay fraction,bulk density, height of vegetation, map of LAI
     
     CALL restput_p (rest_id, 'veget', nbp_glo, nvm, 1, kjit, veget, 'scatter',  nbp_glo, index_g)
@@ -904,4 +950,5 @@
     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)
+    CALL restput_p (rest_id, 'bulk', nbp_glo, 1, 1, kjit, bulk, '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)
@@ -929,5 +976,5 @@
     ! 2.2 Write restart variables managed by STOMATE
     IF ( ok_stomate ) THEN
-       CALL stomate_finalize (kjit, kjpindex, indexLand, clayfraction, assim_param) 
+       CALL stomate_finalize (kjit, kjpindex, indexLand, clayfraction, bulk, assim_param) 
     ENDIF
     
@@ -1057,4 +1104,8 @@
     siltfraction(:)=undef_sechiba
 
+    ALLOCATE (bulk(kjpindex),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'slowproc_init','Problem in allocation of variable bulk','','')
+    bulk(:)=undef_sechiba
+
     ! Allocation of last year vegetation fraction in case of land use change
     ALLOCATE(veget_max_new(kjpindex, nvm), STAT=ier)
@@ -1357,5 +1408,5 @@
           CALL slowproc_soilt(njsc, ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, &
                lalo, neighbours, resolution, contfrac, soilclass, &
-               clayfraction, sandfraction, siltfraction)
+               clayfraction, sandfraction, siltfraction, bulk)
           
           call_slowproc_soilt=.FALSE.
@@ -1536,4 +1587,9 @@
     ENDIF
 
+    var_name= 'bulk' 
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Bulk density in each mesh')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., bulk, "gather", nbp_glo, index_g) 
+
     var_name= 'lai'
     CALL ioconf_setatt_p('UNITS', '-')
@@ -1611,4 +1667,13 @@
        !Config Units = [-]
        CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax)
+
+       !Config Key   = BULK  
+       !Config Desc  = Bulk density (0-dim mode) 
+       !Config Def   = XXX 
+       !Config If    = IMPOSE_VEG and IMPOSE_SOIL 
+       !Config Help  = Determines the bulk density in the grid box.  The bulk density 
+       !Config         is the weight of soil in a given volume. 
+       !Config Units = [-]  
+       CALL setvar_p (bulk, val_exp, 'BULK', bulk_default) 
 
        !Config Key   = SLOWPROC_HEIGHT
@@ -1723,5 +1788,4 @@
     ENDIF
     
-
     !! 4.3 Dynamic irrigation map
     !  If do_irrigation, it will look to the dynamical irrig. map in restart
@@ -1843,4 +1907,5 @@
     IF (ALLOCATED (sandfraction)) DEALLOCATE (sandfraction)
     IF (ALLOCATED (siltfraction)) DEALLOCATE (siltfraction)
+    IF (ALLOCATED (bulk)) DEALLOCATE (bulk) 
     IF (ALLOCATED (laimap)) DEALLOCATE (laimap)
     IF (ALLOCATED (veget_max_new)) DEALLOCATE (veget_max_new)
@@ -2531,5 +2596,5 @@
     END IF 
 
-    ! Assigning the right values and giving a value where information was not found
+    ! Assiggning the right values and giving a value where information was not found
     DO ib=1,nbpt
       IF (alaimap(ib) < min_sechiba) THEN
@@ -2931,6 +2996,6 @@
 !>\BRIEF         looks for nearest grid point on the fine map
 !!
-!! DESCRIPTION  : (definitions, functional, design, flags): 
-!!
+!! DESCRIPTION  : (definitions, functional, design, flags):
+!!           
 !! RECENT CHANGE(S): None
 !!
@@ -3003,4 +3068,5 @@
 !!
 !>\BRIEF         Interpolate the Zobler or Reynolds/USDA soil type map
+!!               Read and interpolate soil bulk from file. 
 !!
 !! DESCRIPTION  : (definitions, functional, design, flags): 
@@ -3010,5 +3076,5 @@
 !!                             and everything needed to read all maps and assign parameter values.  
 !!
-!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction
+!! MAIN OUTPUT VARIABLE(S): ::soiltype, ::clayfraction, sandfraction, siltfraction, bulk
 !!
 !! REFERENCE(S) : Reynold, Jackson, and Rawls (2000). Estimating soil water-holding capacities 
@@ -3019,5 +3085,5 @@
 !! \n
 !_ ================================================================================================================================
-  SUBROUTINE slowproc_soilt(njsc,  ks,  nvan, avan, mcr, mcs, mcfc, mcw, 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, bulk)
 
     USE interpweight
@@ -3061,4 +3127,5 @@
     REAL(r_std), INTENT(out)      :: sandfraction(nbpt)     !! The fraction of sand (for SP-MIP)
     REAL(r_std), INTENT(out)      :: siltfraction(nbpt)     !! The fraction of silt (for SP-MIP)
+    REAL(r_std), INTENT(out)      :: bulk(nbpt)             !! Bulk density  as used by STOMATE 
     !
     !
@@ -3088,5 +3155,5 @@
     CHARACTER(LEN=80)                                    :: spmipexp          !! designing the number of sp-mip experiment
     CHARACTER(LEN=80)                                    :: unif_case               !! designing the model of experiment 4 (sp_mip)
-
+    REAL(r_std), DIMENSION(nbpt)                         :: abulkph          !!Availability of the bulk and ph interpolation 
     REAL(r_std)                                          :: vmin, vmax       !! min/max values to use for the 
 
@@ -3113,4 +3180,5 @@
                                                                              !!   `maskingtype') 
     CHARACTER(LEN=250)                                   :: namemaskvar      !! name of the variable to use to mask 
+    CHARACTER(LEN=80)                                    :: fieldname        !!name of the field read for the bulk density map
     INTEGER(i_std), DIMENSION(:), ALLOCATABLE            :: vecpos
     REAL(r_std)                                          :: sgn              !! sum of fractions excluding glaciers and ocean
@@ -3637,5 +3705,5 @@
                       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) *                &
+                         clayfraction(ib) = clayfraction(ib) + textfrac_table(solt(ilf),3) * &
                               textrefrac(ib,solt(ilf))
                          sandfraction(ib) = sandfraction(ib) + textfrac_table(solt(ilf),2) * &
@@ -3688,4 +3756,75 @@
 
        ENDIF        !      xios_interpolation
+
+    !! 
+    !! Read and interpolate soil bulk and soil ph using IOIPSL or XIOS 
+    !!
+    IF (ok_moyano_soilhumsat) THEN
+       IF (xios_interpolation) THEN 
+         ! Read and interpolate using XIOS 
+    
+           ! Check if the restart file for sechiba is read. 
+           ! Reading of soilbulk and soilph with XIOS is only activated if restname==NONE. 
+           IF (restname_in /= 'NONE') THEN 
+              CALL ipslerr_p(3,'slowproc_soilt','soilbulk and soilph can not be read with XIOS if sechiba restart file exist', &
+              'Remove sechiba restart file and start again','') 
+           END IF 
+     
+           IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt:Read soilbulk and soilph with XIOS' 
+           CALL xios_orchidee_recv_field('soilbulk', bulk) 
+       ELSE 
+       ! Read using IOIPSL and interpolate using aggregate tool in ORCHIDEE 
+           IF (printlev_loc>=3) WRITE (numout,*) 'slowproc_soilt:Read soilbulk and soilph with IOIPSL'      
+           !! Read soilbulk 
+     
+           !Config Key   = SOIL_BULK_FILE 
+           !Config Desc  = Name of file from which soil bulk should be read 
+           !Config Def   = soil_bulk_and_ph.nc 
+           !Config If    =  
+           !Config Help  =  
+           !Config Units = [FILE] 
+
+       ! By default, bulk and ph is stored in the same file but they could be
+       ! separated if needed. 
+           filename = 'soil_bulk_and_ph.nc' 
+           CALL getin_p('SOIL_BULK_FILE',filename) 
+     
+           variablename = 'soilbulk' 
+           ! Name of the longitude and latitude in the input file 
+           lonname = 'nav_lon' 
+           latname = 'nav_lat' 
+           vmin=0  ! not used in interpweight_2Dcont 
+           vmax=0  ! not used in interpweight_2Dcont 
+     
+           ! 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 = '' 
+           ! Type of calculation of cell fractions 
+           fractype = 'default' 
+
+              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, bulk_default,undef_sechiba,                       &
+                bulk, abulkph) 
+           WRITE(numout,*) 'bulk density map is read _______'
+
+
+       ENDIF        !      xios_interpolation
+    ENDIF
+
+       DO ib = 1, nbpt
+          njsc(ib) = MAXLOC(soilclass(ib,:),1)
+       ENDDO
 
        ! End of soil texture reading, for 'maps' and classical behavior
@@ -3818,4 +3957,5 @@
     CALL xios_orchidee_send_field("interp_diag_njsc",REAL(njsc, r_std))
     CALL xios_orchidee_send_field("interp_diag_clayfraction",clayfraction)
+    CALL xios_orchidee_send_field("interp_diag_bulk",bulk)
     CALL xios_orchidee_send_field("interp_diag_sandfraction",sandfraction)
     CALL xios_orchidee_send_field("interp_diag_siltfraction",siltfraction)
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate.f90	(revision 8418)
@@ -644,6 +644,6 @@
          rest_id_stom,   hist_id_stom,      hist_id_stom_IPCC,               &
          index,          lalo,              neighbours,   resolution,        &
-         contfrac,       totfrac_nobio,     clay,         temp_air,          &
-         lai,            veget,             veget_max,                       &
+         contfrac,       totfrac_nobio,     clay,         bulk,              &
+         temp_air,       lai,               veget,        veget_max,         &
          deadleaf_cover,    assim_param,  temp_growth )
 
@@ -665,4 +665,5 @@
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: totfrac_nobio     !! Fraction of grid cell covered by lakes, land ice, cities, ... (unitless) 
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
+    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: bulk              !! Bulk density (kg/m**3)  
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: temp_air          !! Air temperature at first atmospheric model layer (K)
     REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: lai               !! Leaf area inex @tex $(m^2 m^{-2})$ @endtex
@@ -1188,6 +1189,6 @@
        & (kjit, kjpij, kjpindex, &
        &  index, lalo, neighbours, resolution, contfrac, totfrac_nobio, clay, &
-       &  temp_air, temp_sol, stempdiag, &
-       &  humrel, shumdiag, litterhumdiag, precip_rain, precip_snow, &
+       &  bulk, temp_air, temp_sol, stempdiag, &
+       &  humrel, shumdiag,shumdiagSAT, litterhumdiag,litterhumdiagSAT, precip_rain, precip_snow, &
        &  gpp, deadleaf_cover, assim_param, &
        &  lai, frac_age, height, veget, veget_max, &
@@ -1225,4 +1226,5 @@
                                                                          !! ice, cities, ... (unitless) 
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
+    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: bulk              !! Bulk density (kg/m**3)  
     REAL(r_std),DIMENSION(kjpindex,nvm),INTENT(in)  :: humrel            !! Relative humidity ("moisture availability") 
                                                                          !! (0-1, unitless) 
@@ -1231,5 +1233,9 @@
     REAL(r_std),DIMENSION(kjpindex,nslm),INTENT(in) :: stempdiag         !! Soil temperature (K)
     REAL(r_std),DIMENSION(kjpindex,nslm),INTENT(in) :: shumdiag          !! Relative soil moisture (0-1, unitless)
+    REAL(r_std),DIMENSION(kjpindex,nslm),INTENT(in) :: shumdiagSAT       !! Relative soil moisture (0-1, unitless)
+                                                                         !! with respect to(mcs-mcw)
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: litterhumdiag     !! Litter humidity (0-1, unitless)
+    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: litterhumdiagSAT  !! Litter humidity (0-1, unitless)
+                                                                         !! with respect to(tmc_litter_sat-tmc_litter_wilt)
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: precip_rain       !! Rain precipitation  
                                                                          !! @tex $(mm dt_stomate^{-1})$ @endtex 
@@ -1369,5 +1375,6 @@
 
 !_ ================================================================================================================================
-    
+   
+ 
   !! 1. Initialize variables
 
@@ -1487,9 +1494,9 @@
     CALL littercalc (kjpindex, &
          turnover_littercalc, bm_to_littercalc, &
-         veget_cov_max, temp_sol, stempdiag, shumdiag, litterhumdiag, &
+         veget_cov_max, temp_sol, stempdiag, shumdiag, litterhumdiag,shumdiagSAT,litterhumdiagSAT, &
          litterpart, litter, dead_leaves, lignin_struc, &
          deadleaf_cover, resp_hetero_litter, &
          soilcarbon_input_inst, control_temp_inst, control_moist_inst, &
-         matrixA, vectorB)
+         matrixA, vectorB, bulk, clay, carbon)
     
     ! Heterothropic litter respiration during time step ::dt_sechiba @tex $(gC m^{-2})$ @endtex
@@ -1551,5 +1558,6 @@
             &          (kjpindex, dt_days, &
             &           veget_cov, veget_cov_max, &
-            &           humrel_daily, t2m_daily, tsoil_daily, soilhum_daily, lalo, &
+            &           humrel_daily, t2m_daily, tsoil_daily, &
+            &           soilhum_daily, lalo, &
             &           precip_daily, npp_daily, biomass, &
             &           turnover_daily, gpp_daily, when_growthinit, &
@@ -1562,5 +1570,5 @@
             &           maxfpc_lastyear, maxfpc_thisyear, &
             &           humrel_month, humrel_week, t2m_longterm, tau_longterm, &
-            &           t2m_month, t2m_week, tsoil_month, soilhum_month, &
+            &           t2m_month, t2m_week, tsoil_month,soilhum_month, &
             &           npp_longterm, turnover_longterm, gpp_week, &
             &           gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
@@ -1584,5 +1592,5 @@
             &             clay, herbivores, &
             &             tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
-            &             litterhum_daily, soilhum_daily, &
+            &             litterhum_daily, soilhum_daily,  &
             &             maxhumrel_lastyear, minhumrel_lastyear, &
             &             gdd0_lastyear, precip_lastyear, &
@@ -2013,5 +2021,5 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE stomate_finalize (kjit, kjpindex, index, clay, assim_param) 
+  SUBROUTINE stomate_finalize (kjit, kjpindex, index, clay, bulk, assim_param) 
     
     IMPLICIT NONE
@@ -2023,4 +2031,5 @@
     INTEGER(i_std),DIMENSION(kjpindex),INTENT(in)   :: index             !! Indices of the terrestrial pixels only (unitless)
     REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: clay              !! Clay fraction of soil (0-1, unitless)
+    REAL(r_std),DIMENSION(kjpindex),INTENT(in)      :: bulk              !! Bulk density (kg/m**3)
     REAL(r_std),DIMENSION(kjpindex,nvm,npco2),INTENT(in) :: assim_param  !! min+max+opt temperatures (K) & vmax for photosynthesis  
 
@@ -2109,5 +2118,5 @@
          dt_days, days_since_beg, &
          ind, adapted, regenerate, &
-         humrel_daily, gdd_init_date, litterhum_daily, &
+         humrel_daily, gdd_init_date, litterhum_daily,&
          t2m_daily, t2m_min_daily, tsurf_daily, tsoil_daily, &
          soilhum_daily, precip_daily, &
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_io.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_io.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_io.f90	(revision 8418)
@@ -440,4 +440,5 @@
          &              .TRUE., litterhum_daily, 'gather', nbp_glo, index_g)
     IF (ALL(litterhum_daily(:) == val_exp)) litterhum_daily(:) = zero
+
     !-
     t2m_daily(:) = val_exp
@@ -471,4 +472,5 @@
          &                .TRUE., soilhum_daily, 'gather', nbp_glo, index_g)
     IF (ALL(soilhum_daily(:,:) == val_exp)) soilhum_daily(:,:) = zero
+
     !-
     precip_daily(:) = val_exp
@@ -604,4 +606,6 @@
          &              .TRUE., soilhum_month, 'gather', nbp_glo, index_g)
     IF (ALL(soilhum_month(:,:) == val_exp)) soilhum_month(:,:) = zero
+
+
     !-
     ! 6 fire probability
@@ -993,5 +997,9 @@
              CALL restget_p (rest_id_stomate, var_name, nbp_glo, nlitt , 1, itime, &
                   &                     .TRUE., litter(:,:,m,l,k), 'gather', nbp_glo, index_g)
-             IF (ALL(litter(:,:,m,l,k) == val_exp)) litter(:,:,m,l,k) = zero
+             IF (ok_moyano_soilhumsat)THEN
+               IF (ALL(litter(:,:,m,l,k) == val_exp)) litter(:,:,m,l,k) = Litterini_Moyano
+             ELSE
+               IF (ALL(litter(:,:,m,l,k) == val_exp)) litter(:,:,m,l,k) = zero
+             ENDIF
           ENDDO
        ENDDO
@@ -1013,5 +1021,9 @@
        CALL restget_p (rest_id_stomate, var_name, nbp_glo, ncarb , 1, itime, &
             &                   .TRUE., carbon(:,:,m), 'gather', nbp_glo, index_g)
-       IF (ALL(carbon(:,:,m) == val_exp)) carbon(:,:,m) = zero
+       IF (ok_moyano_soilhumsat)THEN 
+          IF (ALL(carbon(:,:,m) == val_exp)) carbon(:,:,m) = Cini_Moyano
+       ELSE
+          IF (ALL(carbon(:,:,m) == val_exp)) carbon(:,:,m) = zero
+       ENDIF
     ENDDO
     !-
@@ -1689,17 +1701,17 @@
     var_name = 't2m_daily'
     CALL restput_p (rest_id_stomate, var_name, nbp_glo,    1, 1, itime, &
-         &                t2m_daily, 'scatter', nbp_glo, index_g)
+                    &                t2m_daily, 'scatter', nbp_glo, index_g)
     !-
     var_name = 't2m_min_daily'
     CALL restput_p (rest_id_stomate, var_name, nbp_glo,    1, 1, itime, &
-         &                t2m_min_daily, 'scatter', nbp_glo, index_g)
+                    &                t2m_min_daily, 'scatter', nbp_glo, index_g)
     !-
     var_name = 'tsurf_daily'
     CALL restput_p (rest_id_stomate, var_name, nbp_glo,    1, 1, itime, &
-         &                tsurf_daily, 'scatter', nbp_glo, index_g)
+                    &                tsurf_daily, 'scatter', nbp_glo, index_g)
     !-
     var_name = 'tsoil_daily'
     CALL restput_p (rest_id_stomate, var_name, nbp_glo, nslm, 1, itime, &
-         &                tsoil_daily, 'scatter', nbp_glo, index_g)
+                    &                tsoil_daily, 'scatter', nbp_glo, index_g)
     !-
     var_name = 'soilhum_daily'
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_litter.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_litter.f90	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_stomate/stomate_litter.f90	(revision 8418)
@@ -26,5 +26,5 @@
 
   ! modules used:
-
+  USE xios_orchidee
   USE ioipsl_para
   USE stomate_data
@@ -180,8 +180,9 @@
        turnover, bm_to_litter, &
        veget_cov_max, tsurf, tsoil, soilhum, litterhum, &
+       soilhumSAT, litterhumSAT, &
        litterpart, litter, dead_leaves, lignin_struc, &
        deadleaf_cover, resp_hetero_litter, &
        soilcarbon_input, control_temp, control_moist, &
-       MatrixA, VectorB)
+       MatrixA, VectorB, bulk, clay, carbon)
 
     !! 0. Variable and parameter declaration
@@ -201,5 +202,11 @@
     REAL(r_std), DIMENSION(npts,nslm), INTENT(in)               :: soilhum            !! Daily soil humidity of each soil layer 
                                                                                       !! (unitless)
+    REAL(r_std), DIMENSION(npts,nslm), INTENT(in)               :: soilhumSAT         !! Daily soil humidity of each soil layer 
+                                                                                      !! (unitless)
     REAL(r_std), DIMENSION(npts), INTENT(in)                    :: litterhum          !! Daily litter humidity (unitless)
+    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: litterhumSAT       !! Daily litter humidity (unitless)
+    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: clay               !! Clay fraction (unitless, 0-1) 
+    REAL(r_std), DIMENSION(npts),INTENT(in)                     :: bulk               !! Bulk density (kg/m**3)  
+    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in)          :: carbon             !! Soil carbon pools: active, slow, or passive, \f$(gC m^{2})$\f
 
     !! 0.2 Output variables
@@ -261,4 +268,6 @@
     REAL(r_std), DIMENSION(npts)                                :: soilhum_decomp     !! Humidity used for decompostition in soil
                                                                                       !! (unitless)
+    REAL(r_std), DIMENSION(npts)                                :: soilhumSAT_decomp  !! Humidity used for decompostition in soil
+                                                                                      !! (unitless) for Moyano et al 2012 
     REAL(r_std), DIMENSION(npts)                                :: fd                 !! Fraction of structural or metabolic litter
                                                                                       !! decomposed (unitless)
@@ -585,6 +594,13 @@
     
     !! 4.1 above the ground: litter humidity
-    control_moist(:,iabove) = control_moist_func (npts, litterhum)
-
+
+    IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
+       control_moist(:,iabove) = control_moist_func (npts,clay,bulk,carbon, litter ,litterhumSAT, veget_cov_max)
+    ELSE !ok_moyano_soilhumsat is false by default
+       control_moist(:,iabove) = control_moist_func (npts,clay,bulk,carbon, litter, litterhum, veget_cov_max)
+    ENDIF
+
+        CALL xios_orchidee_send_field("ControlMoistLitter",control_moist(:,iabove))
+        CALL xios_orchidee_send_field("litterhumSAT",litterhumSAT(:))
     !
     !! 4.2 below: convolution of humidity and decomposer profiles
@@ -596,14 +612,26 @@
     !! 4.2.2 integrate over the nslm levels
     soilhum_decomp(:) = zero
-
-    DO l = 1, nslm !Loop over soil levels
-
-       soilhum_decomp(:) = &
-            soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
-            ( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
-
-    ENDDO
-
-    control_moist(:,ibelow) = control_moist_func (npts, soilhum_decomp)
+    soilhumSAT_decomp(:) = zero
+
+    IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
+      DO l = 1, nslm !Loop over soil levels
+         soilhumSAT_decomp(:) = &
+              soilhumSAT_decomp(:) + soilhumSAT(:,l) * rpc(:) * &
+              ( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
+      ENDDO
+      control_moist(:,ibelow) = control_moist_func (npts,clay,bulk,carbon, litter, soilhumSAT_decomp,veget_cov_max)
+
+     ELSE !ok_moyano_soilhumsat is false by default
+      DO l = 1, nslm !Loop over soil levels
+
+         soilhum_decomp(:) = &
+              soilhum_decomp(:) + soilhum(:,l) * rpc(:) * &
+              ( EXP( -z_soil(l-1)/z_decomp ) - EXP( -z_soil(l)/z_decomp ) )
+      ENDDO
+      control_moist(:,ibelow) = control_moist_func (npts,clay,bulk,carbon, litter, soilhum_decomp,veget_cov_max)
+    ENDIF
+
+        CALL xios_orchidee_send_field("ControlMoistSoil",control_moist(:,ibelow))
+        CALL xios_orchidee_send_field("soilhumSAT",soilhumSAT_decomp(:))
 
   !! 5. fluxes from litter to carbon pools and respiration
@@ -878,13 +906,20 @@
 !_ ================================================================================================================================
   
-  FUNCTION control_moist_func (npts, moist_in) RESULT (moistfunc_result)
+  FUNCTION control_moist_func (npts,clay,bulk,carbon,litter,moist_in,veget_cov_max) RESULT (moistfunc_result)
 
   !! 0. Variable and parameter declaration
     
     !! 0.1 Input variables
-          
+    INTEGER(i_std)                           :: k,m,n,i,l,e         !!indices              
     INTEGER(i_std), INTENT(in)               :: npts                !! Domain size - number of grid pixel (unitless)
     REAL(r_std), DIMENSION(npts), INTENT(in) :: moist_in            !! relative humidity (unitless)
-
+    REAL(r_std), DIMENSION(npts), INTENT(in) :: clay                !! Clay fraction (unitless, 0-1) 
+    REAL(r_std), DIMENSION(npts), INTENT(in) :: bulk                !! Bulk density (kg/m**3)  
+    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(in) :: carbon    !! Soil carbon pools: active, s     low, or passive, \f$(gC m^{2})$\f
+    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs,nelements), INTENT(inout) ::litter   !! Metabolic and structural litter,above and
+                                                                    !! below ground. The unit is given by m^2 of ground
+                                                                    !! @tex $(gC m^{-2})$ @endtex
+    REAL(r_std),DIMENSION(npts,nvm),INTENT(in):: veget_cov_max      !! PFT "Maximal" coverage fraction of a PFT 
+ 
     !! 0.2 Output variables
    
@@ -894,9 +929,134 @@
 
     !! 0.4 Local variables
-
-!_ ================================================================================================================================
-
-    moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) + moist_coeff(2)* moist_in(:) - moist_coeff(3)
-    moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ) )
+    REAL(r_std), DIMENSION(101)              :: Mint                !! range [0,0.01,1] to compute PRSRmax
+    REAL(r_std), DIMENSION(npts,101)         :: PRSR                !! Proportional Response of Soil Respiration  
+    REAL(r_std), DIMENSION(npts,101)         :: SR                  !! Soil respiration compute using Moyano et al 2012 procedure  
+    REAL(r_std)                              :: SRmax               !! maximum value of SR
+    REAL(r_std), DIMENSION(npts,101)         :: SRnorm              !! SR normalized by the SRmax 
+    REAL(r_std)                              :: ind                 !! i index of maximum SRnorm value
+    REAL(r_std)                              :: SRmin               !! minimum value of SRnorm
+    REAL(r_std), DIMENSION(npts,101)         :: SRsc                !! Soil respiration rescale between 0 and 1 following Moyano et al 2012 procedure  
+    REAL(r_std)                              :: SRscmax             !! maximum value of SRsc
+    REAL(r_std)                              :: indmc               !! i index of maximum mois_round values 
+    REAL(r_std)                              :: moist_round         !! moist_in round to 2decimal points
+    REAL(r_std), DIMENSION(npts)             :: carbon_gCgs         !! total Soil carbon pools: active +slow+ passive, (gC/gsoil)
+    REAL(r_std), DIMENSION(npts)             :: carbon_gCm2s        !! total Soil carbon pools: active +slow+ passive     , (gC/gsoil)
+    REAL(r_std), DIMENSION(npts)             :: clay_Moyano         !! Clay fraction borned within the values used by Moyano et al., (unitless, 0-1)
+
+!_ ===============================================================================================================================
+
+     IF (ok_moyano_soilhumsat)THEN !!ok_moyano_soilhumsat is true
+
+     !! In ok_moyano_soilhumsat total soil carbon per grid cells need to be computed and
+     !!  converted from gC/m2soil to gC/gsoil
+     carbon_gCm2s(:)=zero
+     carbon_gCgs(:)=zero
+     DO n = 1, npts ! loop over grids
+      DO m = 1, nvm ! Loop over # PFTs
+       DO k = 1, ncarb ! Loop over carbon pools 
+        carbon_gCm2s(n)=carbon_gCm2s(n)+carbon(n,k,m)*veget_cov_max(n,m)
+       ENDDO
+      ENDDO
+     ENDDO
+
+     IF (ok_orga) THEN
+        DO n = 1, npts ! loop over grids
+           carbon_gCgs(n)=carbon_gCm2s(n) / (bulk(n) * 1E3 * soilheight)
+          carbon_gCgs(n)= MAX(min_carbon_moyano,MIN(max_carbon_moyano,carbon_gCgs(n))) 
+          clay_Moyano(n) = MAX(min_clay_moyano,MIN(max_clay_moyano,clay(n)))
+      ENDDO
+
+        !!1. compute Max(Prsr(M[0,0.01,1])); Prsr:Proportional Response of Soil
+        !!Respiration (PRSR) and Soil respiration(SR): SR(M)=SRo x Prsr(M)/Max(Prsr(M[0,0.01,1]))
+        !! M:soilhumSAT, SRo=1 (arbitrary)
+        DO n = 1,npts
+          Mint(:)= zero
+          SRmax = zero
+          DO i = 1,101
+             Mint(i)=0.01*(i-1)
+             IF (carbon_gCgs(n) .LT. limit_carbon_orga) THEN
+                 PRSR(n,i) = beta1 * Mint(i)  + beta2 * Mint(i)**2.0 + beta3 * Mint(i)**3.0 &
+                            & + beta4 * clay_Moyano(n) + beta5 * clay_Moyano(n) * Mint(i) &
+                            & + beta6 * carbon_gCgs(n) + intercept
+             ELSE 
+                 PRSR(n,i) = beta1_orga * Mint(i)  + beta2_orga * Mint(i)**2.0 + &
+                            & beta3_orga * Mint(i)**3.0 + intercept_orga
+             ENDIF
+
+             IF (i.LT.2) THEN
+               SR(n,i) = PRSR(n,i)
+               SRmax = MAX(SR(n,i),SRmax)
+             ELSE
+               SR(n,i) = PRSR(n,i) * SR(n,i-1)
+               SRmax = MAX(SR(n,i),SRmax)
+             ENDIF
+          ENDDO
+          SRnorm(n,:)= SRo * SR(n,:)/SRmax
+        ENDDO
+     ELSE
+        DO n = 1, npts ! loop over grids
+           carbon_gCgs(n)=carbon_gCm2s(n) / (bulk(n) * 1E3 * soilheight)
+          carbon_gCgs(n)= MAX(min_carbon_moyano,MIN(limit_carbon_orga,carbon_gCgs(n)))
+          clay_Moyano(n) = MAX(min_clay_moyano,MIN(max_clay_moyano,clay(n)))
+      ENDDO
+
+        !!1. compute Max(Prsr(M[0,0.01,1])); Prsr:Proportional Response of Soil
+        !!Respiration (PRSR) and Soil respiration(SR): SR(M)=SRo x
+        !Prsr(M)/Max(Prsr(M[0,0.01,1]))
+        !! M:soilhumSAT, SRo=1 (arbitrary)
+        DO n = 1,npts
+          Mint(:)= zero
+          SRmax = zero
+          DO i = 1,101
+             Mint(i)=0.01*(i-1)
+              PRSR(n,i) = beta1 * Mint(i)  + beta2 * Mint(i)**2.0 + beta3 * Mint(i)**3.0 &
+                         & + beta4 * clay_Moyano(n) + beta5 * clay_Moyano(n) * Mint(i) &
+                         & + beta6 * carbon_gCgs(n) + intercept
+
+             IF (i.LT.2) THEN
+               SR(n,i) = PRSR(n,i)
+               SRmax = MAX(SR(n,i),SRmax)
+             ELSE
+               SR(n,i) = PRSR(n,i) * SR(n,i-1)
+               SRmax = MAX(SR(n,i),SRmax)
+             ENDIF
+          ENDDO
+          SRnorm(n,:)= SRo * SR(n,:)/SRmax
+        ENDDO
+     ENDIF   
+
+     !!2. Rescale SR values between 0 and 1 and defined SR for moist_in
+       DO n = 1,npts
+         ind = 1.0
+         indmc = zero
+         moist_round = zero
+         SRscmax = zero
+         SRsc(n,:) = zero
+         SRmin = 1.0
+         
+         ind = MAXLOC(SRnorm(n,:),1)
+         DO i = 1,ind
+           SRmin = MIN(SRnorm(n,i),SRmin)
+         ENDDO
+         DO i = 1,101
+           SRsc(n,i) = SRnorm(n,i)- SRmin
+           IF (i.LE.ind)THEN
+             SRscmax = MAX(SRsc(n,i),SRscmax)
+           ENDIF
+         ENDDO
+         SRsc(n,:)=SRsc(n,:)/SRscmax
+         
+         moist_round=(NINT(moist_in(n)*100.))/100.
+         indmc=INT(moist_round/0.01) + un
+        
+        moistfunc_result(n) = SRsc(n,indmc)
+        moistfunc_result(n) = MAX( moistcontSAT_min, MIN( un, moistfunc_result(n) ))
+       ENDDO
+
+     ELSE !ok_moyano_soilhumsat is false by default
+       moistfunc_result(:) = -moist_coeff(1) * moist_in(:) * moist_in(:) &
+                             & + moist_coeff(2)* moist_in(:) - moist_coeff(3)
+       moistfunc_result(:) = MAX( moistcont_min, MIN( un, moistfunc_result(:) ))
+     ENDIF
 
   END FUNCTION control_moist_func
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/context_input_orchidee.xml
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/context_input_orchidee.xml	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/context_input_orchidee.xml	(revision 8418)
@@ -30,4 +30,13 @@
       <field id="soilcolor" name="soilcolor" domain_ref="dom_file_soil" operation="instant" freq_offset="1ts" enabled="false"/>
       <field id="soiltext" name="soiltext" domain_ref="dom_file_soil" operation="instant" freq_offset="1ts"/>
+    </file>
+
+    <!-- Default file name: soil_bulk_and_ph.nc for both soilbulk and soilph. The variables can be read from different files, specified in run.def
+                  This file is read in the initialization phase only if there is no restart file.
+         For soilbulk and soilph, set values lower than min_sechiba as masked values because missing values are zero in the file.
+        Renormalization will be done later using the interpolated mask.  -->
+    <file id="soilbulk_file" mode="read" output_freq="1y" enabled="true">
+      <field id="soilbulk_mask_read" name="soilbulk" domain_ref="dom_file_soilbulk" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? 1 : 0 </field>
+      <field id="soilbulk_read" name="soilbulk" domain_ref="dom_file_soilbulk" operation="instant" freq_offset="1ts" > this &gt; $min_sechiba ? this : 0 </field>
     </file>
 
@@ -134,4 +143,11 @@
         <field id="soiltext13" > soiltext == 13 </field>
       </field_group>
+
+      <!-- If soilbulk_mask .gt. min_sechiba, then renormalize by dividing with soilbulk_mask else set bulk_default.
+      soilbulk_mask interpolated to the model grid corresponds to the fraction of the grid cell covered by non-masked
+      values in the input file. Where soilbulk_mask=0, no data where found on the source grid. -->
+      <field id="soilbulk_tmp"  field_ref="soilbulk_read" domain_ref="dom_from_soilbulk" />
+      <field id="soilbulk_mask" field_ref="soilbulk_mask_read" domain_ref="dom_from_soilbulk" read_access="true" enabled="true" />
+      <field id="soilbulk"      field_ref="soilbulk_tmp"  domain_ref="dom_from_soilbulk" read_access="true" enabled="true"  > soilbulk_mask &gt; $min_sechiba ? soilbulk_tmp/(soilbulk_mask+1e-100) : $bulk_default </field>
       
       <field id="frac_veget"      field_ref="maxvegetfrac_read" domain_ref="dom_from_pft" read_access="true" enabled="true" />
@@ -220,4 +236,8 @@
     </domain>
 
+    <domain id="dom_from_soilbulk" domain_ref="domain_landpoints">
+    <interpolate_domain order="1" renormalize="false" />
+    </domain>
+
     <domain id="dom_from_pft" domain_ref="domain_landpoints">
       <interpolate_domain order="1" renormalize="false" />
@@ -258,4 +278,8 @@
       <!-- generate_rectilinear_domain bounds_lon_start="-180" bounds_lon_end="180" bounds_lat_start="-90" bounds_lat_end="90"/-->
       <generate_rectilinear_domain bounds_lon_start="-180" bounds_lon_end="180" bounds_lat_start="90" bounds_lat_end="-90"/>
+    </domain>
+
+    <domain id="dom_file_soilbulk" type="rectilinear">
+    <generate_rectilinear_domain bounds_lon_start="-180" bounds_lon_end="180" bounds_lat_start="90" bounds_lat_end="-90"/>/>
     </domain>
     
@@ -346,4 +370,5 @@
     <variable id="albbg_vis_default" type="double">xxx</variable>
     <variable id="albbg_nir_default" type="double">xxx</variable>
+    <variable id="bulk_default" type="double">xxx</variable>
     
   </variable_definition>
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml	(revision 8418)
@@ -358,5 +358,4 @@
     <field id="gpplut" name="gppLut" long_name="Gross primary productivity on Landuse tiles" unit="kgC m-2 s-1" grid_ref="grid_nlut" />
 
-
     <!-- Varaibles interp_diag_xxx and interp_avail_xxx are written just after reading and interpolating from file and are used for validation -->
     <field id="interp_avail_aveget" name="interp_avail_aveget" long_name="Availability of data for the vegetation fraction interpolation (negative values if no data found)" unit="1" operation="once" freq_offset="0ts"/>
@@ -366,4 +365,5 @@
     <field id="interp_diag_soilclass" name="interp_diag_soilclass"     grid_ref="grid_nscm" long_name="soilclass read from soils_param file" operation="once" freq_offset="0ts"/>
     <field id="interp_diag_njsc" name="interp_diag_njsc" long_name="njsc after reading of soilclass from soils_param file" operation="once" freq_offset="0ts"/>
+    <field id="interp_diag_bulk" name="interp_diag_bulk" long_name="Soil bulk density read from file" unit="kg/m**3" operation="once" freq_offset="0ts"/>
     <field id="interp_diag_clayfraction" name="interp_diag_clayfraction" long_name="clayfraction read from soils_param file" operation="once" freq_offset="0ts" />
     <field id="interp_diag_sandfraction" name="interp_diag_sandfraction" long_name="sandfraction read from soils_param file" operation="once" freq_offset="0ts" />
@@ -636,4 +636,10 @@
   <field id="vegetmax_soil" name="vegetsoil" long_name="PFT fraction for soil tiles" unit="1" grid_ref="grid_nvm_nstm"/>
 
+  <!-- Output variables for Moyano et al. 2012 control humidity function -->
+  <field id="ControlMoistLitter" name="ControlMoistLitter" long_name="Controle moist value for iabove" unit="1" grid_ref="grid_landpoints" />
+  <field id="ControlMoistSoil" name="ControlMoistSoil" long_name="Controle moist value for ibelow" unit="1" grid_ref="grid_landpoints"/>
+  <field id="litterhumSAT" name="litterhumSAT" long_name="litter moisture for the controle moisture fonction" unit="1" grid_ref="grid_landpoints"/>
+  <field id="soilhumSAT" name="soilhumSAT" long_name="soil moisture for the control moisture fonction" unit="1" grid_ref="grid_landpoints"/>
+
 </field_definition>
 
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml	(revision 8417)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml	(revision 8418)
@@ -991,4 +991,10 @@
     <field field_ref="nppOther" level="1"/>
     <field field_ref="nep" level="0"/>
+   <!-- variables for Moyano control moisture -->
+    <field field_ref="ControlMoistLitter" level="1"/>
+    <field field_ref="ControlMoistSoil" level="1"/>
+    <field field_ref="litterhumSAT" level="1"/>
+    <field field_ref="soilhumSAT" level="1"/>
+
    </field_group>
     
