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 '
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 5804)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 5805)
@@ -173,8 +173,4 @@
                                                                          !!  @tex $(kg m^{-2})$ @endtex
 !$OMP THREADPRIVATE(ae_ns)
-  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: evap_bare_lim_ns !! Limitation factor (beta) for bare soil evaporation 
-                                                                         !! per soiltile (used to deconvoluate vevapnu)  
-                                                                         !!  (0-1, unitless) 
-!$OMP THREADPRIVATE(evap_bare_lim_ns)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)    :: free_drain_coef  !! Coefficient for free drainage at bottom
                                                                          !!  (0-1, unitless) 
@@ -458,5 +454,5 @@
                                  humrel,         vegstress, drysoil_frac,                    &
                                  shumdiag_perma,    qsintveg,                        &
-                                 evap_bare_lim,  snow,      snow_age,      snow_nobio,       &
+                                 evap_bare_lim,  evap_bare_lim_ns,  snow,      snow_age,      snow_nobio,       &
                                  snow_nobio_age, snowrho,   snowtemp,      snowgrain,        &
                                  snowdz,         snowheat,  &
@@ -481,4 +477,5 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)   :: qsintveg       !! Water on vegetation due to interception
     REAL(r_std),DIMENSION (kjpindex), INTENT(out)        :: evap_bare_lim  !! Limitation factor for bare soil evaporation
+    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT(out)   :: evap_bare_lim_ns !! Limitation factor for bare soil evaporation
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: snow           !! Snow mass [Kg/m^2]
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)       :: snow_age       !! Snow age
@@ -503,5 +500,5 @@
          humrel, vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
          snowdz, snowgrain, snowrho,    snowtemp,   snowheat, &
-         drysoil_frac, evap_bare_lim)
+         drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
     
     CALL hydrol_var_init (kjpindex, veget, veget_max, &
@@ -565,5 +562,6 @@
        & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,  &
        & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
-       & humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, flood_frac, flood_res, &
+       & humrel, vegstress, drysoil_frac, evapot, evapot_penm, evap_bare_lim, evap_bare_lim_ns, &
+       & flood_frac, flood_res, &
        & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope, rest_id, hist_id, hist2_id,&
        & contfrac, stempdiag, &
@@ -643,5 +641,6 @@
 
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(inout):: qsintveg         !! Water on vegetation due to interception
-    REAL(r_std),DIMENSION (kjpindex), INTENT(inout)    :: evap_bare_lim    !! Limitation factor (beta) for bare soil evaporation    
+    REAL(r_std),DIMENSION (kjpindex), INTENT(inout)    :: evap_bare_lim    !! Limitation factor (beta) for bare soil evaporation 
+    REAL(r_std),DIMENSION (kjpindex,nstm),INTENT(inout):: evap_bare_lim_ns !! Limitation factor (beta) for bare soil evaporation    
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT(inout):: humrel           !! Relative humidity
     REAL(r_std),DIMENSION (kjpindex), INTENT (inout)   :: vevapnu          !! Bare soil evaporation
@@ -782,5 +781,5 @@
          transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
          returnflow, reinfiltration, irrigation, &
-         tot_melt,evap_bare_lim, shumdiag, shumdiag_perma, &
+         tot_melt,evap_bare_lim,evap_bare_lim_ns, shumdiag, shumdiag_perma, &
          k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,&
          stempdiag,snow,snowdz, tot_bare_soil,  u, v, tq_cdrag, &
@@ -937,5 +936,4 @@
 
     CALL xios_orchidee_send_field("frac_bare",frac_bare)
-
     CALL xios_orchidee_send_field("soilwet",soilwet)
     CALL xios_orchidee_send_field("delsoilmoist",delsoilmoist)
@@ -949,5 +947,4 @@
     CALL xios_orchidee_send_field("profil_froz_hydro_ns", profil_froz_hydro_ns)
     CALL xios_orchidee_send_field("kk_moy",kk_moy) ! in mm/d
-
 
     !! Calculate diagnostic variables on Landuse tiles for LUMIP/CMIP6
@@ -1167,5 +1164,5 @@
                               snow_nobio_age, snowrho,    snowtemp, snowdz,     &
                               snowheat,       snowgrain,  &
-                              drysoil_frac, evap_bare_lim)
+                              drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
 
     !! 0. Variable and parameter declaration
