Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90	(revision 7576)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90	(revision 7709)
@@ -12,6 +12,6 @@
 !!\n DESCRIPTION: None
 !!
-!! RECENT CHANGE(S): None
-!!
+!! RECENT CHANGE(S): July 2022: New irrigation scheme. Here new irrigation offer with more information from maps, and
+!!                    calculation of water withdrawal.
 !! REFERENCE(S)	:
 !!
@@ -226,7 +226,22 @@
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: reinfiltration_mean          !! Mean water flow which returns to the grid box (kg/m^2/dt)
 !$OMP THREADPRIVATE(reinfiltration_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigdeficit_mean              !! Mean irrigation deficit.
+                                                                                           !! This is between irrigation and water requirement (kg/m^2/dt)
+!$OMP THREADPRIVATE(irrigdeficit_mean)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigation_mean              !! Mean irrigation flux.
                                                                                              !! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt)
 !$OMP THREADPRIVATE(irrigation_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_gw_source_mean       !! Mean groundwater irrigation flux.
+                                                                                         !! This is the water taken from the GW reservoir only (kg/m^2/dt)
+!$OMP THREADPRIVATE(irrig_gw_source_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_fast_source_mean       !! Mean irrigation flux from fast.
+                                                                                       !! This is the water taken from the fast reservoir only (kg/m^2/dt)
+!$OMP THREADPRIVATE(irrig_fast_source_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_str_source_mean       !! Mean streamflow irrigation flux.
+                                                                                       !! This is the water taken from the streamflow reservoir only (kg/m^2/dt)
+!$OMP THREADPRIVATE(irrig_str_source_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigadduct_mean       !! Irrigation that comes from adduction. It includes water from basins inside the grid cell
+                                                                                        ! and water from nearby grid cells
+!$OMP THREADPRIVATE(irrigadduct_mean)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: riverflow_mean               !! Mean Outflow of the major rivers.
                                                                                              !! The flux will be located on the continental grid but this should be a coastal point (kg/dt)
@@ -291,6 +306,7 @@
                                 neighbours,  resolution,     contfrac,   stempdiag, &
                                 returnflow,  reinfiltration, irrigation, riverflow, &
-                                coastalflow, flood_frac,     flood_res )
-       
+                                coastalflow, flood_frac,     flood_res, soiltile, irrig_frac, &
+                                veget_max, irrigated_next) !
+
     IMPLICIT NONE
     
@@ -309,4 +325,8 @@
     REAL(r_std), INTENT(in)        :: contfrac(nbpt)       !! Fraction of land in each grid box (unitless;0-1)
     REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless)
+    REAL(r_std), INTENT(in)        :: veget_max(nbpt,nvm)  !! Maximal fraction of vegetation (unitless;0-1) !
+    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)  !! Dynamic irrig. area, calculated in slowproc and passed to routing!
+    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)      !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
 
     !! 0.2 Output variables
@@ -323,5 +343,5 @@
     !! 0.3 Local variables
     REAL(r_std), DIMENSION(nbp_glo):: mask_coast_glo       !! Mask with coastal gridcells on global grid (1/0)
-    LOGICAL                        :: init_irrig           !! Logical to initialize the irrigation (true/false)
+    !LOGICAL                        :: init_irrig          !! Logical to initialize the irrigation (true/false)
     LOGICAL                        :: init_flood           !! Logical to initialize the floodplains (true/false)
     LOGICAL                        :: init_swamp           !! Logical to initialize the swamps (true/false)
@@ -385,7 +405,6 @@
     !! Initialisation of flags for irrigated land, flood plains and swamps
     !
-    init_irrig = .FALSE.
     IF ( do_irrigation ) THEN 
-       IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE.
+       irrigated(:) = irrigated_next(:)
     END IF
     
@@ -399,49 +418,19 @@
        IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE.
     END IF
-       
+
     !! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree
     !! base data set to the resolution of the model.
-    
-    IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN
-       CALL routing_irrigmap(nbpt, index, lalo, neighbours, resolution, &
-            contfrac, init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)
+    !IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN
+    IF ( init_flood .OR. init_swamp ) THEN
+       CALL routing_floodmap(nbpt, index, lalo, neighbours, resolution, &
+            contfrac, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)
     ENDIF
-    
-    IF ( do_irrigation ) THEN 
-       CALL xios_orchidee_send_field("irrigmap",irrigated)
-       
-       IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) 
-       IF ( .NOT. almaoutput ) THEN
-          CALL histwrite_p(hist_id, 'irrigmap', 1, irrigated, nbpt, index)
-       ELSE
-          CALL histwrite_p(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index)
-       ENDIF
-       IF ( hist2_id > 0 ) THEN
-          IF ( .NOT. almaoutput ) THEN
-             CALL histwrite_p(hist2_id, 'irrigmap', 1, irrigated, nbpt, index)
-          ELSE
-             CALL histwrite_p(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index)
-          ENDIF
-       ENDIF
+
+    IF ( do_irrigation ) THEN
+
+       IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated)
+       IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrig_frac : ', MINVAL(irrig_frac), MAXVAL(irrig_frac)
     ENDIF
-    
-    IF ( do_floodplains ) THEN
-       CALL xios_orchidee_send_field("floodmap",floodplains)
-       
-       IF (printlev>=3) WRITE(numout,*) 'Verification : range of floodplains : ', MINVAL(floodplains), MAXVAL(floodplains) 
-       IF ( .NOT. almaoutput ) THEN
-          CALL histwrite_p(hist_id, 'floodmap', 1, floodplains, nbpt, index)
-       ELSE
-          CALL histwrite_p(hist_id, 'FloodplainsMap', 1, floodplains, nbpt, index)
-       ENDIF
-       IF ( hist2_id > 0 ) THEN
-          IF ( .NOT. almaoutput ) THEN
-             CALL histwrite_p(hist2_id, 'floodmap', 1, floodplains, nbpt, index)
-          ELSE
-             CALL histwrite_p(hist2_id, 'FloodplainsMap', 1, floodplains, nbpt, index)
-          ENDIF
-       ENDIF
-    ENDIF
-    
+
     IF ( doswamps ) THEN
        CALL xios_orchidee_send_field("swampmap",swamp)
@@ -558,5 +547,6 @@
        & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
        & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
-       & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id)
+       & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id,&
+       soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) !
 
     IMPLICIT NONE
@@ -584,4 +574,9 @@
     REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
     REAL(r_std), INTENT(in)        :: reinf_slope(nbpt)    !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: root_deficit(nbpt)   !! soil water deficit
+    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless)
+    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)     !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
+    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)!! Dynamic irrig. area, calculated in slowproc and passed to routing
+    REAL(r_std), INTENT(in)        :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac
 
     !! 0.2 Output variables
@@ -670,10 +665,28 @@
     ! the growing seasons we will give more weight to these areas.
     !
-    DO jv=2,nvm
-       DO ig=1,nbpt
-          humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing
-          vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
-       ENDDO
-    ENDDO
+    ! New irrigation scheme uses mean of vegtot with jv 1 to nvm
+    ! Old scheme keeps jv 2 to nvm, even if possibly an error
+    IF ( .NOT. old_irrig_scheme ) THEN
+      DO jv=1,nvm
+         DO ig=1,nbpt
+            humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing
+            vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
+         ENDDO
+      ENDDO
+    ELSE
+      DO jv=2,nvm
+         DO ig=1,nbpt
+            humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing
+            vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
+         ENDDO
+      ENDDO
+    ENDIF
+    ! Here updates irrigmap to irrigated_next from slowproc, every timestep
+    !irrigated_next is updated in slowproc when time comes
+    !irrig_frac was also updated in slowproc, here used as input variable
+
+    IF ( do_irrigation .AND. irrig_map_dynamic_flag  ) THEN
+        irrigated(:) = irrigated_next(:)
+    ENDIF
     !
     time_counter = time_counter + dt_sechiba 
@@ -685,9 +698,11 @@
        !! Computes the transport of water in the various reservoirs
        !
-       CALL routing_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, &
+       CALL routing_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, root_deficit, soiltile, &
             & vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, reinf_slope, &
-            & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrigation_mean, riverflow_mean, &
+            & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrig_frac, irrigation_mean, irrigdeficit_mean, &
+            & irrigadduct_mean, irrig_gw_source_mean, & !
+            & irrig_fast_source_mean, irrig_str_source_mean, riverflow_mean, & !
             & coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res, &
-            & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag)
+            & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, fraction_aeirrig_sw)
        !
        !! Responsible for storing the water in lakes
@@ -718,4 +733,9 @@
        reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba
        irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba
+       irrigdeficit_mean(:) = irrigdeficit_mean/dt_routing*dt_sechiba !
+       irrigadduct_mean(:) = irrigadduct_mean(:)/dt_routing*dt_sechiba!
+       irrig_gw_source_mean(:) = irrig_gw_source_mean(:)/dt_routing*dt_sechiba !
+       irrig_fast_source_mean(:) = irrig_fast_source_mean(:)/dt_routing*dt_sechiba !
+       irrig_str_source_mean(:) = irrig_str_source_mean(:)/dt_routing*dt_sechiba !
        irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba
        
@@ -766,4 +786,9 @@
     ! Water fluxes converted from kg/m^2/dt_sechiba into kg/m^2/s 
     CALL xios_orchidee_send_field("irrigation",irrigation/dt_sechiba)
+    CALL xios_orchidee_send_field("irrig_deficit",irrigdeficit_mean/dt_sechiba)!
+    CALL xios_orchidee_send_field("irrig_adduct",irrigadduct_mean/dt_sechiba)!
+    CALL xios_orchidee_send_field("irrig_gw_source",irrig_gw_source_mean/dt_sechiba) !
+    CALL xios_orchidee_send_field("irrig_fast_source",irrig_fast_source_mean/dt_sechiba) !
+    CALL xios_orchidee_send_field("irrig_str_source",irrig_str_source_mean/dt_sechiba) !
     CALL xios_orchidee_send_field("netirrig",irrig_netereq/dt_sechiba)
     CALL xios_orchidee_send_field("riversret",returnflow/dt_sechiba)
@@ -879,5 +904,5 @@
     
     !! 0.2 Local variables
-    REAL(r_std), DIMENSION(1)      :: tmp_day              
+    REAL(r_std), DIMENSION(1)      :: tmp_day
 
 !_ ================================================================================================================================
@@ -931,6 +956,11 @@
 
     IF ( do_irrigation ) THEN
-       CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter',  nbp_glo, index_g)
+       !CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter',  nbp_glo, index_g)
        CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter',  nbp_glo, index_g)
+       CALL restput_p (rest_id, 'irrigdeficit', nbp_glo, 1, 1, kjit, irrigdeficit_mean, 'scatter',  nbp_glo, index_g)
+       CALL restput_p (rest_id, 'irrigadduct', nbp_glo, 1, 1, kjit, irrigadduct_mean, 'scatter',  nbp_glo, index_g)
+       CALL restput_p (rest_id, 'irrig_gw_source', nbp_glo, 1, 1, kjit, irrig_gw_source_mean, 'scatter',  nbp_glo, index_g) !
+       CALL restput_p (rest_id, 'irrig_fast_source', nbp_glo, 1, 1, kjit, irrig_fast_source_mean, 'scatter',  nbp_glo, index_g) !
+       CALL restput_p (rest_id, 'irrig_str_source', nbp_glo, 1, 1, kjit, irrig_str_source_mean, 'scatter',  nbp_glo, index_g) !
     ENDIF
 
@@ -1022,10 +1052,10 @@
     !
     !Config Key   = DO_FLOODINFILT
-    !Config Desc  = Should floodplains reinfiltrate into the soil 
+    !Config Desc  = Should floodplains reinfiltrate into the soil
     !Config If    = RIVER_ROUTING
     !Config Def   = n
     !Config Help  = This parameters allows the user to ask the model
-    !Config         to take into account the flood plains reinfiltration 
-    !Config         into the soil moisture. It then can go 
+    !Config         to take into account the flood plains reinfiltration
+    !Config         into the soil moisture. It then can go
     !Config         back to the slow and fast reservoirs
     !Config Units = [FLAG]
@@ -1039,7 +1069,7 @@
     !Config Def   = n
     !Config Help  = This parameters allows the user to ask the model
-    !Config         to take into account the swamps and return 
-    !Config         the water into the bottom of the soil. It then can go 
-    !Config         back to the atmopshere. This tried to simulate 
+    !Config         to take into account the swamps and return
+    !Config         the water into the bottom of the soil. It then can go
+    !Config         back to the atmopshere. This tried to simulate
     !Config         internal deltas of rivers.
     !Config Units = [FLAG]
@@ -1049,11 +1079,11 @@
     !
     !Config Key   = DO_PONDS
-    !Config Desc  = Should we include ponds 
+    !Config Desc  = Should we include ponds
     !Config If    = RIVER_ROUTING
     !Config Def   = n
     !Config Help  = This parameters allows the user to ask the model
-    !Config         to take into account the ponds and return 
-    !Config         the water into the soil moisture. It then can go 
-    !Config         back to the atmopshere. This tried to simulate 
+    !Config         to take into account the ponds and return
+    !Config         the water into the soil moisture. It then can go
+    !Config         back to the atmopshere. This tried to simulate
     !Config         little ponds especially in West Africa.
     !Config Units = [FLAG]
