Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/diffuco.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/diffuco.f90	(revision 5804)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/diffuco.f90	(revision 5805)
@@ -234,5 +234,5 @@
   SUBROUTINE diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
      & zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, q_cdrag, qsurf, qair, pb, &
-     & evap_bare_lim, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
+     & evap_bare_lim, evap_bare_lim_ns, evapot, evapot_corr, snow, flood_frac, flood_res, frac_nobio, snow_nobio, totfrac_nobio, &
      & swnet, swdown, coszang, ccanopy, humrel, veget, veget_max, lai, qsintveg, qsintmax, assim_param, &
      & vbeta , vbeta1, vbeta2, vbeta3, vbeta3pot, vbeta4, vbeta5, gsmean, rveget, rstruct, cimean, gpp, co2_to_bm, &
@@ -267,5 +267,7 @@
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: flood_res        !! Reservoir in floodplains (estimation to avoid over-evaporation)
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: pb               !! Surface level pressure (hPa)
-    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evap_bare_lim    !! Limit to the bare soil evaporation when the 
+    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)      :: evap_bare_lim    !! Limit to the bare soil evaporation when the 
+                                                                           !! 11-layer hydrology is used (-)
+    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! Limit to the bare soil evaporation when the 
                                                                            !! 11-layer hydrology is used (-)
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evapot           !! Soil Potential Evaporation (mm day^{-1}) 
@@ -389,5 +391,5 @@
     !
 
-    CALL diffuco_bare (kjpindex, evap_bare_lim, vbeta2, vbeta3, vbeta4)
+    CALL diffuco_bare (kjpindex, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)
 
   !! 6. Combine the alpha and beta coefficients
@@ -395,5 +397,6 @@
     ! Ajout qsintmax dans les arguments de la routine.... Nathalie / le 13-03-2006
     CALL diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, snow, &
-       & veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, vbeta, qsintmax)
+       & veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3, vbeta4, &
+       & evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax)
 
     CALL xios_orchidee_send_field("q_cdrag",q_cdrag)
@@ -1191,5 +1194,5 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE diffuco_bare (kjpindex, evap_bare_lim, vbeta2, vbeta3, vbeta4)
+  SUBROUTINE diffuco_bare (kjpindex, evap_bare_lim, evap_bare_lim_ns, vbeta2, vbeta3, vbeta4)
 
     !! 0. Variable and parameter declaration
@@ -1197,5 +1200,4 @@
     !! 0.1 Input variables
     INTEGER(i_std), INTENT(in)                         :: kjpindex       !! Domain size (-)
-    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: evap_bare_lim  !! limiting factor for bare soil evaporation when the 
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbeta2         !! Beta for Interception 
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: vbeta3         !! Beta for Transpiration 
@@ -1204,4 +1206,9 @@
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: vbeta4         !! Beta for bare soil evaporation (-)
     
+    !! 0.3 Modified variables 
+    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: evap_bare_lim  !! limiting factor for bare soil evaporation
+                                                                         !! when the 11-layer hydrology is used (-)
+    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout) :: evap_bare_lim_ns !! limiting factor for bare soil evaporation 
+                                                                         !! when the 11-layer hydrology is used (-)       
     !! 0.4 Local variables
     INTEGER(i_std)                                     :: ji
@@ -1211,4 +1218,9 @@
   !! 1. Calculation of the soil resistance and the beta (beta_4) for bare soil
 
+    ! To use the new version of hydrol_split_soil and ensure the water conservation, we must have throughout sechiba at all time:
+    !        a) evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji))
+    !        b) all the terms (vbeta4, evap_bare_lim, evap_bare_lim_ns) =0 if evap_bare_lim(ji) LE min_sechiba
+    ! This must also be kept true in diffuco_comb
+    
     DO ji = 1, kjpindex
        
@@ -1217,7 +1229,21 @@
        !!     \input{diffucobare3.tex}
        !! \endlatexonly
-       vbeta4(ji) = MIN(evap_bare_lim(ji), un - SUM(vbeta2(ji,:)+vbeta3(ji,:)))
+       
+       IF ( (evap_bare_lim(ji) .GT. min_sechiba) .AND. &  
+            ! in this case we can't have vegtot LE min_sechina, cf hydrol_soil
+            (un - SUM(vbeta2(ji,:)+vbeta3(ji,:)) .GT. min_sechiba) ) THEN 
+          ! eventually, none of the left-hand term is close to zero
+          
+          vbeta4(ji) = MIN(evap_bare_lim(ji), un - SUM(vbeta2(ji,:)+vbeta3(ji,:))) 
+          ! so that vbeta does not exceed 1
+          
+       ELSE ! instead of having very small vbeta4, we set everything to zero
+          vbeta4(ji) = zero
+          evap_bare_lim(ji) = zero
+          evap_bare_lim_ns(ji,:) = zero
+       ENDIF
+       
     END DO