@@ -1188,4 +1185,5 @@
     REAL(r_std),DIMENSION (kjpindex),INTENT(in)          :: drysoil_frac   !! function of litter wetness
     REAL(r_std),DIMENSION (kjpindex),INTENT(in)          :: evap_bare_lim
+    REAL(r_std),DIMENSION (kjpindex,nstm),INTENT(in)     :: evap_bare_lim_ns
 
     !! 0.4 Local variables
@@ -1294,5 +1292,5 @@
        humrel,  vegstress, snow,       snow_age,   snow_nobio, snow_nobio_age, qsintveg, &
        snowdz,  snowgrain, snowrho,    snowtemp,   snowheat, &
-       drysoil_frac, evap_bare_lim)
+       drysoil_frac, evap_bare_lim, evap_bare_lim_ns)
     
 
@@ -1324,4 +1322,5 @@
     REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: drysoil_frac     !! function of litter wetness
     REAL(r_std),DIMENSION (kjpindex),INTENT(out)          :: evap_bare_lim
+    REAL(r_std),DIMENSION (kjpindex,nstm),INTENT(out)     :: evap_bare_lim_ns
 
     !! 0.4 Local variables
@@ -1701,7 +1700,4 @@
     ALLOCATE (ae_ns(kjpindex,nstm),stat=ier) 
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable ae_ns','','')
-
-    ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier) 
-    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable evap_bare_lim_ns','','')
 
     ALLOCATE (rootsink(kjpindex,nslm,nstm),stat=ier) 
@@ -2478,5 +2474,4 @@
     IF (ALLOCATED  (water2infilt)) DEALLOCATE (water2infilt)
     IF (ALLOCATED  (ae_ns)) DEALLOCATE (ae_ns)
-    IF (ALLOCATED  (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
     IF (ALLOCATED  (rootsink)) DEALLOCATE (rootsink)
     IF (ALLOCATED  (subsnowveg)) DEALLOCATE (subsnowveg)
@@ -3743,5 +3738,5 @@
        & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &
        & returnflow, reinfiltration, irrigation, &
-       & tot_melt, evap_bare_lim, shumdiag, shumdiag_perma,&
+       & tot_melt, evap_bare_lim, evap_bare_lim_ns, shumdiag, shumdiag_perma,&
        & k_litt, litterhumdiag, humrel,vegstress, drysoil_frac, &
        & stempdiag,snow, &
@@ -3791,4 +3786,6 @@
                                                                                  !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
     REAL(r_std), DIMENSION (kjpindex), INTENT(out)           :: evap_bare_lim    !! Limitation factor (beta) for bare soil evaporation  
+                                                                                 !! on each soil column (unitless, [0-1])
+    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(out)      :: evap_bare_lim_ns !! Limitation factor (beta) for bare soil evaporation  
                                                                                  !! on each soil column (unitless, [0-1])
     REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: shumdiag         !! Relative soil moisture in each diag soil layer 
@@ -3970,5 +3967,6 @@
     !  At step 7, the reverse transformation is used for the fluxes produced in hydrol_soil
     !    flux_cell(ji)=sum(flux_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji))
-    CALL hydrol_split_soil (kjpindex, veget_max, soiltile, vevapnu, transpir, humrel, evap_bare_lim, tot_bare_soil)
+    CALL hydrol_split_soil (kjpindex, veget_max, soiltile, vevapnu, transpir, humrel, &
+         evap_bare_lim, evap_bare_lim_ns, tot_bare_soil)
     
     !! 0.3 Common variables related to routing, with all return flow applied to the soil surface
@@ -4807,9 +4805,5 @@
                 evap_soil(ji) = evapot_penm(ji)
              ENDIF
-             
-       ! AD16*** et si evap_bare_lim_ns<0 ?? car on suppose que tmcint > tmc(new)
-       ! (water2inflit permet de propager de la ponded water d'un pas de temps a l'autre:
-       ! peut-on s'en servir pour creer des cas d'evapnu potentielle negative ? a gerer dans diffuco ?)
-             
+                          
              flux_top(ji) = evap_soil(ji) * &
                   AINT(frac_bare_ns(ji,jst)+un-min_sechiba)
@@ -4998,4 +4992,6 @@
        r_soil(ji) =  SUM(r_soil_ns(ji,:)*vegtot(ji)*soiltile(ji,:))
     ENDDO
+    ! si vegtot LE min_sechiba, evap_bare_lim_ns et evap_bare_lim valent zero
+
 
     !! 10. XIOS export of local variables, including water conservation checks
@@ -6109,5 +6105,4 @@
 !! 1.1 Throughfall
 !! 1.2 Bare soil evaporation
-!! 1.2.1 vevapnu_old
 !! 1.2.2 ae_ns new
 !! 1.3 transpiration
@@ -6130,5 +6125,6 @@
 !_ hydrol_split_soil
 
-  SUBROUTINE hydrol_split_soil (kjpindex, veget_max, soiltile, vevapnu, transpir, humrel, evap_bare_lim, tot_bare_soil)
+  SUBROUTINE hydrol_split_soil (kjpindex, veget_max, soiltile, vevapnu, transpir, humrel, &
+       evap_bare_lim, evap_bare_lim_ns, tot_bare_soil)
     ! 
     ! interface description
@@ -6145,4 +6141,5 @@
     REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in)       :: humrel           !! Relative humidity
     REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: evap_bare_lim    !!   
