Index: /branches/ORCHIDEE_2_2/ORCHIDEE/orchidee.default
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/orchidee.default	(revision 6101)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/orchidee.default	(revision 6102)
@@ -1693,4 +1693,7 @@
 IRRIGATION_FILE =  floodplains.nc 
 
+# ROUTING_METHOD (character string) :  Choice of routing module to be used 	{RIVER_ROUTING=T}
+ROUTING_METHOD =  standard 
+
 # EPS_CARBON ([%]   ) :  Allowed error on carbon stock 	{SPINUP_ANALYTIC}
 EPS_CARBON =  0.01 
@@ -1720,2 +1723,5 @@
 XIOS_INTERPOLATION =  n 
 
+# XIOS_REMAP_OUTPUT ([FLAG]) :  Actiave remaping of diagnostic output files to regular grid 	{XIOS_ORCHIDEE_OK .AND. grid_type=unstructured}
+XIOS_REMAP_OUTPUT =  True 
+
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90	(revision 6101)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90	(revision 6102)
@@ -58,5 +58,5 @@
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: routing_main, routing_initialize, routing_finalize, routing_clear
+  PUBLIC :: routing_main, routing_initialize, routing_finalize, routing_clear, routing_names
 
 !! PARAMETERS
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90	(revision 6102)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90	(revision 6102)
@@ -0,0 +1,2582 @@
+! =================================================================================================================================
+! MODULE       : routing_simple
+!
+!  CONTACT      : orchidee-help _at_ listes.ipsl.fr
+!
+! LICENCE      : IPSL (2006)
+! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
+!
+!>\BRIEF       This module routes the water over the continents into the oceans and computes the water
+!!             stored in floodplains or taken for irrigation.
+!!
+!!\n DESCRIPTION: The subroutines in this subroutine is only called when ROUTING_METHOD=simple is set in run.def.
+!!                The method can be used for regular latitude-longitude grid or for unstructured grid.
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S)	:
+!!
+!! SVN          :
+!! $HeadURL$
+!! $Date$
+!! $Revision$
+!! \n
+!_ ================================================================================================================================
+
+ 
+
+! Histoire Salee
+!---------------
+! La douce riviere
+! Sortant de son lit
+! S'est jetee ma chere
+! dans les bras mais oui
+! du beau fleuve
+!
+! L'eau coule sous les ponts
+! Et puis les flots s'emeuvent
+! - N'etes vous pas au courant ?
+! Il parait que la riviere 
+! Va devenir mer
+!                       Roland Bacri
+!
+
+
+MODULE routing_simple
+
+  USE ioipsl   
+  USE xios_orchidee
+  USE ioipsl_para 
+  USE constantes
+  USE constantes_soil
+  USE pft_parameters
+  USE sechiba_io_p
+  USE interpol_help
+  USE grid
+  USE mod_orchidee_para
+
+
+  IMPLICIT NONE
+  PRIVATE
+  PUBLIC :: routing_simple_main, routing_simple_xios_initialize
+  PUBLIC :: routing_simple_initialize, routing_simple_finalize, routing_simple_clear
+
+  !! PARAMETERS
+  INTEGER(i_std), PARAMETER                 :: nbasmax=5                   !! The maximum number of basins we wish to have per grid box (truncation of the model) (unitless)
+  INTEGER(i_std), SAVE                      :: nbvmax                      !! The maximum number of basins we can handle at any time during the generation of the maps (unitless)
+  !$OMP THREADPRIVATE(nbvmax)
+  REAL(r_std), SAVE                         :: fast_tcst = 3.0             !! Property of the fast reservoir (day/m)
+  !$OMP THREADPRIVATE(fast_tcst)
+  REAL(r_std), SAVE                         :: slow_tcst = 25.0            !! Property of the slow reservoir (day/m)
+  !$OMP THREADPRIVATE(slow_tcst)
+  REAL(r_std), SAVE                         :: stream_tcst = 0.24          !! Property of the stream reservoir (day/m)
+  !$OMP THREADPRIVATE(stream_tcst)
+  REAL(r_std), SAVE                         :: flood_tcst = 4.0            !! Property of the floodplains reservoir (day/m)
+  !$OMP THREADPRIVATE(flood_tcst)
+  REAL(r_std), SAVE                         :: swamp_cst = 0.2             !! Fraction of the river transport that flows to the swamps (unitless;0-1)
+  !$OMP THREADPRIVATE(swamp_cst)
+  !
+  !  Relation between volume and fraction of floodplains
+  !
+  REAL(r_std), SAVE                         :: beta = 2.0                  !! Parameter to fix the shape of the floodplain (>1 for convex edges, <1 for concave edges) (unitless)
+  !$OMP THREADPRIVATE(beta)
+  REAL(r_std), SAVE                         :: betap = 0.5                 !! Ratio of the basin surface intercepted by ponds and the maximum surface of ponds (unitless;0-1)
+  !$OMP THREADPRIVATE(betap)
+  REAL(r_std), SAVE                         :: floodcri = 2000.0           !! Potential height for which all the basin is flooded (mm)
+  !$OMP THREADPRIVATE(floodcri)
+  !
+  !  Relation between maximum surface of ponds and basin surface, and drainage (mm/j) to the slow_res
+  !
+  REAL(r_std), PARAMETER                    :: pond_bas = 50.0             !! [DISPENSABLE] - not used
+  REAL(r_std), SAVE                         :: pondcri = 2000.0            !! Potential height for which all the basin is a pond (mm)
+  !$OMP THREADPRIVATE(pondcri)
+
+  REAL(r_std), PARAMETER                    :: maxevap_lake = 7.5/86400.   !! Maximum evaporation rate from lakes (kg/m^2/s)
+  REAL(r_std),SAVE                          :: dt_routing                  !! Routing time step (s)
+  !$OMP THREADPRIVATE(dt_routing)
+  INTEGER(i_std), SAVE                      :: diagunit = 87               !! Diagnostic file unit (unitless)
+  !$OMP THREADPRIVATE(diagunit)
+
+  ! Logicals to control model configuration
+  !
+  LOGICAL, SAVE                             :: dofloodinfilt = .FALSE.     !! Logical to choose if floodplains infiltration is activated or not (true/false)
+  !$OMP THREADPRIVATE(dofloodinfilt)
+  LOGICAL, SAVE                             :: doswamps = .FALSE.          !! Logical to choose if swamps are activated or not (true/false)
+  !$OMP THREADPRIVATE(doswamps)
+  LOGICAL, SAVE                             :: doponds = .FALSE.           !! Logical to choose if ponds are activated or not (true/false)
+  !$OMP THREADPRIVATE(doponds)
+
+  INTEGER(i_std), SAVE                                       :: num_largest                 !! Number of largest river basins which should be treated as independently as rivers
+                                                                                            !! (not flow into ocean as diffusion coastal flow) (unitless)
+  !$OMP THREADPRIVATE(num_largest)
+  REAL(r_std), SAVE                                          :: time_counter                !! Time counter (s)
+  !$OMP THREADPRIVATE(time_counter)
+  REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)     :: routing_area_loc            !! Surface of basin (m^2)
+  !$OMP THREADPRIVATE(routing_area_loc)
+  REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)     :: topo_resid_loc              !! Topographic index of the retention time (m)
+  !$OMP THREADPRIVATE(topo_resid_loc)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: route_togrid_loc            !! Grid into which the basin flows (unitless)
+  !$OMP THREADPRIVATE(route_togrid_loc)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: route_tobasin_loc           !! Basin in to which the water goes (unitless)
+  !$OMP THREADPRIVATE(route_tobasin_loc)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: route_nbintobas_loc         !! Number of basin into current one (unitless)
+  !$OMP THREADPRIVATE(route_nbintobas_loc)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: global_basinid_loc          !! ID of basin (unitless)
+  !$OMP THREADPRIVATE(global_basinid_loc)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: hydrodiag_loc               !! Variable to diagnose the hydrographs
+  !$OMP THREADPRIVATE(hydrodiag_loc)
+  REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:)       :: hydroupbasin_loc            !! The area upstream of the gauging station (m^2)
+  !$OMP THREADPRIVATE(hydroupbasin_loc)
+
+  ! parallelism
+  REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)     :: routing_area_glo            !! Surface of basin (m^2)
+  !$OMP THREADPRIVATE(routing_area_glo)
+  REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)     :: topo_resid_glo              !! Topographic index of the retention time (m)
+  !$OMP THREADPRIVATE(topo_resid_glo)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: route_togrid_glo            !! Grid into which the basin flows (unitless)
+  !$OMP THREADPRIVATE(route_togrid_glo)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: route_tobasin_glo           !! Basin in to which the water goes (unitless)
+  !$OMP THREADPRIVATE(route_tobasin_glo)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: route_nbintobas_glo         !! Number of basin into current one (unitless)
+  !$OMP THREADPRIVATE(route_nbintobas_glo)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: global_basinid_glo          !! ID of basin (unitless)
+  !$OMP THREADPRIVATE(global_basinid_glo)
+  INTEGER(i_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:,:)  :: hydrodiag_glo               !! Variable to diagnose the hydrographs
+  !$OMP THREADPRIVATE(hydrodiag_glo)
+  REAL(r_std), SAVE, ALLOCATABLE, TARGET, DIMENSION(:)       :: hydroupbasin_glo            !! The area upstream of the gauging station (m^2)
+  !$OMP THREADPRIVATE(hydroupbasin_glo)
+
+  REAL(r_std), SAVE, POINTER, DIMENSION(:,:)                 :: routing_area                !! Surface of basin (m^2)
+  !$OMP THREADPRIVATE(routing_area)
+  REAL(r_std), SAVE, POINTER, DIMENSION(:,:)                 :: topo_resid                  !! Topographic index of the retention time (m)
+  !$OMP THREADPRIVATE(topo_resid)
+  INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:)              :: route_togrid                !! Grid into which the basin flows (unitless)
+  !$OMP THREADPRIVATE(route_togrid)
+  INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:)              :: route_tobasin               !! Basin in to which the water goes (unitless)
+  !$OMP THREADPRIVATE(route_tobasin)
+  INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:)              :: route_nbintobas             !! Number of basin into current one (unitless)
+  !$OMP THREADPRIVATE(route_nbintobas)
+  INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:)              :: global_basinid              !! ID of basin (unitless)
+  !$OMP THREADPRIVATE(global_basinid)
+  INTEGER(i_std), SAVE, POINTER, DIMENSION(:,:)              :: hydrodiag                   !! Variable to diagnose the hydrographs
+  !$OMP THREADPRIVATE(hydrodiag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: slowflow_diag               !! Diagnostic slow flow hydrographs (kg/dt)
+  !$OMP THREADPRIVATE(slowflow_diag)  
+  REAL(r_std), SAVE, POINTER, DIMENSION(:)                   :: hydroupbasin                !! The area upstream of the gauging station (m^2)
+  !$OMP THREADPRIVATE(hydroupbasin)
+
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrigated                   !! Area equipped for irrigation in each grid box (m^2)
+  !$OMP THREADPRIVATE(irrigated)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: floodplains                 !! Maximal surface which can be inundated in each grid box (m^2)
+  !$OMP THREADPRIVATE(floodplains)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: swamp                       !! Maximal surface of swamps in each grid box (m^2)
+  !$OMP THREADPRIVATE(swamp)
+  !
+  ! The reservoirs, variables kept in the restart file.
+  !
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)             :: fast_reservoir              !! Water amount in the fast reservoir (kg)
+  !$OMP THREADPRIVATE(fast_reservoir)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)             :: slow_reservoir              !! Water amount in the slow reservoir (kg)
+  !$OMP THREADPRIVATE(slow_reservoir)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)             :: stream_reservoir            !! Water amount in the stream reservoir (kg)
+  !$OMP THREADPRIVATE(stream_reservoir)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)             :: flood_reservoir             !! Water amount in the floodplains reservoir (kg)
+  !$OMP THREADPRIVATE(flood_reservoir)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: lake_reservoir              !! Water amount in the lake reservoir (kg)
+  !$OMP THREADPRIVATE(lake_reservoir)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: pond_reservoir              !! Water amount in the pond reservoir (kg)
+  !$OMP THREADPRIVATE(pond_reservoir)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)             :: flood_frac_bas              !! Flooded fraction per basin (unitless;0-1)
+  !$OMP THREADPRIVATE(flood_frac_bas)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: pond_frac                   !! Pond fraction per grid box (unitless;0-1)
+  !$OMP THREADPRIVATE(pond_frac)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: flood_height                !! Floodplain height (mm)
+  !$OMP THREADPRIVATE(flood_height)
+  !
+  ! The accumulated fluxes.
+  !
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: floodout_mean               !! Accumulated flow out of floodplains (kg/m^2/dt)
+  !$OMP THREADPRIVATE(floodout_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: runoff_mean                 !! Accumulated runoff (kg/m^2/dt)
+  !$OMP THREADPRIVATE(runoff_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: drainage_mean               !! Accumulated drainage (kg/m^2/dt)
+  !$OMP THREADPRIVATE(drainage_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: transpot_mean               !! Mean potential transpiration from the plants (kg/m^2/dt)
+  !$OMP THREADPRIVATE(transpot_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: precip_mean                 !! Accumulated precipitation (kg/m^2/dt)
+  !$OMP THREADPRIVATE(precip_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: humrel_mean                 !! Mean soil moisture stress, mean root extraction potential (unitless)
+  !$OMP THREADPRIVATE(humrel_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: totnobio_mean               !! Mean last total fraction of no bio (unitless;0-1)
+  !$OMP THREADPRIVATE(totnobio_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: vegtot_mean                 !! Mean potentially vegetated fraction (unitless;0-1)
+  !$OMP THREADPRIVATE(vegtot_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: k_litt_mean                 !! Mean averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
+  !$OMP THREADPRIVATE(k_litt_mean)
+  !
+  ! The averaged outflow fluxes.
+  !
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: lakeinflow_mean              !! Mean lake inflow (kg/m^2/dt)
+  !$OMP THREADPRIVATE(lakeinflow_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: returnflow_mean              !! Mean water flow from lakes and swamps which returns to the grid box.
+  !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
+  !$OMP THREADPRIVATE(returnflow_mean)
+  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(:)               :: 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(:)               :: 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)
+  !$OMP THREADPRIVATE(riverflow_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: coastalflow_mean             !! Mean outflow on coastal points by small basins.
+  !! This is the water which flows in a disperse way into the ocean (kg/dt)
+  !$OMP THREADPRIVATE(coastalflow_mean)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: floodtemp                    !! Temperature to decide if floodplains work (K)
+  !$OMP THREADPRIVATE(floodtemp)
+  INTEGER(i_std), SAVE                                       :: floodtemp_lev                !! Temperature level to decide if floodplains work (K)
+  !$OMP THREADPRIVATE(floodtemp_lev)
+  !
+  ! Diagnostic variables ... well sort of !
+  !
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: irrig_netereq                !! Irrigation requirement (water requirements by the crop for its optimal growth (kg/m^2/dt)
+  !$OMP THREADPRIVATE(irrig_netereq)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: hydrographs                  !! Hydrographs at the outflow of the grid box for major basins (kg/dt)
+  !$OMP THREADPRIVATE(hydrographs)
+  !
+  ! Diagnostics for the various reservoirs we use (Kg/m^2)
+  !
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: fast_diag                    !! Diagnostic for the fast reservoir (kg/m^2)
+  !$OMP THREADPRIVATE(fast_diag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: slow_diag                    !! Diagnostic for the slow reservoir (kg/m^2)
+  !$OMP THREADPRIVATE(slow_diag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: stream_diag                  !! Diagnostic for the stream reservoir (kg/m^2)
+  !$OMP THREADPRIVATE(stream_diag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: flood_diag                   !! Diagnostic for the floodplain reservoir (kg/m^2)
+  !$OMP THREADPRIVATE(flood_diag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: pond_diag                    !! Diagnostic for the pond reservoir (kg/m^2)
+  !$OMP THREADPRIVATE(pond_diag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: lake_diag                    !! Diagnostic for the lake reservoir (kg/m^2)
+  !$OMP THREADPRIVATE(lake_diag)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: delsurfstor                  !! Diagnostic of the change in surface water storage (flood, pond and lake reservoirs) (kg/m^2)
+  !$OMP THREADPRIVATE(delsurfstor)
+  
+  !
+  ! Specific variables for simple routing
+  !
+  REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:)                !! Topographic index of the retention time (m) index - (local routing grid)
+  !$OMP THREADPRIVATE(topoind_r)
+  INTEGER,SAVE,ALLOCATABLE     :: route_flow_rp1(:)           !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo)   
+  !$OMP THREADPRIVATE(route_flow_rp1)
+  REAL(r_std),SAVE,ALLOCATABLE :: fast_reservoir_r(:)         !! Water amount in the fast reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(fast_reservoir_r)
+  REAL(r_std),SAVE,ALLOCATABLE :: slow_reservoir_r(:)         !! Water amount in the slow reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(slow_reservoir_r)
+  REAL(r_std),SAVE,ALLOCATABLE :: stream_reservoir_r(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(stream_reservoir_r)
+  LOGICAL,SAVE,ALLOCATABLE     :: is_lakeinflow_r(:)          !! is lake inflow point  - (local routing grid)
+  !$OMP THREADPRIVATE(is_lakeinflow_r)
+  LOGICAL,SAVE,ALLOCATABLE     :: is_coastalflow_r(:)         !! is coastal flow point - (local routing grid)
+  !$OMP THREADPRIVATE(is_coastalflow_r)
+  LOGICAL,SAVE,ALLOCATABLE     :: is_riverflow_r(:)           !! is river flow point - (local routing grid) 
+  !$OMP THREADPRIVATE(is_riverflow_r)
+  LOGICAL,SAVE,ALLOCATABLE     :: is_streamflow_r(:)          !! is stream flow point - (local routing grid) 
+  !$OMP THREADPRIVATE(is_streamflow_r)
+  LOGICAL,SAVE,ALLOCATABLE     :: routing_mask_r(:)           !! valid routing point - (local routing grid) 
+  !$OMP THREADPRIVATE(routing_mask_r)
+  LOGICAL,SAVE,ALLOCATABLE     :: coast_mask(:)               !! is a coast point - (local native grid)
+  !$OMP THREADPRIVATE(coast_mask)
+  INTEGER,SAVE                 :: nbpt_r                      !! number of point in local routing grid
+  !$OMP THREADPRIVATE(nbpt_r)
+  INTEGER,SAVE                 :: nbpt_rp1                    !! number of point in local routing grid with halo of 1
+  !$OMP THREADPRIVATE(nbpt_rp1)
+  REAL(r_std),SAVE,ALLOCATABLE :: routing_weight(:)           !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
+  !$OMP THREADPRIVATE(routing_weight)
+  REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_r(:)         !! Weight to transfer lost water into intersecting cell when doing interpolation from local routing grid to local native grid (lakeinflow)
+  !$OMP THREADPRIVATE(routing_weight_r)
+  REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_coast_r(:)   !! Weight to transfer lost water into intersecting cell on coast 
+  !$OMP THREADPRIVATE(routing_weight_coast_r)
+  !! when doing interpolation from local routing grid to local native grid (river flow+coastal flow)
+  REAL(r_std),SAVE,ALLOCATABLE :: basins_extended_r(:)        !! basins riverflow id (local routing grid) 
+  !$OMP THREADPRIVATE(basins_extended_r)
+  INTEGER(i_std)               :: basins_count                !! number of basins (local routing grid) 
+  !$OMP THREADPRIVATE(basins_count)
+  INTEGER(i_std)               :: basins_out                  !! number of basins to output for diag 
+  !$OMP THREADPRIVATE(basins_out)
+
+  INTEGER(i_std)               :: split_routing               !! time spliting for routing
+  !$OMP THREADPRIVATE(split_routing)
+
+
+  INTEGER(i_std), PARAMETER :: nb_stations=14
+  REAL(r_std),PARAMETER  :: station_lon(nb_stations) = &
+       (/ -162.8830, -90.9058, -55.5110, -49.3242, -133.7447, -63.6000,  28.7167, &
+       15.3000,  66.5300,  89.6700,  86.5000,  127.6500,   3.3833, 117.6200 /)
+  REAL(r_std),PARAMETER  :: station_lat(nb_stations) = &
+       (/  61.9340,  32.3150, -1.9470, -5.1281, 67.4583,  8.1500, 45.2167, &
+       -4.3000,  66.5700, 25.1800, 67.4800, 70.7000, 11.8667, 30.7700 /)
+  CHARACTER(LEN=17),PARAMETER  :: station_name(nb_stations) = &
+       (/ "Pilot station    ", "Vicksburg        ", "Obidos           ", &
+       "Itupiranga       ", "Arctic red river ", "Puente Angostura ", &
+       "Ceatal Izmail    ", "Kinshasa         ", "Salekhard        ", &
+       "Bahadurabad      ", "Igarka           ", "Kusur            ", &
+       "Malanville       ", "Datong           " /)
+
+CONTAINS
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_simple_xios_initialize
+!!
+!>\BRIEF	  Initialize xios dependant defintion before closing context defintion
+!!
+!! DESCRIPTION:	  Initialize xios dependant defintion before closing context defintion.
+!!                This subroutine is called before the xios context is closed. 
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+
+  SUBROUTINE routing_simple_xios_initialize
+    USE xios
+    USE routing, ONLY : routing_names
+    IMPLICIT NONE     
+
+    !! 0 Variable and parameter description
+    CHARACTER(LEN=60),ALLOCATABLE :: label(:)
+    LOGICAL :: file_exists
+
+    IF (is_omp_root) THEN   
+       CALL xios_get_axis_attr("basins", n_glo=basins_out)    ! get nb basins to output
+       ALLOCATE(label(basins_out))
+       CALL routing_names(basins_out,label)
+       CALL xios_set_axis_attr("basins", label=label)         ! set riverflow basins name
+       INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
+       IF (file_exists) CALL xios_set_file_attr("routing_start", enabled=.TRUE.)  
+    ENDIF
+
+  END SUBROUTINE routing_simple_xios_initialize
+
+
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_simple_init_2
+!!
+!>\BRIEF	  Privat subroutine
+!!
+!! DESCRIPTION:	  Privat subroutine to the module routing_simple. This subroutine is called in the end of routing_simple_initialize
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+
+  SUBROUTINE routing_simple_init_2(nbpt, contfrac)
+    USE xios
+    USE grid, ONLY : area
+    IMPLICIT NONE
+    INCLUDE "mpif.h"
+
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+    INTEGER,INTENT(in) :: nbpt             !! nb points native grid
+    REAL,INTENT(in)    :: contfrac(nbpt)   !! fraction of land
+
+    !! 0.2 Local variables
+    INTEGER :: ni                                         !! longitude dimension of local routing grid
+    INTEGER :: nj                                         !! latitude dimension of local routing grid
+
+    REAL(r_std),ALLOCATABLE :: trip_rp1(:,: )             !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid + halo (0:ni+1,0:nj+1)
+    REAL(r_std),ALLOCATABLE :: trip_extended_r(:,:)       !! direction of flow (1-8) or river flow (99) or coastal flow (98) or lake inflow (97) - local routing grid (ni,nj)
+    !! routing is artificially computed on sea and endoric basins
+    LOGICAL,ALLOCATABLE     :: routing_mask_rp1(:,:)      !! valid routing point - local routing grid+halo (0:ni+1,0:nj+1) 
+    REAL(r_std),ALLOCATABLE :: mask_native(:)             !! some mask on native grid (nbpt)
+    REAL(r_std),ALLOCATABLE :: frac_routing_r(:)          !! fraction of routing cell intesected by cells of native grid
+    REAL(r_std),ALLOCATABLE :: frac_routing_coast_r(:)    !! fraction of routing cell intesected by coastal cells of native grid
+    INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir
+    REAL(r_std) :: sum1,sum2,sum_frac_routing_r
+    REAL(r_std) :: contfrac_mpi(nbp_mpi)
+    REAL(r_std) :: area_mpi(nbp_mpi)
+    INTEGER :: basins_count_mpi
+    INTEGER :: ierr
+    LOGICAL :: file_exists
+
+!_ ================================================================================================================================
+
+    split_routing=1
+    CALL getin_p("SPLIT_ROUTING",split_routing)
+
+    CALL gather_omp(contfrac,contfrac_mpi)
+    CALL gather_omp(area, area_mpi)
+
+    IF (is_omp_root) THEN
+
+       CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)    ! get routing domain dimension
+
+       nbpt_r= ni*nj                                                
+       nbpt_rp1= (ni+2)*(nj+2)
+
+       ALLOCATE(trip_extended_r(ni,nj))
+       ALLOCATE(topoind_r(ni*nj))       
+       ALLOCATE(is_lakeinflow_r(ni*nj))       
+       ALLOCATE(is_coastalflow_r(ni*nj))       
+       ALLOCATE(is_riverflow_r(ni*nj)) 
+       ALLOCATE(is_streamflow_r(ni*nj)) 
+       ALLOCATE(routing_mask_r(ni*nj)) 
+       ALLOCATE(fast_reservoir_r(ni*nj))       
+       ALLOCATE(slow_reservoir_r(ni*nj))       
+       ALLOCATE(stream_reservoir_r(ni*nj))       
+
+       ALLOCATE(mask_native(nbp_mpi))
+       ALLOCATE(coast_mask(nbp_mpi))
+       ALLOCATE(frac_routing_r(nbpt_r))
+       ALLOCATE(frac_routing_coast_r(nbpt_r))
+       ALLOCATE(routing_weight(nbp_mpi))
+       ALLOCATE(routing_weight_r(nbpt_r))
+       ALLOCATE(routing_weight_coast_r(nbpt_r))
+       ALLOCATE(basins_extended_r(nbpt_r))
+
+       coast_mask=.FALSE.
+       WHERE( contfrac_mpi(:)< 1.-1.e-5) coast_mask(:)=.TRUE.            ! create mask for coastal cells on native grid
+
+       INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
+
+       IF (file_exists) THEN  
+          CALL xios_recv_field("fast_reservoir_start",fast_reservoir_r)
+          CALL xios_recv_field("slow_reservoir_start",slow_reservoir_r)
+          CALL xios_recv_field("stream_reservoir_start",stream_reservoir_r)
+       ELSE
+          fast_reservoir_r(:)=0
+          slow_reservoir_r(:)=0
+          stream_reservoir_r(:)=0
+       ENDIF
+
+       ALLOCATE(trip_rp1(0:ni+1,0:nj+1))
+       trip_rp1(:,:)=1e10
+       ALLOCATE(routing_mask_rp1(0:ni+1,0:nj+1))
+       ALLOCATE(route_flow_rp1((ni+2)*(nj+2)))
+
+       CALL xios_send_field("routing_contfrac",contfrac_mpi)       ! diag
+       CALL xios_recv_field("trip_r",trip_rp1(1:ni,1:nj))          ! recv trip array with halo of 1
+       CALL xios_recv_field("trip_extended_r",trip_extended_r)     ! recv extended trip array from file
+       CALL xios_recv_field("topoind_r",topoind_r)                 ! recv topo index array from file
+       CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file
+
+       is_lakeinflow_r(:) = .FALSE.
+       is_coastalflow_r(:) = .FALSE.
+       is_riverflow_r(:) = .FALSE.
+
+       route_flow_rp1(:)=-1
+
+       routing_mask_rp1=.TRUE.
+       DO j=1,nj
+          DO i=1,ni
+             IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE.         ! if ocean or endoric basins => masked point
+          ENDDO
+       ENDDO
+
+       mask_native(:)=0
+       frac_routing_r(:)=0
+       WHERE( .NOT. coast_mask(:)) mask_native(:)=1
+       CALL xios_send_field("mask_native_lake",mask_native)              ! send full land point to XIOS (native grid)
+       CALL xios_recv_field("frac_routing_lake_r",frac_routing_r)        ! receive fraction of intersected cell by full land, on routing grid 
+
+
+       DO j=1,nj
+          DO i=1,ni
+             ij=i+ni*(j-1)
+             IF (frac_routing_r(ij)>1-1e-5) THEN                                                       ! if full land point (no part of coast)
+                IF ((.NOT. routing_mask_rp1(i,j)) .OR. trip_rp1(i,j)==98 .OR. trip_rp1(i,j)==99) THEN   ! if masked point or river flow or coastal flow
+                   routing_mask_rp1(i,j)=.TRUE.                                                          ! transform to non masked
+                   trip_rp1(i,j)=trip_extended_r(i,j)                                                  ! replace by extended trip value         
+                ENDIF
+             ENDIF
+          ENDDO
+       ENDDO
+
+       mask_native(:)=0
+       frac_routing_coast_r(:)=0
+       WHERE( coast_mask(:)) mask_native(:)=1
+       CALL xios_send_field("mask_native_coast",mask_native)               ! send only coastal point to XIOS (native grid)
+       CALL xios_recv_field("frac_routing_coast_r",frac_routing_coast_r)   ! receive fraction of intersected cell by coastal point, on routing grid 
+
+       DO j=1,nj
+          DO i=1,ni
+             ij=i+ni*(j-1)
+             IF (frac_routing_coast_r(ij) > 0) THEN                        ! fraction of coastal point ?
+                IF (.NOT. routing_mask_rp1(i,j)) THEN                      ! if masked, transform to coastal flow to conserve the water balance
+                   routing_mask_rp1(i,j)=.TRUE.
+                   trip_rp1(i,j)=98     ! transform to coastal flow
+                   topoind_r(ij)=1e-10  ! no residence time
+                ENDIF
+             ENDIF
+          ENDDO
+       ENDDO
+
+
+       mask_native(:)=1
+       frac_routing_r(:)=0
+       CALL xios_send_field("mask_native",mask_native)              ! send only all point to XIOS (native grid)
+       CALL xios_recv_field("frac_routing_r",frac_routing_r)        ! receive fraction of intersected cell, (routing grid)
+       DO j=1,nj
+          DO i=1,ni
+             ij=i+ni*(j-1)
+             IF (frac_routing_r(ij)>0) THEN
+                IF (.NOT. routing_mask_rp1(i,j)) THEN
+                   routing_mask_rp1(i,j)=.TRUE.                         ! may never happen
+                ENDIF
+             ELSE
+                routing_mask_rp1(i,j)=.FALSE.                           ! if no intersection from native cells, mask the point (no routage) 
+                trip_rp1(i,j)=1e10
+             ENDIF
+          ENDDO
+       ENDDO
+
+       CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj))        ! send to xios trip array to update for halo
+       CALL xios_recv_field("trip_rp1",trip_rp1)                        ! recv trip array with halo of 1
+
+
+       routing_mask_rp1=.TRUE.
+       DO j=0,nj+1
+          DO i=0,ni+1
+             IF (trip_rp1(i,j)>99) routing_mask_rp1(i,j)=.FALSE.        !  update mask with halo of 1
+          ENDDO
+       ENDDO
+
+       !! Compute the routing
+       !! Loop on all point point of the local routing grid + halo
+       DO j=0,nj+1                                                 
+          jp1=j+1
+          jm1=j-1
+          DO i=0,ni+1
+             ij_rp1=i+(ni+2)*j+1
+             ij_r=i+ni*(j-1)
+             ip1=i+1
+             im1=i-1
+
+             IF (trip_rp1(i,j) < 100) THEN                    
+
+                ir=-1 ; jr=-1  !  -> -1 for 97,98,99
+                SELECT CASE (NINT(trip_rp1(i,j)))                     ! get the trip value for each points
+                CASE (1)
+                   jr=jm1 ; ir=i    ! north
+                CASE (2)
+                   jr=jm1 ; ir=ip1  ! north-east
+                CASE (3)
+                   jr=j ;   ir=ip1  ! east
+                CASE (4)
+                   jr=jp1 ;  ir=ip1 ! south-east
+                CASE (5)
+                   jr=jp1 ;  ir=i   ! south
+                CASE(6)
+                   jr=jp1 ;  ir=im1 ! south-west
+                CASE (7)
+                   jr=j ;   ir=im1  ! west
+                CASE (8)
+                   jr=jm1 ;  ir=im1 ! north-west
+                CASE (97)
+                   IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1)  THEN         ! if inside my local domain
+                      is_lakeinflow_r(ij_r)=.TRUE.                             ! I am a lakeinflow point and route to myself
+                      jr=j ;   ir=i                                
+                   ENDIF
+                CASE (98)
+                   IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN         ! if inside my local domain
+                      is_coastalflow_r(ij_r)=.TRUE.                           ! I am a coastal flow point and route to myself           
+                      jr=j ;   ir=i
+                   ENDIF
+                CASE (99)
+                   IF ( i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN         ! if inside my local domain
+                      is_riverflow_r(ij_r)=.TRUE.                             ! I am a riverflow point and route to myself
+                      jr=j ;   ir=i
+                   ENDIF
+                END SELECT
+
+                IF (ir<0 .OR. ir>ni+1 .OR. jr<0 .OR. jr>nj+1) THEN
+                   route_flow_rp1(ij_rp1)=-1                                     ! if route outside my local domain+halo, no routing (will be done by other process)
+                ELSE 
+
+                   IF ( .NOT. routing_mask_rp1(ir,jr)) THEN                      ! if routing to a masked point, cut the flow
+                      jr=j ; ir=i ;                                              ! and route to myself
+                      IF (i>0 .AND. i<ni+1 .AND. j>0 .AND. j<nj+1) THEN          ! if some flow come from neighbour cell 
+                         IF (     trip_rp1(i-1,j-1)==4 .OR. trip_rp1(i,j-1)==5 .OR. trip_rp1(i+1,j)==6     &
+                              .OR. trip_rp1(i-1,j)==3   .OR. trip_rp1(i+1,j)==7                             & 
+                              .OR. trip_rp1(i-1,j+1)==2 .OR. trip_rp1(i-1,j+1)==1 .OR. trip_rp1(i+1,j+1)==8) THEN
+                            is_riverflow_r(ij_r)=.TRUE.                          !  => transform to riverflow
+                         ELSE
+                            is_coastalflow_r(ij_r)=.TRUE.                        !  else transform to coastalflow
+                         ENDIF
+                      ENDIF
+                   ENDIF
+
+                   IF (ir<1 .OR. ir>ni .OR. jr<1 .OR. jr>nj) THEN  
+                      route_flow_rp1(ij_rp1)=-1                                  ! if route outside my local domain, no routing (will be done by other process)
+                   ELSE
+                      route_flow_rp1(ij_rp1)=ir+ni*(jr-1)                        ! define the cell where to flow
+                   ENDIF
+                ENDIF
+
+             ENDIF
+          ENDDO
+       ENDDO
+
+       routing_mask_r(:)=reshape(routing_mask_rp1(1:ni,1:nj),(/ni*nj/))
+
+       is_streamflow_r(:)=.NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:)
+
+       DO ij=1,nbp_mpi
+          routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij)  
+       ENDDO
+
+       routing_weight_r(:)=0
+       DO ij=1,nbpt_r
+          IF (frac_routing_r(ij)>0) routing_weight_r(ij)=1./frac_routing_r(ij)
+       ENDDO
+
+       routing_weight_coast_r(:)=0.
+       DO ij=1,nbpt_r
+          IF (frac_routing_coast_r(ij)>0) routing_weight_coast_r(ij)=1./frac_routing_coast_r(ij)
+       ENDDO
+
+       CALL xios_send_field("routing_weight_coast_r",routing_weight_coast_r)
+       ! looking for basins
+       basins_count_mpi=0
+       DO ij=1,nbpt_r
+          IF (basins_extended_r(ij) > basins_count_mpi) basins_count_mpi=basins_extended_r(ij)
+       ENDDO
+       CALL MPI_ALLREDUCE(basins_count_mpi,basins_count,1,MPI_INT_ORCH, MPI_MAX,MPI_COMM_ORCH,ierr)
+
+    ELSE
+
+       nbpt_r=0
+
+    ENDIF
+
+  END SUBROUTINE routing_simple_init_2
+
+
+  !! ================================================================================================================================
+  !! SUBROUTINE 	: routing_simple_flow
+  !!
+  !>\BRIEF         This subroutine computes the transport of water in the various reservoirs
+  !!                (including ponds and floodplains) and the water withdrawals from the reservoirs for irrigation.
+  !!
+  !! 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 
+  !! 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 
+  !! 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.
+  !!
+  !! RECENT CHANGE(S): None
+  !!
+  !! MAIN OUTPUT VARIABLE(S): lakeinflow, returnflow, reinfiltration, irrigation, riverflow, coastalflow, hydrographs, flood_frac, flood_res
+  !!
+  !! REFERENCES   :
+  !! - Ngo-Duc, T., K. Laval, G. Ramillien, J. Polcher, and A. Cazenave (2007)
+  !!   Validation of the land water storage simulated by Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) with Gravity Recovery and Climate Experiment (GRACE) data.
+  !!   Water Resour. Res., 43, W04427, doi:10.1029/2006WR004941.
+  !! * Irrigation:
+  !! - de Rosnay, P., J. Polcher, K. Laval, and M. Sabre (2003)
+  !!   Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula.
+  !!   Geophys. Res. Lett., 30(19), 1986, doi:10.1029/2003GL018024.
+  !! - A.C. Vivant (2003)
+  !!   Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte.
+  !!   , , 51pp.
+  !! - N. Culson (2004)
+  !!   Impact de l'irrigation sur le cycle de l'eau
+  !!   Master thesis, Paris VI University, 55pp.
+  !! - X.-T. Nguyen-Vinh (2005)
+  !!   Analyse de l'impact de l'irrigation en Amerique du Nord - plaine du Mississippi - sur la climatologie regionale
+  !!   Master thesis, Paris VI University, 33pp.
+  !! - M. Guimberteau (2006)
+  !!   Analyse et modifications proposees de la modelisation de l'irrigation dans un modele de surface.
+  !!   Master thesis, Paris VI University, 46pp.
+  !! - Guimberteau M. (2010)
+  !!   Modelisation de l'hydrologie continentale et influences de l'irrigation sur le cycle de l'eau.
+  !!   Ph.D. thesis, Paris VI University, 195pp.
+  !! - Guimberteau M., Laval K., Perrier A. and Polcher J. (2011).
+  !!   Global effect of irrigation and its impact on the onset of the Indian summer monsoon.
+  !!   In press, Climate Dynamics, doi: 10.1007/s00382-011-1252-5.
+  !! * Floodplains:
+  !! - A.C. Vivant (2002)
+  !!   L'ecoulement lateral de l'eau sur les surfaces continentales. Prise en compte des plaines d'inondations dans ORCHIDEE.
+  !!   Master thesis, Paris VI University, 46pp.
+  !! - A.C. Vivant (2003)
+  !!   Les plaines d'inondations et l'irrigation dans ORCHIDEE, impacts de leur prise en compte.
+  !!   , , 51pp.
+  !! - T. d'Orgeval (2006)
+  !!   Impact du changement climatique sur le cycle de l'eau en Afrique de l'Ouest: modelisation et incertitudes.
+  !!   Ph.D. thesis, Paris VI University, 188pp.
+  !! - T. d'Orgeval, J. Polcher, and P. de Rosnay (2008)
+  !!   Sensitivity of the West African hydrological cycle in ORCHIDEE to infiltration processes.
+  !!   Hydrol. Earth Syst. Sci., 12, 1387-1401
+  !! - M. Guimberteau, G. Drapeau, J. Ronchail, B. Sultan, J. Polcher, J.-M. Martinez, C. Prigent, J.-L. Guyot, G. Cochonneau,
+  !!   J. C. Espinoza, N. Filizola, P. Fraizy, W. Lavado, E. De Oliveira, R. Pombosa, L. Noriega, and P. Vauchel (2011)
+  !!   Discharge simulation in the sub-basins of the Amazon using ORCHIDEE forced by new datasets.
+  !!   Hydrol. Earth Syst. Sci. Discuss., 8, 11171-11232, doi:10.5194/hessd-8-11171-2011
+  !!
+  !! FLOWCHART    :None
+  !! \n
+  !_ ================================================================================================================================
+
+  SUBROUTINE routing_simple_flow(nbpt, dt_routing, lalo, floodout, runoff_omp, drainage_omp, &
+                                 vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, &
+                                 lakeinflow_omp, returnflow, reinfiltration, irrigation, riverflow_omp, &
+                                 coastalflow_omp, hydrographs, slowflow_diag, flood_frac, flood_res)
+    
+    USE xios
+    USE grid, ONLY : area
+    IMPLICIT NONE
+    INCLUDE "mpif.h"
+    
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+
+    INTEGER(i_std), INTENT(in)                   :: nbpt                      !! Domain size (unitless)
+    REAL(r_std), INTENT (in)                     :: dt_routing                !! Routing time step (s)
+    REAL(r_std), INTENT(in)                      :: lalo(nbpt,2)              !! Vector of latitude and longitudes
+    REAL(r_std), INTENT(in)                      :: runoff_omp(nbpt)              !! Grid-point runoff (kg/m^2/dt)
+    REAL(r_std), INTENT(in)                      :: floodout(nbpt)            !! Grid-point flow out of floodplains (kg/m^2/dt)
+    REAL(r_std), INTENT(in)                      :: drainage_omp(nbpt)            !! Grid-point drainage (kg/m^2/dt)
+    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
+    REAL(r_std), INTENT(in)                      :: transpot_mean(nbpt)       !! Mean potential transpiration of the vegetation (kg/m^2/dt)
+    REAL(r_std), INTENT(in)                      :: precip(nbpt)              !! Rainfall (kg/m^2/dt)
+    REAL(r_std), INTENT(in)                      :: humrel(nbpt)              !! Soil moisture stress, root extraction potential (unitless)
+    REAL(r_std), INTENT(in)                      :: k_litt(nbpt)              !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
+    REAL(r_std), INTENT(in)                      :: floodtemp(nbpt)           !! Temperature to decide if floodplains work (K)
+    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)
+
+    !! 0.2 Output variables
+
+    REAL(r_std), INTENT(out)                     :: lakeinflow_omp(nbpt)          !! Water inflow to the lakes (kg/dt)
+    REAL(r_std), INTENT(out)                     :: returnflow(nbpt)          !! The water flow from lakes and swamps which returns into the grid box.
+    !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
+    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)
+    REAL(r_std), INTENT(out)                     :: riverflow_omp(nbpt)           !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
+    REAL(r_std), INTENT(out)                     :: coastalflow_omp(nbpt)         !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
+    REAL(r_std), INTENT(out)                     :: hydrographs(nbpt)         !! Hydrographs at the outflow of the grid box for major basins (kg/dt)
+    REAL(r_std), INTENT(out)                     :: slowflow_diag(nbpt)       !! Hydrographs of slow_flow = routed slow_flow for major basins (kg/dt)
+    REAL(r_std), INTENT(out)                     :: flood_frac(nbpt)          !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(out)                     :: flood_res(nbpt)           !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+    !! 0.4 Local variables
+    REAL(r_std)                                  :: runoff(nbp_mpi)              !! Grid-point runoff (kg/dt)
+    REAL(r_std)                                  :: drainage(nbp_mpi)            !! Grid-point drainage (kg/dt)
+    REAL(r_std)                                  :: riverflow(nbp_mpi)             
+    REAL(r_std)                                  :: coastalflow(nbp_mpi)           
+    REAL(r_std)                                  :: lakeinflow(nbp_mpi)            
+    REAL(r_std)                                  :: fast_diag_mpi(nbp_mpi)          
+    REAL(r_std)                                  :: slow_diag_mpi(nbp_mpi)        
+    REAL(r_std)                                  :: stream_diag_mpi(nbp_mpi)         
+
+    ! from input model -> routing_grid
+    REAL(r_std)                      :: runoff_r(nbpt_r)              !! Grid-point runoff (kg/m^2/dt)
+    REAL(r_std)                      :: floodout_r(nbpt_r)            !! Grid-point flow out of floodplains (kg/m^2/dt)
+    REAL(r_std)                      :: drainage_r(nbpt_r)            !! Grid-point drainage (kg/m^2/dt)
+    REAL(r_std)                      :: vegtot_r(nbpt_r)              !! Potentially vegetated fraction (unitless;0-1)
+    REAL(r_std)                      :: totnobio_r(nbpt_r)            !! Other areas which can not have vegetation
+    REAL(r_std)                      :: transpot_mean_r(nbpt_r)       !! Mean potential transpiration of the vegetation (kg/m^2/dt)
+    REAL(r_std)                      :: precip_r(nbpt_r)              !! Rainfall (kg/m^2/dt)
+    REAL(r_std)                      :: humrel_r(nbpt_r)              !! Soil moisture stress, root extraction potential (unitless)
+    REAL(r_std)                      :: k_litt_r(nbpt_r)              !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
+    REAL(r_std)                      :: floodtemp_r(nbpt_r)           !! Temperature to decide if floodplains work (K)
+    REAL(r_std)                      :: reinf_slope_r(nbpt_r)         !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1)
+
+
+
+    REAL(r_std), DIMENSION(nbpt_r)        :: fast_flow_r                 !! Outflow from the fast reservoir (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: slow_flow_r                 !! Outflow from the slow reservoir (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: stream_flow_r               !! Outflow from the stream reservoir (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: hydrographs_r                !! hydrograph (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: flood_flow_r                !! Outflow from the floodplain reservoir (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: pond_inflow_r               !! Inflow to the pond reservoir (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: pond_drainage_r             !! Drainage from pond (kg/m^2/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: flood_drainage_r            !! Drainage from floodplains (kg/m^2/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: return_swamp_r              !! Inflow to the swamp (kg/dt)
+    !
+    ! Irrigation per basin
+    !
+    REAL(r_std), DIMENSION(nbpt_r)        :: irrig_needs_r               !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
+    REAL(r_std), DIMENSION(nbpt_r)        :: irrig_actual_r              !! Possible irrigation according to the water availability in the reservoirs (kg)
+    REAL(r_std), DIMENSION(nbpt_r)        :: irrig_deficit_r             !! Amount of water missing for irrigation (kg)
+    REAL(r_std), DIMENSION(nbpt_r)        :: irrig_adduct_r              !! Amount of water carried over from other basins for irrigation (kg)
+    !
+    REAL(r_std), DIMENSION(nbpt_r)        :: transport_r                 !! Water transport between basins (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: floods_r                    !! Water flow in to the floodplains (kg/dt)
+    REAL(r_std), DIMENSION(nbpt_r)        :: potflood_r                  !! Potential inflow to the swamps (kg/dt)
+    REAL(r_std), DIMENSION(nbpt)                 :: tobeflooded               !! Maximal surface which can be inundated in each grid box (m^2)
+    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) (not used)
+    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)                                  :: 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)
+    REAL(r_std)                                  :: adduction                 !! Importation of water from a stream reservoir of a neighboring grid box (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)                               :: ig, ib, ib2, ig2          !! Indices (unitless)
+    INTEGER(i_std)                               :: rtg, rtb, in              !! Indices (unitless)
+    INTEGER(i_std)                               :: ier                       !! Error handling 
+    INTEGER(i_std)                               :: isplit
+    
+    LOGICAL, PARAMETER                           :: check_reservoir = .TRUE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false)
+    REAL(r_std), DIMENSION(nbpt_r)        :: flood_frac_bas_r   
+    REAL(r_std), DIMENSION(nbpt_r)        :: flood_reservoir_r   
+    REAL(r_std), DIMENSION(nbpt_r)        :: pond_reservoir_r   
+
+    REAL(r_std), DIMENSION(nbpt_r)        :: lakeinflow_r   
+    REAL(r_std), DIMENSION(nbpt_r)        :: coastalflow_r   
+    REAL(r_std), DIMENSION(nbpt_r)        :: riverflow_r   
+
+    REAL(r_std), DIMENSION(nbpt_r)        :: flow_r   
+    REAL(r_std), DIMENSION(nbpt_rp1)      :: flow_rp1   
+    REAL(r_std) :: water_balance_before, water_balance_after
+    REAL(r_std) :: basins_riverflow_mpi(0:basins_count)
+    REAL(r_std) :: basins_riverflow(0:basins_count)
+
+    INTEGER :: ierr
+
+    !_ ================================================================================================================================
+    !
+
+
+    irrig_netereq(:) = zero
+    totarea(:) = zero
+    totflood(:) = zero
+
+    flood_frac_bas_r(:)=zero
+    !
+    !> 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
+    !> given by a 0.5 degree resolution map for each pixel performed from a simplification of Manning's formula
+    !> (Dingman, 1994; Ducharne et al., 2003).
+    !> The resulting product of tcst (in day/m) and topo_resid (in m) represents the time constant (day)
+    !> which is an e-folding time, the time necessary for the water amount
+    !> in the stream reservoir to decrease by a factor e. Hence, it gives an order of
+    !> magnitude of the travel time through this reservoir between
+    !> the sub-basin considered and its downstream neighbor.
+    CALL gather_omp(runoff_omp,runoff)
+    CALL gather_omp(drainage_omp, drainage)
+
+    IF (is_omp_root) THEN
+
+       hydrographs_r(:)=0
+
+       DO isplit=1,split_routing
+
+
+          irrig_needs_r(:) = zero
+          irrig_actual_r(:) = zero
+          irrig_deficit_r(:) = zero
+          irrig_adduct_r(:) = zero
+
+          DO ig=1,nbpt_r
+             IF ( routing_mask_r(ig) ) THEN
+                !
+                ! Each of the fluxes is limited by the water in the reservoir and a small margin 
+                ! (min_reservoir) to avoid rounding errors.
+                !
+                !               IF (fast_reservoir_r(ig)>1e-30) THEN
+                !                 PRINT*,"fast_reservoir > 0  ",ig, fast_reservoir_r(ig), topoind_r(ig)
+                !               ENDIF 
+                flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),&
+                     & fast_reservoir_r(ig)-min_sechiba)
+                fast_flow_r(ig) = MAX(flow, zero)
+                !
+                flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),&
+                     & slow_reservoir_r(ig)-min_sechiba)
+                slow_flow_r(ig) = MAX(flow, zero)
+                !
+                flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)*stream_tcst* & 
+                     & MAX(un-SQRT(flood_frac_bas_r(ig)),min_sechiba)*one_day/(dt_routing/split_routing)),&
+                     & stream_reservoir_r(ig)-min_sechiba)
+                stream_flow_r(ig) = MAX(flow, zero)
+                ! 
+                !
+             ELSE
+                fast_flow_r(ig) = zero
+                slow_flow_r(ig) = zero
+                stream_flow_r(ig) = zero
+             ENDIF
+          ENDDO
+
+          DO ig=1,nbpt_r
+             flood_drainage_r(ig) = zero
+             flood_flow_r(ig) = zero
+             flood_reservoir_r(ig) = zero
+          ENDDO
+
+          DO ig=1,nbpt_r
+             pond_inflow_r(ig) = zero
+             pond_drainage_r(ig) = zero
+             pond_reservoir_r(ig) = zero
+          ENDDO
+
+
+          !-
+          !- Compute the transport
+          !-
+
+          flow_r(:)=fast_flow_r(:) + slow_flow_r(:) + stream_flow_r(:)
+
+          CALL xios_send_field("flow_r",flow_r)       ! transfer halo
+          CALL xios_recv_field("flow_rp1",flow_rp1)
+
+          transport_r(:)=0
+
+          DO ig=1,nbpt_rp1
+             IF ( route_flow_rp1(ig) > 0 ) THEN
+                transport_r(route_flow_rp1(ig))=transport_r(route_flow_rp1(ig))+ flow_rp1(ig)
+             ENDIF
+          ENDDO
+
+
+          !-
+          !- Do the floodings - First initialize
+          !-
+          return_swamp_r(:)=zero
+          floods_r(:)=zero
+
+          !
+          ! Update all reservoirs
+          !> The slow and deep reservoir (slow_reservoir) collect the deep drainage whereas the
+          !> fast_reservoir collects the computed surface runoff. Both discharge into a third reservoir
+          !> (stream_reservoir) of the next sub-basin downstream.
+          !> Water from the floodplains reservoir (flood_reservoir) flows also into the stream_reservoir of the next sub-basin downstream.
+          !> Water that flows into the pond_reservoir is withdrawn from the fast_reservoir.
+          !
+          runoff=runoff*routing_weight
+          CALL xios_send_field("routing_runoff",runoff)   ! interp conservative model -> routing
+          CALL xios_recv_field("routing_runoff_r",runoff_r)
+          drainage=drainage*routing_weight
+          CALL xios_send_field("routing_drainage",drainage)   ! interp conservative model -> routing
+          CALL xios_recv_field("routing_drainage_r",drainage_r)
+
+          CALL xios_send_field("fast_reservoir_r",fast_reservoir_r)
+          CALL xios_send_field("slow_reservoir_r",slow_reservoir_r)
+          CALL xios_send_field("stream_reservoir_r",stream_reservoir_r)
+
+          WHERE (.NOT. routing_mask_r) runoff_r(:)=0
+          WHERE (.NOT. routing_mask_r) drainage_r(:)=0
+
+          CALL MPI_ALLREDUCE(sum(runoff+drainage),water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+          CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+
+          PRINT *,"routing water Balance ;  before : ", water_balance_before," ; after : ",water_balance_after,  &
+               " ; delta : ", 100.*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%"
+
+          CALL xios_recv_field("water_balance_before",water_balance_before)
+
+
+          DO ig=1,nbpt_r
+             IF ( routing_mask_r(ig) ) THEN
+
+                !
+                fast_reservoir_r(ig) =  fast_reservoir_r(ig) + runoff_r(ig) - &
+                     & fast_flow_r(ig) - pond_inflow_r(ig)
+                !
+                slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) - &
+                     & slow_flow_r(ig)
+                !
+
+                stream_reservoir_r(ig) = stream_reservoir_r(ig) + flood_flow_r(ig)  - &
+                     & stream_flow_r(ig) - return_swamp_r(ig) - floods_r(ig)
+
+                IF (is_streamflow_r(ig)) stream_reservoir_r(ig)= stream_reservoir_r(ig) +  transport_r(ig)  
+                !
+                flood_reservoir_r(ig) = flood_reservoir_r(ig) + floods_r(ig) - &
+                     & flood_flow_r(ig) 
+                !
+                pond_reservoir_r(ig) = pond_reservoir_r(ig) + pond_inflow_r(ig) - pond_drainage_r(ig)
+                !
+                IF ( flood_reservoir_r(ig) .LT. zero ) THEN
+                   IF ( check_reservoir ) THEN
+                      WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ". Problem is being corrected."
+                      WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir_r(ig), floods_r(ig), &
+                           & flood_flow_r(ig) 
+                   ENDIF
+                   stream_reservoir_r(ig) = stream_reservoir_r(ib) + flood_reservoir_r(ib)
+                   flood_reservoir_r(ib) = zero
+                ENDIF
+                !
+                IF ( stream_reservoir_r(ig) .LT. zero ) THEN
+                   IF ( check_reservoir ) THEN
+                      WRITE(numout,*) "WARNING : negative stream reservoir at :", ig, ". Problem is being corrected."
+                      WRITE(numout,*) "stream_reservoir, flood_flow, transport : ", stream_reservoir_r(ig), flood_flow_r(ig), &
+                           &  transport_r(ig)
+                      WRITE(numout,*) "stream_flow, return_swamp, floods :", stream_flow_r(ig), return_swamp_r(ig), floods_r(ig)
+                   ENDIF
+                   fast_reservoir_r(ig) =  fast_reservoir_r(ig) + stream_reservoir_r(ig)
+                   stream_reservoir_r(ig) = zero
+                ENDIF
+                !
+                IF ( fast_reservoir_r(ig) .LT. zero ) THEN
+                   IF ( check_reservoir ) THEN
+                      WRITE(numout,*) "WARNING : negative fast reservoir at :", ig, ". Problem is being corrected."
+                      WRITE(numout,*) "fast_reservoir, runoff, fast_flow, ponf_inflow  : ", fast_reservoir_r(ig), &
+                           &runoff_r(ig), fast_flow_r(ig), pond_inflow_r(ig)
+                   ENDIF
+                   slow_reservoir_r(ig) =  slow_reservoir_r(ig) + fast_reservoir_r(ig)
+                   fast_reservoir_r(ig) = zero
+                ENDIF
+
+                IF ( slow_reservoir_r(ig) .LT. - min_sechiba ) THEN
+                   WRITE(numout,*) 'WARNING : There is a negative reservoir at :', ig
+                   WRITE(numout,*) 'WARNING : slowr, slow_flow, drainage', &
+                        & slow_reservoir_r(ig), slow_flow_r(ig), drainage_r(ig)
+                   WRITE(numout,*) 'WARNING : pondr, pond_inflow, pond_drainage', &
+                        & pond_reservoir_r(ig), pond_inflow_r(ig), pond_drainage_r(ig)
+                   CALL ipslerr_p(2, 'routing_simple_flow', 'WARNING negative slow_reservoir.','','')
+                ENDIF
+
+             ENDIF
+          ENDDO
+
+          DO ig=1,nbpt_r
+             IF ( routing_mask_r(ig) ) THEN
+                hydrographs_r(ig)=hydrographs_r(ig)+fast_flow_r(ig)+slow_flow_r(ig)+stream_flow_r(ig)
+             ENDIF
+          ENDDO
+
+       ENDDO ! isplit
+
+       !
+       !
+       !
+       ! The three output types for the routing : endoheric basins,, rivers and 
+       ! diffuse coastal flow.
+       !
+       lakeinflow_r(:)=0
+       coastalflow_r(:)=0
+       riverflow_r(:)=0
+       basins_riverflow_mpi(:)=0
+
+       DO ig=1,nbpt_r
+          IF ( routing_mask_r(ig) ) THEN
+
+             IF (is_lakeinflow_r(ig))  THEN
+                lakeinflow_r(ig) = transport_r(ig)
+                basins_riverflow_mpi(basins_extended_r(ig)) = &
+                     basins_riverflow_mpi(basins_extended_r(ig))+lakeinflow_r(ig)
+             ENDIF
+
+             IF (is_coastalflow_r(ig)) THEN
+                coastalflow_r(ig) = transport_r(ig)
+                basins_riverflow_mpi(basins_extended_r(ig)) = &
+                     basins_riverflow_mpi(basins_extended_r(ig))+coastalflow_r(ig)
+             ENDIF
+
+             IF (is_riverflow_r(ig)) THEN
+                riverflow_r(ig) = transport_r(ig)
+                basins_riverflow_mpi(basins_extended_r(ig)) = &
+                     basins_riverflow_mpi(basins_extended_r(ig))+riverflow_r(ig)
+             ENDIF
+
+          ENDIF
+       ENDDO
+
+
+       CALL MPI_ALLREDUCE(basins_riverflow_mpi,basins_riverflow,basins_count+1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+       CALL xios_send_field("basins_riverflow",basins_riverflow(1:basins_out)/1000./dt_routing)
+
+       !
+       flood_res = flood_diag + pond_diag
+       !
+       CALL xios_send_field("routing_lakeinflow_r"  ,lakeinflow_r*routing_weight_r)
+       CALL xios_recv_field("routing_lakeinflow"    ,lakeinflow)
+       CALL xios_send_field("routing_coastalflow_r" ,coastalflow_r*routing_weight_coast_r)
+       CALL xios_recv_field("routing_coastalflow_temp"   ,coastalflow)
+       WHERE(.NOT. coast_mask(:)) coastalflow(:)=0
+       CALL xios_send_field("routing_coastalflow"   ,coastalflow)
+
+       CALL xios_send_field("routing_riverflow_r"   ,riverflow_r*routing_weight_coast_r)
+       CALL xios_recv_field("routing_riverflow_temp"     ,riverflow)
+       WHERE(.NOT. coast_mask(:)) riverflow(:)=0
+       CALL xios_send_field("routing_riverflow"     ,riverflow)
+
+       CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow)
+       CALL xios_send_field("routing_hydrographs_r", hydrographs_r/1000./dt_routing)
+
+       ! diag
+       CALL xios_send_field("routing_fast_reservoir_r"  , fast_reservoir_r)
+       CALL xios_recv_field("routing_fast_reservoir"  , fast_diag_mpi)
+
+       CALL xios_send_field("routing_slow_reservoir_r"  , slow_reservoir_r)
+       CALL xios_recv_field("routing_slow_reservoir"  ,   slow_diag_mpi)
+
+       CALL xios_send_field("routing_stream_reservoir_r"  , stream_reservoir_r)
+       CALL xios_recv_field("routing_stream_reservoir"    , stream_diag_mpi)
+
+       CALL xios_recv_field("water_balance_after"  ,water_balance_after)
+
+       PRINT *,"routing water Balance ;  before : ", water_balance_before," ; after : ",water_balance_after,  &
+            " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%"
+
+
+
+    ENDIF ! is_omp_root
+
+    CALL scatter_omp(riverflow,riverflow_omp)
+    CALL scatter_omp(coastalflow,coastalflow_omp)
+    CALL scatter_omp(lakeinflow,lakeinflow_omp)
+    CALL scatter_omp(fast_diag_mpi,fast_diag)
+    CALL scatter_omp(slow_diag_mpi,slow_diag)
+    CALL scatter_omp(stream_diag_mpi,stream_diag)
+
+    DO ig=1,nbpt
+       fast_diag(ig)=fast_diag(ig)/area(ig)  
+       slow_diag(ig)=slow_diag(ig)/area(ig)    
+       stream_diag(ig)=stream_diag(ig)/area(ig)  
+    ENDDO
+
+
+    flood_frac(:) = zero
+    flood_height(:) = zero
+    flood_frac_bas_r(:) = zero
+
+    returnflow(:) = zero
+    reinfiltration(:) = zero
+
+    !
+    !
+    ! Compute the fluxes which leave the routing scheme
+    !
+    ! Lakeinflow is in Kg/dt
+    ! returnflow is in Kg/m^2/dt
+    !
+    delsurfstor(:) = -flood_diag(:)-pond_diag(:)-lake_diag(:)
+    hydrographs(:) = zero
+    slowflow_diag(:) = zero
+    flood_diag(:) =  zero
+    pond_diag(:) =  zero
+    irrigation(:) = zero
+
+
+  END SUBROUTINE routing_simple_flow
+
+
+  !!  =============================================================================================================================
+  !! SUBROUTINE:         routing_simple_initialize
+  !!
+  !>\BRIEF	         Initialize the routing_simple module
+  !!
+  !! DESCRIPTION:        Initialize the routing_simple module. Read from restart file or read the routing.nc file to initialize the
+  !!                     routing scheme. 
+  !!
+  !! RECENT CHANGE(S)
+  !!
+  !! REFERENCE(S)
+  !! 
+  !! FLOWCHART   
+  !! \n
+  !_ ==============================================================================================================================
+
+  SUBROUTINE routing_simple_initialize( kjit,       nbpt,           index,                 &
+       rest_id,     hist_id,        hist2_id,   lalo,      &
+       neighbours,  resolution,     contfrac,   stempdiag, &
+       returnflow,  reinfiltration, irrigation, riverflow, &
+       coastalflow, flood_frac,     flood_res )
+
+    IMPLICIT NONE
+
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+    INTEGER(i_std), INTENT(in)     :: kjit                 !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)     :: nbpt                 !! Domain size (unitless)
+    INTEGER(i_std), INTENT(in)     :: index(nbpt)          !! Indices of the points on the map (unitless)
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    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)
+    REAL(r_std), INTENT(in)        :: lalo(nbpt,2)         !! Vector of latitude and longitudes (beware of the order !)
+
+    INTEGER(i_std), INTENT(in)     :: neighbours(nbpt,8)   !! Vector of neighbours for each grid point 
+                                                           !! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
+    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)
+    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+
+    !! 0.2 Output variables
+    REAL(r_std), INTENT(out)       :: returnflow(nbpt)     !! The water flow from lakes and swamps which returns to the grid box.
+                                                           !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
+    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)
+    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)
+
+    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)
+    REAL(r_std), INTENT(out)       :: flood_frac(nbpt)     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(out)       :: flood_res(nbpt)      !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+    !! 0.3 Local variables
+    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)
+
+    !_ ================================================================================================================================
+
+    !
+    ! do initialisation
+    !
+    nbvmax = 440
+    ! Here we will allocate the memory and get the fixed fields from the restart file.
+    ! If the info is not found then we will compute the routing map.
+    !
+    CALL routing_simple_init_1 (kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
+         riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id)
+
+    routing_area => routing_area_loc  
+    topo_resid => topo_resid_loc
+    route_togrid => route_togrid_loc
+    route_tobasin => route_tobasin_loc
+    global_basinid => global_basinid_loc
+    hydrodiag => hydrodiag_loc
+
+    ! This routine computes the routing map if the route_togrid_glo is undefined. This means that the
+    ! map has not been initialized during the restart process..
+    !
+    !! Reads in the map of the basins and flow directions to construct the catchments of each grid box
+    !
+    IF ( COUNT(route_togrid_glo .GE. undef_int) .GT. 0 ) THEN
+       !       CALL routing_basins_p(nbpt, lalo, neighbours, resolution, contfrac)
+    ENDIF
+    !
+    ! Do we have what we need if we want to do irrigation
+    !! 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.
+    END IF
+
+    init_flood = .FALSE.
+    IF ( do_floodplains ) THEN
+       IF (COUNT(floodplains .GE. undef_sechiba-1) > 0) init_flood = .TRUE.
+    END IF
+
+    init_swamp = .FALSE.
+    IF ( doswamps ) THEN
+       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)
+    ENDIF
+
+    IF ( do_irrigation ) THEN 
+       CALL xios_orchidee_send_field("irrigmap",irrigated)
+
+       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
+    ENDIF
+
+    IF ( do_floodplains ) THEN
+       CALL xios_orchidee_send_field("floodmap",floodplains)
+
+       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)
+
+       WRITE(numout,*) 'Verification : range of swamp : ', MINVAL(swamp), MAXVAL(swamp) 
+       IF ( .NOT. almaoutput ) THEN
+          CALL histwrite_p(hist_id, 'swampmap', 1, swamp, nbpt, index)
+       ELSE
+          CALL histwrite_p(hist_id, 'SwampMap', 1, swamp, nbpt, index)
+       ENDIF
+       IF ( hist2_id > 0 ) THEN
+          IF ( .NOT. almaoutput ) THEN
+             CALL histwrite_p(hist2_id, 'swampmap', 1, swamp, nbpt, index)
+          ELSE
+             CALL histwrite_p(hist2_id, 'SwampMap', 1, swamp, nbpt, index)
+          ENDIF
+       ENDIF
+    ENDIF
+
+    !! This routine gives a diagnostic of the basins used.
+    !    CALL routing_diagnostic_p(nbpt, index, lalo, resolution, contfrac, hist_id, hist2_id)
+
+    CALL routing_simple_init_2(nbpt, contfrac)
+
+  END SUBROUTINE routing_simple_initialize
+
+
+  !! ================================================================================================================================
+  !! SUBROUTINE   : routing_simple_main 
+  !!
+  !>\BRIEF          This module routes the water over the continents (runoff and
+  !!                drainage produced by the hydrolc or hydrol module) into the oceans. 
+  !!
+  !! DESCRIPTION (definitions, functional, design, flags):
+  !! The routing scheme (Polcher, 2003) carries the water from the runoff and drainage simulated by SECHIBA
+  !! to the ocean through reservoirs, with some delay. The routing scheme is based on
+  !! a parametrization of the water flow on a global scale (Miller et al., 1994; Hagemann
+  !! and Dumenil, 1998). Given the global map of the main watersheds (Oki et al., 1999;
+  !! Fekete et al., 1999; Vorosmarty et al., 2000) which delineates the boundaries of subbasins
+  !! and gives the eight possible directions of water flow within the pixel, the surface
+  !! runoff and the deep drainage are routed to the ocean. The time-step of the routing is one day.
+  !! The scheme also diagnoses how much water is retained in the foodplains and thus return to soil
+  !! moisture or is taken out of the rivers for irrigation. \n
+  !!
+  !! RECENT CHANGE(S): None
+  !!
+  !! MAIN OUTPUT VARIABLE(S): 
+  !! The result of the routing are 3 fluxes :
+  !! - riverflow   : The water which flows out from the major rivers. The flux will be located
+  !!                 on the continental grid but this should be a coastal point.
+  !! - coastalflow : This is the water which flows in a disperse way into the ocean. Essentially these
+  !!                 are the outflows from all of the small rivers.
+  !! - returnflow  : This is the water which flows into a land-point - typically rivers which end in
+  !!                 the desert. This water will go back into the hydrol module to allow re-evaporation.
+  !! - irrigation  : This is water taken from the reservoir and is being put into the upper 
+  !!                 layers of the soil.
+  !! The two first fluxes are in kg/dt and the last two fluxes are in kg/(m^2dt).\n
+  !!
+  !! REFERENCE(S) :
+  !! - Miller JR, Russell GL, Caliri G (1994)
+  !!   Continental-scale river flow in climate models.
+  !!   J. Clim., 7:914-928
+  !! - Hagemann S and Dumenil L. (1998)
+  !!   A parametrization of the lateral waterflow for the global scale.
+  !!   Clim. Dyn., 14:17-31
+  !! - Oki, T., T. Nishimura, and P. Dirmeyer (1999)
+  !!   Assessment of annual runoff from land surface models using total runoff integrating pathways (TRIP)
+  !!   J. Meteorol. Soc. Jpn., 77, 235-255
+  !! - Fekete BM, Charles V, Grabs W (2000)
+  !!   Global, composite runoff fields based on observed river discharge and simulated water balances.
+  !!   Technical report, UNH/GRDC, Global Runoff Data Centre, Koblenz
+  !! - Vorosmarty, C., B. Fekete, B. Meybeck, and R. Lammers (2000)
+  !!   Global system of rivers: Its role in organizing continental land mass and defining land-to-ocean linkages
+  !!   Global Biogeochem. Cycles, 14, 599-621
+  !! - Vivant, A-C. (?? 2002)
+  !!   DÃ©veloppement du schÃ©ma de routage et des plaines d'inondation, MSc Thesis, Paris VI University
+  !! - J. Polcher (2003)
+  !!   Les processus de surface a l'echelle globale et leurs interactions avec l'atmosphere
+  !!   Habilitation a diriger les recherches, Paris VI University, 67pp.
+  !!
+  !! FLOWCHART    :
+  !! \latexonly 
+  !! \includegraphics[scale=0.75]{routing_main_flowchart.png}
+  !! \endlatexonly
+  !! \n
+  !_ ================================================================================================================================
+
+
+  SUBROUTINE routing_simple_main(kjit, nbpt, index, &
+                                 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)
+
+    IMPLICIT NONE
+
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+
+    INTEGER(i_std), INTENT(in)     :: kjit                 !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)     :: nbpt                 !! Domain size (unitless)
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    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)
+    INTEGER(i_std), INTENT(in)     :: 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 !)
+    INTEGER(i_std), INTENT(in)     :: neighbours(nbpt,8)   !! Vector of neighbours for each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
+    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)
+    REAL(r_std), INTENT(in)        :: totfrac_nobio(nbpt)  !! Total fraction of no-vegetation (continental ice, lakes ...) (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: veget_max(nbpt,nvm)  !! Maximal fraction of vegetation (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: floodout(nbpt)       !! Grid-point flow out of floodplains (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: runoff(nbpt)         !! Grid-point runoff (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: drainage(nbpt)       !! Grid-point drainage (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: transpot(nbpt,nvm)   !! Potential transpiration of the vegetation (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: precip_rain(nbpt)    !! Rainfall (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: k_litt(nbpt)         !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: humrel(nbpt,nvm)     !! Soil moisture stress, root extraction potential (unitless)
+    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)
+
+    !! 0.2 Output variables
+    REAL(r_std), INTENT(out)       :: returnflow(nbpt)     !! The water flow from lakes and swamps which returns to the grid box.
+    !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
+    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)
+    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)
+    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)
+    REAL(r_std), INTENT(out)       :: flood_frac(nbpt)     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(out)       :: flood_res(nbpt)      !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+    !! 0.3 Local variables
+    CHARACTER(LEN=30)              :: var_name             !! To store variables names for I/O (unitless)
+    REAL(r_std), DIMENSION(1)      :: tmp_day              !! 
+    REAL(r_std), DIMENSION(nbpt)   :: return_lakes         !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
+    INTEGER(i_std)                 :: ig, jv               !! Indices (unitless)
+    REAL(r_std), DIMENSION(nbpt)   :: tot_vegfrac_nowoody  !! Total fraction occupied by grass (0-1,unitless)
+    REAL(r_std),PARAMETER :: delta_lon=(360./144.)/2.
+    REAL(r_std),PARAMETER :: delta_lat=(180./144.)/2.
+
+    !_ ================================================================================================================================
+
+    !
+    !! Computes the variables averaged between routing time steps and which will be used in subsequent calculations
+    !
+    floodout_mean(:) = floodout_mean(:) + floodout(:)
+    runoff_mean(:) = runoff_mean(:) + runoff(:)
+    drainage_mean(:) = drainage_mean(:) + drainage(:)
+    floodtemp(:) = stempdiag(:,floodtemp_lev)
+    precip_mean(:) =  precip_mean(:) + precip_rain(:)
+    !
+    !! Computes the total fraction occupied by the grasses and the crops for each grid cell
+    tot_vegfrac_nowoody(:) = zero
+    DO jv  = 1, nvm
+       IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
+          tot_vegfrac_nowoody(:) = tot_vegfrac_nowoody(:) + veget_max(:,jv) 
+       END IF
+    END DO
+
+    DO ig = 1, nbpt
+       IF ( tot_vegfrac_nowoody(ig) .GT. min_sechiba ) THEN
+          DO jv = 1,nvm
+             IF ( (jv /= ibare_sechiba) .AND. .NOT.(is_tree(jv)) ) THEN
+                transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/tot_vegfrac_nowoody(ig)  
+             END IF
+          END DO
+       ELSE
+          IF (MAXVAL(veget_max(ig,2:nvm)) .GT. min_sechiba) THEN
+             DO jv = 2, nvm
+                transpot_mean(ig) = transpot_mean(ig) + transpot(ig,jv) * veget_max(ig,jv)/ SUM(veget_max(ig,2:nvm))
+             ENDDO
+          ENDIF
+       ENDIF
+    ENDDO
+
+    !
+    ! Averaged variables (i.e. *dt_sechiba/dt_routing). This accounts for the difference between the shorter
+    ! timestep dt_sechiba of other parts of the model and the long dt_routing timestep (set to one day at present)
+    !
+    totnobio_mean(:) = totnobio_mean(:) + totfrac_nobio(:)*dt_sechiba/dt_routing
+    k_litt_mean(:) = k_litt_mean(:) + k_litt(:)*dt_sechiba/dt_routing
+    !
+    ! Only potentially vegetated surfaces are taken into account. At the start of
+    ! 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
+    !
+    time_counter = time_counter + dt_sechiba 
+    !
+    ! If the time has come we do the routing.
+    !
+    IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN 
+       !
+       ! Check the water balance if needed
+       !
+       !       IF ( check_waterbal ) THEN
+       !          CALL routing_waterbal(nbpt, .TRUE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, &
+       !               & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean)
+       !       ENDIF
+       !
+       ! Make sure we do not flood north of 49N as there freezing processes start to play a role and they
+       ! are not yet well treated in ORCHIDEE.
+       !
+       DO ig=1,nbpt
+          IF ( lalo(ig,1) > 49.0 ) THEN
+             floodtemp(ig) = tp_00 - un
+          ENDIF
+       ENDDO
+       !
+       !! Computes the transport of water in the various reservoirs
+       !
+       !ym       runoff_mean=0
+       !ym       drainage_mean=0
+       !ym      DO ig=1,nbpt
+       !ym         IF ( lalo(ig,1)-delta_lat < 0 .AND. lalo(ig,1)+delta_lat > 0 .AND. lalo(ig,2)-delta_lon < 32. .AND. lalo(ig,2)+delta_lon > 32.) runoff_mean(ig)=one_day/dt_routing        
+       !ym       ENDDO
+
+       CALL routing_simple_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, &
+            & 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, &
+            & coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res)
+       !
+       !! Responsible for storing the water in lakes
+       !
+       CALL routing_simple_lake(nbpt, dt_routing, lakeinflow_mean, humrel_mean, contfrac, return_lakes)
+       !
+       returnflow_mean(:) = returnflow_mean(:) + return_lakes(:)
+       !
+       !! Check the water balance in the routing scheme
+       !
+       !       IF ( check_waterbal ) THEN
+       !          CALL routing_waterbal(nbpt, .FALSE., floodout_mean, runoff_mean, drainage_mean, returnflow_mean, &
+       !               & reinfiltration_mean, irrigation_mean, riverflow_mean, coastalflow_mean)
+       !       ENDIF
+       !
+       time_counter = zero
+       !
+       floodout_mean(:) = zero
+       runoff_mean(:) = zero
+       drainage_mean(:) = zero
+       transpot_mean(:) = zero
+       precip_mean(:) = zero
+       !
+       humrel_mean(:) = zero
+       totnobio_mean(:) = zero
+       k_litt_mean(:) = zero
+       vegtot_mean(:) = zero
+       !
+       ! Change the units of the routing fluxes from kg/dt_routing into kg/dt_sechiba
+       !                                    and from m^3/dt_routing into m^3/dt_sechiba
+       !
+       returnflow_mean(:) = returnflow_mean(:)/dt_routing*dt_sechiba
+       reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba
+       irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba
+       irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba
+       !
+       ! Change units as above but at the same time transform the kg/dt_sechiba to m^3/dt_sechiba
+       !
+       riverflow_mean(:) = riverflow_mean(:)/dt_routing*dt_sechiba/mille
+       coastalflow_mean(:) = coastalflow_mean(:)/dt_routing*dt_sechiba/mille
+       hydrographs(:) = hydrographs(:)/dt_routing*dt_sechiba/mille
+       slowflow_diag(:) = slowflow_diag(:)/dt_routing*dt_sechiba/mille
+    ENDIF
+
+    !
+    ! Return the fraction of routed water for this time step.
+    !
+    returnflow(:) = returnflow_mean(:)
+    reinfiltration(:) = reinfiltration_mean(:)
+    irrigation(:) = irrigation_mean(:)
+    riverflow(:) = riverflow_mean(:)
+    coastalflow(:) = coastalflow_mean(:)
+
+    !
+    ! Write diagnostics
+    !
+
+    CALL xios_orchidee_send_field("riversret",returnflow*one_day/dt_sechiba)
+    CALL xios_orchidee_send_field("hydrographs",hydrographs/dt_sechiba)
+    CALL xios_orchidee_send_field("slowflow",slowflow_diag/dt_sechiba) ! Qb in m3/s
+    CALL xios_orchidee_send_field("reinfiltration",reinfiltration)
+    CALL xios_orchidee_send_field("fastr",fast_diag)
+    CALL xios_orchidee_send_field("slowr",slow_diag)
+    CALL xios_orchidee_send_field("streamr",stream_diag)
+    CALL xios_orchidee_send_field("laker",lake_diag)
+    CALL xios_orchidee_send_field("pondr",pond_diag)
+    CALL xios_orchidee_send_field("irrigation",irrigation*one_day/dt_sechiba)
+    CALL xios_orchidee_send_field("netirrig",irrig_netereq*one_day/dt_sechiba)
+    CALL xios_orchidee_send_field("floodh",flood_height)
+    CALL xios_orchidee_send_field("floodr",flood_diag)
+    delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag
+    CALL xios_orchidee_send_field("SurfStor",flood_diag+pond_diag+lake_diag)
+
+    ! Transform from kg/dt_sechiba into m^3/s
+    CALL xios_orchidee_send_field("hydrographs",hydrographs/mille/dt_sechiba)
+    CALL xios_orchidee_send_field("slowflow",slowflow_diag/mille/dt_sechiba) ! previous id name: Qb
+    CALL xios_orchidee_send_field("coastalflow",coastalflow/dt_sechiba)
+    CALL xios_orchidee_send_field("riverflow",riverflow/dt_sechiba)
+
+
+    IF ( .NOT. almaoutput ) THEN
+       !
+       CALL histwrite_p(hist_id, 'riversret', kjit, returnflow, nbpt, index)
+       IF (do_floodplains .OR. doponds) THEN
+          CALL histwrite_p(hist_id, 'reinfiltration', kjit, reinfiltration, nbpt, index)
+       ENDIF
+       CALL histwrite_p(hist_id, 'hydrographs', kjit, hydrographs, nbpt, index)
+       !
+       CALL histwrite_p(hist_id, 'fastr', kjit, fast_diag, nbpt, index)
+       CALL histwrite_p(hist_id, 'slowr', kjit, slow_diag, nbpt, index)
+       CALL histwrite_p(hist_id, 'streamr', kjit, stream_diag, nbpt, index)
+       IF ( do_floodplains ) THEN
+          CALL histwrite_p(hist_id, 'floodr', kjit, flood_diag, nbpt, index)
+          CALL histwrite_p(hist_id, 'floodh', kjit, flood_height, nbpt, index)
+       ENDIF
+       CALL histwrite_p(hist_id, 'pondr', kjit, pond_diag, nbpt, index)
+       CALL histwrite_p(hist_id, 'lakevol', kjit, lake_diag, nbpt, index)
+       !
+       IF ( do_irrigation ) THEN
+          CALL histwrite_p(hist_id, 'irrigation', kjit, irrigation, nbpt, index)
+          CALL histwrite_p(hist_id, 'returnflow', kjit, returnflow, nbpt, index)
+          CALL histwrite_p(hist_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
+       ENDIF
+       !
+    ELSE
+       !
+       delsurfstor = delsurfstor + flood_diag + pond_diag + lake_diag
+       CALL histwrite_p(hist_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index)
+       CALL histwrite_p(hist_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index)
+       CALL histwrite_p(hist_id, 'Dis', kjit, hydrographs, nbpt, index)
+       IF ( do_irrigation ) THEN
+          CALL histwrite_p(hist_id, 'Qirrig', kjit, irrigation, nbpt, index)
+          CALL histwrite_p(hist_id, 'Qirrig_req', kjit, irrig_netereq, nbpt, index)
+       ENDIF
+       !
+    ENDIF
+    IF ( hist2_id > 0 ) THEN
+       IF ( .NOT. almaoutput ) THEN
+          !
+          CALL histwrite_p(hist2_id, 'riversret', kjit, returnflow, nbpt, index)
+          IF (do_floodplains .OR. doponds) THEN
+             CALL histwrite_p(hist2_id, 'reinfiltration', kjit, reinfiltration, nbpt, index)
+          ENDIF
+          CALL histwrite_p(hist2_id, 'hydrographs', kjit, hydrographs, nbpt, index)
+          !
+          CALL histwrite_p(hist2_id, 'fastr', kjit, fast_diag, nbpt, index)
+          CALL histwrite_p(hist2_id, 'slowr', kjit, slow_diag, nbpt, index)
+          IF ( do_floodplains ) THEN
+             CALL histwrite_p(hist2_id, 'floodr', kjit, flood_diag, nbpt, index)
+             CALL histwrite_p(hist2_id, 'floodh', kjit, flood_height, nbpt, index)
+          ENDIF
+          CALL histwrite_p(hist2_id, 'pondr', kjit, pond_diag, nbpt, index)
+          CALL histwrite_p(hist2_id, 'streamr', kjit, stream_diag, nbpt, index)
+          CALL histwrite_p(hist2_id, 'lakevol', kjit, lake_diag, nbpt, index)
+          !
+          IF ( do_irrigation ) THEN
+             CALL histwrite_p(hist2_id, 'irrigation', kjit, irrigation, nbpt, index)
+             CALL histwrite_p(hist2_id, 'returnflow', kjit, returnflow, nbpt, index)
+             CALL histwrite_p(hist2_id, 'netirrig', kjit, irrig_netereq, nbpt, index)
+          ENDIF
+          !
+       ELSE
+          !
+          delsurfstor=delsurfstor + flood_diag + pond_diag + lake_diag
+          CALL histwrite_p(hist2_id, 'DelSurfStor', kjit, delsurfstor, nbpt, index)
+          CALL histwrite_p(hist2_id, 'SurfStor', kjit, flood_diag+pond_diag+lake_diag, nbpt, index)
+          CALL histwrite_p(hist2_id, 'Dis', kjit, hydrographs, nbpt, index)
+          !
+       ENDIF
+    ENDIF
+    !
+    !
+  END SUBROUTINE routing_simple_main
+
+  !!  =============================================================================================================================
+  !! SUBROUTINE:         routing_simple_finalize
+  !!
+  !>\BRIEF	         Write to restart file
+  !!
+  !! DESCRIPTION:        Write module variables to restart file
+  !!
+  !! RECENT CHANGE(S)
+  !!
+  !! REFERENCE(S)
+  !! 
+  !! FLOWCHART   
+  !! \n
+  !_ ==============================================================================================================================
+
+  SUBROUTINE routing_simple_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
+    USE xios
+    IMPLICIT NONE
+
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+    INTEGER(i_std), INTENT(in)     :: kjit                 !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)     :: nbpt                 !! Domain size (unitless)
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    REAL(r_std), INTENT(in)        :: flood_frac(nbpt)     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: flood_res(nbpt)      !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+    !! 0.2 Local variables
+    REAL(r_std), DIMENSION(1)      :: tmp_day              
+
+    !_ ================================================================================================================================
+
+    IF (is_omp_root) THEN
+       CALL xios_send_field("fast_reservoir_restart",fast_reservoir_r)
+       CALL xios_send_field("slow_reservoir_restart",slow_reservoir_r)
+       CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r)
+    ENDIF
+
+  END SUBROUTINE routing_simple_finalize
+
+  !! ================================================================================================================================
+  !! SUBROUTINE 	: routing_simple_init_1
+  !!
+  !>\BRIEF         This subroutine allocates the memory and get the fixed fields from the restart file.
+  !!
+  !! DESCRIPTION:	  Privat subroutine to the module routing_simple. This subroutine is called in the begining 
+  !!                      of routing_simple_initialize
+  !!
+  !! RECENT CHANGE(S): None
+  !!
+  !! MAIN OUTPUT VARIABLE(S):
+  !!
+  !! REFERENCES   : None
+  !!
+  !! FLOWCHART    :None
+  !! \n
+  !_ ================================================================================================================================
+
+  SUBROUTINE routing_simple_init_1(kjit, nbpt, index, returnflow, reinfiltration, irrigation, &
+                                   riverflow, coastalflow, flood_frac, flood_res, stempdiag, rest_id)
+    
+    IMPLICIT NONE
+        
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+
+    INTEGER(i_std), INTENT(in)                   :: kjit           !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)                   :: nbpt           !! Domain size (unitless)
+    INTEGER(i_std), DIMENSION (nbpt), INTENT(in) :: index          !! Indices of the points on the map (unitless)
+    REAL(r_std), DIMENSION(nbpt,nslm),INTENT(in) :: stempdiag      !! Temperature profile in soil
+    INTEGER(i_std), INTENT(in)                   :: rest_id        !! Restart file identifier (unitless)
+    
+    !! 0.2 Output variables
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: returnflow     !! The water flow from lakes and swamps which returns into the grid box.
+                                                                   !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: reinfiltration !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt)
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: irrigation     !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil.(kg/m^2/dt)
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: riverflow      !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt)
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: coastalflow    !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt)
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: flood_frac     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), DIMENSION (nbpt),INTENT(out)    :: flood_res      !! Diagnostic of water amount in the floodplains reservoir (kg)
+    
+    !! 0.3 Local variables
+    CHARACTER(LEN=80)                            :: var_name       !! To store variables names for I/O (unitless)
+    REAL(r_std), ALLOCATABLE, DIMENSION(:,:)     :: tmp_real_g     !! A temporary real array for the integers
+    REAL(r_std), DIMENSION(1)                    :: tmp_day        !!
+    REAL(r_std)                                  :: ratio          !! Diagnostic ratio to check that dt_routing is a multiple of dt_sechiba (unitless)
+    REAL(r_std)                                  :: totarea        !! Total area of basin (m^2)
+    INTEGER(i_std)                               :: ier, ig, ib, ipn(1) !! Indices (unitless)
+
+    !_ ================================================================================================================================
+    !
+    !
+    ! These variables will require the configuration infrastructure
+    !
+    !Config Key   = DT_ROUTING 
+    !Config If    = RIVER_ROUTING
+    !Config Desc  = Time step of the routing scheme
+    !Config Def   = one_day
+    !Config Help  = This values gives the time step in seconds of the routing scheme. 
+    !Config         It should be multiple of the main time step of ORCHIDEE. One day
+    !Config         is a good value.
+    !Config Units = [seconds]
+    !
+    dt_routing = one_day
+    CALL getin_p('DT_ROUTING', dt_routing)
+    !
+    !Config Key   = ROUTING_RIVERS
+    !Config If    = RIVER_ROUTING
+    !Config Desc  = Number of rivers 
+    !Config Def   = 50
+    !Config Help  = This parameter chooses the number of largest river basins
+    !Config         which should be treated as independently as rivers and not
+    !Config         flow into the oceans as diffusion coastal flow.
+    !Config Units = [-]
+    num_largest = 50
+    CALL getin_p('ROUTING_RIVERS', num_largest)
+    !
+    !Config Key   = DO_FLOODINFILT
+    !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         back to the slow and fast reservoirs
+    !Config Units = [FLAG]
+    !
+    dofloodinfilt = .FALSE.
+    CALL getin_p('DO_FLOODINFILT', dofloodinfilt)
+    !
+    !Config Key   = DO_SWAMPS
+    !Config Desc  = Should we include swamp parameterization 
+    !Config If    = RIVER_ROUTING
+    !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         internal deltas of rivers.
+    !Config Units = [FLAG]
+    !
+    doswamps = .FALSE.
+    CALL getin_p('DO_SWAMPS', doswamps)
+    !
+    !Config Key   = DO_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         little ponds especially in West Africa.
+    !Config Units = [FLAG]
+    !
+    doponds = .FALSE.
+    CALL getin_p('DO_PONDS', doponds)
+
+    !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         time constant (in days) of the slow reservoir
+    !Config         in order to get better river flows for 
+    !Config         particular regions.
+    !Config Units = [days]
+    !
+    !> A value for property of each reservoir (in day/m) is given to compute a time constant (in day)
+    !> for each reservoir (product of tcst and topo_resid).
+    !> The value of tcst has been calibrated for the three reservoirs over the Senegal river basin only,
+    !> 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. 
+    !> 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.
+    CALL getin_p('SLOW_TCST', slow_tcst)
+    !
+    !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         time constant (in days) of the fast reservoir
+    !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 If    = RIVER_ROUTING
+    !Config Def   = 0.24
+    !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         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 If    = RIVER_ROUTING
+    !Config Def   = 4.0
+    !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         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 If    = RIVER_ROUTING
+    !Config Def   = 0.2
+    !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 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 = [-] 
+    CALL getin_p("FLOOD_BETA", beta)
+
+    !Config Key   = POND_BETAP
+    !Config Desc  = Ratio of the basin surface intercepted by ponds and the maximum surface of ponds
+    !Config If    = RIVER_ROUTING
+    !Config Def   = 0.5
+    !Config Help  = 
+    !Config Units = [-] 
+    CALL getin_p("POND_BETAP", betap)    
+
+    !Config Key   = FLOOD_CRI
+    !Config Desc  = Potential height for which all the basin is flooded
+    !Config If    = DO_FLOODPLAINS or DO_PONDS
+    !Config Def   = 2000.
+    !Config Help  = 
+    !Config Units = [mm] 
+    CALL getin_p("FLOOD_CRI", floodcri)
+
+    !Config Key   = POND_CRI
+    !Config Desc  = Potential height for which all the basin is a pond
+    !Config If    = DO_FLOODPLAINS or DO_PONDS
+    !Config Def   = 2000.
+    !Config Help  = 
+    !Config Units = [mm] 
+    CALL getin_p("POND_CRI", pondcri)
+
+
+    ! In order to simplify the time cascade check that dt_routing
+    ! is a multiple of dt_sechiba
+    !
+    ratio = dt_routing/dt_sechiba
+    IF ( ABS(NINT(ratio) - ratio) .GT. 10*EPSILON(ratio)) THEN
+       WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
+       WRITE(numout,*) "The chosen time step for the routing is not a multiple of the"
+       WRITE(numout,*) "main time step of the model. We will change dt_routing so that"
+       WRITE(numout,*) "this condition os fulfilled"
+       dt_routing = NINT(ratio) * dt_sechiba
+       WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
+    ENDIF
+    !
+    IF ( dt_routing .LT. dt_sechiba) THEN
+       WRITE(numout,*) 'WARNING -- WARNING -- WARNING -- WARNING'
+       WRITE(numout,*) 'The routing timestep can not be smaller than the one'
+       WRITE(numout,*) 'of the model. We reset its value to the model''s timestep.'
+       WRITE(numout,*) 'The old DT_ROUTING is : ', dt_routing
+       dt_routing = dt_sechiba
+       WRITE(numout,*) 'THE NEW DT_ROUTING IS : ', dt_routing
+    ENDIF
+    !
+    var_name ="routingcounter"
+    IF (is_root_prc) THEN
+       CALL ioconf_setatt('UNITS', 's')
+       CALL ioconf_setatt('LONG_NAME','Time counter for the routing scheme')
+       CALL restget (rest_id, var_name, 1, 1, 1, kjit, .TRUE., tmp_day)
+       IF (tmp_day(1) == val_exp) THEN
+          time_counter = zero
+       ELSE
+          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_simple_init_1','Pb in allocate for routing_area_loc','','')
+
+    ALLOCATE (routing_area_glo(nbp_glo,nbasmax))
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for routing_area_glo','','')
+    var_name = 'routingarea'
+    IF (is_root_prc) THEN
+       CALL ioconf_setatt('UNITS', 'm^2')
+       CALL ioconf_setatt('LONG_NAME','Area of basin')
+       CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., routing_area_glo, "gather", nbp_glo, index_g)
+    ENDIF
+    CALL scatter(routing_area_glo,routing_area_loc)
+    routing_area=>routing_area_loc
+
+    ALLOCATE (tmp_real_g(nbp_glo,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for tmp_real_g','','')
+
+    ALLOCATE (route_togrid_loc(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_loc','','')
+    ALLOCATE (route_togrid_glo(nbp_glo,nbasmax), stat=ier)      ! used in global in routing_simple_flow
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_togrid_glo','','')
+
+    IF (is_root_prc) THEN
+       var_name = 'routetogrid'
+       CALL ioconf_setatt('UNITS', '-')
+       CALL ioconf_setatt('LONG_NAME','Grid into which the basin flows')
+       CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
+       route_togrid_glo(:,:) = undef_int
+       WHERE ( tmp_real_g .LT. val_exp )
+          route_togrid_glo = NINT(tmp_real_g)
+       ENDWHERE
+    ENDIF
+    CALL bcast(route_togrid_glo)                      ! used in global in routing_simple_flow
+    CALL scatter(route_togrid_glo,route_togrid_loc)
+    route_togrid=>route_togrid_loc
+    !
+    ALLOCATE (route_tobasin_loc(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_loc','','')
+
+    ALLOCATE (route_tobasin_glo(nbp_glo,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_tobasin_glo','','')
+
+    IF (is_root_prc) THEN
+       var_name = 'routetobasin'
+       CALL ioconf_setatt('UNITS', '-')
+       CALL ioconf_setatt('LONG_NAME','Basin in to which the water goes')
+       CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
+       route_tobasin_glo = undef_int
+       WHERE ( tmp_real_g .LT. val_exp )
+          route_tobasin_glo = NINT(tmp_real_g)
+       ENDWHERE
+    ENDIF
+    CALL scatter(route_tobasin_glo,route_tobasin_loc)
+    route_tobasin=>route_tobasin_loc
+    !
+    ! nbintobasin
+    !
+    ALLOCATE (route_nbintobas_loc(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_loc','','')
+    ALLOCATE (route_nbintobas_glo(nbp_glo,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for route_nbintobas_glo','','')
+
+    IF (is_root_prc) THEN
+       var_name = 'routenbintobas'
+       CALL ioconf_setatt('UNITS', '-')
+       CALL ioconf_setatt('LONG_NAME','Number of basin into current one')
+       CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
+       route_nbintobas_glo = undef_int
+       WHERE ( tmp_real_g .LT. val_exp )
+          route_nbintobas_glo = NINT(tmp_real_g)
+       ENDWHERE
+    ENDIF
+    CALL scatter(route_nbintobas_glo,route_nbintobas_loc)
+    route_nbintobas=>route_nbintobas_loc
+    !
+    ALLOCATE (global_basinid_loc(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_loc','','')
+    ALLOCATE (global_basinid_glo(nbp_glo,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for global_basinid_glo','','')
+
+    IF (is_root_prc) THEN
+       var_name = 'basinid'
+       CALL ioconf_setatt('UNITS', '-')
+       CALL ioconf_setatt('LONG_NAME','ID of basin')
+       CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., tmp_real_g, "gather", nbp_glo, index_g)
+       global_basinid_glo = undef_int
+       WHERE ( tmp_real_g .LT. val_exp )
+          global_basinid_glo = NINT(tmp_real_g)
+       ENDWHERE
+    ENDIF
+    CALL scatter(global_basinid_glo,global_basinid_loc)
+    global_basinid=>global_basinid_loc
+    !
+    ALLOCATE (topo_resid_loc(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_loc','','')
+    ALLOCATE (topo_resid_glo(nbp_glo,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for topo_resid_glo','','')
+
+    IF (is_root_prc) THEN
+       var_name = 'topoindex'
+       CALL ioconf_setatt('UNITS', 'm')
+       CALL ioconf_setatt('LONG_NAME','Topographic index of the residence time')
+       CALL restget (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., topo_resid_glo, "gather", nbp_glo, index_g)
+    ENDIF
+    CALL scatter(topo_resid_glo,topo_resid_loc)
+    topo_resid=>topo_resid_loc
+
+    ALLOCATE (fast_reservoir(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_reservoir','','')
+    var_name = 'fastres'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Water in the fast reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., fast_reservoir, "gather", nbp_glo, index_g)
+    CALL setvar_p (fast_reservoir, val_exp, 'NO_KEYWORD', zero)
+
+    ALLOCATE (slow_reservoir(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for slow_reservoir','','')
+    var_name = 'slowres'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Water in the slow reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., slow_reservoir, "gather", nbp_glo, index_g)
+    CALL setvar_p (slow_reservoir, val_exp, 'NO_KEYWORD', zero)
+
+    ALLOCATE (stream_reservoir(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for stream_reservoir','','')
+    var_name = 'streamres'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Water in the stream reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., stream_reservoir, "gather", nbp_glo, index_g)
+    CALL setvar_p (stream_reservoir, val_exp, 'NO_KEYWORD', zero)
+
+    ALLOCATE (flood_reservoir(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_reservoir','','')
+    var_name = 'floodres'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Water in the flood reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_reservoir, "gather", nbp_glo, index_g)
+    CALL setvar_p (flood_reservoir, val_exp, 'NO_KEYWORD', zero)
+
+    ALLOCATE (flood_frac_bas(nbpt,nbasmax), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_frac_bas','','')
+    var_name = 'flood_frac_bas'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per basin')
+    CALL restget_p (rest_id, var_name, nbp_glo, nbasmax, 1, kjit, .TRUE., flood_frac_bas, "gather", nbp_glo, index_g)
+    CALL setvar_p (flood_frac_bas, val_exp, 'NO_KEYWORD', zero)
+
+    ALLOCATE (flood_height(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for flood_height','','')
+    var_name = 'floodh'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','')
+    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_simple_init_1','Pb in allocate for pond_frac','','')
+    var_name = 'pond_frac'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Pond fraction per grid box')
+    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', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Flooded fraction per grid box')
+    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')
+    CALL ioconf_setatt_p('LONG_NAME','Flooded quantity (estimation)')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_res, "gather", nbp_glo, index_g)
+    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_simple_init_1','Pb in allocate for lake_reservoir','','')
+    var_name = 'lakeres'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Water in the lake reservoir')
+    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_simple_init_1','Pb in allocate for pond_reservoir','','')
+    var_name = 'pondres'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Water in the pond reservoir')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_reservoir, "gather", nbp_glo, index_g)
+    CALL setvar_p (pond_reservoir, val_exp, 'NO_KEYWORD', zero)
+    !
+    ! Map of irrigated areas
+    !
+    IF ( do_irrigation ) THEN
+       ALLOCATE (irrigated(nbpt), stat=ier)
+       IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','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
+
+    IF ( do_floodplains ) THEN
+       ALLOCATE (floodplains(nbpt), stat=ier)
+       IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodplains','','')
+       var_name = 'floodplains'
+       CALL ioconf_setatt_p('UNITS', 'm^2')
+       CALL ioconf_setatt_p('LONG_NAME','Surface which can be flooded')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodplains, "gather", nbp_glo, index_g)
+       CALL setvar_p (floodplains, val_exp, 'NO_KEYWORD', undef_sechiba)
+    ENDIF
+    IF ( doswamps ) THEN
+       ALLOCATE (swamp(nbpt), stat=ier)
+       IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for swamp','','')
+       var_name = 'swamp'
+       CALL ioconf_setatt_p('UNITS', 'm^2')
+       CALL ioconf_setatt_p('LONG_NAME','Surface which can become swamp')
+       CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., swamp, "gather", nbp_glo, index_g)
+       CALL setvar_p (swamp, val_exp, 'NO_KEYWORD', undef_sechiba)
+    ENDIF
+    !
+    ! Put into the restart file the fluxes so that they can be regenerated at restart.
+    !
+    ALLOCATE (lakeinflow_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for lakeinflow_mean','','')
+    var_name = 'lakeinflow'
+    CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+    CALL ioconf_setatt_p('LONG_NAME','Lake inflow')
+    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_simple_init_1','Pb in allocate for returnflow_mean','','')
+    var_name = 'returnflow'
+    CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
+    CALL ioconf_setatt_p('LONG_NAME','Deep return flux')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., returnflow_mean, "gather", nbp_glo, index_g)
+    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_simple_init_1','Pb in allocate for reinfiltration_mean','','')
+    var_name = 'reinfiltration'
+    CALL ioconf_setatt_p('UNITS', 'Kg/m^2/dt')
+    CALL ioconf_setatt_p('LONG_NAME','Top return flux')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., reinfiltration_mean, "gather", nbp_glo, index_g)
+    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_simple_init_1','Pb in allocate for irrigation_mean','','')
+    ALLOCATE (irrig_netereq(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for irrig_netereq','','')
+    irrig_netereq(:) = zero
+
+    IF ( do_irrigation ) THEN
+       var_name = 'irrigation'
+       CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+       CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux')
+       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)
+    ELSE
+       irrigation_mean(:) = zero
+    ENDIF
+    irrigation(:) = irrigation_mean(:) 
+
+    ALLOCATE (riverflow_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for riverflow_mean','','')
+    var_name = 'riverflow'
+    CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+    CALL ioconf_setatt_p('LONG_NAME','River flux into the sea')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., riverflow_mean, "gather", nbp_glo, index_g)
+    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_simple_init_1','Pb in allocate for coastalflow_mean','','')
+    var_name = 'coastalflow'
+    CALL ioconf_setatt_p('UNITS', 'Kg/dt')
+    CALL ioconf_setatt_p('LONG_NAME','Diffuse flux into the sea')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., coastalflow_mean, "gather", nbp_glo, index_g)
+    CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero)
+    coastalflow(:) = coastalflow_mean(:)
+
+    ! Locate it at the 2m level
+    ipn = MINLOC(ABS(diaglev-2))
+    floodtemp_lev = ipn(1)
+    ALLOCATE (floodtemp(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodtemp','','')
+    floodtemp(:) = stempdiag(:,floodtemp_lev)
+
+    ALLOCATE(hydrographs(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrographs','','')
+    var_name = 'hydrographs'
+    CALL ioconf_setatt_p('UNITS', 'm^3/dt')
+    CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid')
+    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_simple_init_1','Pb in allocate for slowflow_diag','','')
+    var_name = 'slowflow_diag'
+    CALL ioconf_setatt_p('UNITS', 'm^3/dt')
+    CALL ioconf_setatt_p('LONG_NAME','Slowflow hydrograph at outlow of grid')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE.,slowflow_diag, "gather", nbp_glo, index_g)
+    CALL setvar_p (slowflow_diag, val_exp, 'NO_KEYWORD', zero)
+
+    !
+    ! The diagnostic variables, they are initialized from the above restart variables.
+    !
+    ALLOCATE(fast_diag(nbpt), slow_diag(nbpt), stream_diag(nbpt), flood_diag(nbpt), &
+         & pond_diag(nbpt), lake_diag(nbpt), delsurfstor(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for fast_diag,..','','')
+
+    fast_diag(:) = zero
+    slow_diag(:) = zero
+    stream_diag(:) = zero
+    flood_diag(:) = zero
+    pond_diag(:) = zero
+    lake_diag(:) = zero
+    delsurfstor(:) = zero
+
+    DO ig=1,nbpt
+       totarea = zero
+       DO ib=1,nbasmax
+          totarea = totarea + routing_area(ig,ib)
+          fast_diag(ig) = fast_diag(ig) + fast_reservoir(ig,ib)
+          slow_diag(ig) = slow_diag(ig) + slow_reservoir(ig,ib)
+          stream_diag(ig) = stream_diag(ig) + stream_reservoir(ig,ib)
+          flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib)
+       ENDDO
+       !
+       fast_diag(ig) = fast_diag(ig)/totarea
+       slow_diag(ig) = slow_diag(ig)/totarea
+       stream_diag(ig) = stream_diag(ig)/totarea
+       flood_diag(ig) = flood_diag(ig)/totarea
+       !
+       ! This is the volume of the lake scaled to the entire grid.
+       ! It would be batter to scale it to the size of the lake
+       ! but this information is not yet available.
+       !
+       lake_diag(ig) = lake_reservoir(ig)/totarea
+       !
+    ENDDO
+    !
+    ! Get from the restart the fluxes we accumulated.
+    !
+    ALLOCATE (floodout_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for floodout_mean','','')
+    var_name = 'floodout_route'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Accumulated flow out of floodplains for routing')
+    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_simple_init_1','Pb in allocate for runoff_mean','','')
+    var_name = 'runoff_route'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Accumulated runoff for routing')
+    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_simple_init_1','Pb in allocate for drainage_mean','','')
+    var_name = 'drainage_route'
+    CALL ioconf_setatt_p('UNITS', 'Kg')
+    CALL ioconf_setatt_p('LONG_NAME','Accumulated drainage for routing')
+    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_simple_init_1','Pb in allocate for transpot_mean','','')
+    var_name = 'transpot_route'
+    CALL ioconf_setatt_p('UNITS', 'Kg/m^2')
+    CALL ioconf_setatt_p('LONG_NAME','Accumulated potential transpiration for routing/irrigation')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., transpot_mean, "gather", nbp_glo, index_g)
+    CALL setvar_p (transpot_mean, val_exp, 'NO_KEYWORD', zero)
+
+    ALLOCATE(precip_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for precip_mean','','')
+    var_name = 'precip_route'
+    CALL ioconf_setatt_p('UNITS', 'Kg/m^2')
+    CALL ioconf_setatt_p('LONG_NAME','Accumulated rain precipitation for irrigation')
+    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_simple_init_1','Pb in allocate for humrel_mean','','')
+    var_name = 'humrel_route'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Mean humrel for irrigation')
+    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_simple_init_1','Pb in allocate for k_litt_mean','','')
+    var_name = 'k_litt_route'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Mean cond. for litter')
+    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_simple_init_1','Pb in allocate for totnobio_mean','','')
+    var_name = 'totnobio_route'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of no bio for irrigation')
+    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_simple_init_1','Pb in allocate for vegtot_mean','','')
+    var_name = 'vegtot_route'
+    CALL ioconf_setatt_p('UNITS', '-')
+    CALL ioconf_setatt_p('LONG_NAME','Last Total fraction of vegetation')
+    CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., vegtot_mean, "gather", nbp_glo, index_g)
+    CALL setvar_p (vegtot_mean, val_exp, 'NO_KEYWORD', un)
+    !
+    !
+    DEALLOCATE(tmp_real_g)
+    !
+    ! Allocate diagnostic variables
+    !
+    ALLOCATE(hydrodiag_loc(nbpt,nbasmax),hydrodiag_glo(nbp_glo,nbasmax),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydrodiag_glo','','')
+    hydrodiag=>hydrodiag_loc
+
+    ALLOCATE(hydroupbasin_loc(nbpt),hydroupbasin_glo(nbp_glo), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_simple_init_1','Pb in allocate for hydroupbasin_glo','','')
+    hydroupbasin=>hydroupbasin_loc
+
+  END SUBROUTINE routing_simple_init_1
+
+
+  !! ================================================================================================================================
+  !! SUBROUTINE 	: routing_simple_clear
+  !!
+  !>\BRIEF         This subroutine deallocates the block memory previously allocated.
+  !!
+  !! DESCRIPTION:  This subroutine deallocates the block memory previously allocated.
+  !!
+  !! RECENT CHANGE(S): None
+  !!
+  !! MAIN OUTPUT VARIABLE(S):
+  !!
+  !! REFERENCES   : None
+  !!
+  !! FLOWCHART    :None
+  !! \n
+  !_ ================================================================================================================================
+
+  SUBROUTINE routing_simple_clear()
+
+    IF (is_omp_root) THEN
+       IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r)
+       IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1)
+       IF (ALLOCATED(fast_reservoir_r)) DEALLOCATE(fast_reservoir_r)
+       IF (ALLOCATED(slow_reservoir_r)) DEALLOCATE(slow_reservoir_r)
+       IF (ALLOCATED(is_lakeinflow_r)) DEALLOCATE(is_lakeinflow_r) 
+       IF (ALLOCATED(is_coastalflow_r)) DEALLOCATE(is_coastalflow_r)    
+       IF (ALLOCATED(is_riverflow_r)) DEALLOCATE(is_riverflow_r)   
+    ENDIF
+
+
+    IF (ALLOCATED(routing_area_loc)) DEALLOCATE(routing_area_loc)
+    IF (ALLOCATED(route_togrid_loc)) DEALLOCATE(route_togrid_loc)
+    IF (ALLOCATED(route_tobasin_loc)) DEALLOCATE(route_tobasin_loc)
+    IF (ALLOCATED(route_nbintobas_loc)) DEALLOCATE(route_nbintobas_loc)
+    IF (ALLOCATED(global_basinid_loc)) DEALLOCATE(global_basinid_loc)
+    IF (ALLOCATED(topo_resid_loc)) DEALLOCATE(topo_resid_loc)
+    IF (ALLOCATED(routing_area_glo)) DEALLOCATE(routing_area_glo)
+    IF (ALLOCATED(route_togrid_glo)) DEALLOCATE(route_togrid_glo)
+    IF (ALLOCATED(route_tobasin_glo)) DEALLOCATE(route_tobasin_glo)
+    IF (ALLOCATED(route_nbintobas_glo)) DEALLOCATE(route_nbintobas_glo)
+    IF (ALLOCATED(global_basinid_glo)) DEALLOCATE(global_basinid_glo)
+    IF (ALLOCATED(topo_resid_glo)) DEALLOCATE(topo_resid_glo)
+    IF (ALLOCATED(fast_reservoir)) DEALLOCATE(fast_reservoir)
+    IF (ALLOCATED(slow_reservoir)) DEALLOCATE(slow_reservoir)
+    IF (ALLOCATED(stream_reservoir)) DEALLOCATE(stream_reservoir)
+    IF (ALLOCATED(flood_reservoir)) DEALLOCATE(flood_reservoir)
+    IF (ALLOCATED(flood_frac_bas)) DEALLOCATE(flood_frac_bas)
+    IF (ALLOCATED(flood_height)) DEALLOCATE(flood_height)
+    IF (ALLOCATED(pond_frac)) DEALLOCATE(pond_frac)
+    IF (ALLOCATED(lake_reservoir)) DEALLOCATE(lake_reservoir)
+    IF (ALLOCATED(pond_reservoir)) DEALLOCATE(pond_reservoir)
+    IF (ALLOCATED(returnflow_mean)) DEALLOCATE(returnflow_mean)
+    IF (ALLOCATED(reinfiltration_mean)) DEALLOCATE(reinfiltration_mean)
+    IF (ALLOCATED(riverflow_mean)) DEALLOCATE(riverflow_mean)
+    IF (ALLOCATED(coastalflow_mean)) DEALLOCATE(coastalflow_mean)
+    IF (ALLOCATED(lakeinflow_mean)) DEALLOCATE(lakeinflow_mean)
+    IF (ALLOCATED(runoff_mean)) DEALLOCATE(runoff_mean)
+    IF (ALLOCATED(floodout_mean)) DEALLOCATE(floodout_mean)
+    IF (ALLOCATED(drainage_mean)) DEALLOCATE(drainage_mean)
+    IF (ALLOCATED(transpot_mean)) DEALLOCATE(transpot_mean)
+    IF (ALLOCATED(precip_mean)) DEALLOCATE(precip_mean)
+    IF (ALLOCATED(humrel_mean)) DEALLOCATE(humrel_mean)
+    IF (ALLOCATED(k_litt_mean)) DEALLOCATE(k_litt_mean)
+    IF (ALLOCATED(totnobio_mean)) DEALLOCATE(totnobio_mean)
+    IF (ALLOCATED(vegtot_mean)) DEALLOCATE(vegtot_mean)
+    IF (ALLOCATED(floodtemp)) DEALLOCATE(floodtemp)
+    IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc)
+    IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo)
+    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(irrigated)) DEALLOCATE(irrigated)
+    IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains)
+    IF (ALLOCATED(swamp)) DEALLOCATE(swamp)
+    IF (ALLOCATED(fast_diag)) DEALLOCATE(fast_diag)
+    IF (ALLOCATED(slow_diag)) DEALLOCATE(slow_diag)
+    IF (ALLOCATED(stream_diag)) DEALLOCATE(stream_diag)
+    IF (ALLOCATED(flood_diag)) DEALLOCATE(flood_diag)
+    IF (ALLOCATED(pond_diag)) DEALLOCATE(pond_diag)
+    IF (ALLOCATED(lake_diag)) DEALLOCATE(lake_diag)
+    IF (ALLOCATED(delsurfstor)) DEALLOCATE(delsurfstor)
+    !
+  END SUBROUTINE routing_simple_clear
+  !
+
+  !! ================================================================================================================================
+  !! SUBROUTINE 	: routing_simple_lake
+  !!
+  !>\BRIEF        : This subroutine stores water in lakes so that it does not cycle through the runoff.
+  !!                For the moment it only works for endoheric lakes but I can be extended in the future.
+  !!
+  !! DESCRIPTION (definitions, functional, design, flags): The return flow to the soil moisture reservoir
+  !! is based on a maximum lake evaporation rate (maxevap_lake). \n
+  !!
+  !! RECENT CHANGE(S): None
+  !!
+  !! MAIN OUTPUT VARIABLE(S):
+  !!
+  !! REFERENCES   : None
+  !!
+  !! FLOWCHART    :None
+  !! \n
+  !_ ================================================================================================================================
+
+  SUBROUTINE routing_simple_lake(nbpt, dt_routing, lakeinflow, humrel, contfrac, return_lakes)
+
+    USE grid, ONLY : area
+    
+    IMPLICIT NONE
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+
+    INTEGER(i_std), INTENT(in) :: nbpt               !! Domain size (unitless)
+    REAL(r_std), INTENT (in)   :: dt_routing         !! Routing time step (s)
+    REAL(r_std), INTENT(in)    :: lakeinflow(nbpt)   !! Water inflow to the lakes (kg/dt)
+    REAL(r_std), INTENT(in)    :: humrel(nbpt)       !! Soil moisture stress, root extraction potential (unitless)
+    REAL(r_std), INTENT(in)    :: contfrac(nbpt)     !! Fraction of land in each grid box (unitless;0-1)
+    
+    !! 0.2 Output variables
+    REAL(r_std), INTENT(out)   :: return_lakes(nbpt) !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
+    
+    !! 0.3 Local variables 
+    INTEGER(i_std)             :: ig                 !! Indices (unitless)
+    REAL(r_std)                :: refill             !!
+    REAL(r_std)                :: total_area         !! Sum of all the surfaces of the basins (m^2)
+
+    !_ ================================================================================================================================
+
+    
+    DO ig=1,nbpt
+       !
+       total_area = area(ig)*contfrac(ig)
+       !
+       lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig)
+       !uptake in Kg/dt
+       refill = MAX(zero, maxevap_lake * (un - humrel(ig)) * dt_routing * total_area)
+       return_lakes(ig) = MIN(refill, lake_reservoir(ig))
+       lake_reservoir(ig) = lake_reservoir(ig) - return_lakes(ig)
+       !Return in Kg/m^2/dt
+       IF ( doswamps ) THEN
+          return_lakes(ig) = return_lakes(ig)/total_area
+       ELSE
+          return_lakes(ig) = zero
+       ENDIF
+       !
+       ! This is the volume of the lake scaled to the entire grid.
+       ! It would be batter to scale it to the size of the lake
+       ! but this information is not yet available.
+       lake_diag(ig) = lake_reservoir(ig)/total_area
+       !
+    ENDDO
+
+  END SUBROUTINE routing_simple_lake
+
+END MODULE routing_simple
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90	(revision 6102)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90	(revision 6102)
@@ -0,0 +1,326 @@
+! ================================================================================================================================
+!  MODULE       : routing_wrapper
+!
+!  CONTACT      : orchidee-help _at_ listes.ipsl.fr
+!
+!  LICENCE      : IPSL (2006)
+!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
+!
+!>\BRIEF          Interfaces to all routing schemes.
+!!
+!!\n DESCRIPTION: This module contains uniformed subroutines called from sechiba. These subroutines make the swich the between 
+!!                the different existing routing modules.
+!!                
+!!                Depending on the key world ROUTING_METHOD set in run.def, this module calls one of the 
+!!                available routing modules: 
+!!                - ROUTING_METOD=standard for the standard routing scheme available in module routing. 
+!!                - ROUTING_METHOD=simple for the routing scheme in module routing_simple. 
+!!
+!! REFERENCE(S) : None
+!!
+!! SVN          :
+!! $HeadURL$
+!! $Date$
+!! $Revision$
+!! \n
+!_ ================================================================================================================================
+
+MODULE routing_wrapper
+
+  USE defprec
+  USE pft_parameters
+  USE grid
+  USE routing
+  USE routing_simple
+
+  IMPLICIT NONE
+
+  CHARACTER(LEN=255), SAVE :: routing_method                      !! 'standard' or 'simple': Character string used to switch between routing modules
+  !$OMP THREADPRIVATE(routing_method) 
+
+  PUBLIC :: routing_wrapper_xios_initialize, routing_wrapper_initialize, &
+            routing_wrapper_main, routing_wrapper_finalize, routing_wrapper_clear 
+  PRIVATE
+
+CONTAINS
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_wrapper_xios_initialize
+!!
+!>\BRIEF	  First initialization phase of the choosen routing module
+!!
+!! DESCRIPTION:	  Read ROUTING_METHOD from run.def and call the xios initialization subroutine from corresponding routing module.
+!!                This subroutine is called before the xios context is closed. 
+!!                It is called from sechiba_initialize only if 1 is activated.
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+  SUBROUTINE routing_wrapper_xios_initialize()
+
+    ! Get ROUTING_METHOD from run.def. Note that this is also done in 
+    ! routing_wrapper_initialize because current subroutine is not alwyas called.
+    routing_method='standard'
+    CALL getin_p("ROUTING_METHOD",routing_method)
+    IF(routing_method=='simple') THEN  
+       CALL routing_simple_xios_initialize
+    ENDIF
+
+  END SUBROUTINE routing_wrapper_xios_initialize
+
+
+
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_wrapper_initialize
+!!
+!>\BRIEF	  Initialize the choosen routing module
+!!
+!! DESCRIPTION:	  Read ROUTING_METHOD from run.def and call the initialization subroutine from corresponding routing module
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+  SUBROUTINE routing_wrapper_initialize( &
+       kjit,        nbpt,           index,                 &
+       rest_id,     hist_id,        hist2_id,   lalo,      &
+       neighbours,  resolution,     contfrac,   stempdiag, &
+       returnflow,  reinfiltration, irrigation, riverflow, &
+       coastalflow, flood_frac,     flood_res )
+
+
+    !! 0.1 Input variables
+    INTEGER(i_std), INTENT(in)     :: kjit                 !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)     :: nbpt                 !! Domain size (unitless)
+    INTEGER(i_std), INTENT(in)     :: index(nbpt)          !! Indices of the points on the map (unitless)
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    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)
+    REAL(r_std), INTENT(in)        :: lalo(nbpt,2)         !! Vector of latitude and longitudes (beware of the order !)
+
+    INTEGER(i_std), INTENT(in)     :: neighbours(nbpt,8)   !! Vector of neighbours for each grid point 
+                                                           !! (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
+    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)
+    REAL(r_std), INTENT(in)        :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile
+
+    !! 0.2 Output variables
+    REAL(r_std), INTENT(out)       :: returnflow(nbpt)     !! The water flow from lakes and swamps which returns to the grid box.
+                                                           !! This water will go back into the hydrol or hydrolc module to allow re-evaporation (kg/m^2/dt)
+    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)
+    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)
+    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)
+    REAL(r_std), INTENT(out)       :: flood_frac(nbpt)     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(out)       :: flood_res(nbpt)      !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+
+    !_ ================================================================================================================================
+
+    !! 1. Get routing_method from run.def
+    !!    This variable will switch between the existing modules for the routing scheme.
+
+    !Config Key   = ROUTING_METHOD
+    !Config Desc  = Choice of routing module to be used
+    !Config If    = RIVER_ROUTING=T
+    !Config Def   = standard
+    !Config Help  = Possible options are standard and simple
+    !Config Units = character string
+
+    routing_method='standard'
+    CALL getin_p("ROUTING_METHOD",routing_method)
+
+
+    !! 2. Initialize the choosen routing module
+    IF (routing_method == 'standard') THEN
+
+       CALL routing_initialize(  kjit,        nbpt,           index,                 &
+                                 rest_id,     hist_id,        hist2_id,   lalo,      &
+                                 neighbours,  resolution,     contfrac,   stempdiag, &
+                                 returnflow,  reinfiltration, irrigation, riverflow, &
+                                 coastalflow, flood_frac,     flood_res )
+
+    ELSE IF(routing_method== 'simple') THEN 
+
+       CALL routing_simple_initialize(    kjit,        nbpt,           index,                 &
+                                          rest_id,     hist_id,        hist2_id,   lalo,      &
+                                          neighbours,  resolution,     contfrac,   stempdiag, &
+                                          returnflow,  reinfiltration, irrigation, riverflow, &
+                                          coastalflow, flood_frac,     flood_res )
+
+       riverflow(:) = zero
+       coastalflow(:) = zero
+       returnflow(:) = zero
+       reinfiltration(:) = zero
+       irrigation(:) = zero
+       flood_frac(:) = zero
+       flood_res(:) = zero
+
+    ELSE
+       ! Bad choice of routing_method. Exit the model now. 
+       WRITE(numout,*) 'Following routing method is not implemented, ROUTING_METHOD=',routing_method
+       CALL ipslerr_p(3,'routing_wrapper_inititalize','ROUTING_METHOD can only be standard or simple','Error in run.def','')
+    ENDIF
+
+  END SUBROUTINE routing_wrapper_initialize
+
+
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_wrapper_main
+!!
+!>\BRIEF	  Call the main subroutine for the choosen routing module
+!!
+!! DESCRIPTION:	  Call the main subroutine for the choosen routing module according to ROUTING_METHOD
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+  SUBROUTINE routing_wrapper_main(kjit, nbpt, index, &
+       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)
+
+
+
+    IMPLICIT NONE
+
+    !! 0.1 Input variables
+    INTEGER(i_std), INTENT(in)     :: kjit                 !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)     :: nbpt                 !! Domain size (unitless)
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    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)
+    INTEGER(i_std), INTENT(in)     :: 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 !)
+    INTEGER(i_std), INTENT(in)     :: neighbours(nbpt,NbNeighb)   !! Vector of neighbours for each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW) (unitless)
+    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)
+    REAL(r_std), INTENT(in)        :: totfrac_nobio(nbpt)  !! Total fraction of no-vegetation (continental ice, lakes ...) (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: veget_max(nbpt,nvm)  !! Maximal fraction of vegetation (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: floodout(nbpt)       !! Grid-point flow out of floodplains (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: runoff(nbpt)         !! Grid-point runoff (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: drainage(nbpt)       !! Grid-point drainage (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: transpot(nbpt,nvm)   !! Potential transpiration of the vegetation (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: precip_rain(nbpt)    !! Rainfall (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: k_litt(nbpt)         !! Averaged conductivity for saturated infiltration in the 'litter' layer (kg/m^2/dt)
+    REAL(r_std), INTENT(in)        :: humrel(nbpt,nvm)     !! Soil moisture stress, root extraction potential (unitless)
+    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)
+
+    !! 0.2 Output variables
+    REAL(r_std), INTENT(out)       :: returnflow(nbpt)     !! The water flow from lakes and swamps which returns to the grid box.
+    !! This water will go back into the hydrol module to allow re-evaporation (kg/m^2/dt)
+    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)
+    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)
+    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)
+    REAL(r_std), INTENT(out)       :: flood_frac(nbpt)     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(out)       :: flood_res(nbpt)      !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+    !_ ================================================================================================================================
+
+    !! 1. Call the main subroutine from the routing module corresponding to the choice of ROUTING_METHOD
+
+    IF (routing_method=='standard') THEN
+
+       CALL routing_main (kjit, nbpt, index, &
+            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)
+
+    ELSE IF(routing_method=='simple') THEN 
+
+       CALL routing_simple_main (kjit, nbpt, index, &
+            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)
+    ENDIF
+
+
+  END SUBROUTINE routing_wrapper_main
+
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_wrapper_finalize
+!!
+!>\BRIEF	  Call the finalization subroutine for the choosen routing module
+!!
+!! DESCRIPTION:	  Call the subroutine for finalization for the choosen routing module according to ROUTING_METHOD
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+  SUBROUTINE routing_wrapper_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
+
+    IMPLICIT NONE
+    !! 0.1 Input variables
+    INTEGER(i_std), INTENT(in)     :: kjit                 !! Time step number (unitless)
+    INTEGER(i_std), INTENT(in)     :: nbpt                 !! Domain size (unitless)
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    REAL(r_std), INTENT(in)        :: flood_frac(nbpt)     !! Flooded fraction of the grid box (unitless;0-1)
+    REAL(r_std), INTENT(in)        :: flood_res(nbpt)      !! Diagnostic of water amount in the floodplains reservoir (kg)
+
+    !_ ================================================================================================================================
+
+    !! 1. Call the finalization subroutine from the routing module corresponding to the choice of ROUTING_METHOD
+
+    IF (routing_method=='standard') THEN
+
+       CALL routing_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
+
+    ELSE IF(routing_method=='simple') THEN 
+
+       CALL routing_simple_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
+
+    ENDIF
+
+  END SUBROUTINE routing_wrapper_finalize
+
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_wrapper_clear
+!!
+!>\BRIEF	  Call the clear subroutine for the choosen routing module
+!!
+!! DESCRIPTION:	  Call the clear subroutine for the choosen routing module according to ROUTING_METHOD
+!!
+!! RECENT CHANGE(S): None
+!!
+!! REFERENCE(S): None
+!! 
+!! FLOWCHART: None
+!! \n
+!_ ==============================================================================================================================
+  SUBROUTINE routing_wrapper_clear
+
+    IF (routing_method=='standard') THEN
+
+       CALL routing_clear
+
+    ELSE IF(routing_method=='simple') THEN 
+
+       CALL routing_simple_clear
+
+    ENDIF
+
+  END SUBROUTINE routing_wrapper_clear
+
+END MODULE routing_wrapper
Index: /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90
===================================================================
--- /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 6101)
+++ /branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/sechiba.f90	(revision 6102)
@@ -45,5 +45,5 @@
   USE sechiba_io_p
   USE slowproc
-  USE routing
+  USE routing_wrapper
   USE ioipsl_para
   USE chemistry
@@ -295,4 +295,5 @@
       CALL chemistry_xios_initialize
       CALL thermosoil_xios_initialize
+      CALL routing_wrapper_xios_initialize
     END IF
     IF (printlev_loc>=3) WRITE(numout,*) 'End sechiba_xios_initialize'
@@ -492,6 +493,8 @@
     !! 1.12 Initialize river routing
     IF ( river_routing .AND. nbp_glo .GT. 1) THEN
+
        !! 1.12.1 Initialize river routing
-       CALL routing_initialize( kjit,        kjpindex,       index,                 &
+       CALL routing_wrapper_initialize( &
+            kjit,        kjpindex,       index,                 &
             rest_id,     hist_id,        hist2_id,   lalo,      &
             neighbours,  resolution,     contfrac,   stempdiag, &
@@ -762,5 +765,5 @@
     IF ( river_routing .AND. nbp_glo .GT. 1) THEN
        !! 8.1 River routing
-       CALL routing_main (kjit, kjpindex, index, &
+       CALL routing_wrapper_main (kjit, kjpindex, index, &
             & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, &
             & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, &
@@ -1364,5 +1367,5 @@
     IF ( river_routing .AND. nbp_glo .GT. 1) THEN
        !! 6.1 Call river routing to write restart files 
-       CALL routing_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res )
+       CALL routing_wrapper_finalize( kjit, kjpindex, rest_id, flood_frac, flood_res )
     ELSE
        !! 6.2 No routing, set variables to zero
@@ -1960,5 +1963,5 @@
     CALL thermosoil_clear
     CALL condveg_clear 
-    CALL routing_clear
+    CALL routing_wrapper_clear
 
   END SUBROUTINE sechiba_clear