@@ -1064,10 +1094,10 @@
 
     !Config Key   = SLOW_TCST
-    !Config Desc  = Time constant for the slow reservoir 
-    !Config If    = RIVER_ROUTING 
-    !Config Def   = 25.0 
-    !Config Help  = This parameters allows the user to fix the 
+    !Config Desc  = Time constant for the slow reservoir
+    !Config If    = RIVER_ROUTING
+    !Config Def   = 25.0
+    !Config Help  = This parameters allows the user to fix the
     !Config         time constant (in days) of the slow reservoir
-    !Config         in order to get better river flows for 
+    !Config         in order to get better river flows for
     !Config         particular regions.
     !Config Units = [days]
@@ -1078,5 +1108,5 @@
 !> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and
 !> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir"
-!> have the highest value in order to simulate the groundwater. 
+!> have the highest value in order to simulate the groundwater.
 !> The "stream reservoir", which represents all the water of the stream, has the lowest value.
 !> Those figures are the same for all the basins of the world.
@@ -1087,53 +1117,53 @@
     !
     !Config Key   = FAST_TCST
-    !Config Desc  = Time constant for the fast reservoir 
-    !Config If    = RIVER_ROUTING 
-    !Config Def   = 3.0 
-    !Config Help  = This parameters allows the user to fix the 
+    !Config Desc  = Time constant for the fast reservoir
+    !Config If    = RIVER_ROUTING
+    !Config Def   = 3.0
+    !Config Help  = This parameters allows the user to fix the
     !Config         time constant (in days) of the fast reservoir
-    !Config         in order to get better river flows for 
+    !Config         in order to get better river flows for
     !Config         particular regions.
     !Config Units = [days]
     CALL getin_p('FAST_TCST', fast_tcst)
-    
+
     !Config Key   = STREAM_TCST
-    !Config Desc  = Time constant for the stream reservoir 
+    !Config Desc  = Time constant for the stream reservoir
     !Config If    = RIVER_ROUTING
     !Config Def   = 0.24
-    !Config Help  = This parameters allows the user to fix the 
+    !Config Help  = This parameters allows the user to fix the
     !Config         time constant (in days) of the stream reservoir
-    !Config         in order to get better river flows for 
+    !Config         in order to get better river flows for
     !Config         particular regions.
     !Config Units = [days]
     CALL getin_p('STREAM_TCST', stream_tcst)
-    
+
     !Config Key   = FLOOD_TCST
-    !Config Desc  = Time constant for the flood reservoir 
+    !Config Desc  = Time constant for the flood reservoir
     !Config If    = RIVER_ROUTING
     !Config Def   = 4.0
-    !Config Help  = This parameters allows the user to fix the 
+    !Config Help  = This parameters allows the user to fix the
     !Config         time constant (in days) of the flood reservoir
-    !Config         in order to get better river flows for 
+    !Config         in order to get better river flows for
     !Config         particular regions.
     !Config Units = [days]
     CALL getin_p('FLOOD_TCST', flood_tcst)
-    
+
     !Config Key   = SWAMP_CST
-    !Config Desc  = Fraction of the river that flows back to swamps 
+    !Config Desc  = Fraction of the river that flows back to swamps
     !Config If    = RIVER_ROUTING
     !Config Def   = 0.2
-    !Config Help  = This parameters allows the user to fix the 
+    !Config Help  = This parameters allows the user to fix the
     !Config         fraction of the river transport
     !Config         that flows to swamps
     !Config Units = [-]
     CALL getin_p('SWAMP_CST', swamp_cst)
-    
+
     !Config Key   = FLOOD_BETA
-    !Config Desc  = Parameter to fix the shape of the floodplain  
+    !Config Desc  = Parameter to fix the shape of the floodplain
     !Config If    = RIVER_ROUTING
     !Config Def   = 2.0
     !Config Help  = Parameter to fix the shape of the floodplain
     !Config         (>1 for convex edges, <1 for concave edges)
-    !Config Units = [-] 
+    !Config Units = [-]
     CALL getin_p("FLOOD_BETA", beta)
     !
@@ -1142,7 +1172,7 @@
     !Config If    = RIVER_ROUTING
     !Config Def   = 0.5
-    !Config Help  = 
-    !Config Units = [-] 
-    CALL getin_p("POND_BETAP", betap)    
+    !Config Help  =
+    !Config Units = [-]
+    CALL getin_p("POND_BETAP", betap)
     !
     !Config Key   = FLOOD_CRI
@@ -1150,6 +1180,6 @@
     !Config If    = DO_FLOODPLAINS or DO_PONDS
     !Config Def   = 2000.
-    !Config Help  = 
-    !Config Units = [mm] 
+    !Config Help  =
+    !Config Units = [mm]
     CALL getin_p("FLOOD_CRI", floodcri)
     !
@@ -1158,6 +1188,6 @@
     !Config If    = DO_FLOODPLAINS or DO_PONDS
     !Config Def   = 2000.
-    !Config Help  = 
-    !Config Units = [mm] 
+    !Config Help  =
+    !Config Units = [mm]
     CALL getin_p("POND_CRI", pondcri)
 
@@ -1166,6 +1196,6 @@
     !Config If    = RIVER_ROUTING
     !Config Def   = 7000
-    !Config Help  = 
-    !Config Units = [kg/m2(routing area)] 
+    !Config Help  =
+    !Config Units = [kg/m2(routing area)]
     max_lake_reservoir = 7000
     CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir)
@@ -1205,10 +1235,10 @@
        ELSE
           ! Take the value from restart file
-          time_counter = tmp_day(1) 
+          time_counter = tmp_day(1)
        ENDIF
     ENDIF
     CALL bcast(time_counter)
 
-    
+
     ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for routing_area_loc','','')
@@ -1365,5 +1395,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g)
     CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE (pond_frac(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for pond_frac','','')
@@ -1373,5 +1403,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g)
     CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero)
-    
+
     var_name = 'flood_frac'
     CALL ioconf_setatt_p('UNITS', '-')
@@ -1379,5 +1409,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g)
     CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero)
-    
+
     var_name = 'flood_res'
     CALL ioconf_setatt_p('UNITS','mm')
@@ -1386,5 +1416,5 @@
     CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero)
 !    flood_res = zero
-    
+
     ALLOCATE (lake_reservoir(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for lake_reservoir','','')
@@ -1394,5 +1424,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g)
     CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE (pond_reservoir(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for pond_reservoir','','')
@@ -1404,15 +1434,13 @@
     !
     ! Map of irrigated areas
-    !
-    IF ( do_irrigation ) THEN
-       ALLOCATE (irrigated(nbpt), stat=ier)
-       IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigated','','')
-       var_name = 'irrigated'
-       CALL ioconf_setatt_p('UNITS', 'm^2')
-       CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area')
-       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g)
-       CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba)
-    ENDIF
-    
+    ! irrigated equal to irrigated_next from slowproc. Here, we initialize to zero.
+    ! Values from slowproc given in routing_main
+    ALLOCATE (irrigated(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigated','','')
+    var_name = 'irrigated'
+    CALL ioconf_setatt_p('UNITS', 'm^2')
+    CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area')
+    irrigated(:) = zero
+
     IF ( do_floodplains ) THEN
        ALLOCATE (floodplains(nbpt), stat=ier)
@@ -1443,5 +1471,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE (returnflow_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for returnflow_mean','','')
@@ -1452,5 +1480,5 @@
     CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero)
     returnflow(:) = returnflow_mean(:)
-    
+
     ALLOCATE (reinfiltration_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for reinfiltration_mean','','')
@@ -1461,11 +1489,21 @@
     CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero)
     reinfiltration(:) = reinfiltration_mean(:)
-    
+
     ALLOCATE (irrigation_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigation_mean','','')
+    ALLOCATE (irrigdeficit_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigdeficit_mean','','')
+    ALLOCATE (irrigadduct_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigadduct_mean','','')
     ALLOCATE (irrig_netereq(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_netereq','','')
+    ALLOCATE (irrig_gw_source_mean(nbpt), stat=ier) !
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_gw_source_mean','','')!
+    ALLOCATE (irrig_fast_source_mean(nbpt), stat=ier) !
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_fast_source_mean','','')!
+    ALLOCATE (irrig_str_source_mean(nbpt), stat=ier) !
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_str_source_mean','','')!
     irrig_netereq(:) = zero
-    
+
     IF ( do_irrigation ) THEN
        var_name = 'irrigation'
@@ -1474,9 +1512,39 @@
        CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g)
        CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero)
+       var_name = 'irrigdeficit'
+       CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+       CALL ioconf_setatt_p('LONG_NAME','Irrigation deficit')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigdeficit_mean, "gather", nbp_glo, index_g)
+       CALL setvar_p (irrigdeficit_mean, val_exp, 'NO_KEYWORD', zero)
+       var_name = 'irrigadduct'
+       CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+       CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux from adduction')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigadduct_mean, "gather", nbp_glo, index_g)
+       CALL setvar_p (irrigadduct_mean, val_exp, 'NO_KEYWORD', zero)
+       var_name = 'irrig_gw_source'!
+       CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+       CALL ioconf_setatt_p('LONG_NAME','Irrigation from GW reservoir')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_gw_source_mean, "gather", nbp_glo, index_g)
+       CALL setvar_p (irrig_gw_source_mean, val_exp, 'NO_KEYWORD', zero)
+       var_name = 'irrig_fast_source'!
+       CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+       CALL ioconf_setatt_p('LONG_NAME','Irrigation from fast reservoir')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_fast_source_mean, "gather", nbp_glo, index_g)
+       CALL setvar_p (irrig_fast_source_mean, val_exp, 'NO_KEYWORD', zero)
+       var_name = 'irrig_str_source'!
+       CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+       CALL ioconf_setatt_p('LONG_NAME','Irrigation from stream reservoir')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_str_source_mean, "gather", nbp_glo, index_g)
+       CALL setvar_p (irrig_str_source_mean, val_exp, 'NO_KEYWORD', zero)
     ELSE
        irrigation_mean(:) = zero
+       irrig_gw_source_mean(:) = zero
+       irrig_fast_source_mean(:) = zero
+       irrig_str_source_mean(:) = zero
+       irrigdeficit_mean(:) = zero
+       irrigadduct_mean(:) = zero
     ENDIF
-    irrigation(:) = irrigation_mean(:) 
-    
+    irrigation(:) = irrigation_mean(:)
+
     ALLOCATE (riverflow_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for riverflow_mean','','')
@@ -1487,5 +1555,5 @@
     CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero)
     riverflow(:) = riverflow_mean(:)
-    
+
     ALLOCATE (coastalflow_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for coastalflow_mean','','')
@@ -1496,5 +1564,5 @@
     CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
     coastalflow(:) = coastalflow_mean(:)
-    
+
     ! Locate it at the 2m level
     ipn = MINLOC(ABS(diaglev-2))
@@ -1503,5 +1571,5 @@
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for floodtemp','','')
     floodtemp(:) = stempdiag(:,floodtemp_lev)
-    
+
     ALLOCATE(hydrographs(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for hydrographs','','')
@@ -1511,5 +1579,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g)
     CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero)
- 
+
     ALLOCATE(slowflow_diag(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for slowflow_diag','','')
@@ -1526,5 +1594,5 @@
          & pond_diag(nbpt), lake_diag(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for fast_diag,..','','')
-    
+
     fast_diag(:) = zero
     slow_diag(:) = zero
@@ -1533,5 +1601,5 @@
     pond_diag(:) = zero
     lake_diag(:) = zero
-    
+
     DO ig=1,nbpt
        totarea = zero
@@ -1566,5 +1634,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE (runoff_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for runoff_mean','','')
@@ -1574,5 +1642,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE(drainage_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for drainage_mean','','')
@@ -1582,5 +1650,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE(transpot_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for transpot_mean','','')
@@ -1598,5 +1666,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE(humrel_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for humrel_mean','','')
@@ -1606,5 +1674,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un)
-    
+
     ALLOCATE(k_litt_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for k_litt_mean','','')
@@ -1614,5 +1682,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE(totnobio_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for totnobio_mean','','')
@@ -1622,5 +1690,5 @@
     CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g)
     CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero)
-    
+
     ALLOCATE(vegtot_mean(nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for vegtot_mean','','')
@@ -1693,9 +1761,14 @@
     IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc)
     IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo)
-    IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc)    
+    IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc)
     IF (ALLOCATED(hydroupbasin_glo)) DEALLOCATE(hydroupbasin_glo)
     IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs)
     IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag)
     IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean)
+    IF (ALLOCATED(irrigdeficit_mean)) DEALLOCATE(irrigdeficit_mean)!
+    IF (ALLOCATED(irrigadduct_mean)) DEALLOCATE(irrigadduct_mean)!
+    IF (ALLOCATED(irrig_gw_source_mean)) DEALLOCATE(irrig_gw_source_mean) !
+    IF (ALLOCATED(irrig_fast_source_mean)) DEALLOCATE(irrig_fast_source_mean) !
+    IF (ALLOCATED(irrig_str_source_mean)) DEALLOCATE(irrig_str_source_mean) !
     IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated)
     IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains)
@@ -1718,9 +1791,9 @@
 !!
 !! DESCRIPTION (definitions, functional, design, flags) :
-!! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an 
-!! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes 
-!! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of 
+!! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an
+!! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes
+!! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of
 !! the up-stream basin and adding it to the down-stream one.
-!! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once 
+!! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once
 !! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them
 !! the hydrographs of the largest rivers we have chosen to monitor.
@@ -1778,9 +1851,11 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE routing_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, &
+  SUBROUTINE routing_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, root_deficit, soiltile, &
        &                  vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, &
-       &                  lakeinflow, returnflow, reinfiltration, irrigation, riverflow, &
+       &                  lakeinflow, returnflow, reinfiltration, irrig_frac, irrigation, irrigdeficit, &
+       &                  irrigadduct, irrig_gw_source, &
+       &                  irrig_fast_source, irrig_str_source, riverflow, & !
        &                  coastalflow, hydrographs, slowflow_diag, flood_frac, flood_res, &
-                          netflow_stream_diag, netflow_fast_diag, netflow_slow_diag)
+                          netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, fraction_aeirrig_sw) !
     !
     IMPLICIT NONE