+    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT(in)       :: evap_bare_lim_ns !!   
     REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: tot_bare_soil    !! Total evaporating bare soil fraction 
 
@@ -6150,9 +6147,8 @@
 
     INTEGER(i_std)                                :: ji, jv, jsl, jst
-    REAL(r_std), DIMENSION (kjpindex)             :: vevapnu_old
     REAL(r_std), DIMENSION (kjpindex)             :: tmp_check1
     REAL(r_std), DIMENSION (kjpindex)             :: tmp_check2
     REAL(r_std), DIMENSION (kjpindex,nstm)        :: tmp_check3
-
+    LOGICAL                                       :: error
 !_ ================================================================================================================================
     
@@ -6185,59 +6181,14 @@
     END DO
     
-    !! 1.2 Bare soil evaporation
-    !! 1.2.1 vevapnu_old
-! AD16*** vevapnu_old ne sert que pour le split suivant de vevapnu (issu de enerbil) en ae_ns pour hydrol_soil
-!           mais il ne semble y avoir aucune bonne raison de contraindre ae_ns en fonction de vevapnu_old
-    vevapnu_old(:)=zero
+    !! 1.2 Bare soil evaporation and ae_ns
+    ae_ns(:,:)=zero
     DO jst=1,nstm
        DO ji=1,kjpindex
-          IF ( vegtot(ji) .GT. min_sechiba) THEN
-             vevapnu_old(ji)=vevapnu_old(ji)+ &
-                  & ae_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji)
+          IF(evap_bare_lim(ji).GT.min_sechiba) THEN       
+             ae_ns(ji,jst) = vevapnu(ji) * evap_bare_lim_ns(ji,jst)/evap_bare_lim(ji)
           ENDIF
-       END DO
-    END DO
-    
-    !! 1.2.2 ae_ns new
-! AD16*** les lignes ci-dessous sont excessivement compliquees et ne garantissent pas que ae_ns = 0 si evap_bare_lim=0 
-!           c'est notamment le cas pour les 3emes et 6emes conditions
-    DO jst=1,nstm
-       DO ji=1,kjpindex
-          IF (vevapnu_old(ji).GT.min_sechiba) THEN   
-             IF(evap_bare_lim(ji).GT.min_sechiba) THEN       
-                ae_ns(ji,jst) = vevapnu(ji) * evap_bare_lim_ns(ji,jst)/evap_bare_lim(ji) 
-             ELSE
-                IF(vevapnu_old(ji).GT.min_sechiba) THEN  
-                   ae_ns(ji,jst)=ae_ns(ji,jst) * vevapnu(ji)/vevapnu_old(ji) ! 3Ãšme condition
-                ELSE
-                   ae_ns(ji,jst)=zero
-                ENDIF
-             ENDIF
-          ELSEIF(frac_bare_ns(ji,jst).GT.min_sechiba) THEN
-             IF(evap_bare_lim(ji).GT.min_sechiba) THEN  
-                ae_ns(ji,jst) = vevapnu(ji) * evap_bare_lim_ns(ji,jst)/evap_bare_lim(ji)
-             ELSE
-                IF(tot_bare_soil(ji).GT.min_sechiba) THEN  
-                   ae_ns(ji,jst) = vevapnu(ji) * frac_bare_ns(ji,jst)/tot_bare_soil(ji) ! 6Ãšme condition
-                ELSE
-                   ae_ns(ji,jst) = zero
-                ENDIF
-             ENDIF
-          ENDIF
-       END DO
-    END DO
-! ADNV27072016: we believe the following block should be used (tests needed before committ, since AD16*** had pb with it)    
-!!$    ! given the definition of evap_bare_lim, it leads to sum(ae_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji))=vevapnu(ji)
-!!$    ae_ns(:,:)=zero
-!!$    DO jst=1,nstm
-!!$       DO ji=1,kjpindex
-!!$          IF(evap_bare_lim(ji).GT.min_sechiba) THEN       
-!!$             ae_ns(ji,jst) = vevapnu(ji) * evap_bare_lim_ns(ji,jst)/evap_bare_lim(ji)
-!            ELSE
-!               ae_ns(ji,jst) = zero
-!!$          ENDIF
-!!$       ENDDO
-!!$    ENDDO
-    
+       ENDDO
+    ENDDO
+
     !! 1.3 transpiration
     ! Transformation from transpir (flux from PFT jv in m2 of grid-mesh)