-        
+    
     IF (printlev>=3) WRITE (numout,*) ' diffuco_bare done '
     
@@ -2350,5 +2376,6 @@
 
   SUBROUTINE diffuco_comb (kjpindex, humrel, rau, u, v, q_cdrag, pb, qair, temp_sol, temp_air, &
-       & snow, veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, vbeta, qsintmax)    
+       & snow, veget, lai, tot_bare_soil, vbeta1, vbeta2, vbeta3 , vbeta4, &
+       & evap_bare_lim, evap_bare_lim_ns, veget_max, vbeta, qsintmax)    
     
     ! Ajout qsintmax dans les arguments de la routine Nathalie / le 13-03-2006
@@ -2372,4 +2399,5 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: qsintmax   !! Maximum water on vegetation (kg m^{-2})
     REAL(r_std), DIMENSION (kjpindex), INTENT(in)        :: tot_bare_soil!! Total evaporating bare soil fraction 
+    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)    :: veget_max   !! Max. fraction of vegetation type (LAI->infty)
 
     !! 0.2 Output variables
@@ -2384,5 +2412,8 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta2     !! Beta for interception loss (-)
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (inout) :: vbeta3     !! Beta for Transpiration (-)
-    
+    REAL(r_std),DIMENSION (kjpindex), INTENT (inout)     :: evap_bare_lim  !! limiting factor for bare soil evaporation  
+                                                                           !! when the 11-layer hydrology is used (-)
+    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (inout):: evap_bare_lim_ns !! limiting factor for bare soil evaporation 
+                                                                             !! when the 11-layer hydrology is used (-)    
     !! 0.4 Local variables
     
@@ -2392,4 +2423,5 @@
     LOGICAL, DIMENSION(kjpindex)                         :: toveg, tosnow
     REAL(r_std)                                          :: coeff_dew_veg
+    REAL(r_std), DIMENSION(kjpindex)                     :: vegtot
 !_ ================================================================================================================================
     
@@ -2427,13 +2459,44 @@
     !!     \input{diffucocomb6.tex}
     !! \endlatexonly
+
+    ! We need to keep consistency between evap_bare_lim, evap_bare_lim_ns and vbeta4 (thus vevapnu)
+    ! or we have a water conservation issue in hydrol_split_soil
+
     DO ji = 1, kjpindex
-      IF ( toveg(ji) ) THEN
-        vbeta1(ji) = zero
-        vbeta4(ji) = tot_bare_soil(ji)
-      ENDIF
-      IF ( tosnow(ji) ) THEN
-        vbeta1(ji) = un
-        vbeta4(ji) = zero
-      ENDIF
+
+       IF ( toveg(ji) ) THEN
+
+          vbeta1(ji) = zero
+          vegtot(ji) = SUM(veget_max(ji,:))
+
+          IF ( (tot_bare_soil(ji) .GT. min_sechiba) .AND. (vegtot(ji).GT. min_sechiba) ) THEN
+             
+             vbeta4(ji) = tot_bare_soil(ji)
+             
+             ! We now have to redefine evap_bare_lim & evap_bare_lim_ns
+             IF (evap_bare_lim(ji) .GT. min_sechiba) THEN               
+                evap_bare_lim_ns(ji,:) = evap_bare_lim_ns(ji,:) * vbeta4(ji) / evap_bare_lim(ji)
+             ELSE ! we must re-invent evap_bare_lim_ns => uniform across soiltiles             
+                evap_bare_lim_ns(ji,:) = tot_bare_soil(ji)/vegtot(ji)               
+             ENDIF
+
+             evap_bare_lim(ji) = vbeta4(ji)
+             ! consistent with evap_bare_lim(ji) = SUM(evap_bare_lim_ns(ji,:)*soiltile(ji,:)*vegtot(ji))
+             ! as SUM(soiltile(ji,:)) = 1
+
+          ELSE          
+             vbeta4(ji) = zero
+             evap_bare_lim_ns(ji,:) = zero
+             evap_bare_lim(ji) = zero
+          ENDIF
+       ENDIF
+
+       IF ( tosnow(ji) ) THEN
+          vbeta1(ji) = un
+          vbeta4(ji) = zero
+          evap_bare_lim_ns(ji,:) = zero
+          evap_bare_lim(ji) = zero
+       ENDIF
+
     ENDDO
 
@@ -2518,6 +2581,11 @@
              vbeta2(ji,:)= zero
              vbeta3(ji,:)= zero
+             evap_bare_lim_ns(ji,:) = zero
+             evap_bare_lim(ji) = zero
           END IF
     ENDDO 
+
+    CALL xios_orchidee_send_field("evap_bare_lim",evap_bare_lim) 
+    CALL xios_orchidee_send_field("evap_bare_lim_ns",evap_bare_lim_ns)
 
     IF (printlev>=3) WRITE (numout,*) ' diffuco_comb done '