@@ -1793,4 +1868,5 @@
     REAL(r_std), INTENT(in)                      :: floodout(nbpt)            !! Grid-point flow out of floodplains (kg/m^2/dt)
     REAL(r_std), INTENT(in)                      :: drainage(nbpt)            !! Grid-point drainage (kg/m^2/dt)
+    REAL(r_std), INTENT(in)                      :: root_deficit(nbpt)        !! soil water deficit
     REAL(r_std), INTENT(in)                      :: vegtot(nbpt)              !! Potentially vegetated fraction (unitless;0-1)
     REAL(r_std), INTENT(in)                      :: totnobio(nbpt)            !! Other areas which can not have vegetation
@@ -1802,4 +1878,8 @@
     REAL(r_std), INTENT(in)                      :: reinf_slope(nbpt)         !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
     REAL(r_std), INTENT(out)                     :: lakeinflow(nbpt)          !! Water inflow to the lakes (kg/dt)
+    REAL(r_std), INTENT(in)                      :: soiltile(nbpt,nstm)       !! Fraction of each soil tile within vegtot (0-1, unitless)
+    REAL(r_std), INTENT(in)                      :: irrig_frac(nbpt)          !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
+    REAL(r_std), INTENT(in)                      :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac
+
     !
 !! OUTPUT VARIABLES
@@ -1808,4 +1888,7 @@
     REAL(r_std), INTENT(out)                     :: reinfiltration(nbpt)      !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
     REAL(r_std), INTENT(out)                     :: irrigation(nbpt)          !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt_routing)
+    REAL(r_std), INTENT(out)                     :: irrigdeficit(nbpt)          !! Irrigation deficit. Difference btw irrig. demand and irrigation
+    REAL(r_std), INTENT(out)                     :: irrigadduct(nbpt)          !! Irrigation from adducted water. Included in irrigation
+
     REAL(r_std), INTENT(out)                     :: riverflow(nbpt)           !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt_routing)
     REAL(r_std), INTENT(out)                     :: coastalflow(nbpt)         !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt_routing)
@@ -1818,4 +1901,8 @@
     REAL(r_std), INTENT(out)                     :: netflow_fast_diag(nbpt)   !! Input - Output flow to fast reservoir
     REAL(r_std), INTENT(out)                     :: netflow_slow_diag(nbpt)   !! Input - Output flow to slow reservoir
+
+    REAL(r_std), INTENT(out)                     :: irrig_gw_source(nbpt)     !! Irrigation from  groundwater per cell
+    REAL(r_std), INTENT(out)                     :: irrig_fast_source(nbpt)     !! Irrigation from  fast reservoir per cell
+    REAL(r_std), INTENT(out)                     :: irrig_str_source(nbpt)     !! Irrigation from  stramflow per cell
     !
 !! LOCAL VARIABLES
@@ -1843,8 +1930,8 @@
     REAL(r_std), DIMENSION(nbpt)                 :: totarea                   !! Total area of basin (m^2)
     REAL(r_std), DIMENSION(nbpt)                 :: totflood                  !! Total amount of water in the floodplains reservoir (kg)
-    REAL(r_std), DIMENSION(nbasmax)              :: pond_excessflow           !! 
+    REAL(r_std), DIMENSION(nbasmax)              :: pond_excessflow           !!
     REAL(r_std)                                  :: flow                      !! Outflow computation for the reservoirs (kg/dt)
     REAL(r_std)                                  :: floodindex                !! Fraction of grid box area inundated (unitless;0-1)
-    REAL(r_std)                                  :: pondex                    !! 
+    REAL(r_std)                                  :: pondex                    !!
     REAL(r_std)                                  :: flood_frac_pot            !! Total fraction of the grid box which is flooded at optimum repartition (unitless;0-1)
     REAL(r_std)                                  :: stream_tot                !! Total water amount in the stream reservoirs (kg)
@@ -1855,10 +1942,16 @@
     REAL(r_std)                                  :: total_lake_overflow       !! Sum of lake_overflow over full grid (kg)
     REAL(r_std), DIMENSION(8,nbasmax)            :: streams_around            !! Stream reservoirs of the neighboring grid boxes (kg)
-    INTEGER(i_std), DIMENSION(8)                 :: igrd                      !! 
-    INTEGER(i_std), DIMENSION(2)                 :: ff                        !! 
-    INTEGER(i_std), DIMENSION(1)                 :: fi                        !! 
+    INTEGER(i_std), DIMENSION(8)                 :: igrd                      !!
+    INTEGER(i_std), DIMENSION(2)                 :: ff                        !!
+    INTEGER(i_std), DIMENSION(1)                 :: fi                        !!
     INTEGER(i_std)                               :: ig, ib, ib2, ig2          !! Indices (unitless)
     INTEGER(i_std)                               :: rtg, rtb, in              !! Indices (unitless)
-    INTEGER(i_std)                               :: ier                       !! Error handling 
+    INTEGER(i_std)                               :: ier                       !! Error handling
+    REAL(r_std), DIMENSION(nbpt)              :: Count_failure_slow        !! Counter times slow reserv. does not fit irrigation demand
+    REAL(r_std), DIMENSION(nbpt)              :: Count_failure_fast        !! Counter times fast reserv. does not fit irrigation demand
+    REAL(r_std), DIMENSION(nbpt)              :: Count_failure_stre        !! Counter times stream reserv. does not fit irrigation demand
+    LOGICAL              :: IsFail_slow        !! Logical to ask if slow reserv. does not fit irrigation demand
+    LOGICAL              :: IsFail_fast        !! Logical to ask if fast reserv. does not fit irrigation demand
+    LOGICAL              :: IsFail_stre        !! Logical to ask if stream reserv. does not fit irrigation demand
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: fast_flow_g               !! Outflow from the fast reservoir (kg/dt)
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: slow_flow_g               !! Outflow from the slow reservoir (kg/dt)
@@ -1872,5 +1965,7 @@
     REAL(r_std), DIMENSION(nbpt, nbasmax)        :: netflow_slow              !! Input - Output flow to slow reservoir
 
-
+    REAL(r_std)                :: pcent_vol_irrig               !! Percentage of surface water for irrigation from total available water
+    REAL(r_std)                :: slow_wdr_dummy,fast_wdr_dummy,stre_wdr_dummy !! Dummy variables, real abstraction in each reservoir for irrigation
+    REAL(r_std)                :: pot_slow_wdr_dummy,pot_fast_wdr_dummy,pot_stre_wdr_dummy !! Dummy variables, potential abstraction in each reservoir for irrigatio
     !! PARAMETERS
     LOGICAL, PARAMETER                           :: check_reservoir = .FALSE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false)
@@ -1886,4 +1981,10 @@
     totarea(:) = zero
     totflood(:) = zero
+    irrig_gw_source(:) = zero !
+    irrig_fast_source(:) = zero !
+    irrig_str_source(:) = zero !
+    Count_failure_slow(:) = zero !
+    Count_failure_fast(:) = zero !
+    Count_failure_stre(:) = zero !
     !
     ! Compute all the fluxes
@@ -1897,5 +1998,5 @@
     ENDDO
           !
-!> The outflow fluxes from the three reservoirs are computed. 
+!> The outflow fluxes from the three reservoirs are computed.
 !> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume.
 !> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid
@@ -1912,5 +2013,5 @@
           IF ( route_tobasin(ig,ib) .GT. 0 ) THEN
              !
-             ! Each of the fluxes is limited by the water in the reservoir and a small margin 
+             ! Each of the fluxes is limited by the water in the reservoir and a small margin
              ! (min_reservoir) to avoid rounding errors.
              !
@@ -1923,5 +2024,5 @@
              slow_flow(ig,ib) = MAX(flow, zero)
 
-             flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst* & 
+             flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst* &
                   & MAX(un-SQRT(flood_frac_bas(ig,ib)),min_sechiba)*one_day/dt_routing),&
                   & stream_reservoir(ig,ib)-min_sechiba)
@@ -1944,5 +2045,5 @@
              flow = MIN(floodout(ig)*totarea(ig)*pond_frac(ig)/flood_frac(ig), pond_reservoir(ig)+totflood(ig))
              pondex = MAX(flow - pond_reservoir(ig), zero)
-             pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex) 
+             pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex)
              !
              ! If demand was over reservoir size, we will take it out from floodplains
@@ -1981,5 +2082,5 @@
     !-
     !- Computing the drainage and outflow from floodplains
-!> Drainage from floodplains is depending on a averaged conductivity (k_litt) 
+!> Drainage from floodplains is depending on a averaged conductivity (k_litt)
 !> for saturated infiltration in the 'litter' layer. Flood_drainage will be
 !> a component of the total reinfiltration that leaves the routing scheme.
@@ -1997,5 +2098,5 @@
           DO ib=1,nbasmax
              DO ig=1,nbpt
-                flood_drainage(ig,ib) = zero 
+                flood_drainage(ig,ib) = zero
              ENDDO
           ENDDO
@@ -2045,5 +2146,5 @@
              pond_drainage(ig,ib) = MIN(pond_reservoir(ig)*routing_area(ig,ib)/totarea(ig), &
                   & pond_frac(ig)*routing_area(ig,ib)*k_litt(ig)*dt_routing/one_day)
-             fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib) 
+             fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib)
           ENDDO
        ENDDO
@@ -2073,5 +2174,5 @@
     ENDIF
     IF (ier /= 0) CALL ipslerr_p(3,'routing_flow','Pb in allocate for fast_flow_g','','')
-       
+
     CALL gather(fast_flow,fast_flow_g)
     CALL gather(slow_flow,slow_flow_g)
@@ -2092,5 +2193,5 @@
 
     DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g )
-   
+
     CALL scatter(transport_glo,transport)
 
@@ -2112,5 +2213,5 @@
        DO ib=1,nbasmax
           DO ig=1,nbpt
-             potflood(ig,ib) = transport(ig,ib) 
+             potflood(ig,ib) = transport(ig,ib)
              !
              IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > 0. .AND. floodtemp(ig) > tp_00 ) THEN
@@ -2123,5 +2224,5 @@
                 return_swamp(ig,ib) = swamp_cst * potflood(ig,ib) * floodindex
                 !
-                tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib) 
+                tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib)
                 !
              ENDIF
@@ -2136,5 +2237,5 @@
           IF (floodplains(ig) .GT. min_sechiba .AND. floodtemp(ig) .GT. tp_00) THEN
              DO ib=1,nbasmax
-                floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib) 
+                floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib)
              ENDDO
           ENDIF
@@ -2162,5 +2263,5 @@
           !
           flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + floods(ig,ib) - &
-               & flood_flow(ig,ib) 
+               & flood_flow(ig,ib)
           !
           pond_reservoir(ig) = pond_reservoir(ig) + pond_inflow(ig,ib) - pond_drainage(ig,ib)
@@ -2170,5 +2271,5 @@
                 WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ib, ". Problem is being corrected."
                 WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir(ig,ib), floods(ig,ib), &
-                     & flood_flow(ig,ib) 
+                     & flood_flow(ig,ib)
              ENDIF
              stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_reservoir(ig,ib)
@@ -2229,5 +2330,5 @@
              flood_frac_pot = (totflood(ig) / (totarea(ig)*floodcri/(beta+un)))**(beta/(beta+un))
              flood_frac(ig) = MIN(floodplains(ig) / totarea(ig), flood_frac_pot)
-             ! Then we diagnose the fraction for each basin with the size of its flood_reservoir 
+             ! Then we diagnose the fraction for each basin with the size of its flood_reservoir
              ! (flood_frac_bas may be > 1)
              DO ib=1,nbasmax
@@ -2238,7 +2339,7 @@
              ENDDO
              ! We diagnose the maximum height of floodplain
-             flood_height(ig) = (beta/(beta+1))*floodcri*(flood_frac(ig))**(un/beta) + totflood(ig)/(totarea(ig)*flood_frac(ig)) 
+             flood_height(ig) = (beta/(beta+1))*floodcri*(flood_frac(ig))**(un/beta) + totflood(ig)/(totarea(ig)*flood_frac(ig))
              ! And finally add the pond surface
-             pond_frac(ig) = MIN(un-flood_frac(ig), ((betap+1)*pond_reservoir(ig) / (pondcri*totarea(ig)))**(betap/(betap+1)) ) 
+             pond_frac(ig) = MIN(un-flood_frac(ig), ((betap+1)*pond_reservoir(ig) / (pondcri*totarea(ig)))**(betap/(betap+1)) )
              flood_frac(ig) = flood_frac(ig) + pond_frac(ig)
              !