@@ -6276,5 +6227,150 @@
 
 
-    !!! ADNV270716 *** we are here
+    !! 2. Verification: Check if the deconvolution is correct and conserves the fluxes (grid-cell average)
+
+    IF (check_cwrr) THEN
+
+       error=.FALSE.
+
+       !! 2.1 precisol 
+
+       tmp_check1(:)=zero
+       DO jst=1,nstm
+          DO ji=1,kjpindex
+             tmp_check1(ji)=tmp_check1(ji) + precisol_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji)
+          END DO
+       END DO
+       
+       tmp_check2(:)=zero  
+       DO jv=1,nvm
+          DO ji=1,kjpindex
+             tmp_check2(ji)=tmp_check2(ji) + precisol(ji,jv)
+          END DO
+       END DO
+
+       DO ji=1,kjpindex   
+          IF(ABS(tmp_check1(ji) - tmp_check2(ji)).GT.allowed_err) THEN
+             WRITE(numout,*) 'PRECISOL SPLIT FALSE:ji=',ji,tmp_check1(ji),tmp_check2(ji)
+             WRITE(numout,*) 'err',ABS(tmp_check1(ji)- tmp_check2(ji))
+             WRITE(numout,*) 'vegtot',vegtot(ji)
+             DO jv=1,nvm
+                WRITE(numout,'(a,i2.2,"|",F13.4,"|",F13.4,"|",3(F9.6))') &
+                     'jv,veget_max, precisol, vegetmax_soil ', &
+                     jv,veget_max(ji,jv),precisol(ji,jv),vegetmax_soil(ji,jv,:)
+             END DO
+             DO jst=1,nstm
+                WRITE(numout,*) 'jst,precisol_ns',jst,precisol_ns(ji,jst)
+                WRITE(numout,*) 'soiltile', soiltile(ji,jst)
+             END DO
+             error=.TRUE.
+             CALL ipslerr_p(2, 'hydrol_split_soil', 'We will STOP in the end of this subroutine.',&
+                  & 'check_CWRR','PRECISOL SPLIT FALSE')
+          ENDIF
+       END DO
+       
+       !! 2.2 ae_ns and evapnu
+
+       tmp_check1(:)=zero
+       DO jst=1,nstm
+          DO ji=1,kjpindex
+             tmp_check1(ji)=tmp_check1(ji) + ae_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji)
+          END DO
+       END DO
+
+       DO ji=1,kjpindex   
+
+          IF(ABS(tmp_check1(ji) - vevapnu(ji)).GT.allowed_err) THEN
+             WRITE(numout,*) 'VEVAPNU SPLIT FALSE:ji, Sum(ae_ns), vevapnu =',ji,tmp_check1(ji),vevapnu(ji)
+             WRITE(numout,*) 'err',ABS(tmp_check1(ji)- vevapnu(ji))
+             WRITE(numout,*) 'ae_ns',ae_ns(ji,:)
+             WRITE(numout,*) 'vegtot',vegtot(ji)
+             WRITE(numout,*) 'evap_bare_lim, evap_bare_lim_ns',evap_bare_lim(ji), evap_bare_lim_ns(ji,:)
+             DO jst=1,nstm
+                WRITE(numout,*) 'jst,ae_ns',jst,ae_ns(ji,jst)
+                WRITE(numout,*) 'soiltile', soiltile(ji,jst)
+             END DO
+             error=.TRUE.
+             CALL ipslerr_p(2, 'hydrol_split_soil', 'We will STOP in the end of this subroutine.',&
+                  & 'check_CWRR','VEVAPNU SPLIT FALSE')
+          ENDIF
+       ENDDO
+
+    !! 2.3 transpiration
+
+       tmp_check1(:)=zero
+       DO jst=1,nstm
+          DO ji=1,kjpindex
+             tmp_check1(ji)=tmp_check1(ji) + tr_ns(ji,jst)*soiltile(ji,jst)*vegtot(ji)
+          END DO
+       END DO
+       
+       tmp_check2(:)=zero  
+       DO jv=1,nvm
+          DO ji=1,kjpindex
+             tmp_check2(ji)=tmp_check2(ji) + transpir(ji,jv)
+          END DO
+       END DO
+
+       DO ji=1,kjpindex   
+          IF(ABS(tmp_check1(ji)- tmp_check2(ji)).GT.allowed_err) THEN
+             WRITE(numout,*) 'TRANSPIR SPLIT FALSE:ji=',ji,tmp_check1(ji),tmp_check2(ji)
+             WRITE(numout,*) 'err',ABS(tmp_check1(ji)- tmp_check2(ji))
+             WRITE(numout,*) 'vegtot',vegtot(ji)
+             DO jv=1,nvm
+                WRITE(numout,*) 'jv,veget_max, transpir',jv,veget_max(ji,jv),transpir(ji,jv)
+                DO jst=1,nstm
+                   WRITE(numout,*) 'vegetmax_soil:ji,jv,jst',ji,jv,jst,vegetmax_soil(ji,jv,jst)
+                END DO
+             END DO
+             DO jst=1,nstm
+                WRITE(numout,*) 'jst,tr_ns',jst,tr_ns(ji,jst)
+                WRITE(numout,*) 'soiltile', soiltile(ji,jst)
+             END DO
+             error=.TRUE.
+             CALL ipslerr_p(2, 'hydrol_split_soil', 'We will STOP in the end of this subroutine.',&
+                  & 'check_CWRR','TRANSPIR SPLIT FALSE')
+          ENDIF
+
+       END DO
+
+    !! 2.4 root sink
+
+       tmp_check3(:,:)=zero
+       DO jst=1,nstm
+          DO jsl=1,nslm
+             DO ji=1,kjpindex
+                tmp_check3(ji,jst)=tmp_check3(ji,jst) + rootsink(ji,jsl,jst)
+             END DO
+          END DO
+       ENDDO
+
+       DO jst=1,nstm
+          DO ji=1,kjpindex
+             IF(ABS(tmp_check3(ji,jst) - tr_ns(ji,jst)).GT.allowed_err) THEN
+                WRITE(numout,*) 'ROOTSINK SPLIT FALSE:ji,jst=', ji,jst,&
+                     & tmp_check3(ji,jst),tr_ns(ji,jst)
+                WRITE(numout,*) 'err',ABS(tmp_check3(ji,jst)- tr_ns(ji,jst))
+                WRITE(numout,*) 'HUMREL(jv=1:13)',humrel(ji,:)
+                WRITE(numout,*) 'TRANSPIR',transpir(ji,:)
+                DO jv=1,nvm 
+                   WRITE(numout,*) 'jv=',jv,'us=',us(ji,jv,jst,:)
+                ENDDO
+                error=.TRUE.
+                CALL ipslerr_p(2, 'hydrol_split_soil', 'We will STOP in the end of this subroutine.',&
+                  & 'check_CWRR','ROOTSINK SPLIT FALSE')
+             ENDIF
+          END DO
+       END DO
+
+
+       !! Exit if error was found previously in this subroutine
+       IF ( error ) THEN
+          WRITE(numout,*) 'One or more errors have been detected in hydrol_split_soil. Model stops.'
+          CALL ipslerr_p(3, 'hydrol_split_soil', 'We will STOP now.',&
+               & 'One or several fatal errors were found previously.','')
+       END IF
+
+    ENDIF ! end of check_cwrr
+
 
   END SUBROUTINE hydrol_split_soil
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 5804)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 5805)
@@ -208,5 +208,6 @@
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: evap_bare_lim  !! Bare soil stress
 !$OMP THREADPRIVATE(evap_bare_lim)
