Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90	(revision 8423)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/slowproc.f90	(revision 8462)
@@ -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,6 @@
     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)
+    IF (ok_moyano_soilhumsat) 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 +977,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 +1105,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)
@@ -1211,4 +1263,18 @@
     END IF
 
+
+    IF (ok_moyano_soilhumsat) THEN
+       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) 
+       IF ( ALL(bulk(:) .EQ. val_exp ) ) THEN
+          ! bulk is not in restart file
+          call_slowproc_soilt=.TRUE.
+       END IF
+    END IF
+
+    
+    
     ! Calculate siltfraction not needed to be in restart file
     IF ( ALL( sandfraction(:) .EQ. val_exp) ) THEN
@@ -1308,4 +1374,13 @@
        CALL setvar_p (ks, val_exp, 'KS_IMP', ks_default)
 
+       !Config Key   = BULK  
+       !Config Desc  = Bulk density (0-dim mode) 
+       !Config Def   = 1000.0
+       !Config If    = 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) 
+       
        ! By default, we calculate mcf and mcw from the above values, as in slowproc_soilt,
        ! but they can be overruled by values from run.def
@@ -1357,5 +1432,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 +1611,5 @@
     ENDIF
 
+
     var_name= 'lai'
     CALL ioconf_setatt_p('UNITS', '-')
@@ -1611,4 +1687,5 @@
        !Config Units = [-]
        CALL setvar_p (lai, val_exp, 'SECHIBA_LAI', llaimax)
+ 
 
        !Config Key   = SLOWPROC_HEIGHT
@@ -1723,5 +1800,4 @@
     ENDIF
     
-
     !! 4.3 Dynamic irrigation map
     !  If do_irrigation, it will look to the dynamical irrig. map in restart
@@ -1843,4 +1919,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 +2608,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 +3008,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 +3080,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 +3088,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 +3097,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 +3139,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 +3167,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 
 
@@ -3637,5 +3716,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 +3767,71 @@
 
        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
 
        ! End of soil texture reading, for 'maps' and classical behavior
@@ -3818,4 +3964,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)