@@ -2265,5 +2366,5 @@
           DO ig=1,nbpt
              returnflow(ig) =  returnflow(ig) + return_swamp(ig,ib)
-             reinfiltration(ig) =  reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib) 
+             reinfiltration(ig) =  reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib)
           ENDDO
        ENDDO
@@ -2304,34 +2405,237 @@
 !
     IF ( do_irrigation ) THEN
-       DO ig=1,nbpt
-          !
-          IF ((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. &
-               & (runoff(ig) .LT. min_sechiba) ) THEN
-             
-             irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - &
+      DO ig=1,nbpt
+            !It enters to the new irrigation module only if there is an irrigated fraction, if not irrig_netereq = zero for that cell
+            IF ((irrig_frac(ig) .GT. min_sechiba) .AND. .NOT. old_irrig_scheme ) THEN
+
+                  irrig_netereq(ig) = irrig_netereq(ig) + MIN( irrig_dosmax/3600*dt_routing, &
+                                      root_deficit(ig) ) * soiltile(ig, irrig_st) * vegtot(ig)
+                  ! By definition, irrig_dosmax is in kg/m2 of soil tile/hour,dividing by 3600(seconds/hour) * DT_ROUTING  !
+                  ! = kg/m2 of soil tile/(routing timestep)
+                  ! irrig_netereq(kg/m2 of grid cell / routing timstep ) is equal to
+                  ! root_deficit (kg/m2 of soil tile) * soiltile*vegtot (fraction of soil tile at cell level) = kg/m2 of grid cell
+
+                  IF (.NOT. irrigated_soiltile .AND. ( soiltile(ig,irrig_st) .GT. min_sechiba ) .AND. (vegtot(ig) .GT. min_sechiba) ) THEN
+                      ! Irrigated_soiltile asks if there is an independent soil tile for irrigated crops. If not,
+                      ! actual volume calculated for irrig_netereq assumed that the whole SOILTILE was irrigated, but in this case
+                      ! just a fraction of the irrig_st (irrigated soil tile, by default = 3) is actually irrigated,
+                      ! and irrig_netereq needs to be reduced by irrig_frac/( soiltile * vegtot ) (note that it is max = 1 thanks to irrig_frac calculation in l. 424)
+                      ! Demand(ST3)*irrig_frac/(soiltile(3)*vegtot) = irrig_netereq_In_ST3, then
+                      !irrig_netereq_In_ST3 * (soiltile(3)*vegtot) = irrig_netereq at grid scale = Demand(ST3)*irrig_frac.
+                      irrig_netereq(ig) = irrig_netereq(ig) * irrig_frac(ig) / ( soiltile(ig,irrig_st) * vegtot(ig) )
+                      !irrig_netereq = kg/m2 of grid cell
+
+                  ENDIF
+            !Old irrigation scheme as in tag 2.0
+            ELSE IF((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. &
+                    & (runoff(ig) .LT. min_sechiba) .AND.  old_irrig_scheme) THEN
+
+                  irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - &
                   & (precip(ig)+reinfiltration(ig)) )
-             
-          ENDIF
-          !
-          DO ib=1,nbasmax
-             IF ( routing_area(ig,ib) .GT. 0 ) THEN
-             
-                irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib)
-
-                irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),&
-                     &   stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) )
-                
-                slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + &
-                     & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))))
-
-                fast_reservoir(ig,ib) = MAX( zero, &
-                     &  fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))
-
-                stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) )
-
-                irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib)
-
-             ENDIF
-          ENDDO
+
+            ENDIF
+
+            DO ib=1,nbasmax
+                IF ( routing_area(ig,ib) .GT. 0 ) THEN
+
+                    IF (.NOT. old_irrig_scheme .AND. select_source_irrig) THEN
+
+                      ! For   irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
+                      ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
+
+                      ! The new priorization scheme takes into account irrig. infrastructur according to GMIA map
+                      ! It also withdraw water according to availability, it means that it wont seek for all the water in then
+                      ! stream reservoir, even if this one could respond to the demand by itself
+
+                      pot_slow_wdr_dummy = ( 1 - fraction_aeirrig_sw(ig)) * avail_reserve(3)*slow_reservoir(ig,ib)
+                      pot_fast_wdr_dummy = fraction_aeirrig_sw(ig) * avail_reserve(2)*fast_reservoir(ig,ib)
+                      pot_stre_wdr_dummy = fraction_aeirrig_sw(ig) * avail_reserve(1) * stream_reservoir(ig,ib)
+                      pcent_vol_irrig = zero
+                      IsFail_slow = .FALSE. !
+                      IsFail_fast = .FALSE. !
+                      IsFail_stre = .FALSE. !
+                      irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib)
+
+                      irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),&
+                            pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)
+
+                       !!   additional IF to calculate pcent_vol_irrig, in the case the total avail.
+                       !! water is zero, I.E. when there is no water in surface and fraction_ae = 0,
+                       !! so GW is not taken into account
+                       !! Note on pcent_vol_irrig: It correspond to the fraction of available water in surface,
+                       !! considering environmental needs and irrigation equipement by source from map. It controls
+                       !! how the source of water withdrawl, especially when requirements < available water
+
+                      IF (  (pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)  .GT. min_sechiba ) THEN
+
+                          pcent_vol_irrig = ( pot_stre_wdr_dummy + pot_fast_wdr_dummy ) / &
+                                ( pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
+
+                          !Irrig_actual set to zero, because there is no available water.
+                          !Put to avoid negative values due to problems in the Min function
+                          irrig_actual(ig,ib) = MAX(irrig_actual(ig,ib), zero)
+                      !Already zero because pcent_vol_irrig initialized to zero
+                      !Put here to readability but not necessary
+                      !ELSE
+                      !    pcent_vol_irrig = zero
+
+                      ENDIF
+
+                      !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
+                      !Volume used for irrigation that comes from GW
+                      ! Idem for irrig_fast_source and irrig_str_source
+
+                      slow_wdr_dummy = slow_reservoir(ig,ib)
+                      slow_reservoir(ig,ib) = MAX( (un - ( un - fraction_aeirrig_sw(ig) ) * avail_reserve(3) ) * &
+                                              slow_reservoir(ig,ib), slow_reservoir(ig,ib) + &
+                                              MIN( - irrig_actual(ig,ib) * (un - pcent_vol_irrig ), &
+                                              avail_reserve(2) * fraction_aeirrig_sw(ig) * fast_reservoir(ig,ib) + &
+                                              MIN(zero, avail_reserve(1) * fraction_aeirrig_sw(ig) * stream_reservoir(ig,ib)  - &
+                                              pcent_vol_irrig * irrig_actual(ig,ib) ) ) )
+
+                      slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib)
+                      irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy
+
+                      fast_wdr_dummy = fast_reservoir(ig,ib)
+                      fast_reservoir(ig,ib) = MAX( (un - avail_reserve(2) * fraction_aeirrig_sw(ig) ) * fast_reservoir(ig,ib) , &
+                                              fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1) * fraction_aeirrig_sw(ig) * stream_reservoir(ig,ib)  - &
+                                              pcent_vol_irrig * irrig_actual(ig,ib) ) )
+                      fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib)
+                      irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy
+
+                      stre_wdr_dummy = stream_reservoir(ig,ib)
+                      stream_reservoir(ig,ib) = MAX((un - avail_reserve(1)* fraction_aeirrig_sw(ig) )*stream_reservoir(ig,ib), &
+                                                stream_reservoir(ig,ib)  - &
+                                                pcent_vol_irrig * irrig_actual(ig,ib) )
+                      stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib)
+                      irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy
+
+                      irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib)
+
+                      !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
+                      !We assume that the pot. requirement = Needs * fraction of area equipped for SW/GW
+                      !In the case of surface. we also sustract the withdrawal from Fast/Stream, because both are
+                      !  considered as surface water
+                      IsFail_slow = ( ( irrig_needs(ig,ib)*(un - fraction_aeirrig_sw(ig)) ) > pot_slow_wdr_dummy )
+                      IsFail_fast = ( ( irrig_needs(ig,ib)*fraction_aeirrig_sw(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
+                      IsFail_stre = ( ( irrig_needs(ig,ib)*fraction_aeirrig_sw(ig) - fast_wdr_dummy ) > pot_stre_wdr_dummy )
+
+                      IF( IsFail_stre ) THEN
+                        Count_failure_stre(ig) = un
+                      ENDIF
+                      IF( IsFail_fast ) THEN
+                        Count_failure_fast(ig) = un
+                      ENDIF
+                      IF( IsFail_slow ) THEN
+                        Count_failure_slow(ig) = un
+                      ENDIF
+
+                    ELSE IF (.NOT. old_irrig_scheme .AND. .NOT. select_source_irrig) THEN
+                        ! For   irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
+                        ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
+                        irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib)
+
+                        pot_slow_wdr_dummy = avail_reserve(3)*slow_reservoir(ig,ib)
+                        pot_fast_wdr_dummy = avail_reserve(2)*fast_reservoir(ig,ib)
+                        pot_stre_wdr_dummy = avail_reserve(1)*stream_reservoir(ig,ib)
+                        IsFail_slow = .FALSE. !
+                        IsFail_fast = .FALSE. !
+                        IsFail_stre = .FALSE. !
+
+                        irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),&
+                             & pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
+
+                        !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
+                        !Volume used for irrigation that comes from GW
+                        ! Idem for irrig_fast_source and irrig_str_source
+                        slow_wdr_dummy = slow_reservoir(ig,ib)
+                        slow_reservoir(ig,ib) = MAX( (1-avail_reserve(3) )*slow_reservoir(ig,ib), slow_reservoir(ig,ib) + &
+                             & MIN(zero, avail_reserve(2)*fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1)*stream_reservoir(ig,ib)-irrig_actual(ig,ib))))
+                        slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib)
+                        irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy
+
+                        fast_wdr_dummy = fast_reservoir(ig,ib)
+                        fast_reservoir(ig,ib) = MAX(  (1-avail_reserve(2) )*fast_reservoir(ig,ib) , &
+                             fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1)*stream_reservoir(ig,ib)-irrig_actual(ig,ib)))
+                        fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib)
+                        irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy
+
+                        stre_wdr_dummy = stream_reservoir(ig,ib)
+                        stream_reservoir(ig,ib) = MAX( (1-avail_reserve(1) )*stream_reservoir(ig,ib), stream_reservoir(ig,ib)-irrig_actual(ig,ib) )
+                        stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib)
+                        irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy
+
+                        irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib)
+                        !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
+                        ! Because it follows the old scheme, we do not separate between surface/gw, but consider that
+                        ! priority is given in this order: River, Fast and Slow reservoir.
+                        IsFail_slow = ( ( irrig_needs(ig,ib) - stre_wdr_dummy - fast_wdr_dummy  ) > pot_slow_wdr_dummy )
+                        IsFail_fast = ( ( irrig_needs(ig,ib) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
+                        IsFail_stre = ( irrig_needs(ig,ib) > pot_stre_wdr_dummy )
+
+                        IF( IsFail_stre ) THEN
+                          Count_failure_stre(ig) = un
+                        ENDIF
+                        IF( IsFail_fast ) THEN
+                          Count_failure_fast(ig) = un
+                        ENDIF
+                        IF( IsFail_slow ) THEN
+                          Count_failure_slow(ig) = un
+                        ENDIF
+
+                    ELSE !Old irrigation scheme as in tag 2.0
+                      !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
+                      !Volume used for irrigation that comes from GW
+                      ! Idem for irrig_fast_source and irrig_str_source
+                        irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib)
+
+                        pot_slow_wdr_dummy = slow_reservoir(ig,ib)
+                        pot_fast_wdr_dummy = fast_reservoir(ig,ib)
+                        pot_stre_wdr_dummy = stream_reservoir(ig,ib)
+                        IsFail_slow = .FALSE. !
+                        IsFail_fast = .FALSE. !
+                        IsFail_stre = .FALSE. !
+                        irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),&
+                             &   stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) )
+
+                        slow_wdr_dummy = slow_reservoir(ig,ib)
+                        slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + &
+                             & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))))
+                        slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib)
+                        irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy
+
+                        fast_wdr_dummy = fast_reservoir(ig,ib)
+                        fast_reservoir(ig,ib) = MAX( zero, &
+                             &  fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))
+                        fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib)
+                        irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy
+
+                        stre_wdr_dummy = stream_reservoir(ig,ib)
+                        stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) )
+                        stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib)
+                        irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy
+
+                        irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib)
+                        !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
+                        ! Because it follows the old scheme, we do not separate between surface/gw, but consider that
+                        ! priority is given in this order: River, Fast and Slow reservoir.
+                        IsFail_slow = ( ( irrig_needs(ig,ib) - stre_wdr_dummy - fast_wdr_dummy  ) > pot_slow_wdr_dummy )
+                        IsFail_fast = ( ( irrig_needs(ig,ib) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
+                        IsFail_stre = ( irrig_needs(ig,ib) > pot_stre_wdr_dummy )
+
+                        IF( IsFail_stre ) THEN
+                          Count_failure_stre(ig) = un
+                        ENDIF
+                        IF( IsFail_fast ) THEN
+                          Count_failure_fast(ig) = un
+                        ENDIF
+                        IF( IsFail_slow ) THEN
+                          Count_failure_slow(ig) = un
+                        ENDIF
+
+                    ENDIF
+
+                ENDIF
+            ENDDO
           !
           ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only).
@@ -2342,32 +2646,32 @@
 !> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction
 !> from that subbasin to the one where we need it for irrigation.
-!> 
-          DO ib=1,nbasmax
-
-             stream_tot = SUM(stream_reservoir(ig,:))
-
-             DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)
-                
-                fi = MAXLOC(stream_reservoir(ig,:))
-                ib2 = fi(1)
-
-                irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), stream_reservoir(ig,ib2))
-                stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)
-                irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)
-             
-                stream_tot = SUM(stream_reservoir(ig,:))
-                
-             ENDDO
-             
-          ENDDO
-          !
-       ENDDO
-       !
-       ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams
-       ! which can feed irrigation
-!
-!> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes
-!> to the one where we need it for irrigation.
-       !
+!>
+            DO ib=1,nbasmax
+
+               stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
+
+                 DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)
+
+                      fi = MAXLOC(stream_reservoir(ig,:))
+                      ib2 = fi(1)
+
+                      irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), a_stream_adduction * stream_reservoir(ig,ib2))
+                      stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)
+                      irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)
+
+                      stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
+
+                 ENDDO
+
+            ENDDO
+          !
+      ENDDO
+      !
+      ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams
+      ! which can feed irrigation
+      !
+      !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes
+      !> to the one where we need it for irrigation.
+      !
        IF (is_root_prc) THEN
           ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), &