-
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: evap_bare_lim_ns !! Bare soil stress
+!$OMP THREADPRIVATE(evap_bare_lim_ns)
   REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:)  :: co2_flux       !! CO2 flux (gC/m**2 of average ground/s)
 !$OMP THREADPRIVATE(co2_flux)
@@ -463,5 +464,5 @@
          humrel,         vegstress, drysoil_frac,                    &
          shumdiag_perma,    qsintveg,                        &
-         evap_bare_lim,  snow,      snow_age,      snow_nobio,       &
+         evap_bare_lim, evap_bare_lim_ns,  snow,      snow_age,      snow_nobio,       &
          snow_nobio_age, snowrho,   snowtemp,      snowgrain,        &
          snowdz,         snowheat,  &
@@ -703,5 +704,5 @@
     CALL diffuco_main (kjit, kjpindex, index, indexveg, indexlai, u, v, &
          & zlev, z0m, z0h, roughheight, temp_sol, temp_air, temp_growth, rau, tq_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, &
@@ -726,5 +727,5 @@
          & qsintmax, qsintveg, vevapnu, vevapsno, vevapflo, snow, snow_age, snow_nobio, snow_nobio_age,  &
          & tot_melt, transpir, precip_rain, precip_snow, returnflow, reinfiltration, irrigation, &
-         & humrel, vegstress, drysoil_frac, evapot, evapot_corr, evap_bare_lim, flood_frac, flood_res, &
+         & 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,&
          & rest_id, hist_id, hist2_id,&
@@ -1349,5 +1350,5 @@
          snow_nobio_age, snowrho,  snowtemp, snowdz,     &
          snowheat,       snowgrain,    &
-         drysoil_frac, evap_bare_lim)
+         drysoil_frac,   evap_bare_lim, evap_bare_lim_ns)
     
     !! 4. Call condveg to write surface variables to restart files
