Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 7476)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/hydrol.f90	(revision 7709)
@@ -24,4 +24,6 @@
 !!                    We can also impose kfact_root=1 in all soil layers to cancel the effect of
 !!                    roots on ks profile (keyword KFACT_ROOT_CONST).
+!!                    July 2022: New irrigation scheme. Here new irrigation demand based in
+!!                    soil moisture deficit, and irrigation application.
 !!                 
 !! REFERENCE(S) :
@@ -366,4 +368,10 @@
                                                                          !!  @tex $(m^{3} m^{-3})$ @endtex
 !$OMP THREADPRIVATE(mc)
+  REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:)    :: root_mc_fc         !! Max Field Capacity moisture content, for layers with roots, in soil tile of irrig_st
+                                                                         !!  @tex $(kg m^{-2})$ @endtex
+!$OMP THREADPRIVATE(root_mc_fc)
+  INTEGER(i_std), ALLOCATABLE, SAVE, DIMENSION (:) :: nslm_root          !! max. layers defining the root zone
+                                                                         !!  @tex $(layer)$ @endtex
+!$OMP THREADPRIVATE(nslm_root)
 
    REAL(r_std), ALLOCATABLE, SAVE, DIMENSION (:,:,:) :: mc_read_prev       !! Soil moisture from file at previous timestep in the file
@@ -565,5 +573,5 @@
        & 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,&
+       & shumdiag,shumdiag_perma, k_litt, litterhumdiag, soilcap, soiltile, fraclut, reinf_slope_soil, rest_id, hist_id, hist2_id,&
        & contfrac, stempdiag, &
        & temp_air, pb, u, v, tq_cdrag, swnet, pgflux, &
@@ -571,5 +579,5 @@
        & grndflux,gtemp,tot_bare_soil, &
        & lambda_snow,cgrnd_snow,dgrnd_snow,frac_snow_veg,temp_sol_add, &
-       & mc_layh, mcl_layh, soilmoist_out )
+       & mc_layh, mcl_layh, soilmoist_out, root_deficit )
 
     !! 0. Variable and parameter declaration
@@ -607,5 +615,5 @@
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: qsintmax         !! Maximum water on vegetation for interception
     REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)  :: transpir         !! Transpiration
-    REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: reinf_slope      !! Slope coef
+    REAL(r_std),DIMENSION (kjpindex,nstm), INTENT (in) :: reinf_slope_soil !! Slope coef per soil tile
 
     REAL(r_std),DIMENSION (kjpindex), INTENT (in)      :: ks               !! Hydraulic conductivity at saturation (mm {-1})
@@ -647,5 +655,6 @@
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: tot_melt         !! Total melt    
     REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: floodout         !! Flux out of floodplains
-    
+    REAL(r_std),DIMENSION (kjpindex), INTENT (out)     :: root_deficit     !! water deficit to reach field capacity of soil
+
     !! 0.3 Modified variables
 
@@ -794,5 +803,5 @@
     !! 3.5 computes soil hydrology ==>hydrol_soil
 
-    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope,  &
+    CALL hydrol_soil(ks, nvan, avan, mcr, mcs, mcfc, mcw, kjpindex, veget_max, soiltile, njsc, reinf_slope_soil,  &
          transpir, vevapnu, evapot, evapot_penm, runoff, drainage, & 
          returnflow, reinfiltration, irrigation, &
@@ -800,9 +809,13 @@
          k_litt, litterhumdiag, humrel, vegstress, drysoil_frac,&
          stempdiag,snow,snowdz, tot_bare_soil,  u, v, tq_cdrag, &
-         mc_layh, mcl_layh)
+         mc_layh, mcl_layh, root_deficit, veget)
 
     ! The update fluxes come from hydrol_vegupd
     drainage(:) =  drainage(:) +  drain_upd(:)
     runoff(:) =  runoff(:) +  runoff_upd(:)