@@ -2398,5 +2702,5 @@
                          ig2 = neighbours_g(ig,in)
                          IF (ig2 .GT. 0 ) THEN
-                            streams_around(in,:) = stream_reservoir_glo(ig2,:)
+                            streams_around(in,:) = a_stream_adduction * stream_reservoir_glo(ig2,:)
                             igrd(in) = ig2
                          ENDIF
@@ -2409,6 +2713,6 @@
                          ib2=ff(2)
                          !
-                         IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. stream_reservoir_glo(ig2,ib2) > zero ) THEN
-                            adduction = MIN(irrig_deficit_glo(ig,ib), stream_reservoir_glo(ig2,ib2))
+                         IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. a_stream_adduction * stream_reservoir_glo(ig2,ib2) > zero ) THEN
+                            adduction = MIN(irrig_deficit_glo(ig,ib), a_stream_adduction * stream_reservoir_glo(ig2,ib2))
                             stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction
                             irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction
@@ -2476,4 +2780,6 @@
     pond_diag(:) =  zero
     irrigation(:) = zero
+    irrigdeficit(:) = zero !
+    irrigadduct(:) = zero !
     !
     !
@@ -2482,6 +2788,6 @@
        DO ig=1,nbpt
           IF (hydrodiag(ig,ib) > 0 ) THEN
-             hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + & 
-                  &  stream_flow(ig,ib) 
+             hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + &
+                  &  stream_flow(ig,ib)
              slowflow_diag(ig) = slowflow_diag(ig) + slow_flow(ig,ib)
           ENDIF
@@ -2491,4 +2797,6 @@
           flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib)
           irrigation (ig) = irrigation (ig) + irrig_actual(ig,ib) + irrig_adduct(ig,ib)
+          irrigdeficit(ig) = irrigdeficit(ig) + irrig_deficit(ig,ib)!
+          irrigadduct (ig) = irrigadduct (ig) + irrig_adduct(ig,ib)!
        ENDDO
     ENDDO
@@ -2502,6 +2810,12 @@
        !
        irrigation(ig) = irrigation(ig)/totarea(ig)
-       !
-       ! The three output types for the routing : endoheric basins,, rivers and 
+       irrigdeficit(ig) = irrigdeficit(ig)/totarea(ig)!
+       irrigadduct(ig) = irrigadduct(ig)/totarea(ig)!
+       irrig_gw_source(ig) =  irrig_gw_source(ig)/totarea(ig)!
+       irrig_fast_source(ig) =  irrig_fast_source(ig)/totarea(ig)!
+       irrig_str_source(ig) =  irrig_str_source(ig)/totarea(ig)!
+
+       !
+       ! The three output types for the routing : endoheric basins,, rivers and
        ! diffuse coastal flow.
        !
@@ -2513,9 +2827,9 @@
     !
     flood_res = flood_diag + pond_diag
-    
-
-    !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it 
+
+
+    !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it
     !! uniformly over all possible the coastflow gridcells
-    
+
     ! Calculate lake_overflow and remove it from lake_reservoir
     DO ig=1,nbpt
@@ -2541,5 +2855,10 @@
     ! Transform from kg/grid-cell/dt_routing into m^3/grid-cell/s to match output unit of coastalflow
     CALL xios_orchidee_send_field("lake_overflow_coast",lake_overflow_coast/mille/dt_routing)
-   
+
+    ! Counter of reservoir failure to assure irrigation
+    CALL xios_orchidee_send_field("Count_failure_slow", Count_failure_slow)
+    CALL xios_orchidee_send_field("Count_failure_fast", Count_failure_fast)
+    CALL xios_orchidee_send_field("Count_failure_stre", Count_failure_stre)
+
 
   END SUBROUTINE routing_flow
@@ -2590,5 +2909,5 @@
        !
        lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig)
-       
+
        IF ( doswamps ) THEN
           ! Calculate a return flow that will be extracted from the lake reservoir and reinserted in the soil in hydrol
@@ -2635,5 +2954,5 @@
     !
     IMPLICIT NONE
-    
+
 !! INPUT VARIABLES
     INTEGER(i_std), INTENT(in)      :: nbpt               !! Domain size (unitless)
@@ -2652,5 +2971,5 @@
 
 !_ ================================================================================================================================
-    routing_area => routing_area_glo  
+    routing_area => routing_area_glo
     topo_resid => topo_resid_glo
     route_togrid => route_togrid_glo
@@ -2660,8 +2979,8 @@
     hydrodiag=>hydrodiag_glo
     hydroupbasin=>hydroupbasin_glo
-    
+
     IF (is_root_prc) CALL routing_diagnostic(nbp_glo, index_g, lalo_g, resolution_g, contfrac_g, nbrivers_g,basinmap_g)
 
-    routing_area => routing_area_loc  
+    routing_area => routing_area_loc
     topo_resid => topo_resid_loc
     route_togrid => route_togrid_loc
@@ -2671,10 +2990,10 @@
     hydrodiag=>hydrodiag_loc
     hydroupbasin=>hydroupbasin_loc
-    
+
     CALL scatter(nbrivers_g,nbrivers)
     CALL scatter(basinmap_g,basinmap)
     CALL scatter(hydrodiag_glo,hydrodiag_loc)
     CALL scatter(hydroupbasin_glo,hydroupbasin_loc)
-       
+
     CALL xios_orchidee_send_field("basinmap",basinmap)
     CALL xios_orchidee_send_field("nbrivers",nbrivers)
@@ -2692,6 +3011,6 @@
        ENDIF
     ENDIF
-    
-        
+
+
   END SUBROUTINE routing_diagnostic_p
 
@@ -2699,11 +3018,11 @@
 !! SUBROUTINE 	: routing_diagnostic
 !!
-!>\BRIEF         This non-parallelized subroutine gives a diagnostic of the basins used. This produces some information 
+!>\BRIEF         This non-parallelized subroutine gives a diagnostic of the basins used. This produces some information
 !!               on the rivers which are being diagnosed.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : As not all rivers can be monitored in the model, we will only
 !! archive num_largest rivers. In this routine we will diagnose the num_largest largest rivers and print to the standard
-!! output the names of these basins and their area. The list of names of these largest rivers are taken from a list coded in the 
-!! routine routing_names. As this standard output is not sufficient, we will also write it to a netCDF file with the routine 
+!! output the names of these basins and their area. The list of names of these largest rivers are taken from a list coded in the
+!! routine routing_names. As this standard output is not sufficient, we will also write it to a netCDF file with the routine
 !! routing_diagncfile. It is important to keep for diagnostic the fraction of the largest basins in each grid box and keep information
 !! how they are linked one to the other.
@@ -2726,5 +3045,5 @@
     INTEGER(i_std), INTENT(in)                   :: nbpt                !! Domain size  (unitless)
     INTEGER(i_std), INTENT(in)                   :: l_index(nbpt)       !! Indices of the points on the map (unitless)
-    REAL(r_std), INTENT(in)                      :: lalo(nbpt,2)        !! Vector of latitude and longitudes (beware of the order !) 
+    REAL(r_std), INTENT(in)                      :: lalo(nbpt,2)        !! Vector of latitude and longitudes (beware of the order !)
     REAL(r_std), INTENT(in)                      :: resolution(nbpt,2)  !! The size of each grid box in X and Y (m)
     REAL(r_std), INTENT(in)                      :: contfrac(nbpt)      !! Fraction of land in each grid box (unitless;0-1)
@@ -2741,5 +3060,5 @@
     INTEGER(i_std), ALLOCATABLE, DIMENSION(:)    :: nb_pts              !! Number of points in the basin (unitless)
     REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: totarea             !! Total area of basin (m^2)
-    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: tmparea             !! 
+    REAL(r_std), ALLOCATABLE, DIMENSION(:)       :: tmparea             !!
     INTEGER(i_std), ALLOCATABLE, DIMENSION(:)    :: topids              !! The IDs of the first num_largest basins (unitless)
     CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: basin_names         !! Names of the rivers (unitless)
@@ -2818,8 +3137,8 @@
     ALLOCATE(streams_resid(num_largest), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for streams_resid','','')
-    
+
     ALLOCATE(lbasin_area(num_largest,nbpt), lbasin_uparea(num_largest,nbpt), lrivercode(num_largest,nbpt), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for lbasin_area','','')
-    
+
     IF ( .NOT. is_root_prc) THEN
        WRITE(numout,*) "routing_diagnostic is not suitable for running in parallel"
@@ -2828,11 +3147,11 @@
        CALL ipslerr_p(3,'routing_diagnostic','This routine is not suitable for running in parallel','','')
     ENDIF
-    
-    
+
+
     !Config Key   = RIVER_DESC
     !Config Desc  = Writes out a description of the rivers
     !Config If    = RIVER_ROUTING
     !Config Def   = n
-    !Config Help  = This flag allows to write out a file containing the list of 
+    !Config Help  = This flag allows to write out a file containing the list of
     !Config         rivers which are beeing simulated. It provides location of outflow
     !Config         drainage area, name and ID.
@@ -2880,7 +3199,7 @@
        ENDDO
     ENDDO
-    
+
     nb_small = MIN(nb_small, 349)
-    
+
     ALLOCATE(basin_names(nb_small), stat=ier)
     IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for basins_names','','')
@@ -2922,5 +3241,5 @@
              longest_river = MAX(longest_river, ic)
              !
-             ! Now that we have an outflow check if it is one of the num_largest rivers. 
+             ! Now that we have an outflow check if it is one of the num_largest rivers.
              ! In this case we keeps the location so we diagnose it.
              !
@@ -3123,5 +3442,5 @@
        tuparea(1:ii) = lbasin_uparea(idbas,1:ii)
        !
-       CALL routing_diagcode(ii, tpts, tptbas, tuparea, tslen, MAXVAL(tslen), allstreams, upstreamchange, tcode)  
+       CALL routing_diagcode(ii, tpts, tptbas, tuparea, tslen, MAXVAL(tslen), allstreams, upstreamchange, tcode)
        !
        lrivercode(idbas,:) = 0
@@ -3178,5 +3497,5 @@
           ENDIF
        ENDDO
-       
+
        IF ( name_found > 1 ) THEN
           DO ic=num_largest,1,-1
@@ -3199,5 +3518,5 @@
           ENDDO
        ENDIF
-       
+
     ENDDO
     !
@@ -3281,5 +3600,5 @@
        WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. zero)
     END IF
-    
+
     DEALLOCATE(pts)
     DEALLOCATE(outpt)
@@ -3301,5 +3620,5 @@
        WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. 0.)
     END IF
-    
+
   END SUBROUTINE routing_diagnostic
   !
@@ -3308,5 +3627,5 @@
 !!
 !>\BRIEF       This subroutine determines the code in the Pfafstetter system for all points
-!!              within the given catchment.  
+!!              within the given catchment.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -3322,5 +3641,5 @@
 !_ ================================================================================================================================
 
-  SUBROUTINE routing_diagcode(ip, tpts, tpbas, tuparea, tslen, ls, allstreams, upstreamchange, streamcode)  
+  SUBROUTINE routing_diagcode(ip, tpts, tpbas, tuparea, tslen, ls, allstreams, upstreamchange, streamcode)
     !
     IMPLICIT NONE
@@ -3388,5 +3707,5 @@
              IF ( streamcode(ic) == indsubbas(ib) ) THEN
                 it =it+1
-                iw(it)=ic 
+                iw(it)=ic
              ENDIF
           ENDDO
@@ -3407,7 +3726,7 @@
              ENDDO
           ELSE
-             ! 
+             !
              ! Else do the Pfafstetter numbering
-             ! 
+             !
              !
              ! Where do we have the 4 largest change in upstream area on this stream.