@@ -1499,4 +1500,7 @@
     ALLOCATE (evap_bare_lim(kjpindex),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim','','')
+
+    ALLOCATE (evap_bare_lim_ns(kjpindex,nstm),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'sechiba_init','Pb in alloc for evap_bare_lim_ns','','')
 
     ALLOCATE (evapot(kjpindex),stat=ier)
@@ -1867,4 +1871,5 @@
     IF ( ALLOCATED (drysoil_frac)) DEALLOCATE (drysoil_frac)
     IF ( ALLOCATED (evap_bare_lim)) DEALLOCATE (evap_bare_lim)
+    IF ( ALLOCATED (evap_bare_lim_ns)) DEALLOCATE (evap_bare_lim_ns)
     IF ( ALLOCATED (evapot)) DEALLOCATE (evapot)
     IF ( ALLOCATED (evapot_corr)) DEALLOCATE (evapot_corr)
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml	(revision 5804)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/field_def_orchidee.xml	(revision 5805)
@@ -71,4 +71,6 @@
     <field id="drainage" name="drainage_sec" long_name="Deep drainage" unit="mm/s"/>
     <field id="mrro" name="mrro" long_name="Total Runoff" unit="kg m-2 s-1" > runoff + drainage </field>
+    <field id="evap_bare_lim" name="evap_bare_lim" long_name="Bare soil resistance (hydrol)" unit="0-1" />
+    <field id="evap_bare_lim_ns" name="evap_bare_lim_ns" long_name="Bare soil resistance (hydrol) for soil type" unit="0-1" grid_ref="grid_nstm" />
     <field id="evapnu_soil" name="evapnu_soil" long_name="Bare soil evap for soil type" unit="mm/d" grid_ref="grid_nstm" > this*86400 </field>
     <field id="drainage_soil" name="drainage_soil" long_name="Drainage for soil type" unit="mm/d" grid_ref="grid_nstm" > this*86400 </field>
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml	(revision 5804)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_xml/file_def_orchidee.xml	(revision 5805)
@@ -70,4 +70,6 @@
     <field field_ref="hydrographs" level="2"/>
     <field field_ref="vevapnu" name="evapnu" level="0" unit="mm/d" > @vevapnu*86400  </field>
+    <field field_ref="evap_bare_lim" level="2"/>
+    <field field_ref="evap_bare_lim_ns" grid_ref="grid_nstm_out" level="2"/>
     <field field_ref="evapnu_soil" grid_ref="grid_nstm_out" level="2"/>
     <field field_ref="drainage_soil" grid_ref="grid_nstm_out" level="2"/>