+
+   ! Output irrigation related variables
+    CALL xios_orchidee_send_field("root_deficit", root_deficit)
+    CALL xios_orchidee_send_field("root_mc_fc", root_mc_fc)
 
 
@@ -1691,4 +1704,10 @@
     ALLOCATE (mc(kjpindex,nslm,nstm),stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable mc','','')
+
+    ALLOCATE (root_mc_fc(kjpindex),stat=ier) !
+    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable root_mc_fc','','') !
+
+    ALLOCATE (nslm_root(kjpindex),stat=ier) !
+    IF (ier /= 0) CALL ipslerr_p(3,'hydrol_init','Problem in allocate of variable nslm_root','','') !
 
 
@@ -2365,4 +2384,6 @@
     IF (ALLOCATED  (vegetmax_soil)) DEALLOCATE (vegetmax_soil)
     IF (ALLOCATED  (mc)) DEALLOCATE (mc)
+    IF (ALLOCATED  (root_mc_fc)) DEALLOCATE (root_mc_fc)
+    IF (ALLOCATED  (nslm_root)) DEALLOCATE (nslm_root)
     IF (ALLOCATED  (soilmoist)) DEALLOCATE (soilmoist)
     IF (ALLOCATED  (soilmoist_liquid)) DEALLOCATE (soilmoist_liquid)
@@ -2728,4 +2749,8 @@
     INTEGER(i_std)                                      :: iiref           !! To identify the mc_lins where k_lin and d_lin
                                                                            !! need special treatment
+    REAL(r_std), DIMENSION (kjpindex,nvm)               :: cum_nroot       !! local variable to acummulate nroot
+    INTEGER(i_std), DIMENSION (kjpindex,nvm)            :: nslm_pft_root   !! Deeper layer with cum_nroot > cum_nroot_thr per PFT,
+    REAL(r_std), DIMENSION (kjpindex,nslm)              :: smf             !! Soil moisture of each layer at field capacity
+                                                                           !!  @tex $(kg m^{-2})$ @endtex
 
 !_ ================================================================================================================================
@@ -2880,12 +2905,20 @@
     !! 2 Compute the root density profile if not ok_dynroot
     !!   For the case with ok_dynroot, the calculations are done at each time step in hydrol_soil
+    cum_nroot(:,:) = zero
+    nslm_pft_root(:,:) = nslm
+    nslm_root(:) =  nslm
     IF (.NOT. ok_dynroot) THEN
-       DO ji=1, kjpindex
-          !-
-          !! The three following equations concerning nroot computation are derived from the integrals 
-          !! of equations C9 to C11 of De Rosnay's (1999) PhD thesis (page 158).
-          !! The occasional absence of minus sign before humcste parameter is correct.
+       !! Calculation of nroot and of new nslm_root for irrigation 
+       !! The three following equations concerning nroot computation are derived from the integrals
+       !! of equations C9 to C11 of De Rosnay's (1999) PhD thesis (page 158).
+       !! The occasional absence of minus sign before humcste parameter is correct.
+       ! First layer
+       nroot(:,:,1) = zero
+       !From 2 to nslm-1 layers
+       DO jsl = 2, nslm-1
           DO jv = 1,nvm
-             DO jsl = 2, nslm-1
+             DO ji=1, kjpindex
+                
+                
                 nroot(ji,jv,jsl) = (EXP(-humcste(jv)*zz(jsl)/mille)) * &
                      & (EXP(humcste(jv)*dz(jsl)/mille/deux) - &
@@ -2893,7 +2926,20 @@
                      & (EXP(-humcste(jv)*dz(2)/mille/deux) &
                      & -EXP(-humcste(jv)*zz(nslm)/mille))
+                !We acum. nroot, and change nslm_pft_root if necessary
+                cum_nroot(ji,jv) = cum_nroot(ji,jv) + nroot(ji,jv,jsl) !
+                
+                IF(cum_nroot(ji,jv) > cum_nroot_thr .AND. nslm_pft_root(ji,jv) .GE. jsl ) THEN
+                   nslm_pft_root(ji,jv) =  jsl
+                ENDIF
              ENDDO
-             nroot(ji,jv,1) = zero
-
+             
+          ENDDO
+          
+       ENDDO
+       !Last layer
+       !No need to put and IF here, if it is the case, nslm_pft_root is already
+       ! equal to nslm,
+       DO jv = 1,nvm
+          DO ji=1, kjpindex
              nroot(ji,jv,nslm) = (EXP(humcste(jv)*dz(nslm)/mille/deux) -un) * &
                   & EXP(-humcste(jv)*zz(nslm)/mille) / &
@@ -2902,8 +2948,34 @@
           ENDDO
        ENDDO
+       !New loop to compute min. soil layer per cell, using just PFT that
+       !are not natural and exist inside the cell
+       DO jv=1, nvm
+          IF(.NOT. natural(jv)) THEN
+             DO ji=1,kjpindex
+                IF(veget_max(ji, jv) > min_sechiba) THEN
+                   nslm_root(ji) = MIN(nslm_root(ji), nslm_pft_root(ji,jv))
+                ENDIF
+             ENDDO
+          ENDIF
+       ENDDO
+       
     END IF
 
- 
-
+    ! Calculates field capacity soil moisture per soil layers
+    ! then calculate field capacity soil moisture over root zone
+    smf(:,:) = zero
+    root_mc_fc(:) = zero
+    smf(:,1) = dz(2) * (quatre*mcfc(:))/huit
+    
+    DO jsl = 2,nslm-1
+       smf(:,jsl) = dz(jsl) * ( quatre*mcfc(:) )/huit &
+            + dz(jsl+1) * ( quatre*mcfc(:) )/huit
+    ENDDO
+    
+    smf(:,nslm) = dz(nslm) * (quatre*mcfc(:))/huit
+    DO ji = 1,kjpindex
+       root_mc_fc(ji) = SUM(smf(ji,1:nslm_root(ji) ))
+    ENDDO
+    
     !-
     !! 3 Compute the profile for a and n
@@ -3592,5 +3664,5 @@
 !_ hydrol_soil
   SUBROUTINE hydrol_soil (ks, nvan, avan, mcr, mcs, mcfc, mcw, &
-       kjpindex, veget_max, soiltile, njsc, reinf_slope, &
+       kjpindex, veget_max, soiltile, njsc, reinf_slope_soil, &
        & transpir, vevapnu, evapot, evapot_penm, runoff, drainage, &
        & returnflow, reinfiltration, irrigation, &
@@ -3598,5 +3670,5 @@
        & k_litt, litterhumdiag, humrel,vegstress, drysoil_frac, &
        & stempdiag,snow, &
-       & snowdz, tot_bare_soil, u, v, tq_cdrag, mc_layh, mcl_layh)
+       & snowdz, tot_bare_soil, u, v, tq_cdrag, mc_layh, mcl_layh, root_deficit, veget)
     ! 
     ! interface description
@@ -3608,4 +3680,5 @@
     INTEGER(i_std), INTENT(in)                               :: kjpindex 
     REAL(r_std), DIMENSION (kjpindex,nvm), INTENT (in)       :: veget_max        !! Map of max vegetation types [-]
+    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: veget            !! Fraction of vegetation type
     INTEGER(i_std),DIMENSION (kjpindex), INTENT (in)         :: njsc             !! Index of the dominant soil textural class 
                                                                                  !!   in the grid cell (1-nscm, unitless)
@@ -3622,5 +3695,5 @@
     REAL(r_std), DIMENSION (kjpindex,nvm), INTENT(in)        :: transpir         !! Transpiration  
                                                                                  !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
-    REAL(r_std), DIMENSION (kjpindex), INTENT (in)           :: reinf_slope      !! Fraction of surface runoff that reinfiltrates
+    REAL(r_std), DIMENSION (kjpindex,nstm), INTENT (in)      :: reinf_slope_soil !! Fraction of surface runoff that reinfiltrates per soil tile
                                                                                  !!  (unitless, [0-1])
     REAL(r_std), DIMENSION (kjpindex), INTENT(in)            :: returnflow       !! Water returning to the soil from the bottom
@@ -3672,5 +3745,7 @@
     REAL(r_std), DIMENSION (kjpindex,nslm), INTENT (out)     :: mcl_layh         !! Volumetric liquid water content for each soil layer 
                                                                                  !! averaged over the mesh (for thermosoil)
-                                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex 
+                                                                                 !!  @tex $(m^{3} m^{-3})$ @endtex     
+    REAL(r_std),DIMENSION (kjpindex), INTENT(out)            :: root_deficit     !! water deficit to reach SM target of soil column, for irrigation demand
+    
     !! 0.3 Modified variables
 
@@ -3698,5 +3773,5 @@
     REAL(r_std), DIMENSION(kjpindex)               :: reinfiltration_soil        !! Water from the routing back to the top of the 
                                                                                  !! soil @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
-    REAL(r_std), DIMENSION(kjpindex)               :: irrigation_soil            !! Water from irrigation returning to soil moisture 
+    REAL(r_std), DIMENSION(kjpindex,nstm)          :: irrigation_soil            !! Water from irrigation returning to soil moisture per soil tile
                                                                                  !!  @tex $(kg m^{-2} dt\_sechiba^{-1})$ @endtex
     REAL(r_std), DIMENSION(kjpindex)               :: flux_infilt                !! Water to infiltrate
@@ -3766,5 +3841,6 @@
     REAL(r_std), DIMENSION (kjpindex)              :: check_over                  !! Water conservation diagnostic at routine scale 
     REAL(r_std), DIMENSION (kjpindex)              :: check_under                 !! Water conservation diagnostic at routine scale
-
+    ! For irrigation triggering
+    INTEGER(i_std), DIMENSION (kjpindex)           :: lai_irrig_trig              !! Number of PFT per cell with LAI> LAI_IRRIG_MIN -
     ! Diagnostic of the vertical soil water fluxes  
     REAL(r_std), DIMENSION (kjpindex,nslm)         :: qflux                       !! Local upward flux into soil layer 
@@ -3796,5 +3872,5 @@
     returnflow_soil(:) = zero
     reinfiltration_soil(:) = zero
-    irrigation_soil(:) = zero
+    irrigation_soil(:,:) = zero
     qflux_ns(:,:,:) = zero
     mc_layh(:,:) = zero ! for thermosoil
@@ -3838,17 +3914,42 @@
     ! AD16*** The transformation in 0.2 and 0.3 is likely to induce conservation problems
     !         when tot_frac_nobio NE 0, since sum(soiltile) NE vegtot in this case
-    
-    DO ji=1,kjpindex
-       IF(vegtot(ji).GT.min_sechiba) THEN
-          ! returnflow_soil is assumed to enter from the bottom, but it is not possible with CWRR
-          returnflow_soil(ji) = zero
-          reinfiltration_soil(ji) = (returnflow(ji) + reinfiltration(ji))/vegtot(ji)
-          irrigation_soil(ji) = irrigation(ji)/vegtot(ji)
-       ELSE
-          returnflow_soil(ji) = zero
-          reinfiltration_soil(ji) = zero
-          irrigation_soil(ji) = zero
+    IF (.NOT. old_irrig_scheme) THEN !
+       IF (.NOT. irrigated_soiltile) THEN
+          DO ji=1,kjpindex
+             IF(vegtot(ji).GT.min_sechiba ) THEN
+                returnflow_soil(ji) = zero
+                reinfiltration_soil(ji) = (returnflow(ji) + reinfiltration(ji))/vegtot(ji)
+                IF(soiltile(ji, irrig_st).GT.min_sechiba) THEN
+                   !irrigation_soil(ji, 1:2) = 0, if irrig_st = 3. Not put because Values
+                   !are already zero, due to initialization
+                   irrigation_soil(ji, irrig_st) = irrigation(ji) / (soiltile(ji, irrig_st) * vegtot(ji) )
+                   !Irrigation is kg/m2 of grid cell. Here, all that water is put on
+                   !irrig_st (irrigated soil tile), by default = 3, for the others
+                   !soil tiles irrigation = zero
+                ENDIF
+             ENDIF
+          ENDDO
        ENDIF
-    ENDDO       
+    ELSE  !
+       DO ji=1,kjpindex
+          IF(vegtot(ji).GT.min_sechiba) THEN
+             ! returnflow_soil is assumed to enter from the bottom, but it is not possible with CWRR
+             returnflow_soil(ji) = zero
+             reinfiltration_soil(ji) = (returnflow(ji) + reinfiltration(ji))/vegtot(ji)
+             irrigation_soil(ji,:) = irrigation(ji)/vegtot(ji)
+             ! irrigation_soil(ji) = irrigation(ji)/vegtot(ji)
+             ! Computed for all the grid cell. New way is equivalent, and coherent
+             ! with irrigation_soil new dimensions (cells, soil tiles)
+             ! Irrigation is kg/m2 of grid cell. For the old irrig. scheme,
+             ! irrigation soil is the same for every soil tile
+             ! Next lines are in tag 2.0, deleted because values are already init to zero
+             ! ELSE
+             ! returnflow_soil(ji) = zero
+             ! reinfiltration_soil(ji) = zero
+             ! irrigation_soil(ji) = zero
+             ! ENDIF
+          ENDIF
+       ENDDO
+    ENDIF
     
     !! -- START MAIN LOOP (prognostic loop to update mc and mcl) OVER SOILTILES
@@ -3901,5 +4002,5 @@
        ! Here, temp is used as a temporary variable to store the min of water to infiltrate vs evaporate
        DO ji = 1, kjpindex
-          temp(ji) = MIN(water2infilt(ji,jst) + irrigation_soil(ji) + reinfiltration_soil(ji) &
+          temp(ji) = MIN(water2infilt(ji,jst) + irrigation_soil(ji,jst) + reinfiltration_soil(ji) &
                          - MIN(ae_ns(ji,jst),zero) - MIN(subsinksoil(ji),zero) + precisol_ns(ji,jst), &
                            MAX(ae_ns(ji,jst),zero) + MAX(subsinksoil(ji),zero) )
@@ -3912,5 +4013,6 @@
        !   - eventually, water2infilt holds all fluxes to the soil surface except precisol (reduced by water2extract)
        DO ji = 1, kjpindex
-          water2infilt(ji,jst) = water2infilt(ji,jst) + irrigation_soil(ji) + reinfiltration_soil(ji) &
+         !Note that in tag 2.0, irrigation_soil(ji), changed to be coherent with new variable dimension
+          water2infilt(ji,jst) = water2infilt(ji,jst) + irrigation_soil(ji, jst) + reinfiltration_soil(ji) &
                 - MIN(ae_ns(ji,jst),zero) - MIN(subsinksoil(ji),zero) + precisol_ns(ji,jst) &
                 - temp(ji) 
@@ -3955,5 +4057,5 @@
        IF ( .NOT. doponds ) THEN ! this is the general case...
           DO ji = 1, kjpindex
-             water2infilt(ji,jst) = reinf_slope(ji) * ru_ns(ji,jst)
+             water2infilt(ji,jst) = reinf_slope_soil(ji,jst) * ru_ns(ji,jst)
           ENDDO
        ELSE
@@ -4341,5 +4443,47 @@
           tmcs_litter(:) = tmcs_litter(:) + sms(:,jsl)
        END DO
-             
+
+       ! Here we compute root zone deficit, to have an estimate of water demand in irrigated soil column (i.e. crop and grass)
+       IF(jst .EQ. irrig_st ) THEN
+           !It computes water deficit only on the root zone, and only on the layers where
+           !there is actually a deficit. If there is not deficit, it does not take into account that layer
+            DO ji = 1,kjpindex
+
+                root_deficit(ji) = SUM( MAX(zero, beta_irrig*smf(ji,1:nslm_root(ji) ) &
+                - sm(ji,1:nslm_root(ji)  ) )) - water2infilt(ji,jst)
+
+                root_deficit(ji) = MAX( root_deficit(ji) ,  zero)
+
+            ENDDO
+            !It COUNTS the number of pft with LAI > lai_irrig_min, inside the soiltile
+            !It compares veget, but it is the same as they are related by a function
+            lai_irrig_trig(:) = 0
+
+            DO jv = 1, nvm
+              IF( .NOT. natural(jv) ) THEN
+                DO ji = 1,kjpindex
+
+                  IF(veget(ji, jv) > veget_max(ji, jv) * ( un - &
+                  exp( - lai_irrig_min * ext_coeff_vegetfrac(jv) ) ) ) THEN
+
+                    lai_irrig_trig(ji) = lai_irrig_trig(ji) + 1
+
+                  ENDIF
+
+                ENDDO
+              ENDIF
+
+            ENDDO
+            !If any of the PFT inside the soil tile have LAI >  lai_irrig_min (I.E. lai_irrig_trig(ji) = 0 )
+            !The root deficit is set to zero, and irrigation is not triggered
+            DO ji = 1,kjpindex
+
+              IF( lai_irrig_trig(ji) < 1 ) THEN
+                root_deficit(ji) = zero
+              ENDIF
+
+            ENDDO
+       ENDIF
+
        ! Soil wetness profiles (W-Ww)/(Ws-Ww)
        ! soil_wet_ns is the ratio of available soil moisture to max available soil moisture