@@ -3427,5 +3746,5 @@
              ENDDO
              !
-             ! Find all streams which are identical up to that junction and increase their code accordingly 
+             ! Find all streams which are identical up to that junction and increase their code accordingly
              !
              DO i=1,it
@@ -3474,5 +3793,5 @@
 !!
 !>\BRIEF         This subroutine creates a netCDF file containing all the informations
-!!                on the largest rivers which can be used for a refined analysis. 
+!!                on the largest rivers which can be used for a refined analysis.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -3571,5 +3890,5 @@
                &         TRIM(river_file_name),'(Solution ?)')
        ENDIF
-       
+
        iret = NF90_DEF_DIM(fid, 'y', jjm_g, dims(2))
        IF (iret /= NF90_NOERR) THEN
@@ -3624,5 +3943,5 @@
     ! 1.3 Add attributes to the coordinate variables
     !
-    iret = NF90_PUT_ATT(fid, nlonid, 'units', "degrees_east") 
+    iret = NF90_PUT_ATT(fid, nlonid, 'units', "degrees_east")
     IF (iret /= NF90_NOERR) THEN
        CALL ipslerr_p (3,'routing_diagncfile', 'Could not add attribut to variable '//lon_name//' for the file :', &
@@ -3754,9 +4073,9 @@
                 ELSE
                    basinfrac(i,j) = MIN(basinfrac(i,j), un)
-                ENDIF 
+                ENDIF
              ENDDO
           ENDDO
           !
-          ! 
+          !
           ! 2.2 Define the variables in the netCDF file
           !
@@ -3834,5 +4153,5 @@
           ENDIF
           !
-          ! Longitude of outflow point 
+          ! Longitude of outflow point
           att_str='Longitude_of_outflow_point'
           iret = NF90_PUT_ATT(fid, varid, att_str, lalo(outpt(ib,1),2))
@@ -3877,5 +4196,5 @@
           ENDIF
           !
-          ! Average number of hops to go to the ocean for any stream 
+          ! Average number of hops to go to the ocean for any stream
           att_str='Average_number_of_hops_to_ocean_for_any_stream'
           iret = NF90_PUT_ATT(fid, varid, att_str, streams_avehops(ib))
@@ -3885,5 +4204,5 @@
           ENDIF
           !
-          ! Maximum number of hops to go to the ocean for any stream 
+          ! Maximum number of hops to go to the ocean for any stream
           att_str='Maximum_number_of_hops_to_ocean_for_any_stream'
           iret = NF90_PUT_ATT(fid, varid, att_str, streams_maxhops(ib))
@@ -4027,5 +4346,5 @@
 !! SUBROUTINE 	: routing_basins_p
 !!
-!>\BRIEF        This parallelized subroutine computes the routing map if needed. 
+!>\BRIEF        This parallelized subroutine computes the routing map if needed.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -4056,5 +4375,5 @@
 !    INTEGER(i_std)    :: neighbours_tmp(nbpt,8)
 !    INTEGER(i_std) :: i,j
-    
+
 !    DO i=1,nbp_loc
 !      DO j=1,NbNeighb
@@ -4063,9 +4382,9 @@
 !	ELSE
 !	  neighbours_tmp(i,j)=neighbours(i,j)+nbp_para_begin(mpi_rank)-1
-!	ENDIF  
+!	ENDIF
 !      ENDDO
 !    ENDDO
 
-    routing_area => routing_area_glo  
+    routing_area => routing_area_glo
     topo_resid => topo_resid_glo
     route_togrid => route_togrid_glo
@@ -4073,8 +4392,8 @@
     route_nbintobas => route_nbintobas_glo
     global_basinid => global_basinid_glo
- 
+
     IF (is_root_prc) CALL routing_basins(nbp_glo,lalo_g, neighbours_g, resolution_g, contfrac_g)
 
-    routing_area => routing_area_loc  
+    routing_area => routing_area_loc
     topo_resid => topo_resid_loc
     route_togrid => route_togrid_loc
@@ -4089,16 +4408,16 @@
     CALL scatter(route_nbintobas_glo,route_nbintobas_loc)
     CALL scatter(global_basinid_glo,global_basinid_loc)
-    
+
   END SUBROUTINE routing_basins_p
-  ! 
- 
+  !
+
 !! ================================================================================================================================
 !! SUBROUTINE 	: routing_basins
 !!
 !>\BRIEF        This non-parallelized subroutine reads in the map of basins and flow direction to construct
-!!              the catchments of each grid box. 
+!!              the catchments of each grid box.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) :
-!! The work is done in a number of steps which are performed locally on the 
+!! The work is done in a number of steps which are performed locally on the
 !! GCM grid:
 !!  1) First we find the grid-points of the high resolution routing grid which are
@@ -4130,5 +4449,5 @@
     INTEGER(i_std), INTENT(in)                    :: nbpt                  !! Domain size (unitless)
     REAL(r_std), INTENT(in)                       :: lalo(nbpt,2)          !! Vector of latitude and longitudes (beware of the order !)
-    INTEGER(i_std), INTENT(in)                    :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 
+    INTEGER(i_std), INTENT(in)                    :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
                                                                            !! (1=North and then cloxkwise)
     REAL(r_std), INTENT(in)                       :: resolution(nbpt,2)    !! The size of each grid box in X and Y (m)
@@ -4212,8 +4531,8 @@
     !
     IF (debug) THEN
-       IF (ANY(diagbox .GT. nbpt)) THEN 
+       IF (ANY(diagbox .GT. nbpt)) THEN
           WRITE(numout,*) "Debug diganostics : nbpt, diagbox", nbpt, diagbox
           call ipslerr_p(3,'routing_basin', &
-               &      'Problem with diagbox in debug mode.', & 
+               &      'Problem with diagbox in debug mode.', &
                &      'diagbox values can''t be greater than land points number.', &
                &      '(decrease diagbox wrong value)')
@@ -4230,5 +4549,5 @@
     !Config Def   = routing.nc
     !Config Help  = The file provided here should alow the routing module to
-    !Config         read the high resolution grid of basins and the flow direction 
+    !Config         read the high resolution grid of basins and the flow direction
     !Config         from one mesh to the other.
     !Config Units = [FILE]
@@ -4280,7 +4599,7 @@
     !! Basins : Provides a uniqe ID for each basin. These IDs are also used to get
     !! the name of the basin from the table in routine routing_names.
-    !! 
+    !!
     !! Topoind :  is the topographic index for the retention time of the water in the
-    !! grid box. It has been computed with the following formula : 1000 x sqrt(d^3/Dz) 
+    !! grid box. It has been computed with the following formula : 1000 x sqrt(d^3/Dz)
     !! where d is the distance of the river from the current grid box to the next one
     !! as indicated by the variable trip.
@@ -4333,5 +4652,5 @@
           ! Resolution in longitude
           !
-          coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos )     
+          coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos )
           IF ( ip .EQ. 1 ) THEN
              resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip,jp) ) * pi/180. * R_Earth * coslat
@@ -4361,5 +4680,5 @@
     callsign = "routing_basins"
     ok_interpol = .FALSE.
-    !  
+    !
     nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2
     njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2
@@ -4454,5 +4773,5 @@
     ALLOCATE (nbcoastal(nbpt), coastal_basin(nbpt,nwbas), stat=ALLOC_ERR)
     IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_basins','Pb in allocate for nbcoastal','','')
-    
+
     !    Order all sub points in each grid_box and find the sub basins
     !
@@ -4527,5 +4846,5 @@
             & nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,&
             & basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,&
-            & nbcoastal, coastal_basin) 
+            & nbcoastal, coastal_basin)
        !
        !
@@ -4548,5 +4867,5 @@
          & basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, &
          & nbcoastal, coastal_basin, invented_basins)
-    ! 
+    !
     !
     IF (printlev>=1) WRITE(numout,*) 'The maximum number of basins in any grid :', MAXVAL(basin_count)
@@ -4780,5 +5099,5 @@
     !  We simplify the outflow. We only need the direction normal to the
     !     box boundary and the 4 corners.
-    ! 
+    !
     ! Northern border
     IF ( trip_bx(1,1) .EQ. 102 ) trip_bx(1,1) = 101
@@ -4804,5 +5123,5 @@
     DO jp=2,nbj-1
        IF ( trip_bx(1,jp) .EQ. 106 .OR. trip_bx(1,jp) .EQ. 108 ) trip_bx(1,jp) = 107
-    ENDDO       
+    ENDDO
     !
     !
@@ -4854,5 +5173,5 @@
        IF ( coords(ipos+1) /= undef_sechiba) THEN
          IF ( COUNT(coords(ipos:nb_out) == coords(ipos)) > 1 ) THEN
-            coords(ipos:nb_out-1) = coords(ipos+1:nb_out) 
+            coords(ipos:nb_out-1) = coords(ipos+1:nb_out)
             coords(nb_out:nb_in) = undef_sechiba
             nb_out = nb_out - 1
@@ -4890,5 +5209,5 @@
        DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1)
           ll = MINLOC(coords_tmp(:), coords_tmp /= undef_sechiba)
-          ind(ipos) = ll(1) 
+          ind(ipos) = ll(1)
           coords_tmp(ll(1)) = undef_sechiba
           ipos = ipos + 1
@@ -4897,5 +5216,5 @@
        DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1)
           ll = MAXLOC(coords_tmp(:), coords_tmp /= undef_sechiba)
-          ind(ipos) = ll(1) 
+          ind(ipos) = ll(1)
           coords_tmp(ll(1)) = undef_sechiba
           ipos = ipos + 1
@@ -5066,6 +5385,6 @@
           coast_pts(sz(ipb)) = bname(trans(ip))
           sz(trans(ip)) = 0
-          pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1) 
-          pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2) 
+          pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1)
+          pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2)
        ENDDO
     ENDIF
@@ -5116,5 +5435,5 @@
           ENDIF
           !
-          ! Third issue : we have more than one outflow from the boxes. This could be 
+          ! Third issue : we have more than one outflow from the boxes. This could be
           !             - flow into 2 or more neighboring GCM grids
           !             - flow into a neighboring GCM grids and into the ocean or be a return flow (=97. =98, =99)
@@ -5181,5 +5500,5 @@
     !
     !
-    ! Sort the output by size of the various basins. 
+    ! Sort the output by size of the various basins.
     !
     nb_basin = COUNT(sz(1:nbb) .GT. 0)
@@ -5218,5 +5537,5 @@
        CALL ipslerr_p(3,'routing_findbasins','Probleme less outflow points than basins','','')
     ENDIF
-    
+
   END SUBROUTINE routing_findbasins
   !
@@ -5224,8 +5543,8 @@
 !! SUBROUTINE 	: routing_simplify
 !!
-!>\BRIEF         This subroutine symplifies the routing out of each basin by taking 
+!>\BRIEF         This subroutine symplifies the routing out of each basin by taking
 !!               out redundancies at the borders of the GCM box.
 !!               The aim is to have only one outflow point per basin and grid box.
-!!               But here we will not change the direction of the outflow.  
+!!               But here we will not change the direction of the outflow.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -5353,9 +5672,9 @@
                 IF ( trip_flow(ip,jp,1) .EQ. outflow(ismall,1) .AND.&
                      & trip_flow(ip,jp,2) .EQ. outflow(ismall,2) ) THEN
-                   ! Now that we have found a point of the smallest sub-basin we 
+                   ! Now that we have found a point of the smallest sub-basin we
                    ! look around for another sub-basin
                    ib = 1
                    not_found = .TRUE.
-                   DO WHILE ( not_found .AND. ib .LE. itodo ) 
+                   DO WHILE ( not_found .AND. ib .LE. itodo )
                       IF ( ib .NE. ill(1) ) THEN
                          ibas = todopt(ib)
@@ -5415,5 +5734,5 @@
 !!
 !>\BRIEF        This subroutine cuts the original basin which has more than one outflow
-!!              into as many subbasins as outflow directions.  
+!!              into as many subbasins as outflow directions.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -5499,5 +5818,5 @@
     !
     CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz)
-    ! 
+    !
 !    IF ( debug ) THEN
 !       DO ib = nb_in+1,nb
@@ -5516,5 +5835,5 @@
     !  Take out the small sub-basins. That is those which have only one grid box
     !  This is only done if we need to save space in the number of basins. Else we
-    !  can take it easy and keep diverging sub-basins for the moment. 
+    !  can take it easy and keep diverging sub-basins for the moment.
     !
     IF ( nbbasins .GE. nbasmax ) THEN
@@ -5543,5 +5862,5 @@
                       DO ibb=1,nbout
                          IF ( iip .EQ. outflow(ibb,1) .AND. jjp .EQ. outflow(ibb,2) ) THEN
-                            outsz(ibb) = outsz(ibb)+1 
+                            outsz(ibb) = outsz(ibb)+1
                             trip_flow(ip,jp,1) = outflow(ibb,1)
                             trip_flow(ip,jp,2) = outflow(ibb,2)
@@ -5584,5 +5903,5 @@
           ENDIF
        ENDDO
-       ! A short verification 
+       ! A short verification
        IF ( SUM(sz(nb_in+1:nb)) .NE. basin_sz) THEN
           WRITE(numout,*) 'Lost some points while spliting the basin'
@@ -5598,5 +5917,5 @@
           CALL ipslerr_p(3,'routing_cutbasin','Lost some points while spliting the basin','','')
        ENDIF
-       
+
        IF ( debug ) THEN
           DO ib = nb_in+1,nb
@@ -5623,5 +5942,5 @@
 !!
 !>\BRIEF        This subroutine finds, for each point, the distance to the outflow
-!!               point along the flowlines of the basin. 
+!!               point along the flowlines of the basin.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -5691,5 +6010,5 @@
              ENDIF
              !
-             DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml) 
+             DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml)
                 cnt = cnt + 1
                 ntripi = ntripi + inc(trp,1)
@@ -5707,5 +6026,5 @@
                 trp = NINT(trip(ntripi, ntripj))
              ENDDO
-             
+
              IF ( cnt .EQ. iml*jml) THEN
                 WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp
@@ -5717,7 +6036,7 @@
                 CALL ipslerr_p(3,'routing_hierarchy','We could not route point','','')
              ENDIF
-             
+
              hierarchy(ip,jp) = topohier
-             
+
           ENDIF
        ENDDO
@@ -5731,5 +6050,5 @@
 !!
 !>\BRIEF        This subroutine simply computes the route to each outflow point
-!!              and returns the outflow point for each point in the basin.  
+!!              and returns the outflow point for each point in the basin.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -5759,5 +6078,5 @@
     INTEGER(i_std), DIMENSION(nbvmax,2), INTENT(out)        :: outflow   !!
     INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2), INTENT(out) :: trip_flow !!
-    INTEGER(i_std), INTENT(out)                             :: nbout     !! 
+    INTEGER(i_std), INTENT(out)                             :: nbout     !!
     INTEGER(i_std), DIMENSION(nbvmax), INTENT(out)          :: outsz     !!
     !
@@ -5790,5 +6109,5 @@
     !
     !
-    !  Get the outflows and determine for each point to which outflow point it belong 
+    !  Get the outflows and determine for each point to which outflow point it belong
     !
     nbout = 0
@@ -5800,5 +6119,5 @@
              outflow(nbout,1) = ip
              outflow(nbout,2) = jp
-          ENDIF 
+          ENDIF
           IF ( trip(ip,jp) .GT. 0) THEN
              trip_flow(ip,jp,1) = ip
@@ -5815,5 +6134,5 @@
              trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2))
              cnt = 0
-             DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj) 
+             DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj)
                 cnt = cnt + 1
                 trip_flow(ip,jp,1) = trip_flow(ip,jp,1) + inc(trp,1)
@@ -5862,5 +6181,5 @@
 !!
 !>\BRIEF        This subroutine puts the basins found for grid box in the global map.
-!!               Connection can only be made later when all information is together. 
+!!               Connection can only be made later when all information is together.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -5888,7 +6207,7 @@
     INTEGER(i_std), INTENT (in)                :: nbpt                   !! Domain size (unitless)
     INTEGER(i_std), INTENT (in)                :: ib                     !! Current basin (unitless)
-    INTEGER(i_std), INTENT(in)                 :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 
+    INTEGER(i_std), INTENT(in)                 :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
                                                                          !! (1=North and then clockwise)
-!! LOCAL VARIABLES 
+!! LOCAL VARIABLES
     REAL(r_std), DIMENSION(nbvmax,nbvmax)      :: area_bx                !! Area of each small box in the grid box (m^2)
     INTEGER(i_std), DIMENSION(nbvmax,nbvmax)   :: trip_bx                !! The trip field for each of the smaller boxes (unitless)
@@ -5924,11 +6243,11 @@
        !
        basin_count(ib) = basin_count(ib)+1
-       if (basin_count(ib) > nwbas) then 
+       if (basin_count(ib) > nwbas) then
           WRITE(numout,*) 'ib=',ib
           call ipslerr_p(3,'routing_globalize', &
-               &      'Problem with basin_count : ', & 
+               &      'Problem with basin_count : ', &
                &      'It is greater than number of allocated basin nwbas.', &
                &      '(stop to count basins)')
-       endif 
+       endif
        basin_id(ib,basin_count(ib)) = basin_inbxid(ij)
        !
@@ -6025,5 +6344,5 @@
 !! SUBROUTINE 	: routing_linkup
 !!
-!>\BRIEF         This subroutine makes the connections between the basins and ensure global coherence. 
+!>\BRIEF         This subroutine makes the connections between the basins and ensure global coherence.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) :
@@ -6078,5 +6397,5 @@
     INTEGER(i_std)                                 :: ff(1)                 !!
     !
-    ! ERRORS 
+    ! ERRORS
     LOGICAL                                        :: error1, error2, error3, error4, error5 !! (true/false)
     !
@@ -6084,5 +6403,5 @@
     LOGICAL, PARAMETER                             :: check = .TRUE.       !! (true/false)
 
-!_ ================================================================================================================================    
+!_ ================================================================================================================================
     error1=.FALSE.
     error2=.FALSE.
@@ -6133,5 +6452,5 @@
                 !
              ENDDO
-             ! 
+             !
              !  If we have more than one basin we will take the one which is lowest
              !  in the hierarchy.
@@ -6146,5 +6465,5 @@
                       bop = sbl
                    ELSE
-                      ! The same hierarchy is allowed if both grids flow in about 
+                      ! The same hierarchy is allowed if both grids flow in about
                       ! the same direction :
                       IF ( ( MOD(basin_flowdir(inp,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. &
@@ -6179,5 +6498,5 @@
           !
           ! In case the outflow point was ocean or we did not find the correct basin we start to look
-          ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding 
+          ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding
           ! basin index (bp1 & bm1).
           !
@@ -6202,5 +6521,5 @@
                          bp1 = sbl
                       ELSE
-                         ! The same hierarchy is allowed if both grids flow in about 
+                         ! The same hierarchy is allowed if both grids flow in about
                          ! the same direction :
                          angle=MODULO(basin_flowdir(dp1,sbl)-basin_flowdir(sp,sb)+8,8)
@@ -6223,6 +6542,6 @@
                       IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dm1,sbl) ) THEN
                          bm1 = sbl
-                      ELSE                         
-                         ! The same hierarchy is allowed if both grids flow in about 
+                      ELSE
+                         ! The same hierarchy is allowed if both grids flow in about
                          ! the same direction :
                          angle=MODULO(basin_flowdir(dm1,sbl)-basin_flowdir(sp,sb)+8,8)
@@ -6240,5 +6559,5 @@
              ! First deal with the case on land.
              !
-             ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1 
+             ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1
              ! and not return to our current grid. If it is the current grid
              ! then we can not do anything with that neighbour. Thus we set the
@@ -6250,5 +6569,5 @@
                 IF (basin_flowdir(dp1,bp1) .GT. 0) THEN
                    outdp1 = neighbours(dp1,basin_flowdir(dp1,bp1))
-                   IF ( outdp1 .EQ. sp ) outdp1 = undef_int 
+                   IF ( outdp1 .EQ. sp ) outdp1 = undef_int
                 ELSE
                    outdp1 = nbpt + 1
@@ -6270,5 +6589,5 @@
              bop = undef_int
              IF ( outdp1 .LT. undef_int .AND. outdm1 .LT. undef_int) THEN
-                ! 
+                !
                 ! In this case we let the current basin flow into the smaller one
                 !
@@ -6292,5 +6611,5 @@
              ELSE
                 !
-                ! Now we are at the point where none of the neighboring points is suitable 
+                ! Now we are at the point where none of the neighboring points is suitable
                 ! or we have a coastal point.
                 !
@@ -6301,7 +6620,7 @@
                 ELSE
                    !
-                   ! If we are on a land point with only land neighbors but no one suitable to let the 
+                   ! If we are on a land point with only land neighbors but no one suitable to let the
                    ! water flow into we have to look for a solution in the current grid box.
-                   ! 
+                   !
                    !
                    IF ( bp1 .LT. 0 .AND. bm1 .LT. 0 ) THEN
@@ -6402,5 +6721,5 @@
                & .AND. basin_flowdir(sp,sb) .GT. 0) THEN
              !
-             
+
              IF (check) &
                   WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb
@@ -6451,5 +6770,5 @@
     IF (error2) THEN
        CALL ipslerr_p(3,'routing_linkup', &
-            &      'In the routine which make connections between the basins and ensure global coherence,', & 
+            &      'In the routine which make connections between the basins and ensure global coherence,', &
             &      'there is a problem with outflow linkup without any valid direction. Try with check=.TRUE.', &
             &      '(Perhaps there is a problem with the grid.)')
@@ -6461,6 +6780,6 @@
     ENDIF
     IF (error4) THEN
-       WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", & 
-            &  " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, & 
+       WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", &
+            &  " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, &
             &  basin_hierarchy(sp,sb),basin_hierarchy(sp,sbl)
        CALL ipslerr_p(3,'routing_linkup', &
@@ -6471,5 +6790,5 @@
        WRITE(numout,*) 'We could not find the basin into which we need to flow'
        WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb
-       WRITE(numout,*) 'Explored neighbours :', dm1, dp1 
+       WRITE(numout,*) 'Explored neighbours :', dm1, dp1
        WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb)
        WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb)
@@ -6507,5 +6826,5 @@
 !!               will know how much area is upstream. It will help decide how to procede with the
 !!               the truncation later and allow to set correctly in outflow_grid the distinction
-!!               between coastal and river flow. 
+!!               between coastal and river flow.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -6613,6 +6932,6 @@
     WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin)
     !
-    ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow 
-    ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow 
+    ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow
+    ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow
     ! (i.e. outflow_grid = -2). The return flow is not touched.
     !
@@ -6654,25 +6973,25 @@
 !>\BRIEF         This subroutine reduces the number of basins per grid to the value chosen by the user.
 !!               It also computes the final field which will be used to route the water at the
-!!               requested truncation.  
-!!
-!! DESCRIPTION (definitions, functional, design, flags) : 
+!!               requested truncation.
+!!
+!! DESCRIPTION (definitions, functional, design, flags) :
 !! Truncate if needed and find the path closest to the high resolution data.
 !!
 !! The algorithm :
-!! 
+!!
 !! We only go through this procedure only as many times as there are basins to take out at most.
 !! This is important as it allows the simplifications to spread from one grid to the other.
-!! The for each step of the iteration and at each grid point we check the following options for 
+!! The for each step of the iteration and at each grid point we check the following options for
 !! simplifying the pathways of water :
 !! 1) If the basin of a grid flows into another basin of the same grid. Kill the one which only
 !!    served as a transition
-!! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow. 
-!!    We kill the smallest one and put into the largest basin. There is no need to manage many 
-!!    basins going into the ocean as coastal flows. 
-!! 3) If we have streams run in parallel from one gird box to the others (that is these are 
+!! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow.
+!!    We kill the smallest one and put into the largest basin. There is no need to manage many
+!!    basins going into the ocean as coastal flows.
+!! 3) If we have streams run in parallel from one gird box to the others (that is these are
 !!    different basins) we will put the smaller one in the larger one. This may hapen at any
 !!    level of the flow but in theory it should propagate downstream.
 !! 4) If we have two basins with the same ID but flow into different grid boxes we sacrifice
-!!    the smallest one and route it through the largest. 
+!!    the smallest one and route it through the largest.
 !!
 !! Obviously if any of the options find something then we skip the rest and take out the basin.:\n
@@ -6758,5 +7077,5 @@
        ENDDO
        !
-       ! Go through the basins which need to be truncated.       
+       ! Go through the basins which need to be truncated.
        !
        DO ibt=1,nbtruncate
@@ -6782,6 +7101,6 @@
           ENDDO
           !
-          ! 2) Merge two basins which flow into the ocean as coastal or return flow 
-          ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and 
+          ! 2) Merge two basins which flow into the ocean as coastal or return flow
+          ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and
           ! have not found anything yet!
           !
@@ -6818,5 +7137,5 @@
           ENDIF
           !
-          !   3) If we have basins flowing into the same grid but different basins then we put them 
+          !   3) If we have basins flowing into the same grid but different basins then we put them
           !   together. Obviously we first work with the grid which has most streams running into it
           !   and putting the smallest in the largests catchments.
@@ -6947,5 +7266,5 @@
        ENDDO
        !
-       !      
+       !
     ENDDO
     !
@@ -7103,6 +7422,6 @@
                 WRITE(numout,*) 'We should not be here as the basin flows into the pampa'
                 WRITE(numout,*) 'Last correct point :', pold, bold
-                WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold) 
-                WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold) 
+                WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold)
+                WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold)
                 WRITE(numout,*) 'Where we ended up :', igrif,ibasf
                 CALL ipslerr_p(3,'routing_truncate','Problem with routing..','','')
@@ -7139,5 +7458,5 @@
     DO ib=1, pickmax
        ff = MAXLOC(floflo)
-       ! tdo - To take into account rivers that do not flow to the oceans 
+       ! tdo - To take into account rivers that do not flow to the oceans
        IF ( route_tobasin(ff(1), ff(2)) .GT. nbasmax ) THEN
 !       IF ( route_tobasin(ff(1), ff(2)) .EQ. nbasmax + 2) THEN
@@ -7169,5 +7488,5 @@
 !!
 !>\BRIEF        The aim of this subroutine is to kill a basin (that is put into another larger one).
-!!              When we do this we need to be careful and change all associated variables.  
+!!              When we do this we need to be careful and change all associated variables.
 !!
 !! DESCRIPTION (definitions, functional, design, flags) : None
@@ -7232,5 +7551,5 @@
     inf = 0
     DO WHILE (igrif .GT. 0)
-       fetch_basin(igrif,ibasf) =  fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill) 
+       fetch_basin(igrif,ibasf) =  fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill)
        it = outflow_grid(igrif, ibasf)
        ibasf = outflow_basin(igrif, ibasf)
@@ -7249,5 +7568,5 @@
        ibasf = outflow_basin(igrif, ibasf)
        igrif = it
-    ENDDO   
+    ENDDO
     !
     !  Redirect the flows which went into the basin to be killed before we change everything
@@ -7291,5 +7610,5 @@
        ENDIF
        inflow_number(ing, inb) = inflow_number(ing, inb) - 1
-       
+
     ENDIF
     !
@@ -7340,5 +7659,5 @@
     basin_count(ib) = basin_count(ib) - 1
     !
-  END SUBROUTINE routing_killbas 
+  END SUBROUTINE routing_killbas
   !
 !! ================================================================================================================================
@@ -7741,5 +8060,5 @@
   !
 !! ================================================================================================================================
-!! SUBROUTINE 	: routing_irrigmap
+!! SUBROUTINE 	: routing_floodmap
 !!
 !>\BRIEF         This  subroutine interpolates the 0.5x0.5 degree based map of irrigated areas to the resolution of the model.
@@ -7747,5 +8066,5 @@
 !! DESCRIPTION (definitions, functional, design, flags) : None
 !!
-!! RECENT CHANGE(S): None
+!! RECENT CHANGE(S): Irrigated is interpolated from slowproc as irrigated_next
 !!
 !! MAIN OUTPUT VARIABLE(S):
@@ -7757,6 +8076,6 @@
 !_ ================================================================================================================================
 
-SUBROUTINE routing_irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, &
-       &                       init_irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)
+SUBROUTINE routing_floodmap (nbpt, index, lalo, neighbours, resolution, contfrac, &
+        &                       init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)
     !
     IMPLICIT NONE
@@ -7780,10 +8099,10 @@
     INTEGER(i_std), INTENT(in)                     :: hist_id               !! Access to history file (unitless)
     INTEGER(i_std), INTENT(in)                     :: hist2_id              !! Access to history file 2 (unitless)
-    LOGICAL, INTENT(in)                            :: init_irrig            !! Logical to initialize the irrigation (true/false)
+
     LOGICAL, INTENT(in)                            :: init_flood            !! Logical to initialize the floodplains (true/false)
     LOGICAL, INTENT(in)                            :: init_swamp            !! Logical to initialize the swamps (true/false)
     !
 !! OUTPUT VARIABLES
-    REAL(r_std), INTENT(out)                       :: irrigated(:)          !! Irrigated surface in each grid box (m^2)
+
     REAL(r_std), INTENT(out)                       :: floodplains(:)        !! Surface which can be inundated in each grid box (m^2)
     REAL(r_std), INTENT(out)                       :: swamp(:)              !! Surface which can be swamp in each grid box (m^2)
@@ -7791,5 +8110,5 @@
 !! LOCAL VARIABLES
     ! Interpolation variables
-    ! 
+    !
     INTEGER(i_std)                                 :: nbpmax, nix, njx, fopt !!
     CHARACTER(LEN=30)                              :: callsign              !!
@@ -7807,7 +8126,7 @@
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: latrel                !! Latitude
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: lonrel                !! Longitude
-    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)       :: irrigated_frac        !! Irrigated fraction of the grid box (unitless;0-1)
+
     REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:)     :: flood_fracmax         !! Maximal flooded fraction of the grid box (unitless;0-1)
-    REAL(r_std)                                    :: area_irrig            !! Irrigated surface in the grid box (m^2)
+
     REAL(r_std)                                    :: area_flood(ntype)     !! Flooded surface in the grid box (m^2)
     REAL(r_std)                                    :: resolution_1          !! temporary variable
@@ -7826,5 +8145,5 @@
     !Config         with the area in m^2 of the area irrigated within each
     !Config         0.5 0.5 deg grid box. The map currently used is the one
-    !Config         developed by the Center for Environmental Systems Research 
+    !Config         developed by the Center for Environmental Systems Research
     !Config         in Kassel (1995).
     !Config Units = [FILE]
@@ -7851,14 +8170,11 @@
     !
     ALLOCATE (latrel(iml,jml), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for latrel','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for latrel','','')
 
     ALLOCATE (lonrel(iml,jml), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for lonrel','','')
-
-    ALLOCATE (irrigated_frac(iml,jml), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrigated_frac','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for lonrel','','')
 
     ALLOCATE (flood_fracmax(iml,jml,ntype), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for flood_fracmax','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for flood_fracmax','','')
 
     IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lonrel, latrel, lev, tml, itau, date, dt, fid)
@@ -7867,6 +8183,4 @@
     CALL bcast(latrel)
     !
-    IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 1, 1, irrigated_frac)
-    CALL bcast(irrigated_frac)
     IF (is_root_prc) CALL flinget(fid, 'lake', iml, jml, lml, tml, 1, 1, flood_fracmax(:,:,ilake))
     IF (is_root_prc) CALL flinget(fid, 'dam', iml, jml, lml, tml, 1, 1, flood_fracmax(:,:,idam))
@@ -7883,10 +8197,5 @@
     DO ip=1,iml
        DO jp=1,jml
-          !
-          IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN
-             irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100.
-             IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero
-          ENDIF
-          !
+
           DO itype=1,ntype
              IF ( flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN
@@ -7898,10 +8207,9 @@
        ENDDO
     ENDDO
-    
+
     IF (printlev>=2) THEN
        WRITE(numout,*) 'lonrel : ', MAXVAL(lonrel), MINVAL(lonrel)
        WRITE(numout,*) 'latrel : ', MAXVAL(latrel), MINVAL(latrel)
-       WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), &
-            MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba)
+
        WRITE(numout,*) 'flood_fracmax : ', MINVAL(flood_fracmax, MASK=flood_fracmax .GT. 0), &
             MAXVAL(flood_fracmax, MASK=flood_fracmax .LT. undef_sechiba)
@@ -7911,8 +8219,8 @@
     !
     ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for resol_lu','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for resol_lu','','')
 
     ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for mask','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for mask','','')
     mask(:,:) = 0
 
@@ -7930,5 +8238,5 @@
           ! Resolution in longitude
           !
-          coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos )     
+          coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos )
           IF ( ip .EQ. 1 ) THEN
              resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip,jp) ) * pi/180. * R_Earth * coslat
@@ -7955,5 +8263,5 @@
     ! Some lmargin is taken.
     !
-    callsign = 'Irrigation map'
+    callsign = 'Flood map'
     ok_interpol = .FALSE.
     IF (is_root_prc) THEN
@@ -7969,9 +8277,9 @@
 
     ALLOCATE(irrsub_index(nbpt, nbpmax, 2), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrsub_index','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for irrsub_index','','')
     irrsub_index(:,:,:)=0
 
     ALLOCATE(irrsub_area(nbpt, nbpmax), STAT=ALLOC_ERR)
-    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrsub_area','','')
+    IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for irrsub_area','','')
     irrsub_area(:,:)=zero
 
@@ -7982,10 +8290,9 @@
     !
     WHERE (irrsub_area < 0) irrsub_area=zero
-    !  
+    !
     ! Test here if not all sub_area are larger than 0 if so, then we need to increase nbpmax
     !
     DO ib=1,nbpt
        !
-       area_irrig = 0.0
        area_flood = 0.0
        !
@@ -7994,8 +8301,4 @@
           ip = irrsub_index(ib, fopt, 1)
           jp = irrsub_index(ib, fopt, 2)
-          !
-          IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN
-             area_irrig = area_irrig + irrsub_area(ib,fopt)*irrigated_frac(ip,jp)
-          ENDIF
           !
           DO itype=1,ntype
@@ -8006,33 +8309,10 @@
        ENDDO
        !
-       ! Put the total irrigated and flooded areas in the output variables
-       !
-       IF ( init_irrig ) THEN
+       ! Put the total flooded areas in the output variables
+       !
+       !
+       IF ( init_flood ) THEN
           ! if we are at the poles resolution(ib,1) = 0
           IF (resolution(ib,1) == 0) THEN
-             ! use pi*resolution(ib,2) to get the disc area
-             resolution_1 = pi*resolution(ib,2)
-          ELSE
-             resolution_1 = resolution(ib,1)
-          END IF
-          irrigated(ib) = MIN(area_irrig, resolution_1*resolution(ib,2)*contfrac(ib))
-          IF ( irrigated(ib) < 0 ) THEN
-             WRITE(numout,*) 'We have a problem here : ', irrigated(ib) 
-             WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)
-             WRITE(numout,*) area_irrig
-             CALL ipslerr_p(3,'routing_irrigmap','Problem with irrigated...','','')
-          ENDIF
-!!$          ! Compute a diagnostic of the map.
-!!$          IF(contfrac(ib).GT.zero) THEN
-!!$             irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) )
-!!$          ELSE
-!!$             irrigmap (ib) = zero
-!!$          ENDIF
-          !
-       ENDIF
-       !
-       IF ( init_flood ) THEN
-          ! if we are at the poles resolution(ib,1) = 0
-          IF (resolution(ib,1) == 0) THEN 
              ! use pi*resolution(ib,2) to get the disc area
              resolution_1 = pi*resolution(ib,2)
@@ -8043,8 +8323,8 @@
                & resolution_1*resolution(ib,2)*contfrac(ib))
           IF ( floodplains(ib) < 0 ) THEN
-             WRITE(numout,*) 'We have a problem here : ', floodplains(ib) 
+             WRITE(numout,*) 'We have a problem here : ', floodplains(ib)
              WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)
              WRITE(numout,*) area_flood
-             CALL ipslerr_p(3,'routing_irrigmap','Problem with floodplains..','','')
+             CALL ipslerr_p(3,'routing_floodmap','Problem with floodplains..','','')
           ENDIF
 !!$          ! Compute a diagnostic of the map.
@@ -8058,6 +8338,6 @@
        IF ( init_swamp ) THEN
           ! if we are at the poles resolution(ib,1) = 0
-          IF (resolution(ib,1) == 0) THEN 
-             ! use pi*resolution(ib,2) to get the disc area 
+          IF (resolution(ib,1) == 0) THEN
+             ! use pi*resolution(ib,2) to get the disc area
              resolution_1 = pi*resolution(ib,2)
           ELSE
@@ -8066,8 +8346,8 @@
           swamp(ib) = MIN(area_flood(iswamp), resolution_1*resolution(ib,2)*contfrac(ib))
           IF ( swamp(ib) < 0 ) THEN
-             WRITE(numout,*) 'We have a problem here : ', swamp(ib) 
+             WRITE(numout,*) 'We have a problem here : ', swamp(ib)
              WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)
              WRITE(numout,*) area_flood
-             CALL ipslerr_p(3,'routing_irrigmap','Problem with swamp...','','')
+             CALL ipslerr_p(3,'routing_floodmap','Problem with swamp...','','')
           ENDIF
 !!$          ! Compute a diagnostic of the map.
@@ -8083,7 +8363,6 @@
     !
     !
-    
+
     IF (printlev>=1) THEN
-       IF ( init_irrig ) WRITE(numout,*) "Diagnostics irrigated :", MINVAL(irrigated), MAXVAL(irrigated)
        IF ( init_flood ) WRITE(numout,*) "Diagnostics floodplains :", MINVAL(floodplains), MAXVAL(floodplains)
        IF ( init_swamp ) WRITE(numout,*) "Diagnostics swamp :", MINVAL(swamp), MAXVAL(swamp)
@@ -8112,5 +8391,5 @@
     DEALLOCATE (latrel)
     !
-  END SUBROUTINE routing_irrigmap
+  END SUBROUTINE routing_floodmap
   !
 !! ================================================================================================================================
@@ -8167,15 +8446,15 @@
     REAL(r_std), SAVE          :: totw_out             !! Sum of the water flow out to the routing scheme
 !$OMP THREADPRIVATE(totw_out)
-    REAL(r_std), SAVE          :: totw_return          !! 
+    REAL(r_std), SAVE          :: totw_return          !!
 !$OMP THREADPRIVATE(totw_return)
-    REAL(r_std), SAVE          :: totw_irrig           !! 
+    REAL(r_std), SAVE          :: totw_irrig           !!
 !$OMP THREADPRIVATE(totw_irrig)
-    REAL(r_std), SAVE          :: totw_river           !! 
+    REAL(r_std), SAVE          :: totw_river           !!
 !$OMP THREADPRIVATE(totw_river)
-    REAL(r_std), SAVE          :: totw_coastal         !! 
+    REAL(r_std), SAVE          :: totw_coastal         !!
 !$OMP THREADPRIVATE(totw_coastal)
     REAL(r_std)                :: totarea              !! Total area of basin (m^2)
     REAL(r_std)                :: area                 !! Total area of routing (m^2)
-    INTEGER(i_std)             :: ig                   !! 
+    INTEGER(i_std)             :: ig                   !!
     !
     ! Just to make sure we do not get too large numbers !
@@ -8194,5 +8473,5 @@
        totw_slow = zero
        totw_lake = zero
-       totw_pond = zero 
+       totw_pond = zero
        totw_in = zero
        !
