Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native.f90	(revision 8227)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native.f90	(revision 8227)
@@ -0,0 +1,543 @@
+! =================================================================================================================================
+! 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: svn://forge.ipsl.jussieu.fr/orchidee/perso/pedro.arboleda/ORCHIDEE/src_sechiba/routing_simple.f90 $
+!! $Date: 2022-06-07 19:03:47 +0200 (Tue, 07 Jun 2022) $
+!! $Revision: 7644 $
+!! \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_native_mod
+
+  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_native_main, routing_native_xios_initialize
+  PUBLIC :: routing_native_initialize, routing_native_finalize, routing_native_clear
+
+  REAL,SAVE :: dt_routing
+  !$OMP THREADPRIVATE(dt_routing) 
+
+  INTEGER,SAVE :: nbpt
+  !$OMP THREADPRIVATE(nbpt) 
+  
+  INTEGER,SAVE :: nbpt_r
+  !$OMP THREADPRIVATE(nbpt_r) 
+  
+  LOGICAL,SAVE :: dofloodinfilt
+  !$OMP THREADPRIVATE(dofloodinfilt)   
+
+  LOGICAL,SAVE :: doswamps
+  !$OMP THREADPRIVATE(doswamps)  
+
+  LOGICAL,SAVE :: doponds
+  !$OMP THREADPRIVATE(doponds)  
+
+
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:)               :: returnflow_mean              !! Mean water flow from lakes and swamps which returns to the grid box.
+ !$OMP THREADPRIVATE(returnflow_mean)
+
+ REAL, SAVE :: time_counter = 0
+ !$OMP THREADPRIVATE(time_counter)  
+   
+CONTAINS
+
+!!  =============================================================================================================================
+!! SUBROUTINE:    routing_native_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_native_xios_initialize
+    USE xios
+    USE routing_native_flow_mod
+    IMPLICIT NONE     
+
+     CALL routing_flow_xios_initialize
+     
+   END SUBROUTINE routing_native_xios_initialize
+
+
+
+
+
+
+
+  !!  =============================================================================================================================
+  !! 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_native_initialize( kjit,       nbpt,           index,                 &
+       rest_id,     hist_id,        hist2_id,   lalo,      &
+       neighbours,  resolution,     contfrac,   stempdiag, &
+       returnflow,  reinfiltration, irrigation, riverflow, &
+       coastalflow, flood_frac,     flood_res , irrigated_next)
+    
+    USE routing_native_flow_mod, ONLY : routing_flow_initialize
+    USE routing_native_irrig_mod, ONLY : irrigation_initialize
+    USE routing_native_lake_mod, ONLY : routing_lake_initialize
+    USE routing_native_para, ONLY: routing_para_initialize
+    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)
+    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)  !! Dynamic irrig. area, calculated in slowproc and passed to routing!
+
+    !_ ================================================================================================================================
+    
+    
+      CALL routing_para_initialize 
+      CALL routing_native_init_local(kjit, rest_id, nbpt, contfrac)
+      CALL routing_flow_initialize(kjit, rest_id, nbpt, dt_routing, contfrac, nbpt_r, riverflow, coastalflow)
+      CALL routing_lake_initialize(kjit, rest_id, nbpt, contfrac)
+      CALL irrigation_initialize(kjit, rest_id, nbpt, nbpt_r,irrigated_next )
+      
+      reinfiltration(:)=0
+      irrigation(:)=0
+      flood_frac(:)=0
+      flood_res(:)=0
+      returnflow(:)=returnflow_mean(:)
+
+    END SUBROUTINE routing_native_initialize
+    
+    
+  !! ================================================================================================================================
+  !! SUBROUTINE 	: routing_simple_initialize
+  !!
+  !>\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_native_init_local(kjit, rest_id, nbpt_, contfrac)
+  USE time, ONLY : dt_sechiba 
+  USE routing_native_flow_mod   
+  IMPLICIT NONE
+    INTEGER(i_std),INTENT(in)      :: kjit
+    INTEGER(i_std),INTENT(in)      :: rest_id              !! Restart file identifier (unitless)
+    INTEGER,INTENT(in) :: nbpt_             !! nb points native grid
+    REAL,INTENT(in)    :: contfrac(nbpt)   !! fraction of land
+
+    INTEGER :: ier    
+    CHARACTER(LEN=80)   :: var_name       !! To store variables names for I/O (unitless)
+  
+    REAL(r_std) ::ratio
+    
+    nbpt = nbpt_
+
+    ALLOCATE (returnflow_mean(nbpt), stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','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)
+    
+    !_ ================================================================================================================================
+    !
+    !
+    ! 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   = 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)
+
+    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
+    
+  END SUBROUTINE routing_native_init_local
+  
+
+
+  !! ================================================================================================================================
+  !! 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_native_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 ,&
+                                 soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) 
+
+    USE routing_native_flow_mod,ONLY : routing_flow_make_mean, routing_flow_main, riverflow_mean, coastalflow_mean, lakeinflow_mean, &
+                                routing_flow_get, routing_flow_set
+    USE routing_native_irrig_mod, ONLY: irrigation_mean_make, irrigation_main, irrigation_get
+    USE routing_native_lake_mod, ONLY: routing_lake_mean_make, routing_lake_route_coast, routing_lake_main
+    USE routing_native_para, ONLY: routing_para_initialize
+    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)
+    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)
+    REAL(r_std), INTENT(in)        :: root_deficit(nbpt)   !! soil water deficit
+    REAL(r_std), INTENT(in)        :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless)
+    REAL(r_std), INTENT(in)        :: irrig_frac(nbpt)     !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
+    REAL(r_std), INTENT(in)        :: irrigated_next (nbpt)!! Dynamic irrig. area, calculated in slowproc and passed to routing
+    REAL(r_std), INTENT(in)        :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac
+
+    !! 0.2 Output variables
+    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
+    REAL(r_std), DIMENSION(nbpt)   :: return_lakes         !! Water from lakes flowing back into soil moisture (kg/m^2/dt)
+    REAL(r_std)                    :: lakeinflow(nbpt)
+  
+  
+    CALL routing_flow_make_mean(runoff, drainage)
+    CALL irrigation_mean_make(dt_routing, veget_max, humrel, transpot, runoff , precip_rain)
+    CALL routing_lake_mean_make(dt_routing, humrel, veget_max)
+       
+    time_counter = time_counter + dt_sechiba
+    
+    reinfiltration(:)= 0 ! for now
+    
+    ! If the time has come we do the routing.
+    IF ( NINT(time_counter) .GE. NINT(dt_routing) ) THEN
+       returnflow_mean(:)=0
+ 
+       CALL routing_flow_main(dt_routing)
+
+       CALL routing_flow_get(coastalflow_mean=coastalflow)
+       CALL routing_lake_route_coast(contfrac,coastalflow)
+       CALL routing_flow_set(coastalflow_mean=coastalflow)
+      
+       CALL irrigation_main( dt_routing, reinfiltration, irrigated_next, irrig_frac, root_deficit, soiltile, &
+                             fraction_aeirrig_sw )
+       
+       CALL routing_flow_get(lakeinflow_mean=lakeinflow)
+       CALL routing_lake_main(dt_routing, contfrac, lakeinflow, return_lakes)
+       CALL routing_flow_set(lakeinflow_mean=lakeinflow)
+       
+       returnflow_mean(:)=returnflow_mean(:)+return_lakes(:)
+       
+       time_counter = 0
+    ENDIF
+
+    CALL routing_flow_get(riverflow_mean=riverflow, coastalflow_mean=coastalflow, lakeinflow_mean=lakeinflow)
+    CALL irrigation_get(irrigation_mean=irrigation)
+
+
+    returnflow  = returnflow/dt_routing*dt_sechiba
+    riverflow   = riverflow/dt_routing*dt_sechiba
+    coastalflow = coastalflow/dt_routing*dt_sechiba
+    irrigation = irrigation/dt_routing*dt_sechiba
+    
+    CALL xios_send_field("irrigation", irrigation/dt_sechiba)
+    
+    returnflow(:) = 0
+    reinfiltration(:) =0
+    flood_frac(:) = 0
+    flood_res(:) = 0
+ 
+  END SUBROUTINE routing_native_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_native_finalize(kjit, nbpt_, rest_id, flood_frac, flood_res )
+    USE routing_native_flow_mod
+    USE routing_native_lake_mod
+    USE routing_native_irrig_mod
+    IMPLICIT NONE
+    INTEGER, INTENT(IN)          :: kjit
+    INTEGER,INTENT(IN)           :: nbpt_
+    INTEGER, INTENT(IN)          :: rest_id
+    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)
+
+    CALL routing_native_local_finalize(kjit, rest_id)
+    CALL routing_flow_finalize(kjit, rest_id)
+    CALL routing_lake_finalize(kjit, rest_id)
+    CALL irrigation_finalize(kjit, rest_id)
+
+  END SUBROUTINE routing_native_finalize
+
+  SUBROUTINE routing_native_local_finalize(kjit,rest_id)
+    USE ioipsl_para
+    USE grid
+    IMPLICIT NONE
+    INTEGER, INTENT(IN) :: kjit
+    INTEGER, INTENT(IN) :: rest_id
+
+    CALL restput_p (rest_id, 'returnflow', nbp_glo, 1, 1, kjit, returnflow_mean, 'scatter',  nbp_glo, index_g)
+    DEALLOCATE(returnflow_mean)
+
+  END SUBROUTINE routing_native_local_finalize
+
+
+  SUBROUTINE routing_native_clear
+    USE xios
+    IMPLICIT NONE
+
+
+
+  END SUBROUTINE routing_native_clear
+  
+END MODULE routing_native_mod
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flood_pond.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flood_pond.f90	(revision 8227)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flood_pond.f90	(revision 8227)
@@ -0,0 +1,11 @@
+MODULE routing_native_flood_pond_mod
+
+
+
+
+
+CONTAINS
+
+
+
+END MODULE routing_native_flood_pond_mod
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90	(revision 8227)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90	(revision 8227)
@@ -0,0 +1,1764 @@
+MODULE routing_native_flow_mod
+
+  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_flow_xios_initialize, routing_flow_set, routing_flow_get, routing_flow_main
+  PUBLIC :: routing_flow_initialize, routing_flow_finalize, routing_flow_clear, routing_flow_make_mean
+  PUBLIC :: compute_coastline
+    
+  INTEGER,SAVE :: nbpt
+  !$OMP THREADPRIVATE(nbpt)
+   
+  !! PARAMETERS
+  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,ALLOCATABLE :: runoff_mean(:)               
+  !$OMP THREADPRIVATE(runoff_mean)
+
+  REAL(r_std),SAVE,ALLOCATABLE :: drainage_mean(:)               
+  !$OMP THREADPRIVATE(drainage_mean)
+
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: fast_reservoir_r(:)         !! Water amount in the fast reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(fast_reservoir_r)
+
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:)         !! Water amount in the slow reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(slow_reservoir_r)
+
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_reservoir_r(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(stream_reservoir_r) 
+  
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(riverflow_mean) 
+
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: coastalflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(coastalflow_mean) 
+
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  !$OMP THREADPRIVATE(lakeinflow_mean) 
+
+  !
+  ! 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)
+  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_coastline(:)           !! is coastline point - (local native grid) 
+  !$OMP THREADPRIVATE(is_coastline)
+  LOGICAL,SAVE,ALLOCATABLE     :: is_streamflow_r(:)          !! is stream flow point - (local routing grid) 
+  !$OMP THREADPRIVATE(is_streamflow_r)
+  LOGICAL,SAVE,ALLOCATABLE,PUBLIC,PROTECTED    :: 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                 :: total_coast_points        !! global number of coast point - (local native grid)                   
+  !$OMP THREADPRIVATE(total_coast_points)
+  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_in(:)           !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
+  !$OMP THREADPRIVATE(routing_weight_in)
+  REAL(r_std),SAVE,ALLOCATABLE :: unrouted_weight(:)           !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid)
+  !$OMP THREADPRIVATE(unrouted_weight)
+
+  REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_coast_r(:)   
+  !$OMP THREADPRIVATE(weight_coast_to_coast_r)
+  
+  REAL(r_std),SAVE,ALLOCATABLE :: weight_coast_to_lake_r(:)   
+  !$OMP THREADPRIVATE(weight_coast_to_lake_r)
+  
+  REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_coast_r(:)   
+  !$OMP THREADPRIVATE(weight_lake_to_coast_r)
+  
+  REAL(r_std),SAVE,ALLOCATABLE :: weight_lake_to_lake_r(:)   
+  !$OMP THREADPRIVATE(weight_lake_to_lake_r)
+
+  REAL(r_std),SAVE,ALLOCATABLE :: basins_area_r(:)        
+  !$OMP THREADPRIVATE(basins_area_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 :: nb_station
+  CHARACTER(LEN=60),ALLOCATABLE :: station(:)  
+  REAL,ALLOCATABLE :: station_lonlat(:,:)  
+  INTEGER,ALLOCATABLE :: station_index(:)
+  INTEGER :: station_ts = 0
+
+
+!  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_flow_get(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, &
+                              lakeinflow_mean)
+  IMPLICIT NONE
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean(:)
+
+    CALL routing_flow_get_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean)
+  END SUBROUTINE routing_flow_get
+
+  SUBROUTINE routing_flow_get_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, &
+                               coastalflow_mean_, lakeinflow_mean_)
+  IMPLICIT NONE
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: slow_reservoir_r_(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: fast_reservoir_r_(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: stream_reservoir_r_(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: riverflow_mean_(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: coastalflow_mean_(:)
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: lakeinflow_mean_(:)
+
+    IF (PRESENT(slow_reservoir_r_))   slow_reservoir_r_ = slow_reservoir_r
+    IF (PRESENT(fast_reservoir_r_))   fast_reservoir_r_ = fast_reservoir_r
+    IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r_ = stream_reservoir_r
+    IF (PRESENT(riverflow_mean_))     riverflow_mean_ = riverflow_mean
+    IF (PRESENT(coastalflow_mean_))   coastalflow_mean_ = coastalflow_mean
+    IF (PRESENT(lakeinflow_mean_))      lakeinflow_mean_ = lakeinflow_mean
+
+  END SUBROUTINE routing_flow_get_
+
+  SUBROUTINE routing_flow_set(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, &
+                             lakeinflow_mean)
+  IMPLICIT NONE
+    REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean(:)
+    
+    CALL routing_flow_set_(slow_reservoir_r, fast_reservoir_r, stream_reservoir_r, riverflow_mean, coastalflow_mean, lakeinflow_mean)
+
+  END SUBROUTINE routing_flow_set
+
+  SUBROUTINE routing_flow_set_(slow_reservoir_r_, fast_reservoir_r_, stream_reservoir_r_, riverflow_mean_, &
+                               coastalflow_mean_, lakeinflow_mean_)
+  IMPLICIT NONE
+    REAL(r_std),OPTIONAL, INTENT(IN) :: slow_reservoir_r_(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: fast_reservoir_r_(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: stream_reservoir_r_(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: riverflow_mean_(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: coastalflow_mean_(:)
+    REAL(r_std),OPTIONAL, INTENT(IN) :: lakeinflow_mean_(:)
+
+    IF (PRESENT(slow_reservoir_r_))   slow_reservoir_r = slow_reservoir_r_
+    IF (PRESENT(fast_reservoir_r_))   fast_reservoir_r = fast_reservoir_r_
+    IF (PRESENT(stream_reservoir_r_)) stream_reservoir_r = stream_reservoir_r_
+    IF (PRESENT(riverflow_mean_))     riverflow_mean = riverflow_mean_
+    IF (PRESENT(coastalflow_mean_))   coastalflow_mean = coastalflow_mean_
+    IF (PRESENT(lakeinflow_mean_))     lakeinflow_mean = lakeinflow_mean_
+
+  END SUBROUTINE routing_flow_set_  
+!!  =============================================================================================================================
+!! 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_flow_xios_initialize
+    USE xios
+    USE routing, ONLY : routing_names
+    IMPLICIT NONE     
+
+    INTEGER(i_std) ::ib
+    INTEGER :: nbasmax=1
+    !! 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
+
+    !! Define XIOS axis size needed for the model output
+    ! Add axis for homogeneity between all routing schemes, these dimensions are currently not used in this scheme
+    CALL xios_orchidee_addaxis("nbhtu", nbasmax, (/(REAL(ib,r_std),ib=1,nbasmax)/))
+    CALL xios_orchidee_addaxis("nbasmon", 1, (/(REAL(ib,r_std),ib=1,1)/))
+
+  END SUBROUTINE routing_flow_xios_initialize
+
+
+  SUBROUTINE routing_flow_initialize(kjit, rest_id, nbpt_, dt_routing, contfrac, nbpt_r_, riverflow, coastalflow)
+  IMPLICIT NONE
+    INTEGER ,INTENT(IN)                     :: kjit           
+    INTEGER ,INTENT(IN)                     :: rest_id           
+    INTEGER, INTENT(IN)                     :: nbpt_             !! nb points orchidee grid
+    REAL(r_std), INTENT(IN)                 :: dt_routing                       
+    REAL(r_std), INTENT(IN)                 :: contfrac(nbpt_)   !! fraction of land
+    INTEGER, INTENT(OUT)                    :: nbpt_r_             !! nb points routing native grid
+    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)
+    INTEGER :: ier    
+    CHARACTER(LEN=80)   :: var_name       !! To store variables names for I/O (unitless)
+     
+    nbpt = nbpt_ 
+    CALL routing_flow_init_local(contfrac, nbpt_r_)
+    nbpt_r = nbpt_r_
+    CALL routing_flow_init_mean(kjit, rest_id)
+    CALL initialize_stations(dt_routing)
+    !
+    ! 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 (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)
+
+    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)
+    
+    riverflow(:)   =  riverflow_mean(:)
+    coastalflow(:) =  coastalflow_mean(:)
+
+ END SUBROUTINE routing_flow_initialize
+
+
+  SUBROUTINE compute_coastline(contfrac, coastline)
+  USE mod_orchidee_para
+  USE grid
+  IMPLICIT NONE
+    REAL(r_std),INTENT(IN)     :: contfrac(nbpt)    ! INPUT   : fraction of continent (unitless)
+    LOGICAL,INTENT(OUT)        :: coastline(nbpt)   ! OUTPUT  : coastline mask (true : on coastaline, else false
+    
+    REAL(r_std) :: contfrac_glo(nbp_glo)
+    REAL(r_std) :: contfrac2D_glo(iim_g*jjm_g)
+    LOGICAL     :: coastline_glo(nbp_glo)
+    LOGICAL     :: coastline2D_glo(iim_g*jjm_g)
+    INTEGER(i_std) :: i,j,ij,next_i,next_j,next_ij,m,k
+    REAL(r_std), PARAMETER :: epsilon=1e-5
+    LOGICAL :: is_periodic
+    
+    is_periodic = global
+    CALL gather(contfrac,contfrac_glo)
+    contfrac2D_glo(:)=0
+    DO i=1,nbp_glo
+      contfrac2D_glo(index_g(i)) = contfrac_glo(i)
+    ENDDO
+    
+    IF (is_mpi_root .AND. is_omp_root) THEN
+      
+      coastline2D_glo(:)=.FALSE.
+      DO j=1,jjm_g
+        DO i=1,iim_g
+          ij=(j-1)*iim_g+i
+          IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN
+            coastline2D_glo(ij)=.TRUE.
+          ELSE IF (contfrac2D_glo(ij)>=1-epsilon .AND. grid_type/=unstructured ) THEN
+            DO k=-1,1
+              DO m=-1,1
+                IF (k==0 .AND. m==0) CYCLE
+                next_i=i+k
+                next_j=j+m
+                
+                IF (next_i==0) THEN ! manage periodicity
+                  IF (is_periodic) THEN
+                    next_i=iim_g
+                  ELSE 
+                    coastline2D_glo(ij)=.TRUE.
+                    CYCLE
+                  ENDIF
+                ENDIF 
+
+                IF (next_i==iim_g+1) THEN! manage periodicity
+                  IF (is_periodic) THEN
+                    next_i=1 !! periodic but not true for limited area 
+                  ELSE
+                    coastline2D_glo(ij)=.TRUE.
+                    CYCLE
+                  ENDIF
+                ENDIF
+                IF (next_j==0 .OR. next_j==jjm_g+1) THEN
+                  IF (.NOT. is_periodic) coastline2D_glo(ij)=.TRUE.
+                  CYCLE
+                ENDIF
+                next_ij =  (next_j-1)*iim_g+next_i    
+                IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE.
+              ENDDO
+            ENDDO
+          ENDIF    
+        ENDDO
+      ENDDO
+    
+      DO i=1,nbp_glo
+        coastline_glo(i) = coastline2D_glo(index_g(i))
+      ENDDO
+
+    ENDIF
+    
+    CALL scatter(coastline_glo, coastline)
+    
+  END SUBROUTINE compute_coastline
+
+  
+     
+  SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r)
+  USE xios
+  USE grid, ONLY : global
+  IMPLICIT NONE
+    INCLUDE "mpif.h"
+    INTEGER, INTENT(IN)        :: ni                      ! INPUT : size i (longitude) of the local routing domain
+    INTEGER, INTENT(IN)        :: nj                      ! INPUT : size i (latitude)  of the local routing domain
+    REAL(r_std), INTENT(IN)    :: contfrac(nbp_mpi)       ! INPUT : continental fraction (unittless)
+    LOGICAL, INTENT(IN)        :: coastline(nbp_mpi)      ! INPUT/OUTPUT : coastline mask
+    REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj)           ! INPUT/OUTPUT : diection of flow, which will be modified by the routine
+    REAL(r_std), INTENT(INOUT)    :: trip_extended_r(ni,nj)  ! INPUT  : direction of flow extended to ocean points
+    REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj)        ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean
+    
+    REAL(r_std)                :: contfrac_r(ni,nj)
+    REAL(r_std)                :: mask_coast(nbp_mpi)
+    REAL(r_std)                :: frac_coast_r(ni,nj)
+    REAL(r_std)                :: trip_extended_rp1(0:ni+1,0:nj+1)
+    REAL(r_std)                :: trip_rp1(0:ni+1,0:nj+1)
+    REAL(r_std)                :: state_r(ni,nj)
+    REAL(r_std)                :: state_rp1(0:ni+1,0:nj+1)
+    
+    INTEGER, PARAMETER ::       is_ter=0
+    INTEGER, PARAMETER ::       is_coast=1
+    INTEGER, PARAMETER ::       is_oce=2
+    INTEGER            :: updated
+    INTEGER            :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr
+    INTEGER            :: ibegin, jbegin, ni_glo, nj_glo
+    REAL(r_std)        :: lonvalue_1d(ni), latvalue_1d(nj)
+    REAL(r_std),PARAMETER :: epsilon=1e-3
+ 
+    TYPE(xios_duration) :: ts 
+    TYPE(xios_domain) :: domain_hdl
+    TYPE(xios_domaingroup) :: domain_def_hdl
+    TYPE(xios_field) :: field_hdl
+    TYPE(xios_fieldgroup) :: field_def_hdl
+    TYPE(xios_file) :: file_hdl
+    TYPE(xios_filegroup) :: file_def_hdl
+    TYPE(xios_expand_domain) :: domain_expand_hdl
+    TYPE(xios_date) :: start_date, time_origin
+    CHARACTER(LEN=20)  :: calendar_type
+    LOGICAL :: true=.TRUE.
+
+
+      CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo)    ! get routing domain dimension
+      CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d)    ! get routing domain dimension
+      
+      ! manage non periodic boundary in case of LAM      
+      IF (.NOT. global) THEN
+        IF (ibegin==0) THEN
+          i=1 ;
+          DO j=1,nj
+            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
+          ENDDO
+        ENDIF
+
+        IF (ibegin+ni==ni_glo) THEN
+          i=ni ;
+          DO j=1,nj
+            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
+          ENDDO
+        ENDIF
+        
+        IF (jbegin==0) THEN
+          j=1 ;
+          DO i=1,ni
+            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
+          ENDDO
+        ENDIF
+
+        IF (jbegin+nj==nj_glo) THEN
+          j=nj ;
+          DO i=1,ni
+            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = 97
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = 97
+          ENDDO
+        ENDIF
+      ENDIF
+
+
+      CALL xios_send_field("routing_contfrac",contfrac) 
+      CALL xios_recv_field("routing_contfrac_r", contfrac_r)
+
+      CALL xios_send_field("trip_ext_r",trip_extended_r)
+      CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1)
+
+      mask_coast=0
+      WHERE(coastline) mask_coast=1
+      CALL xios_send_field("mask_coastline",mask_coast)
+      CALL xios_recv_field("frac_coastline_r",frac_coast_r)
+      
+      DO j=1,nj
+        DO i=1,ni
+          IF (contfrac_r(i,j)<epsilon .OR. contfrac_r(i,j)> 2 ) THEN
+            state_r(i,j)=is_oce
+          ELSE IF (contfrac_r(i,j)>1-epsilon) THEN
+            state_r(i,j)=is_ter
+          ELSE 
+            state_r(i,j)=is_coast
+          ENDIF
+          
+          IF (frac_coast_r(i,j)>epsilon .AND. frac_coast_r(i,j)<=1+epsilon ) state_r(i,j)=is_coast
+        ENDDO
+      ENDDO
+
+
+      DO j=1,nj
+        DO i=1,ni
+          IF (trip_r(i,j)>=100) THEN
+            IF (state_r(i,j) /= is_oce) trip_r(i,j)=trip_extended_r(i,j)
+          ENDIF
+        ENDDO
+      ENDDO          
+
+! create a new XIOS context to solve iterative process of correction
+! so we can use filter to update halo
+!      CALL xios_get_calendar_type(calendar_type)
+!      CALL xios_get_start_date(start_date)
+!      CALL xios_get_time_origin(time_origin)
+!
+!      CALL xios_context_initialize("orchidee_routing_correction",MPI_COMM_ORCH)
+!      CALL xios_orchidee_change_context("orchidee_routing_correction")
+!     ts.second=1
+!     calendar_type="gregorian" 
+!     CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts )
+
+!     CALL xios_get_handle("domain_definition",domain_def_hdl)
+!     CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain") 
+!     CALL xios_set_attr(domain_hdl, type="rectilinear", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo)    ! set routing domain dimension
+!     CALL xios_set_attr(domain_hdl, lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d)                           ! set routing domain dimension
+!     CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain_expand") 
+!     CALL xios_set_attr(domain_hdl, domain_ref="routing_domain")
+!     CALL xios_add_child(domain_hdl,domain_expand_hdl)
+!     CALL xios_set_attr(domain_expand_hdl, type="edge", i_periodic=.TRUE., j_periodic=.TRUE.)
+!  
+!     CALL xios_get_handle("field_definition",field_def_hdl)
+!     CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_init")
+!     CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
+!     CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_final")
+!     CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
+
+!     CALL xios_add_child(field_def_hdl,field_hdl,"trip_r")
+!     CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
+!     CALL xios_add_child(field_def_hdl,field_hdl,"trip_rp1")
+!     CALL xios_set_attr(field_hdl, field_ref="trip_r", domain_ref="routing_domain_expand", read_access=.TRUE.)
+
+!     CALL xios_add_child(field_def_hdl,field_hdl,"state_r")
+!     CALL xios_set_attr(field_hdl,domain_ref="routing_domain")
+!     CALL xios_add_child(field_def_hdl,field_hdl,"state_rp1")
+!     CALL xios_set_attr(field_hdl, field_ref="state_r", domain_ref="routing_domain_expand", read_access=.TRUE.)
+!   
+!     CALL xios_get_handle("file_definition",file_def_hdl)
+!     CALL xios_add_child(file_def_hdl,file_hdl,"routing_correction")
+!     CALL xios_set_attr(file_hdl,type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.)
+!     CALL xios_add_child(file_hdl, field_hdl)
+!     CALL xios_set_attr(field_hdl, field_ref="trip_r", operation="instant")
+!     CALL xios_add_child(file_hdl, field_hdl)
+!     CALL xios_set_attr(field_hdl, field_ref="state_r", operation="instant")
+!     CALL xios_add_child(file_hdl, field_hdl)
+!     CALL xios_set_attr(field_hdl, field_ref="trip_r_init", operation="once")
+!     CALL xios_add_child(file_hdl, field_hdl)
+!     CALL xios_set_attr(field_hdl, field_ref="state_r_init", operation="once")
+
+
+      CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH)
+      CALL xios_orchidee_change_context("orchidee_init_routing")
+      CALL xios_set_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo)    ! get routing domain dimension
+      CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d)                           ! get routing domain dimension
+      CALL xios_close_context_definition()
+      CALL xios_update_calendar(1)
+
+      CALL xios_send_field("trip_r_init",trip_r)
+      
+      CALL xios_send_field("trip_r",trip_r)
+      CALL xios_recv_field("trip_rp1",trip_rp1)
+     
+      CALL xios_send_field("state_r",state_r)
+      CALL xios_recv_field("state_rp1",state_rp1)
+      
+      DO j=1,nj
+        DO i=1,ni
+          IF (trip_rp1(i,j)<=8) THEN
+            CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj)
+            IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow
+          ENDIF
+        ENDDO
+      ENDDO          
+
+        
+      state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
+      trip_r(1:ni,1:nj)  = trip_rp1(1:ni,1:nj)
+     
+      updated=1
+      it=2
+      DO WHILE(updated>0 .AND. it<MAX(ni_glo,nj_glo))
+        updated=0
+
+        CALL xios_update_calendar(it)
+        
+        CALL xios_send_field("trip_r", trip_r)
+        CALL xios_recv_field("trip_rp1",trip_rp1)
+     
+        CALL xios_send_field("state_r",state_r)
+        CALL xios_recv_field("state_rp1",state_rp1)
+
+        state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj)
+        trip_r(1:ni,1:nj)  = trip_rp1(1:ni,1:nj)
+
+      
+         
+        DO j=0,nj+1
+          DO i=0,ni+1
+         
+            IF (trip_rp1(i,j)==97 .OR. trip_rp1(i,j)==98 .OR. trip_rp1(i,j)==99) THEN
+            
+              IF (state_rp1(i,j)==is_coast) THEN ! ok
+                
+                IF (trip_rp1(i,j)==97) THEN
+                  IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) trip_r(i,j)=98  ! => lakinflow will becomes coastalflow
+                ELSE
+                 ! ok for riverflow
+                ENDIF
+                 
+              ELSE IF (state_rp1(i,j) == is_ter) THEN ! too short
+              
+                IF (trip_rp1(i,j)==97) THEN
+                  !ok for lakinflow
+                ELSE  ! => for riverflow/coastalflow
+                  IF (trip_extended_rp1(i,j)==97) THEN ! definitively a lake that cannot be routed further away
+                    IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
+                      trip_r(i,j)=97
+                      updated=updated+1
+                    ENDIF
+                  ELSE
+                    CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj)
+                    IF (nexti/=-1 .AND. nextj/=-1) THEN
+                      IF (nexti==i .AND. nextj==j) THEN 
+                        !! routed on itself in midle of land => lake
+                         IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) trip_r(i,j)=97
+                      ELSE
+                    
+                        IF (state_rp1(nexti,nextj)==is_coast .OR. state_rp1(nexti,nextj)==is_ter) THEN
+                
+                          IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
+                            trip_r(i,j)=trip_extended_rp1(i,j)
+                            state_r(i,j)= is_ter
+                            topoind_r(i,j) = 1e-10   ! flow instantanesly 
+                            updated=updated+1
+                          ENDIF
+ 
+                          IF (nexti>=1 .AND. nexti<=ni .AND. nextj>=1 .AND. nextj<=nj) THEN
+                            IF (trip_r(nexti,nextj)/=trip_rp1(i,j)) THEN
+                              trip_r(nexti,nextj)=trip_rp1(i,j)
+                              updated=updated+1
+                            ENDIF
+                          ENDIF
+    
+                        ELSE
+                          IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
+                            STOP 'not possible'
+                          ENDIF
+                        ENDIF      
+                      ENDIF
+                    ENDIF
+                  ENDIF
+                ENDIF 
+              
+              ELSE IF (state_rp1(i,j) == is_oce) THEN ! too long
+                 
+                DO k=-1,1
+                  DO m=-1,1
+                    IF (k==0 .AND. m==0) CYCLE
+                    prevj = j+k ; previ = i+m
+                    CALL next_trip(trip_rp1, ni, nj, previ, prevj, nexti, nextj)
+!                    IF (nexti==i .AND. nextj==j .AND. previ>=1 .AND. previ<=ni .AND. prevj>=1 .AND. prevj<=nj) trip_r(previ,prevj)=trip_rp1(i,j)
+                    IF (nexti==i .AND. nextj==j .AND. previ>=1 .AND. previ<=ni .AND. prevj>=1 .AND. prevj<=nj) trip_r(previ,prevj)=trip_rp1(i,j)
+                 ENDDO
+                ENDDO
+                
+                IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj) THEN
+                  trip_r(i,j)=1000
+                  updated=updated+1
+                ENDIF
+              ENDIF
+
+            ENDIF
+          ENDDO
+        ENDDO                      
+      
+        CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
+        it = it +1  
+      ENDDO 
+     
+      CALL xios_send_field("trip_r_final",trip_r)
+      
+      CALL xios_context_finalize
+      CALL xios_orchidee_change_context("orchidee")
+    
+  CONTAINS
+      
+    SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj)
+    IMPLICIT NONE
+      REAL, INTENT(IN)     :: trip(0:ni+1,0:nj+1)
+      INTEGER, INTENT(IN)  :: ni
+      INTEGER, INTENT(IN)  :: nj
+      INTEGER, INTENT(IN)  :: i
+      INTEGER, INTENT(IN)  :: j
+      INTEGER, INTENT(OUT) :: nexti
+      INTEGER, INTENT(OUT) :: nextj
+   
+      IF (i<0 .OR. i>ni+1 .OR. j<0 .OR. j>nj+1) THEN
+        nexti=-1 
+        nextj=-1
+        RETURN
+      ENDIF 
+   
+      SELECT CASE(NINT(trip(i,j)))
+        CASE (1)
+           nextj=j-1 ; nexti=i    ! north
+        CASE (2)
+           nextj=j-1 ; nexti=i+1  ! north-east
+        CASE (3)
+           nextj=j ;   nexti=i+1  ! east
+        CASE (4)
+           nextj=j+1 ;  nexti=i+1 ! south-east
+        CASE (5)
+           nextj=j+1 ;  nexti=i   ! south
+        CASE(6)
+           nextj=j+1 ;  nexti=i-1 ! south-west
+        CASE (7)
+           nextj=j ;   nexti=i-1  ! west
+        CASE (8)
+           nextj=j-1 ;  nexti=i-1 ! north-west
+        CASE DEFAULT
+           nextj=-1 ; nexti=-1 ! undefined
+       END SELECT
+   
+       IF (nexti<0 .OR. nexti>ni+1) nexti=-1
+       IF (nextj<0 .OR. nextj>nj+1) nextj=-1
+
+     END SUBROUTINE next_trip
+     
+  END SUBROUTINE routing_flow_correct_riverflow
+      
+
+
+
+
+
+  
+  SUBROUTINE routing_flow_init_mean(kjit, rest_id)
+  USE ioipsl   
+  IMPLICIT NONE
+  INTEGER(i_std), INTENT(in)    :: kjit
+  INTEGER(i_std), INTENT(in)    :: rest_id        !! Restart file identifier (unitless)
+ 
+  INTEGER :: ier    
+  CHARACTER(LEN=80)   :: var_name       !! To store variables names for I/O (unitless)
+    
+
+    ! Get from the restart the fluxes we accumulated.
+    !
+  
+    ALLOCATE (runoff_mean(nbpt), stat=ier)
+    runoff_mean(:) = 0
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','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)
+    drainage_mean(:) = 0
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_init_mean','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)
+    
+  END SUBROUTINE routing_flow_init_mean
+
+  SUBROUTINE routing_flow_make_mean(runoff, drainage)
+  IMPLICIT NONE
+    REAL(r_std),INTENT(IN) :: runoff(:)
+    REAL(r_std),INTENT(IN) :: drainage(:)
+    
+    runoff_mean(:) = runoff_mean(:) + runoff
+    drainage_mean(:) =  drainage_mean(:) + drainage
+    
+  END SUBROUTINE routing_flow_make_mean
+  
+
+  SUBROUTINE routing_flow_reset_mean
+  IMPLICIT NONE
+    runoff_mean(:) = 0
+    drainage_mean(:) =  0
+    
+  END SUBROUTINE routing_flow_reset_mean
+
+
+  SUBROUTINE routing_flow_finalize_mean(kjit, rest_id)
+  USE ioipsl
+  IMPLICIT NONE
+  INTEGER,INTENT(IN) :: kjit
+  INTEGER,INTENT(IN) :: rest_id
+ 
+    CALL restput_p (rest_id, 'runoff_route', nbp_glo, 1, 1, kjit, runoff_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'drainage_route', nbp_glo, 1, 1, kjit, drainage_mean, 'scatter',  nbp_glo, index_g)
+     
+    DEALLOCATE(runoff_mean)
+    DEALLOCATE(drainage_mean)
+    
+  END SUBROUTINE routing_flow_finalize_mean
+
+
+  SUBROUTINE routing_flow_init_local(contfrac, nbpt_r)
+    USE xios
+    USE grid, ONLY : area
+    IMPLICIT NONE
+    INCLUDE "mpif.h"
+
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+    REAL(r_std),INTENT(IN)                 :: contfrac(nbpt)   !! fraction of land
+    INTEGER,INTENT(OUT)                    :: nbpt_r             !! nb points routing native grid
+
+    !! 0.2 Local variables
+    INTEGER :: ni                                         !! longitude dimension of local routing grid
+    INTEGER :: nj                                         !! latitude dimension of local routing grid
+    REAL(r_std)             :: contfrac_mpi(nbp_mpi)
+    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
+    REAL(r_std),ALLOCATABLE :: diag_r(:)    !! fraction of routing cell intesected by coastal cells of native grid
+    REAL(r_std),ALLOCATABLE :: frac_lake_r(:)          !! fraction of routing cell intesected by cells of native grid
+    REAL(r_std),ALLOCATABLE :: frac_coast_r(:)    
+  
+    REAL(r_std)             :: diag(nbp_mpi)    !! fraction of routing cell intesected by coastal cells of native grid
+    
+    LOGICAL :: coastline(nbpt)
+    LOGICAL :: coastline_mpi(nbp_mpi)
+    REAL(r_std) :: area_mpi(nbp_mpi)
+
+    INTEGER :: ij, ij_r, ij_rp1, i,j,jp1,jm1,ip1,im1,jr,ir
+    INTEGER :: basins_count_mpi
+    INTEGER :: nb_coast_points
+    INTEGER :: ierr
+    LOGICAL :: file_exists
+    REAL(r_std) :: epsilon = 1e-5
+
+!_ ================================================================================================================================
+    !
+    !> 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]
+
+    split_routing=1
+    CALL getin_p("SPLIT_ROUTING",split_routing)
+
+    CALL compute_coastline(contfrac, coastline)
+    
+    CALL gather_omp(contfrac, contfrac_mpi)
+    CALL gather_omp(area, area_mpi)
+    CALL gather_omp(coastline, coastline_mpi)
+    
+    ALLOCATE(is_coastline(nbp_mpi))
+    is_coastline=coastline_mpi
+
+    IF (is_omp_root) THEN
+
+       WHERE(contfrac_mpi<=epsilon) contfrac_mpi=0
+       WHERE(contfrac_mpi>=epsilon) contfrac_mpi=1
+       
+       diag=0
+       WHERE(coastline_mpi) diag=1
+       CALL xios_send_field("is_coastline", diag)
+       
+       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 module variable     
+       ALLOCATE(fast_reservoir_r(ni*nj))       
+       ALLOCATE(slow_reservoir_r(ni*nj))       
+       ALLOCATE(stream_reservoir_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(topoind_r(ni*nj))       
+       ALLOCATE(routing_mask_r(ni*nj)) 
+       ALLOCATE(basins_extended_r(nbpt_r))
+       ALLOCATE(route_flow_rp1((ni+2)*(nj+2)))
+
+       ALLOCATE(routing_weight(nbp_mpi))
+       ALLOCATE(routing_weight_in(nbp_mpi))
+       ALLOCATE(unrouted_weight(nbp_mpi))
+      
+       ALLOCATE(diag_r(nbpt_r)) ! for diags
+       ALLOCATE(frac_coast_r(nbpt_r))
+       ALLOCATE(frac_lake_r(nbpt_r))
+       ALLOCATE(weight_coast_to_coast_r(nbpt_r))
+       ALLOCATE(weight_coast_to_lake_r(nbpt_r))
+       ALLOCATE(weight_lake_to_coast_r(nbpt_r))
+       ALLOCATE(weight_lake_to_lake_r(nbpt_r))
+
+       ALLOCATE(trip_extended_r(ni,nj))
+
+       is_lakeinflow_r(:) = .FALSE.
+       is_coastalflow_r(:) = .FALSE.
+       is_riverflow_r(:) = .FALSE.
+   
+       
+       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
+       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
+
+       CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r) 
+   
+       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
+       
+
+
+       !! Compute the routing
+       !! Loop on all point point of the local routing grid + halo
+  
+       route_flow_rp1(:)=-1
+
+       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                    
+                
+                IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj)   routing_mask_r(ij_r)=.TRUE.
+                
+                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 (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
+                      IF (trip_rp1(ir,jr)>99) STOP 'Pb point not routed to outflow'
+                   ENDIF
+                ENDIF
+             ELSE
+               IF (i>=1 .AND. i<=ni .AND. j>=1 .AND. j<=nj)  routing_mask_r(ij_r)=.FALSE.             
+             ENDIF
+          ENDDO
+       ENDDO
+       is_streamflow_r(:) = .NOT. (is_lakeinflow_r(:) .OR. is_coastalflow_r(:) .OR. is_riverflow_r(:)) .AND. routing_mask_r(:)
+       
+       diag_r(:)=0
+       WHERE(routing_mask_r) diag_r=1
+       CALL xios_send_field("routing_mask_r", diag_r)
+       CALL xios_recv_field("frac_routing", diag)
+       routing_weight_in (:) = 0
+       unrouted_weight (:) = 0
+       WHERE (diag > 0) routing_weight_in = contfrac_mpi*area_mpi/diag
+       WHERE (diag == 0 ) unrouted_weight = contfrac_mpi*area_mpi
+       
+
+      ! compute the number of coast cells to distribute lakeinflow onto
+       diag=0
+       WHERE (coastline_mpi) diag=1
+       nb_coast_points=SUM(diag)
+       CALL reduce_sum_mpi(nb_coast_points, total_coast_points)
+       CALL bcast_mpi(total_coast_points)
+
+
+       DO ij=1,nbp_mpi
+          routing_weight(ij)=contfrac_mpi(ij)*area_mpi(ij)  
+       ENDDO
+
+       ! 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)
+
+       diag(:) = 0
+       WHERE (coastline_mpi) diag=1
+       CALL xios_send_field("mask_coastal",diag)
+       CALL xios_recv_field("frac_coastal_r",diag_r)
+       WHERE (diag_r > 100.) ! missing_value
+         diag_r=0
+       ELSE WHERE (diag_r < epsilon) 
+         diag_r=0
+       ELSEWHERE (diag_r > 1-epsilon) 
+         diag_r=1
+       END WHERE
+       frac_coast_r=diag_r ;
+
+       diag(:) = 0
+       WHERE (.NOT. coastline_mpi) diag=1
+       CALL xios_send_field("mask_lake",diag)
+       CALL xios_recv_field("frac_lake_r",diag_r)
+       WHERE (diag_r > 100.) ! missing_value
+         diag_r=0.
+       ELSEWHERE (diag_r < epsilon) 
+         diag_r=0.
+       ELSEWHERE (diag_r > 1-epsilon) 
+         diag_r=1.
+       END WHERE
+       frac_lake_r = diag_r
+
+       weight_coast_to_coast_r(:)=0
+       weight_coast_to_lake_r(:)=0
+       weight_lake_to_coast_r(:)=0
+       weight_lake_to_lake_r(:)=0
+
+       DO ij=1,nbpt_r
+         IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) ) THEN
+           IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN
+             PRINT*,"riverflow or costalflow point on routing grid can not be interpolated on orchidee grid, this migh not happen"
+             STOP
+           ELSE IF (frac_coast_r(ij)==0 .AND. frac_lake_r(ij)>0) THEN
+             weight_coast_to_lake_r(ij) = 1./frac_lake_r(ij)
+             weight_coast_to_coast_r(ij) = 0
+           ELSE IF (frac_coast_r(ij)>0 ) THEN
+             weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij))
+             weight_coast_to_lake_r(ij) = 0
+           ENDIF
+         ELSE IF (is_lakeinflow_r(ij)) THEN
+           IF ( frac_coast_r(ij)==0 .AND. frac_lake_r(ij)==0) THEN
+             PRINT*,"lakeinflow on routing grid can not be interpolated on orchidee grid, this might not happen"
+             STOP
+           ELSE IF (frac_lake_r(ij)==0 .AND. frac_coast_r(ij)>0) THEN
+             weight_lake_to_coast_r(ij) = 1./frac_coast_r(ij)
+             weight_lake_to_lake_r(ij) = 0
+           ELSE IF (frac_lake_r(ij)>0 ) THEN
+             weight_lake_to_lake_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij))
+             weight_lake_to_coast_r(ij) = 0
+           ENDIF
+         ENDIF
+       ENDDO
+       
+       CALL compute_basins_area(nbpt_r,nbpt_rp1)
+
+    ELSE
+
+       nbpt_r=0
+
+    ENDIF !is_omp_root
+    
+    
+  END SUBROUTINE routing_flow_init_local     
+  
+  
+  SUBROUTINE compute_basins_area(nbpt_r,nbpt_rp1 )
+  USE xios
+  USE routing_native_para
+  IMPLICIT NONE
+    INCLUDE 'mpif.h'
+    INTEGER, INTENT(IN) :: nbpt_r
+    INTEGER, INTENT(IN) :: nbpt_rp1
+    REAL(r_std),PARAMETER :: PI=atan(1.)*4
+    REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter
+
+    INTEGER ::  ni, nj
+    REAL(r_std),ALLOCATABLE :: lon(:), lat(:)
+    REAL(r_std)  :: area_rp1(nbpt_rp1)
+    REAL(r_std)  :: send_area_rp1(nbpt_rp1)
+    REAL(r_std)  :: recv_area_r(nbpt_rp1)
+    REAL :: delta_lon, delta_lat
+    INTEGER :: i,j,ij,it,ig,ierr
+    REAL(r_std) :: sum_area
+   
+    ALLOCATE(basins_area_r(nbpt_r))
+
+    CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)
+    ALLOCATE(lon(ni),lat(nj))
+    CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat)
+    
+   ! compute routing cell area
+    delta_lon=ABS(lon(2)-lon(1))*Pi/180.
+    delta_lat=ABS(lat(2)-lat(1))*Pi/180.
+
+    DO j=1,nj 
+      DO i=1,ni
+        ig=(j+1-1)*(ni+2)+i+1 
+        area_rp1(ig) = delta_lon*delta_lat*cos(lat(j)*Pi/180.)*earth_radius*earth_radius  ! pretty good approximation
+        ij = (j-1)*ni+i
+        basins_area_r(ij) = delta_lon*delta_lat*cos(lat(j)*Pi/180.)*earth_radius*earth_radius  ! pretty good approximation
+      ENDDO
+    ENDDO
+
+    CALL update_halo(area_rp1)
+
+      
+    send_area_rp1(:) = area_rp1(:)
+    sum_area=1.
+    DO WHILE (sum_area/=0)
+
+      recv_area_r(:) = 0
+      DO ig=1,nbpt_rp1
+        IF ( route_flow_rp1(ig) > 0 ) THEN
+          basins_area_r(route_flow_rp1(ig)) = basins_area_r(route_flow_rp1(ig)) + send_area_rp1(ig) 
+          recv_area_r(route_flow_rp1(ig)) = recv_area_r(route_flow_rp1(ig)) + send_area_rp1(ig)
+        ENDIF
+      ENDDO
+
+      DO j=1,nj 
+        DO i=1,ni
+          ij=((j-1)*ni)+i
+          ig=(j+1-1)*(ni+2)+i+1 
+          send_area_rp1(ig) = recv_area_r(ij)
+          IF (is_riverflow_r(ij) .OR. is_coastalflow_r(ij) .OR. is_lakeinflow_r(ij)) send_area_rp1(ig)=0
+        ENDDO
+      ENDDO
+      
+      CALL update_halo(send_area_rp1)
+      sum_area = sum(send_area_rp1)
+      CALL MPI_ALLREDUCE(MPI_IN_PLACE, sum_area, 1, MPI_REAL_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
+
+    ENDDO
+
+  END SUBROUTINE compute_basins_area
+
+  
+  SUBROUTINE initialize_stations(dt_routing)
+  USE xios
+  IMPLICIT NONE
+    INCLUDE 'mpif.h'
+    REAL(r_std), INTENT (in)                     :: dt_routing                !! Routing time step (s)
+
+    REAL(r_std),PARAMETER :: PI=atan(1.)*4
+    REAL(r_std),PARAMETER :: earth_radius=6.371009E6 ! in meter
+    INTEGER ::  ni, nj
+    REAL(r_std),ALLOCATABLE :: lon(:), lat(:)
+
+    INTEGER :: nb_mouth
+    CHARACTER(LEN=60),ALLOCATABLE :: mouth(:)  
+    INTEGER,ALLOCATABLE :: mouth_id(:)  
+    INTEGER,ALLOCATABLE :: mouth_index(:)  
+  
+    REAL,ALLOCATABLE :: station_lonlat(:,:)  
+    INTEGER,ALLOCATABLE :: station_prec(:)  
+
+    INTEGER ::i,j,r,k
+  
+    TYPE(xios_duration) :: ts 
+    TYPE(xios_scalar) :: scalar_hdl
+    TYPE(xios_scalargroup) :: scalar_def_hdl
+    TYPE(xios_field) :: field_hdl
+    TYPE(xios_fieldgroup) :: field_def_hdl, fieldgroup_hdl
+    TYPE(xios_file) :: file_hdl
+    TYPE(xios_filegroup) :: file_def_hdl
+    TYPE(xios_date) :: start_date, time_origin
+    CHARACTER(LEN=20)  :: calendar_type
+
+    CHARACTER(LEN=256) :: str_station_ind
+    REAL :: lon_a,lat_a,lon_b,lat_b,dist,max_area, max_max_area, min_dist
+    INTEGER :: min_dist_index
+    INTEGER :: ierr
+
+    
+    CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)
+    ALLOCATE(lon(ni),lat(nj))
+    CALL xios_get_domain_attr("routing_domain",lonvalue_1d=lon, latvalue_1d=lat)
+
+    nb_station=0
+    CALL getin("nb_station",nb_station)
+    ALLOCATE(station(nb_station)) 
+    ALLOCATE(station_lonlat(nb_station,2)) 
+    ALLOCATE(station_index(nb_station)) 
+    ALLOCATE(station_prec(nb_station)) 
+    
+    DO k=1,nb_station
+      WRITE(str_station_ind,*) k
+      str_station_ind=ADJUSTL(str_station_ind)
+      CALL getin("station"//TRIM(str_station_ind)//"_id",station(k))
+      CALL getin("station"//TRIM(str_station_ind)//"_coor",station_lonlat(k,:))
+      station_prec(k)=50000
+      CALL getin("station"//TRIM(str_station_ind)//"_prec",station_prec(k))
+      
+      lon_a = station_lonlat(k,1)*Pi/180.
+      lat_a = station_lonlat(k,2)*Pi/180.
+      
+      max_area=0
+      min_dist=HUGE(min_dist)
+      DO j=1,nj 
+        DO i=1,ni
+          r=((j-1)*ni)+i
+          IF (routing_mask_r(r)) THEN
+            lon_b=lon(i)*Pi/180.
+            lat_b=lat(j)*Pi/180.
+            dist=earth_radius*acos(sin(lon_a)*sin(lon_b)+cos(lon_a)*cos(lon_b)*cos(lat_b-lat_a))    
+            IF (dist<station_prec(k)) THEN
+              IF (dist<min_dist) THEN
+                min_dist=dist
+                min_dist_index = r
+              ENDIF
+          
+              IF (basins_area_r(r) > max_area) THEN
+                max_area = basins_area_r(r)
+                station_index(k) = r
+              ENDIF
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO      
+
+      CALL MPI_ALLREDUCE(max_area, max_max_area, 1, MPI_REAL_ORCH, MPI_MAX, MPI_COMM_ORCH, ierr)
+      IF (max_area /= max_max_area) station_index(k)=-1
+      IF (max_area == 0) station_index(k)=-1
+      
+!      PRINT*,"Station ",TRIM(station(k))," old coor :",lon(min_dist_index),lat(min_dist_index)," basin_area (km2)", basins_area_r(min_dist_index)/1000/1000 
+!      PRINT*,"Station ",TRIM(station(k)),"new coor ",lon(station_index(k)),lat(station_index(k))," basin_area(km2) ", basins_area_r(station_index(k))/1000/1000
+    ENDDO
+    
+    CALL xios_get_calendar_type(calendar_type)
+    CALL xios_get_start_date(start_date)
+    CALL xios_get_time_origin(time_origin)
+
+    CALL xios_context_initialize("orchidee_routing_out",MPI_COMM_ORCH)
+    CALL xios_orchidee_change_context("orchidee_routing_out")
+    ts.second=dt_routing
+    calendar_type="gregorian" 
+    CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts )
+
+    CALL xios_get_handle("scalar_definition",scalar_def_hdl)
+
+!    CALL xios_get_handle("mouths",field_def_hdl)
+!    DO k=1,nb_mouth
+!      CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(mouth(k))//"_scalar")
+!      CALL xios_set_attr(scalar_hdl,label=TRIM(mouth(k)))
+!      CALL xios_add_child(field_def_hdl,field_hdl, TRIM(mouth(k)))
+!      CALL xios_set_attr(field_hdl,scalar_ref=TRIM(mouth(k))//"_scalar")
+!      CALL xios_set_attr(field_hdl, freq_op = freq_op*xios_second)
+!    ENDDO
+
+    CALL xios_get_handle("field_definition",field_def_hdl)
+    CALL xios_add_child(field_def_hdl, fieldgroup_hdl, "stations")
+    DO k=1,nb_station
+      CALL xios_add_child(scalar_def_hdl,scalar_hdl,TRIM(station(k))//"_scalar")
+      CALL xios_set_attr(scalar_hdl,label=TRIM(station(k)))
+      CALL xios_add_child(fieldgroup_hdl,field_hdl, TRIM(station(k)))
+      CALL xios_set_attr(field_hdl,scalar_ref=TRIM(station(k))//"_scalar")
+    ENDDO
+     
+    CALL xios_get_handle("file_definition",file_def_hdl)
+    CALL xios_add_child(file_def_hdl, file_hdl, "stations")
+    CALL xios_set_attr(file_hdl, type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.)
+    CALL xios_add_child(file_hdl, fieldgroup_hdl)
+    CALL xios_set_attr(fieldgroup_hdl, group_ref="stations", operation="average")
+    CALL xios_close_context_definition()
+    CALL xios_orchidee_change_context("orchidee")
+  END SUBROUTINE initialize_stations
+
+  SUBROUTINE routing_flow_main(dt_routing) !
+    
+    USE xios
+    USE grid, ONLY : area 
+    USE routing_native_para
+    IMPLICIT NONE
+    INCLUDE "mpif.h"
+    
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+    REAL(r_std), INTENT (in)                     :: dt_routing                !! Routing time step (s)
+
+    !! 0.4 Local variables
+    REAL(r_std)                                  :: runoff(nbp_mpi)              !! Grid-point runoff (kg/dt)
+    REAL(r_std)                                  :: runoff_in(nbp_mpi)              !! Grid-point runoff (kg/dt)
+    REAL(r_std)                                  :: drainage(nbp_mpi)            !! Grid-point drainage (kg/dt)
+    REAL(r_std)                                  :: drainage_in(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)         
+    REAL(r_std)                                  :: area_mpi(nbp_mpi) ! cell area         
+    REAL(r_std)                                  :: flow_coast(nbp_mpi)          
+    REAL(r_std)                                  :: flow_lake(nbp_mpi)          
+
+    ! from input model -> routing_grid
+    REAL(r_std)                      :: runoff_r(nbpt_r)              !! Grid-point runoff (kg/m^2/dt)
+    REAL(r_std)                      :: drainage_r(nbpt_r)            !! Grid-point drainage (kg/m^2/dt)
+
+    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)        :: transport_r                 !! Water transport between basins (kg/dt)
+
+    INTEGER(i_std)                               :: ig        !! Indices (unitless)
+    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)        :: 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)                           :: flow                      !! Outflow computation for the reservoirs (kg/dt)
+    REAL(r_std) :: basins_riverflow_mpi(0:basins_count)
+    REAL(r_std) :: basins_riverflow(0:basins_count)
+    REAL(r_std) :: water_balance_before, water_balance_after
+    REAL(r_std) :: sum_water_before, sum_water_after
+    REAL(r_std)  :: value
+    INTEGER(i_std) :: k
+    INTEGER :: ierr
+
+    !_ ================================================================================================================================
+    !> 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_mean,runoff)
+    CALL gather_omp(drainage_mean, drainage)
+    CALL gather_omp(area, area_mpi)
+
+
+    IF (is_omp_root) THEN
+      CALL xios_send_field("routing_basins_area",basins_area_r)
+
+      hydrographs_r(:)=0
+      ! water balance before
+
+      sum_water_before = sum((runoff(:)+drainage(:))*routing_weight(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
+      CALL MPI_ALLREDUCE(sum_water_before, water_balance_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+
+
+      runoff_in=runoff*routing_weight_in
+      CALL xios_send_field("routing_runoff",runoff_in)   ! interp conservative model -> routing
+      CALL xios_recv_field("routing_runoff_r",runoff_r)
+      WHERE(.NOT. routing_mask_r) runoff_r = 0
+
+      drainage_in=drainage*routing_weight_in
+      CALL xios_send_field("routing_drainage",drainage_in)   ! interp conservative model -> routing
+      CALL xios_recv_field("routing_drainage_r",drainage_r)
+      WHERE(.NOT. routing_mask_r) drainage_r = 0
+
+      CALL MPI_ALLREDUCE(sum(runoff*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(runoff_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      IF (sum_water_after/=0) THEN
+        IF (is_mpi_root) PRINT *,"runoff fixer : ", sum_water_before/sum_water_after
+        DO ig=1,nbpt_r
+          runoff_r(ig) =  runoff_r(ig) * sum_water_before/sum_water_after
+        ENDDO
+      ENDIF
+
+      CALL MPI_ALLREDUCE(sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      IF (sum_water_after/=0) THEN
+        IF (is_mpi_root) PRINT *,"drainage fixer : ", sum_water_before/sum_water_after
+        DO ig=1,nbpt_r
+          drainage_r(ig) =  drainage_r(ig) * sum_water_before/sum_water_after
+        ENDDO
+      ENDIF
+          
+
+      CALL MPI_ALLREDUCE(sum(runoff*routing_weight)+sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(runoff_r) +sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+
+      IF (is_mpi_root) PRINT *,"Loose water by interpolation ;  before : ", sum_water_before," ; after : ",sum_water_after,  &
+                                 " ; delta : ", 100.*(sum_water_after-sum_water_before)/(0.5*(sum_water_after+sum_water_before)),"%"
+
+
+      runoff_in(:) = runoff(:)*unrouted_weight(:)
+      drainage_in(:) = drainage(:)*unrouted_weight(:)
+      runoff_r(:) = runoff_r(:) / split_routing
+      drainage_r(:) = drainage_r(:) / split_routing
+      hydrographs_r(:) = 0
+      transport_r(:)=0
+
+      DO isplit=1,split_routing
+
+        DO ig=1,nbpt_r
+          IF ( routing_mask_r(ig) ) THEN
+
+            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*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
+
+
+        !-
+        !- Compute the transport
+        !-
+
+        DO ig=1,nbpt_r
+          flow_rp1(r_to_rp1(ig))=fast_flow_r(ig) + slow_flow_r(ig) + stream_flow_r(ig)
+        ENDDO
+
+        CALL update_halo(flow_rp1) ! transfert halo
+
+        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 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)
+            slow_reservoir_r(ig) = slow_reservoir_r(ig) + drainage_r(ig) -  slow_flow_r(ig)
+            stream_reservoir_r(ig) = stream_reservoir_r(ig) - stream_flow_r(ig)
+            IF (is_streamflow_r(ig)) THEN
+              stream_reservoir_r(ig)= stream_reservoir_r(ig) +  transport_r(ig)
+              transport_r(ig) = 0
+            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, transport, stream_flow,: ", &
+                                   stream_reservoir_r(ig), transport_r(ig), stream_flow_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 : ", fast_reservoir_r(ig), &
+                                 &runoff_r(ig), fast_flow_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)
+              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
+
+      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
+
+      ! now send riverflow, coastalflow and lakeinflow on orchidee grid
+
+      coastalflow(:)=0.
+      riverflow(:)=0.
+      lakeinflow(:)=0.
+
+      CALL xios_send_field("routing_coastalflow_to_coast_r" ,coastalflow_r*weight_coast_to_coast_r)
+      CALL xios_recv_field("routing_coastalflow_to_coast" ,flow_coast)
+      WHERE(.NOT.is_coastline) flow_coast=0
+
+      CALL xios_send_field("routing_coastalflow_to_lake_r" ,coastalflow_r*weight_coast_to_lake_r)
+      CALL xios_recv_field("routing_coastalflow_to_lake" ,flow_lake)
+      WHERE(is_coastline) flow_lake=0.
+
+      CALL MPI_ALLREDUCE(sum(coastalflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      IF (sum_water_after/=0) THEN
+        IF (is_mpi_root) PRINT *,"coastalflow fixer : ", sum_water_before/sum_water_after
+        flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
+        flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
+      ENDIF
+       
+      coastalflow=coastalflow+flow_coast
+      lakeinflow=lakeinflow+flow_lake
+
+      CALL xios_send_field("routing_riverflow_to_coast_r" ,riverflow_r*weight_coast_to_coast_r)
+      CALL xios_recv_field("routing_riverflow_to_coast" ,flow_coast)
+      WHERE(.NOT.is_coastline) flow_coast=0
+
+      CALL xios_send_field("routing_riverflow_to_lake_r" ,riverflow_r*weight_coast_to_lake_r)
+      CALL xios_recv_field("routing_riverflow_to_lake" ,flow_lake)
+      WHERE(is_coastline) flow_lake=0.
+      
+      CALL MPI_ALLREDUCE(sum(riverflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      IF (sum_water_after/=0) THEN
+        IF (is_mpi_root) PRINT *,"riverflow fixer : ", sum_water_before/sum_water_after
+        flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
+        flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
+      ENDIF   
+       
+      riverflow=riverflow+flow_coast
+      lakeinflow=lakeinflow+flow_lake
+
+      CALL xios_send_field("routing_lakeinflow_to_coast_r" ,lakeinflow_r*weight_lake_to_coast_r)
+      CALL xios_recv_field("routing_lakeinflow_to_coast" ,flow_coast)
+      WHERE(.NOT.is_coastline) flow_coast=0
+
+      CALL xios_send_field("routing_lakeinflow_to_lake_r" ,lakeinflow_r*weight_lake_to_lake_r)
+      CALL xios_recv_field("routing_lakeinflow_to_lake" ,flow_lake)
+      WHERE(is_coastline) flow_lake=0.
+
+      CALL MPI_ALLREDUCE(sum(lakeinflow_r),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(flow_coast+flow_lake),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      IF (sum_water_after/=0) THEN
+        IF (is_mpi_root) PRINT *,"lakeinflow fixer : ", sum_water_before/sum_water_after
+        flow_coast(:)=flow_coast(:)*(sum_water_before/sum_water_after)
+        flow_lake(:)=flow_lake(:)*(sum_water_before/sum_water_after)
+      ENDIF   
+       
+      coastalflow=coastalflow+flow_coast
+      lakeinflow=lakeinflow+flow_lake
+
+
+      WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in
+      WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in
+
+      sum_water_after = sum(coastalflow(:) + riverflow(:) + lakeinflow(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
+      CALL MPI_ALLREDUCE(sum_water_after, water_balance_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+
+      IF (is_mpi_root) 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)),"%"
+  ! diag
+      CALL xios_send_field("routing_fast_reservoir_r"  , fast_reservoir_r)
+      CALL xios_send_field("routing_slow_reservoir_r"  , slow_reservoir_r)
+      CALL xios_send_field("routing_stream_reservoir_r"  , stream_reservoir_r)
+      CALL xios_send_field("routing_riverflow_r"  , riverflow_r)
+      CALL xios_send_field("routing_coastalflow_r"  , coastalflow_r)
+      CALL xios_send_field("routing_lakeinflow_r"  , lakeinflow_r)
+      CALL xios_send_field("out_flow",lakeinflow+coastalflow+riverflow)
+      CALL xios_send_field("routing_hydrographs_r", (hydrographs_r+lakeinflow_r+coastalflow_r+riverflow_r)/1000./dt_routing)
+      CALL xios_send_field("routing_riverflow"  , riverflow)
+      CALL xios_send_field("routing_coastalflow"  , coastalflow)
+      CALL xios_send_field("routing_lakeinflow"  , lakeinflow)
+!
+!     stations
+      CALL xios_orchidee_change_context("orchidee_routing_out")
+      station_ts = station_ts +1
+      CALL xios_update_calendar(station_ts)
+      DO k=1,nb_station
+        IF (station_index(k) /=-1) THEN
+          value = hydrographs_r(station_index(k))/1000./dt_routing
+        ELSE
+          value = 0
+        ENDIF
+        CALL MPI_ALLREDUCE(MPI_IN_PLACE,value,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+        CALL xios_send_field(TRIM(station(k)),value)
+      ENDDO
+      CALL xios_orchidee_change_context("orchidee")
+    ENDIF ! is_omp_root
+
+    CALL scatter_omp(riverflow,riverflow_mean)
+    CALL scatter_omp(coastalflow,coastalflow_mean)
+    CALL scatter_omp(lakeinflow,lakeinflow_mean)
+
+    CALL routing_flow_reset_mean
+    
+  END SUBROUTINE routing_flow_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_flow_finalize(kjit, rest_id)
+    USE xios
+    USE ioipsl
+    IMPLICIT NONE
+    INTEGER, INTENT(IN) :: kjit
+    INTEGER, INTENT(IN) :: rest_id
+    !_ ================================================================================================================================
+
+    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
+
+    CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter',  nbp_glo, index_g)
+    CALL routing_flow_finalize_mean(kjit, rest_id)
+
+  END SUBROUTINE routing_flow_finalize
+
+
+  !! ================================================================================================================================
+  !! 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_flow_clear
+  IMPLICIT NONE
+
+    IF (is_omp_root) THEN
+       IF (ALLOCATED(topoind_r)) DEALLOCATE(topoind_r)
+       IF (ALLOCATED(route_flow_rp1)) DEALLOCATE(route_flow_rp1)
+       IF (ALLOCATED(routing_mask_r)) DEALLOCATE(routing_mask_r)
+       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
+
+  END SUBROUTINE routing_flow_clear
+
+
+END MODULE routing_native_flow_mod
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_irrig.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_irrig.f90	(revision 8227)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_irrig.f90	(revision 8227)
@@ -0,0 +1,772 @@
+MODULE routing_native_irrig_mod
+  USE constantes
+
+  
+  PRIVATE
+  
+  INTEGER, SAVE :: nbpt
+  !$OMP THREADPRIVATE(nbpt)
+  
+  INTEGER, SAVE :: nbpt_r
+  !$OMP THREADPRIVATE(nbpt_r)
+
+  REAL(r_std), ALLOCATABLE, SAVE :: irrig_gw_source_r(:)    ! diag
+  REAL(r_std), ALLOCATABLE, SAVE :: irrig_fast_source_r(:)  ! diag
+  REAL(r_std), ALLOCATABLE, SAVE :: irrig_str_source_r(:)   ! diag
+ 
+  REAL(r_std), ALLOCATABLE, SAVE :: vegtot_mean(:)
+  !$OMP THREADPRIVATE(vegtot_mean)
+  
+  REAL(r_std), ALLOCATABLE, SAVE :: humrel_mean(:)
+  !$OMP THREADPRIVATE(humrel_mean)
+  
+  REAL(r_std), ALLOCATABLE, SAVE :: transpot_mean(:)
+  !$OMP THREADPRIVATE(transpot_mean)
+  
+  REAL(r_std), ALLOCATABLE, SAVE :: runoff_mean(:)
+  !$OMP THREADPRIVATE(runoff_mean)
+  
+  REAL(r_std), ALLOCATABLE, SAVE :: precip_mean(:)
+  !$OMP THREADPRIVATE(precip_mean)
+
+  REAL(r_std), ALLOCATABLE, SAVE :: irrigation_mean(:)
+  !$OMP THREADPRIVATE(irrigation_mean)
+ 
+  REAL(r_std), ALLOCATABLE, SAVE :: irrigated(:)
+  !$OMP THREADPRIVATE(irrigated)
+ 
+  PUBLIC irrigation_initialize, irrigation_main, irrigation_mean_make, irrigation_get, irrigation_finalize
+
+CONTAINS
+
+SUBROUTINE irrigation_get(irrigation_mean)
+  IMPLICIT NONE
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: irrigation_mean(:)
+
+    CALL irrigation_get_(irrigation_mean)
+  END SUBROUTINE irrigation_get
+
+  SUBROUTINE irrigation_get_(irrigation_mean_)
+  IMPLICIT NONE
+    REAL(r_std),OPTIONAL, INTENT(OUT) :: irrigation_mean_(:)
+    IF (PRESENT(irrigation_mean_))      irrigation_mean_ = irrigation_mean
+
+  END SUBROUTINE irrigation_get_
+
+
+  SUBROUTINE irrigation_initialize(kjit, rest_id, nbpt_, nbpt_r_, irrigated_next )
+  USE constantes
+  IMPLICIT NONE
+    INTEGER, INTENT(IN)         :: kjit
+    INTEGER, INTENT(IN)         :: rest_id
+    INTEGER, INTENT(IN)         :: nbpt_
+    INTEGER, INTENT(IN)         :: nbpt_r_
+    REAL(r_std), INTENT(IN)     :: irrigated_next(nbpt_)
+  
+    CALL irrigation_local_init(kjit, rest_id, nbpt_, nbpt_r_,irrigated_next)
+    CALL irrigation_mean_init(kjit, rest_id)
+    
+  END SUBROUTINE irrigation_initialize
+
+
+  SUBROUTINE irrigation_finalize(kjit, rest_id)
+  IMPLICIT NONE
+    INTEGER, INTENT(IN)     :: kjit
+    INTEGER, INTENT(IN)     :: rest_id
+    
+    CALL irrigation_mean_finalize(kjit, rest_id)
+    
+  END SUBROUTINE irrigation_finalize
+  
+  
+  SUBROUTINE irrigation_local_init(kjit, rest_id, nbpt_, nbpt_r_, irrigated_next)
+  IMPLICIT NONE
+    INTEGER, INTENT(IN)         :: kjit
+    INTEGER, INTENT(IN)         :: rest_id
+    INTEGER, INTENT(IN)         :: nbpt_
+    INTEGER, INTENT(IN)         :: nbpt_r_
+    REAL(r_std), INTENT(IN)     :: irrigated_next(nbpt_)
+
+     nbpt = nbpt_
+     nbpt_r = nbpt_r_
+     IF (do_irrigation) THEN
+       ALLOCATE(irrigated(nbpt))
+       irrigated(:)=irrigated_next(:)
+       ALLOCATE(irrig_gw_source_r(nbpt_r))
+       ALLOCATE(irrig_fast_source_r(nbpt_r))
+       ALLOCATE(irrig_str_source_r(nbpt_r))
+     ENDIF
+     
+  END SUBROUTINE irrigation_local_init
+
+
+  
+  SUBROUTINE irrigation_mean_init(kjit, rest_id)
+   USE ioipsl_para
+   USE grid
+   USE sechiba_io_p
+   IMPLICIT NONE
+    INTEGER, INTENT(IN)     :: kjit
+    INTEGER, INTENT(IN)     :: rest_id
+    CHARACTER(LEN=80)                            :: var_name       !! To store variables names for I/O (unitless)
+    INTEGER(i_std)                               :: ier
+      
+    ALLOCATE(irrigation_mean(nbpt), stat=ier)
+    irrigation_mean(:) = 0
+    IF (do_irrigation) THEN
+      IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigation_mean','','')
+      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)
+    
+      ALLOCATE(vegtot_mean(nbpt), stat=ier)
+      ALLOCATE(humrel_mean(nbpt), stat=ier)
+      ALLOCATE(transpot_mean(nbpt), stat=ier)
+      ALLOCATE (runoff_mean(nbpt), stat=ier)
+      ALLOCATE(precip_mean(nbpt), stat=ier)
+      CALL irrigation_mean_reset
+    ENDIF
+
+  END SUBROUTINE irrigation_mean_init
+
+  
+  SUBROUTINE irrigation_mean_reset
+  IMPLICIT NONE
+
+    vegtot_mean(:) = 0
+    humrel_mean(:) = 0
+    transpot_mean(:) = 0
+    runoff_mean(:) = 0
+    precip_mean(:) = 0
+
+  END SUBROUTINE irrigation_mean_reset
+  
+  SUBROUTINE irrigation_mean_make(dt_routing, veget_max, humrel, transpot, runoff ,precip )
+  USE pft_parameters
+  IMPLICIT NONE
+    REAL(r_std),INTENT(IN) :: dt_routing
+    REAL(r_std),INTENT(IN) :: veget_max(:,:)
+    REAL(r_std),INTENT(IN) :: humrel(:,:)
+    REAL(r_std),INTENT(IN) :: transpot(:,:)
+    REAL(r_std),INTENT(IN) :: runoff(:)
+    REAL(r_std),INTENT(IN) :: precip(:)
+
+    REAL(r_std), DIMENSION(nbpt)   :: tot_vegfrac_nowoody  !! Total fraction occupied by grass (0-1,unitless)
+    INTEGER :: jv, ig
+
+    IF (do_irrigation) THEN
+      runoff_mean(:) = runoff_mean(:) + runoff
+      precip_mean(:) = precip_mean(:) + precip
+
+      !! 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
+
+      ! New irrigation scheme uses mean of vegtot with jv 1 to nvm
+      ! Old scheme keeps jv 2 to nvm, even if possibly an error
+      IF ( .NOT. old_irrig_scheme ) THEN
+        DO jv=1,nvm
+           DO ig=1,nbpt
+              humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing
+              vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
+           ENDDO
+        ENDDO
+      ELSE
+        DO jv=2,nvm
+           DO ig=1,nbpt
+              humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing
+              vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing
+           ENDDO
+        ENDDO
+      ENDIF
+    ENDIF !do_irrigation
+  END SUBROUTINE irrigation_mean_make
+ 
+  SUBROUTINE irrigation_mean_finalize(kjit, rest_id)
+  USE grid
+  USE ioipsl_para
+  IMPLICIT NONE
+    INTEGER, INTENT(IN)     :: kjit
+    INTEGER, INTENT(IN)     :: rest_id
+
+    IF (do_irrigation) THEN
+      CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter',  nbp_glo, index_g)  
+    ENDIF
+
+  END SUBROUTINE irrigation_mean_finalize
+  
+
+  SUBROUTINE irrigation_main(dt_routing, reinfiltration, irrigated_next, irrig_frac, root_deficit, soiltile, &
+                             fraction_aeirrig_sw )
+  USE constantes
+  USE constantes_soil
+  USE grid, ONLY : area
+  USE xios
+  IMPLICIT NONE 
+    REAL(r_std), INTENT(IN)       :: dt_routing               
+    REAL(r_std), INTENT(IN)       :: reinfiltration(nbpt)               
+    REAL(r_std), INTENT(IN)       :: irrigated_next(nbpt)           
+    REAL(r_std), INTENT(IN)       :: irrig_frac(nbpt)     !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
+    REAL(r_std), INTENT(IN)       :: root_deficit(nbpt)   !! soil water deficit
+    REAL(r_std), INTENT(IN)       :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless)
+    REAL(r_std), INTENT(IN)       :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac
+
+    REAL(r_std)                   :: irrig_needs_r(nbpt_r)   !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
+    REAL(r_std)                   :: irrig_actual_r(nbpt_r)  !! Possible irrigation according to the water availability in the reservoirs (kg)
+    REAL(r_std)                   :: irrig_actual(nbpt)  !! Possible irrigation according to the water availability in the reservoirs (kg)
+
+    IF (do_irrigation) THEN
+      IF (irrig_map_dynamic_flag  ) irrigated=irrigated_next
+   
+      IF (old_irrig_scheme) THEN
+        CALL irrigation_compute_requested_old(reinfiltration, irrigated, irrig_needs_r)
+        CALL irrigation_old(irrig_needs_r, irrig_actual_r)
+      ELSE
+        CALL irrigation_compute_requested_new(dt_routing, root_deficit, soiltile, irrig_frac, irrig_needs_r)
+        CALL irrigation_new(fraction_aeirrig_sw, irrig_needs_r, irrig_actual_r)
+      ENDIF
+      
+      IF (is_omp_root) THEN
+        CALL xios_send_field("routing_irrigation_r",irrig_actual_r)
+        CALL xios_recv_field("routing_irrigation",irrig_actual)
+      ENDIF
+      CALL scatter_omp(irrig_actual,irrigation_mean)
+      irrigation_mean(:)=irrigation_mean(:)/area(:)
+      CALL irrigation_mean_reset()
+    ELSE
+      irrigation_mean(:) = 0
+    ENDIF
+    
+  END SUBROUTINE irrigation_main
+  
+  
+ 
+
+
+  SUBROUTINE irrigation_new(fraction_aeirrig_sw, irrig_needs_r, irrig_actual_r)
+  USE constantes
+  USE routing_native_flow_mod, ONLY : routing_mask_r, routing_flow_get, routing_flow_set 
+  USE mod_orchidee_para
+  USE xios
+  IMPLICIT NONE
+ 
+    REAL(r_std),INTENT(IN)        :: fraction_aeirrig_sw(nbpt)           !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
+    REAL(r_std),INTENT(IN)        :: irrig_needs_r(nbpt_r)               !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
+    REAL(r_std),INTENT(OUT)       :: irrig_actual_r(nbpt_r)             !! Possible irrigation according to the water availability in the reservoirs (kg)
+
+    REAL(r_std)    :: irrig_deficit_r(nbpt_r)             !! Amount of water missing for irrigation (kg)
+
+    LOGICAL        :: IsFail_slow        !! Logical to ask if slow reserv. does not fit irrigation demand
+    LOGICAL        :: IsFail_fast        !! Logical to ask if fast reserv. does not fit irrigation demand
+    LOGICAL        :: IsFail_stre        !! Logical to ask if stream reserv. does not fit irrigation demand
+
+    REAL(r_std)    :: Count_failure_slow(nbpt_r)       !! Counter times slow reserv. does not fit irrigation demand
+    REAL(r_std)    :: Count_failure_fast(nbpt_r)        !! Counter times fast reserv. does not fit irrigation demand
+    REAL(r_std)    :: Count_failure_stre(nbpt_r)        !! Counter times stream reserv. does not fit irrigation demand
+    
+    REAL(r_std)    :: pot_slow_wdr_dummy, pot_fast_wdr_dummy, pot_stre_wdr_dummy
+    REAL(r_std)    :: slow_wdr_dummy, fast_wdr_dummy, stre_wdr_dummy
+
+    REAL(r_std)    :: slow_reservoir_r(nbpt_r)
+    REAL(r_std)    :: fast_reservoir_r(nbpt_r)
+    REAL(r_std)    :: stream_reservoir_r(nbpt_r)
+    REAL(r_std)    :: fraction_aeirrig_sw_mpi(nbp_mpi)
+    REAL(r_std)    :: fraction_aeirrig_sw_r(nbpt_r)
+    REAL(r_std)    :: pcent_vol_irrig
+    INTEGER :: ig
+
+
+    CALL gather_omp(fraction_aeirrig_sw, fraction_aeirrig_sw_mpi)
+      
+    IF (is_omp_root) THEN
+
+      CALL xios_send_field("fraction_aeirrig_sw", fraction_aeirrig_sw_mpi)
+      CALL xios_recv_field("fraction_aeirrig_sw_r", fraction_aeirrig_sw_r) ! ==> need conservative quantity interp
+
+
+      Count_failure_slow(:) = zero !
+      Count_failure_fast(:) = zero !
+      Count_failure_stre(:) = zero !
+
+      CALL routing_flow_get(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)      
+
+      DO ig=1,nbpt_r
+
+        IF (routing_mask_r(ig)) THEN
+        
+          IF (select_source_irrig) THEN
+
+            ! For   irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
+            ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
+
+            ! The new priorization scheme takes into account irrig. infrastructur according to GMIA map
+            ! It also withdraw water according to availability, it means that it wont seek for all the water in then
+            ! stream reservoir, even if this one could respond to the demand by itself
+
+            pot_slow_wdr_dummy = ( 1 - fraction_aeirrig_sw_r(ig)) * avail_reserve(3)*slow_reservoir_r(ig)
+            pot_fast_wdr_dummy = fraction_aeirrig_sw_r(ig) * avail_reserve(2) * fast_reservoir_r(ig)
+            pot_stre_wdr_dummy = fraction_aeirrig_sw_r(ig) * avail_reserve(1) * stream_reservoir_r(ig)
+            pcent_vol_irrig = zero
+            IsFail_slow = .FALSE. !
+            IsFail_fast = .FALSE. !
+            IsFail_stre = .FALSE. !
+
+            irrig_actual_r(ig) = MIN(irrig_needs_r(ig), pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)
+
+            !!   additional IF to calculate pcent_vol_irrig, in the case the total avail.
+            !! water is zero, I.E. when there is no water in surface and fraction_ae = 0,
+            !! so GW is not taken into account
+            !! Note on pcent_vol_irrig: It correspond to the fraction of available water in surface,
+            !! considering environmental needs and irrigation equipement by source from map. It controls
+            !! how the source of water withdrawl, especially when requirements < available water
+
+            IF (  (pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy)  .GT. min_sechiba ) THEN
+
+              pcent_vol_irrig = ( pot_stre_wdr_dummy + pot_fast_wdr_dummy ) / &
+                                ( pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
+
+              !Irrig_actual set to zero, because there is no available water.
+              !Put to avoid negative values due to problems in the Min function
+            
+              irrig_actual_r(ig) = MAX(irrig_actual_r(ig), zero)
+            
+              !Already zero because pcent_vol_irrig initialized to zero
+              !Put here to readability but not necessary
+              !ELSE
+              !    pcent_vol_irrig = zero
+
+            ENDIF
+
+            !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
+            !Volume used for irrigation that comes from GW
+            ! Idem for irrig_fast_source and irrig_str_source
+
+            slow_wdr_dummy = slow_reservoir_r(ig)
+            slow_reservoir_r(ig) = MAX( (un - ( un - fraction_aeirrig_sw_r(ig) ) * avail_reserve(3) ) * &
+                                      slow_reservoir_r(ig), slow_reservoir_r(ig) + &
+                                      MIN( - irrig_actual_r(ig) * (un - pcent_vol_irrig ), &
+                                      avail_reserve(2) * fraction_aeirrig_sw_r(ig) * fast_reservoir_r(ig) + &
+                                      MIN(zero, avail_reserve(1) * fraction_aeirrig_sw_r(ig) * stream_reservoir_r(ig)  - &
+                                      pcent_vol_irrig * irrig_actual_r(ig) ) ) )
+
+            slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
+            irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
+
+            fast_wdr_dummy = fast_reservoir_r(ig)
+            fast_reservoir_r(ig) = MAX( (un - avail_reserve(2) * fraction_aeirrig_sw_r(ig) ) * fast_reservoir_r(ig) , &
+                                    fast_reservoir_r(ig) + MIN(zero, avail_reserve(1) * fraction_aeirrig_sw_r(ig) * stream_reservoir_r(ig)  - &
+                                     pcent_vol_irrig * irrig_actual_r(ig) ) )
+            fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
+            irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
+
+            stre_wdr_dummy = stream_reservoir_r(ig)
+            stream_reservoir_r(ig) = MAX((un - avail_reserve(1)* fraction_aeirrig_sw_r(ig) )*stream_reservoir_r(ig), &
+                                      stream_reservoir_r(ig)  -  pcent_vol_irrig * irrig_actual_r(ig) )
+            stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
+            irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
+
+            irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
+
+            !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
+            !We assume that the pot. requirement = Needs * fraction of area equipped for SW/GW
+            !In the case of surface. we also sustract the withdrawal from Fast/Stream, because both are
+            !  considered as surface water
+
+            IsFail_slow = ( ( irrig_needs_r(ig)*(un - fraction_aeirrig_sw_r(ig)) ) > pot_slow_wdr_dummy )
+            IsFail_fast = ( ( irrig_needs_r(ig)*fraction_aeirrig_sw_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
+            IsFail_stre = ( ( irrig_needs_r(ig)*fraction_aeirrig_sw_r(ig) - fast_wdr_dummy ) > pot_stre_wdr_dummy )
+
+            IF( IsFail_stre ) Count_failure_stre(ig) = un
+            IF( IsFail_fast ) Count_failure_fast(ig) = un
+            IF( IsFail_slow ) Count_failure_slow(ig) = un
+
+          ELSE IF (.NOT. select_source_irrig) THEN
+
+            ! For   irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir
+            ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir
+
+            pot_slow_wdr_dummy = avail_reserve(3)*slow_reservoir_r(ig)
+            pot_fast_wdr_dummy = avail_reserve(2)*fast_reservoir_r(ig)
+            pot_stre_wdr_dummy = avail_reserve(1)*stream_reservoir_r(ig)
+            IsFail_slow = .FALSE. !
+            IsFail_fast = .FALSE. !
+            IsFail_stre = .FALSE. !
+
+            irrig_actual_r(ig) = MIN(irrig_needs_r(ig), pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy )
+
+            !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
+            !Volume used for irrigation that comes from GW
+            ! Idem for irrig_fast_source and irrig_str_source
+
+            slow_wdr_dummy = slow_reservoir_r(ig)
+            slow_reservoir_r(ig) = MAX( (1-avail_reserve(3) )*slow_reservoir_r(ig), slow_reservoir_r(ig)       &
+                                   + MIN(zero, avail_reserve(2)*fast_reservoir_r(ig)                             &
+                                   + MIN(zero, avail_reserve(1)*stream_reservoir_r(ig)-irrig_actual_r(ig))))
+            slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
+            irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
+
+            fast_wdr_dummy = fast_reservoir_r(ig)
+            fast_reservoir_r(ig) = MAX(  (1-avail_reserve(2) )*fast_reservoir_r(ig) , fast_reservoir_r(ig) &
+                                  + MIN(zero, avail_reserve(1)*stream_reservoir_r(ig)-irrig_actual_r(ig)))
+            fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
+            irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
+
+            stre_wdr_dummy = stream_reservoir_r(ig)
+            stream_reservoir_r(ig) = MAX( (1-avail_reserve(1) )*stream_reservoir_r(ig), stream_reservoir_r(ig)-irrig_actual_r(ig) )
+            stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
+            irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
+
+            irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
+
+            !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
+            ! Because it follows the old scheme, we do not separate between surface/gw, but consider that
+            ! priority is given in this order: River, Fast and Slow reservoir.
+
+            IsFail_slow = ( ( irrig_needs_r(ig) - stre_wdr_dummy - fast_wdr_dummy  ) > pot_slow_wdr_dummy )
+            IsFail_fast = ( ( irrig_needs_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
+            IsFail_stre = ( irrig_needs_r(ig) > pot_stre_wdr_dummy )
+
+            IF( IsFail_stre ) Count_failure_stre(ig) = un
+            IF( IsFail_fast ) Count_failure_fast(ig) = un
+            IF( IsFail_slow ) Count_failure_slow(ig) = un
+
+          ENDIF
+
+        ENDIF
+      ENDDO
+  
+      CALL routing_flow_set(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
+
+    ENDIF 
+   
+  END SUBROUTINE irrigation_new
+
+
+  SUBROUTINE irrigation_old(irrig_needs_r, irrig_actual_r)
+  USE constantes
+  USE mod_orchidee_para
+  USE routing_native_flow_mod, ONLY : routing_mask_r, routing_flow_get, routing_flow_set 
+  IMPLICIT NONE
+    REAL(r_std),INTENT(IN)        :: irrig_needs_r(nbpt_r)               !! Total irrigation requirement (water requirements by the crop for its optimal growth) (kg)
+    REAL(r_std),INTENT(OUT)       :: irrig_actual_r(nbpt_r)             !! Possible irrigation according to the water availability in the reservoirs (kg)
+ 
+    REAL(r_std)    :: irrig_deficit_r(nbpt_r)             !! Amount of water missing for irrigation (kg)
+
+    LOGICAL        :: IsFail_slow        !! Logical to ask if slow reserv. does not fit irrigation demand
+    LOGICAL        :: IsFail_fast        !! Logical to ask if fast reserv. does not fit irrigation demand
+    LOGICAL        :: IsFail_stre        !! Logical to ask if stream reserv. does not fit irrigation demand
+
+    REAL(r_std)    :: Count_failure_slow(nbpt_r)       !! Counter times slow reserv. does not fit irrigation demand
+    REAL(r_std)    :: Count_failure_fast(nbpt_r)        !! Counter times fast reserv. does not fit irrigation demand
+    REAL(r_std)    :: Count_failure_stre(nbpt_r)        !! Counter times stream reserv. does not fit irrigation demand
+    
+    REAL(r_std)    :: pot_slow_wdr_dummy, pot_fast_wdr_dummy, pot_stre_wdr_dummy
+    REAL(r_std)    :: slow_wdr_dummy, fast_wdr_dummy, stre_wdr_dummy
+
+    REAL(r_std)    :: slow_reservoir_r(nbpt_r)
+    REAL(r_std)    :: fast_reservoir_r(nbpt_r)
+    REAL(r_std)    :: stream_reservoir_r(nbpt_r)
+
+    INTEGER :: ig
+
+    IF (is_omp_root) THEN
+
+      Count_failure_slow(:) = zero !
+      Count_failure_fast(:) = zero !
+      Count_failure_stre(:) = zero !
+
+      CALL routing_flow_get(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)      
+
+      DO ig=1,nbpt_r
+
+        IF (routing_mask_r(ig)) THEN
+
+          ! Old irrigation scheme as in tag 2.0
+          ! Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the
+          ! Volume used for irrigation that comes from GW
+          ! Idem for irrig_fast_source and irrig_str_source
+
+          pot_slow_wdr_dummy = slow_reservoir_r(ig)
+          pot_fast_wdr_dummy = fast_reservoir_r(ig)
+          pot_stre_wdr_dummy = stream_reservoir_r(ig)
+          IsFail_slow = .FALSE. !
+          IsFail_fast = .FALSE. !
+          IsFail_stre = .FALSE. !
+           
+          irrig_actual_r(ig) = MIN(irrig_needs_r(ig), stream_reservoir_r(ig) + fast_reservoir_r(ig) + slow_reservoir_r(ig) )
+
+          slow_wdr_dummy = slow_reservoir_r(ig)
+          slow_reservoir_r(ig) = MAX(zero, slow_reservoir_r(ig) + MIN(zero, fast_reservoir_r(ig) &
+                                 + MIN(zero, stream_reservoir_r(ig)-irrig_actual_r(ig))))
+          slow_wdr_dummy = slow_wdr_dummy - slow_reservoir_r(ig)
+          irrig_gw_source_r(ig) = irrig_gw_source_r(ig) + slow_wdr_dummy
+
+          fast_wdr_dummy = fast_reservoir_r(ig)
+          fast_reservoir_r(ig) = MAX( zero, fast_reservoir_r(ig) + MIN(zero, stream_reservoir_r(ig)-irrig_actual_r(ig)))
+          fast_wdr_dummy = fast_wdr_dummy - fast_reservoir_r(ig)
+          irrig_fast_source_r(ig) = irrig_fast_source_r(ig) + fast_wdr_dummy
+
+          stre_wdr_dummy = stream_reservoir_r(ig)
+          stream_reservoir_r(ig) = MAX(zero, stream_reservoir_r(ig)-irrig_actual_r(ig) )
+          stre_wdr_dummy = stre_wdr_dummy - stream_reservoir_r(ig)
+          irrig_str_source_r(ig) = irrig_str_source_r(ig) + stre_wdr_dummy
+
+          irrig_deficit_r(ig) = irrig_needs_r(ig)-irrig_actual_r(ig)
+
+          ! A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal
+          ! Because it follows the old scheme, we do not separate between surface/gw, but consider that
+          ! priority is given in this order: River, Fast and Slow reservoir.
+
+          IsFail_slow = ( ( irrig_needs_r(ig) - stre_wdr_dummy - fast_wdr_dummy  ) > pot_slow_wdr_dummy )
+          IsFail_fast = ( ( irrig_needs_r(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy )
+          IsFail_stre = ( irrig_needs_r(ig) > pot_stre_wdr_dummy )
+
+          IF( IsFail_stre ) Count_failure_stre(ig) = un
+          IF( IsFail_fast ) Count_failure_fast(ig) = un
+          IF( IsFail_slow ) Count_failure_slow(ig) = un
+        ELSE
+          irrig_actual_r(ig) = 0
+        ENDIF
+      ENDDO
+  
+      CALL routing_flow_set(slow_reservoir_r=slow_reservoir_r, fast_reservoir_r=fast_reservoir_r, stream_reservoir_r=stream_reservoir_r)
+
+    ENDIF
+  
+  END SUBROUTINE irrigation_old  
+    
+ 
+  SUBROUTINE irrigation_compute_requested_new(dt_routing, root_deficit, soiltile, irrig_frac, irrig_netereq_r)
+  USE constantes
+  USE constantes_soil
+  USE xios
+  USE mod_orchidee_para
+  USE grid, ONLY : area
+  IMPLICIT NONE
+    REAL(r_std),INTENT(IN)   :: dt_routing          
+    REAL(r_std),INTENT(IN)   :: root_deficit(nbpt)               
+    REAL(r_std),INTENT(IN)   :: soiltile(nbpt,nstm)  !! Fraction of each soil tile within vegtot (0-1, unitless)
+    REAL(r_std),INTENT(IN)   :: irrig_frac(nbpt)      !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE.
+    REAL(r_std),INTENT(OUT)  :: irrig_netereq_r(nbpt_r)
+
+    REAL(r_std)              :: irrig_netereq(nbpt)
+    REAL(r_std)              :: irrig_netereq_mpi(nbp_mpi)
+    REAL(r_std)              :: area_mpi(nbp_mpi)
+
+    INTEGER  :: ig
+   
+    irrig_netereq(:)=0
+   
+    DO ig=1,nbpt
+      !It enters to the new irrigation module only if there is an irrigated fraction, if not irrig_netereq = zero for that cell
+      IF ( irrig_frac(ig) .GT. min_sechiba ) THEN
+
+        irrig_netereq(ig) = irrig_netereq(ig) + MIN( irrig_dosmax/3600*dt_routing, &
+                            root_deficit(ig) ) * soiltile(ig, irrig_st) * vegtot_mean(ig)
+       ! By definition, irrig_dosmax is in kg/m2 of soil tile/hour,dividing by 3600(seconds/hour) * DT_ROUTING  !
+       ! = kg/m2 of soil tile/(routing timestep)
+       ! irrig_netereq(kg/m2 of grid cell / routing timstep ) is equal to
+       ! root_deficit (kg/m2 of soil tile) * soiltile*vegtot (fraction of soil tile at cell level) = kg/m2 of grid cell
+
+       IF (.NOT. irrigated_soiltile .AND. ( soiltile(ig,irrig_st) .GT. min_sechiba ) .AND. (vegtot_mean(ig) .GT. min_sechiba) ) THEN
+         ! Irrigated_soiltile asks if there is an independent soil tile for irrigated crops. If not,
+         ! actual volume calculated for irrig_netereq assumed that the whole SOILTILE was irrigated, but in this case
+         ! just a fraction of the irrig_st (irrigated soil tile, by default = 3) is actually irrigated,
+         ! and irrig_netereq needs to be reduced by irrig_frac/( soiltile * vegtot ) (note that it is max = 1 thanks to irrig_frac calculation in l. 424)
+         ! Demand(ST3)*irrig_frac/(soiltile(3)*vegtot) = irrig_netereq_In_ST3, then
+         !irrig_netereq_In_ST3 * (soiltile(3)*vegtot) = irrig_netereq at grid scale = Demand(ST3)*irrig_frac.
+         irrig_netereq(ig) = irrig_netereq(ig) * irrig_frac(ig) / ( soiltile(ig,irrig_st) * vegtot_mean(ig) )
+          !irrig_netereq = kg/m2 of grid cell
+        ENDIF
+      ENDIF
+    ENDDO
+
+    CALL gather_omp(irrig_netereq, irrig_netereq_mpi)
+    CALL gather_omp(area, area_mpi)
+      
+    IF (is_omp_root) THEN
+      CALL xios_send_field("irrig_netereq",irrig_netereq_mpi*area_mpi) !! contfrac ?
+      CALL xios_recv_field("irrig_netereq_r",irrig_netereq_r) ! ==> need conservative quantity interp
+    ENDIF
+      
+  END SUBROUTINE irrigation_compute_requested_new
+  
+ 
+  SUBROUTINE irrigation_compute_requested_old(reinfiltration, irrigated, irrig_need_r)
+  USE constantes
+  USE xios
+  USE mod_orchidee_para
+  USE grid, ONLY : area
+  IMPLICIT NONE
+    REAL(r_std),INTENT(IN)   :: reinfiltration(nbpt)               
+    REAL(r_std),INTENT(IN)   :: irrigated(nbpt)               
+    REAL(r_std),INTENT(OUT)  :: irrig_need_r(nbpt_r)
+
+    REAL(r_std)              :: irrig_netereq(nbpt)
+    REAL(r_std)              :: irrig_netereq_mpi(nbp_mpi)
+    REAL(r_std)              :: area_mpi(nbp_mpi)
+
+    INTEGER  :: ig, ir
+  
+      irrig_netereq(:)=0
+  
+      DO ig=1,nbpt
+        IF ((vegtot_mean(ig) .GT. min_sechiba) .AND. (humrel_mean(ig) .LT. un-min_sechiba) .AND. (runoff_mean(ig) .LT. min_sechiba) ) THEN
+          irrig_netereq(ig) = (irrigated(ig)/area(ig)) * MAX(zero, transpot_mean(ig) - (precip_mean(ig)+reinfiltration(ig))) ! ==> ok kg
+        ELSE
+          irrig_netereq(ig) = 0
+        ENDIF
+      ENDDO
+    
+      CALL gather_omp(irrig_netereq, irrig_netereq_mpi)
+      CALL gather_omp(area, area_mpi)
+
+      IF (is_omp_root) THEN
+        CALL xios_send_field("irrig_netereq",irrig_netereq_mpi*area_mpi) ! contfrac ??
+        CALL xios_recv_field("irrig_netereq_r",irrig_need_r) ! ==> need conservative quantity interp
+      ENDIF
+      
+  END SUBROUTINE irrigation_compute_requested_old
+ 
+ SUBROUTINE test
+ USE routing_native_flow_mod, ONLY : routing_mask_r
+ 
+ END SUBROUTINE TEST
+ 
+ 
+
+
+
+
+!for now  SUBROUTINE abduction
+!for now    REAL(r_std), DIMENSION(nbpt)        :: irrig_adduct              !! Amount of water carried over from other basins for irrigation (kg)
+!for now  
+!for now    DO ig=1,nbpt
+
+!for now      !
+!for now      ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only).
+!for now      ! If we find that then we create some adduction from that subbasin to the one where we need it for
+!for now      ! irrigation.
+!for now      !
+
+!for now      !> If crops water requirements have not been supplied (irrig_deficit>0), we check if we cannot find the missing water
+!for now      !> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction
+!for now      !> from that subbasin to the one where we need it for irrigation.
+!for now      !>
+!for now    
+!for now      DO ib=1,nbasmax
+!for now        stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
+
+!for now        DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)
+!for now          fi = MAXLOC(stream_reservoir(ig,:))
+!for now          ib2 = fi(1)
+
+!for now          irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), a_stream_adduction * stream_reservoir(ig,ib2))
+!for now          stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)
+!for now          irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)
+!for now          stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:))
+
+!for now        ENDDO
+!for now      ENDDO
+!for now
+!for now    ENDDO
+
+!for now    !
+!for now    ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams
+!for now    ! which can feed irrigation
+!for now    !
+!for now    !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes
+!for now    !> to the one where we need it for irrigation.
+!for now    !
+!for now    IF (is_root_prc) THEN
+!for now      ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), &
+!for now               irrig_adduct_glo(nbp_glo, nbasmax), stat=ier)
+!for now    ELSE
+!for now      ALLOCATE(irrig_deficit_glo(0, 0), stream_reservoir_glo(0, 0), irrig_adduct_glo(0, 0), stat=ier)
+!for now    ENDIF
+!for now    IF (ier /= 0) CALL ipslerr_p(3,'routing_flow','Pb in allocate for irrig_deficit_glo, stream_reservoir_glo,...','','')
+
+!for now    CALL gather(irrig_deficit, irrig_deficit_glo)
+!for now    CALL gather(stream_reservoir,  stream_reservoir_glo)
+!for now    CALL gather(irrig_adduct, irrig_adduct_glo)
+
+!for now    IF (is_root_prc) THEN
+!for now      !
+!for now      DO ig=1,nbp_glo
+!for now        ! Only work if the grid box is smaller than 100x100km. Else the piplines we build
+!for now        ! here would be too long to be reasonable.
+!for now        IF ( resolution_g(ig,1) < 100000. .AND. resolution_g(ig,2) < 100000. ) THEN
+!for now          
+!for now          DO ib=1,nbasmax
+!for now          !
+!for now            IF ( irrig_deficit_glo(ig,ib)  > min_sechiba ) THEN
+!for now            !
+!for now              streams_around(:,:) = zero
+!for now              !
+!for now              DO in=1,NbNeighb
+!for now                ig2 = neighbours_g(ig,in)
+!for now                IF (ig2 .GT. 0 ) THEN
+!for now                  streams_around(in,:) = a_stream_adduction * stream_reservoir_glo(ig2,:)
+!for now                  igrd(in) = ig2
+!for now                ENDIF
+!for now              ENDDO
+!for now              !
+!for now              IF ( MAXVAL(streams_around) .GT. zero ) THEN
+!for now              !
+!for now                ff=MAXLOC(streams_around)
+!for now                ig2=igrd(ff(1))
+!for now                ib2=ff(2)
+!for now                !
+!for now                IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. a_stream_adduction * stream_reservoir_glo(ig2,ib2) > zero ) THEN
+!for now                  adduction = MIN(irrig_deficit_glo(ig,ib), a_stream_adduction * stream_reservoir_glo(ig2,ib2))
+!for now                  stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction
+!for now                  irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction
+!for now                  irrig_adduct_glo(ig,ib) = irrig_adduct_glo(ig,ib) + adduction
+!for now                ENDIF
+!for now                !
+!for now              ENDIF
+!for now              !
+!for now            ENDIF
+!for now            !
+!for now          ENDDO
+!for now      
+!for now        ENDIF
+
+!for now      ENDDO
+!for now          !
+!for now    ENDIF
+
+!for now    CALL scatter(irrig_deficit_glo, irrig_deficit)
+!for now    CALL scatter(stream_reservoir_glo,  stream_reservoir)
+!for now    CALL scatter(irrig_adduct_glo, irrig_adduct)
+
+!for now    DEALLOCATE(irrig_deficit_glo, stream_reservoir_glo, irrig_adduct_glo)  
+!for now  
+!for now  
+!for now  END SUBROUTINE abduction
+  
+
+
+END MODULE routing_native_irrig_mod
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_lake.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_lake.f90	(revision 8227)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_lake.f90	(revision 8227)
@@ -0,0 +1,291 @@
+MODULE routing_native_lake_mod
+  USE constantes
+   
+  PRIVATE
+  
+  PUBLIC :: routing_lake_initialize, routing_lake_mean_make, routing_lake_main, routing_lake_route_coast, routing_lake_finalize
+  
+  REAL(r_std), SAVE, ALLOCATABLE :: lake_reservoir(:)
+!$OMP THREADPRIVATE(lake_reservoir)
+
+  REAL(r_std), SAVE, ALLOCATABLE :: humrel_mean(:)
+!$OMP THREADPRIVATE(humrel_mean)
+  
+  INTEGER,SAVE  :: nbpt
+!$OMP THREADPRIVATE( nbpt)
+ 
+  LOGICAL,SAVE :: do_swamps
+!$OMP THREADPRIVATE(do_swamps)
+
+  INTEGER,SAVE :: nb_coast_cells
+!$OMP THREADPRIVATE(nb_coast_cells)
+
+REAL(r_std), PARAMETER   :: maxevap_lake = 7.5/86400.   !! Maximum evaporation rate from lakes (kg/m^2/s)
+
+REAL(r_std), SAVE    :: max_lake_reservoir           !! Maximum limit of water in lake_reservoir [kg/m2]
+!$OMP THREADPRIVATE(max_lake_reservoir)
+  
+LOGICAL, ALLOCATABLE :: is_coastline(:)
+!$OMP THREADPRIVATE(is_coastline)  
+
+CONTAINS
+ 
+ SUBROUTINE routing_lake_initialize(kjit, rest_id, nbpt_, contfrac)
+ IMPLICIT NONE
+   INTEGER, INTENT(IN)     :: kjit
+   INTEGER, INTENT(IN)     :: rest_id
+   INTEGER,INTENT(IN)      :: nbpt_
+   REAL(r_std), INTENT(IN)    :: contfrac(nbpt)     !! Fraction of land in each grid box (unitless;0-1)
+
+
+   CALL routing_lake_init_local(kjit, rest_id, nbpt_, contfrac)
+   CALL routing_lake_mean_init(kjit, rest_id)
+
+ END SUBROUTINE routing_lake_initialize
+
+
+ SUBROUTINE routing_lake_finalize(kjit, rest_id)
+ IMPLICIT NONE
+   INTEGER, INTENT(IN)     :: kjit
+   INTEGER, INTENT(IN)     :: rest_id
+
+   CALL routing_lake_finalize_local(kjit, rest_id)
+   CALL routing_lake_mean_finalize(kjit, rest_id)
+
+ END SUBROUTINE routing_lake_finalize
+
+ 
+ SUBROUTINE routing_lake_init_local(kjit, rest_id, nbpt_, contfrac)
+ USE mod_orchidee_para, ONLY : reduce_sum, bcast
+ USE ioipsl_para
+ USE grid, ONLY : nbp_glo, index_g
+ USE sechiba_io_p
+ USE routing_native_flow_mod, ONLY : compute_coastline
+ IMPLICIT NONE
+   INTEGER, INTENT(IN) :: kjit
+   INTEGER, INTENT(IN) :: rest_id
+   INTEGER,INTENT(IN) :: nbpt_
+   REAL(r_std), INTENT(IN)    :: contfrac(nbpt)     !! Fraction of land in each grid box (unitless;0-1)
+   
+   INTEGER :: ier    
+   CHARACTER(LEN=80)   :: var_name       !! To store variables names for I/O (unitless)
+ 
+   INTEGER  :: ig
+   INTEGER  :: nb_cells
+   
+   nbpt=nbpt_
+   
+   
+    ALLOCATE(lake_reservoir(nbpt),stat=ier)
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','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)
+       
+    !
+    !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]
+    !
+    do_swamps = .FALSE.
+    CALL getin_p('DO_SWAMPS', do_swamps)
+
+    !Config Key   = MAX_LAKE_RESERVOIR
+    !Config Desc  = Maximum limit of water in lake_reservoir
+    !Config If    = RIVER_ROUTING
+    !Config Def   = 7000
+    !Config Help  =
+    !Config Units = [kg/m2(routing area)]
+    max_lake_reservoir = 7000
+    CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir)
+        
+    
+    ! compute  number of coast cells
+    ! WARNING : contfrac is fraction of ter, but 1-contfrac is not fraction of ocean because of landice
+    ! => to be corrected...
+    ALLOCATE(is_coastline(nbpt))
+    CALL compute_coastline(contfrac, is_coastline)
+    nb_cells=0
+    
+    DO ig=1,nbpt
+      IF (is_coastline(ig)) nb_cells=nb_cells+1
+    ENDDO
+    
+    CALL reduce_sum(nb_cells, nb_coast_cells)
+    CALL bcast(nb_coast_cells)
+  
+  END SUBROUTINE routing_lake_init_local
+
+
+  SUBROUTINE routing_lake_finalize_local(kjit, rest_id)
+  USE ioipsl_para
+  USE grid, ONLY : nbp_glo, index_g
+  IMPLICIT NONE
+    INTEGER, INTENT(IN) :: kjit
+    INTEGER, INTENT(IN) :: rest_id
+ 
+    CALL restput_p (rest_id, 'lakeres', nbp_glo, 1, 1, kjit, lake_reservoir, 'scatter',  nbp_glo, index_g)
+    DEALLOCATE(lake_reservoir)
+
+  END SUBROUTINE  routing_lake_finalize_local
+  
+  
+
+ SUBROUTINE routing_lake_mean_init(kjit, rest_id)
+ USE ioipsl_para
+ USE grid
+ USE sechiba_io_p
+ IMPLICIT NONE
+    INTEGER, INTENT(IN) :: kjit
+    INTEGER, INTENT(IN) :: rest_id
+
+    INTEGER :: ier    
+    CHARACTER(LEN=80)   :: var_name       !! To store variables names for I/O (unitless)
+  
+    ALLOCATE(humrel_mean(nbpt), stat=ier)
+
+    IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for humrel_mean','','')
+    var_name = 'humrel_lake'
+    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)
+   
+ END SUBROUTINE  routing_lake_mean_init
+ 
+ SUBROUTINE routing_lake_mean_make(dt_routing, humrel, veget_max)
+ USE constantes
+ USE pft_parameters
+ IMPLICIT NONE
+   REAL(r_std), INTENT(IN)    :: dt_routing 
+   REAL(r_std), INTENT(IN)    :: humrel(:,:)       !! Soil moisture stress, root extraction potential (unitless)
+   REAL(r_std), INTENT(IN)    :: veget_max(:,:)       !! Soil moisture stress, root extraction potential (unitless)
+   INTEGER :: jv
+   
+   IF ( .NOT. old_irrig_scheme ) THEN
+     DO jv=1,nvm
+       humrel_mean(:) = humrel_mean(:) + humrel(:,jv)*veget_max(:,jv)*dt_sechiba/dt_routing
+     ENDDO
+   ELSE
+     DO jv=2,nvm
+       humrel_mean(:) = humrel_mean(:) + humrel(:,jv)*veget_max(:,jv)*dt_sechiba/dt_routing
+     ENDDO
+   ENDIF
+    
+ END SUBROUTINE routing_lake_mean_make  
+
+  
+ SUBROUTINE routing_lake_mean_reset
+ IMPLICIT NONE
+   humrel_mean(:) = 0
+ END SUBROUTINE routing_lake_mean_reset 
+
+ 
+ SUBROUTINE routing_lake_mean_finalize(kjit, rest_id)
+ USE ioipsl_para
+ USE grid
+ IMPLICIT NONE
+   INTEGER, INTENT(IN) :: kjit
+   INTEGER, INTENT(IN) :: rest_id
+   
+   CALL restput_p (rest_id, 'humrel_lake', nbp_glo, 1, 1, kjit, humrel_mean, 'scatter',  nbp_glo, index_g)
+   DEALLOCATE(humrel_mean)
+ 
+ END SUBROUTINE routing_lake_mean_finalize
+ 
+ 
+ SUBROUTINE routing_lake_route_coast(contfrac, coastalflow)
+ USE grid, ONLY : area
+ USE mod_orchidee_para, ONLY : reduce_sum, bcast
+
+ IMPLICIT NONE
+   REAL(r_std), INTENT(IN)    :: contfrac(nbpt)     !! Fraction of land in each grid box (unitless;0-1)
+   REAL(r_std), INTENT(INOUT) :: coastalflow(nbpt)   !! Water inflow to the lakes (kg/dt)
+   
+   REAL(r_std) :: sum_lake_overflow, lake_overflow, total_lake_overflow
+   INTEGER :: ig
+    
+   !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it 
+   !! uniformly over all possible the coastflow gridcells
+    
+   ! Calculate lake_overflow and remove it from lake_reservoir
+    sum_lake_overflow=0
+    DO ig=1,nbpt
+      lake_overflow = MAX(0., lake_reservoir(ig) - max_lake_reservoir*area(ig)*contfrac(ig))
+      lake_reservoir(ig) = lake_reservoir(ig) - lake_overflow
+      sum_lake_overflow = sum_lake_overflow + lake_overflow
+    END DO
+
+   ! Calculate the sum of the lake_overflow and distribute it uniformly over all gridboxes
+    CALL reduce_sum(sum_lake_overflow, total_lake_overflow)
+    CALL bcast(total_lake_overflow)
+
+    WHERE (is_coastline) coastalflow = coastalflow + total_lake_overflow/nb_coast_cells 
+ 
+ END SUBROUTINE routing_lake_route_coast
+
+
+
+ 
+ SUBROUTINE routing_lake_main(dt_routing, contfrac, lakeinflow, return_lakes)
+
+    USE grid, ONLY : area
+    
+    IMPLICIT NONE
+    !! 0 Variable and parameter description
+    !! 0.1 Input variables
+
+    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)    :: 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
+       IF ( do_swamps ) THEN
+         ! Calculate a return flow that will be extracted from the lake reservoir and reinserted in the soil in hydrol
+         ! Uptake in Kg/dt
+         refill = MAX(zero, maxevap_lake * (un - humrel_mean(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
+          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.
+!ym for now       lake_diag(ig) = lake_reservoir(ig)/total_area
+       !
+    ENDDO
+
+    CALL routing_lake_mean_reset
+ 
+  END SUBROUTINE routing_lake_main
+
+END MODULE routing_native_lake_mod
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_para.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_para.f90	(revision 8227)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_para.f90	(revision 8227)
@@ -0,0 +1,270 @@
+MODULE routing_native_para
+
+ PRIVATE
+
+ INTEGER,SAVE  :: ni
+ INTEGER,SAVE  :: nj
+
+ INTEGER,SAVE  :: ni_glo
+ INTEGER,SAVE  :: nj_glo
+ INTEGER,SAVE  :: nbpt_r 
+ INTEGER,SAVE  :: nbpt_rp1 
+ INTEGER,ALLOCATABLE :: r_to_rp1(:)
+ INTEGER,ALLOCATABLE :: rp1_to_r(:)
+
+ TYPE mpi_neighbours
+   INTEGER             :: rank
+   INTEGER             :: n
+   INTEGER,ALLOCATABLE :: index(:)
+   INTEGER,ALLOCATABLE :: ibuffer(:)
+   REAL,ALLOCATABLE :: rbuffer(:)
+ END TYPE
+ 
+ TYPE(mpi_neighbours),SAVE, ALLOCATABLE :: send_neighbours(:)
+ TYPE(mpi_neighbours),SAVE, ALLOCATABLE :: recv_neighbours(:)
+ 
+ INTERFACE update_halo
+   PROCEDURE update_halo_integer, update_halo_real
+ END INTERFACE
+
+ PUBLIC routing_para_initialize, update_halo, r_to_rp1, rp1_to_r
+
+CONTAINS
+
+  SUBROUTINE routing_para_initialize
+  USE mod_orchidee_para, ONLY : mpi_rank,mpi_size, MPI_COMM_ORCH, MPI_INT_ORCH
+  USE xios
+  IMPLICIT NONE
+    INCLUDE 'mpif.h'
+    REAL, ALLOCATABLE :: my_rank(:)
+    REAL, ALLOCATABLE :: global_ranks(:)
+    REAL, ALLOCATABLE :: my_index(:)
+    REAL, ALLOCATABLE :: global_index(:)
+    INTEGER :: nb_send_neighbours
+    INTEGER :: nb_recv_neighbours
+    INTEGER :: nb_points_to_recv
+    INTEGER, ALLOCATABLE :: list_rank_to_recv(:)
+     INTEGER :: i,j,k,ij,ind,r,rp1
+    INTEGER :: rank, pos_rank
+    INTEGER,ALLOCATABLE :: send_request(:)
+    INTEGER,ALLOCATABLE :: recv_request(:)
+    INTEGER,ALLOCATABLE :: send_status(:,:)
+    INTEGER,ALLOCATABLE :: recv_status(:,:)
+    INTEGER :: nb_points_ranks(0:mpi_size-1)
+    INTEGER :: ierr
+    INTEGER :: test
+
+    CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo)    ! get routing domain dimension
+
+    nbpt_r= ni*nj                                                
+    nbpt_rp1= (ni+2)*(nj+2)
+    
+    ALLOCATE(my_rank(nbpt_r))
+    ALLOCATE(global_ranks(nbpt_rp1))
+    ALLOCATE(my_index(nbpt_r))
+    ALLOCATE(global_index(nbpt_rp1))
+    ALLOCATE(r_to_rp1(nbpt_r))
+    ALLOCATE(rp1_to_r(nbpt_rp1))
+
+    my_rank(:)=mpi_rank
+    rp1_to_r(:)=-1
+
+    DO j=1,nj
+      DO i=1,ni
+        r=((j-1)*ni)+i
+        rp1=(j+1-1)*(ni+2)+i+1 
+        r_to_rp1(r)=rp1
+        rp1_to_r(rp1)=r
+        my_index(r) = rp1
+      ENDDO
+    ENDDO
+
+    CALL xios_send_field("my_rank",my_rank)
+    CALL xios_recv_field("global_ranks",global_ranks)
+    CALL xios_send_field("my_index",my_index)
+    CALL xios_recv_field("global_index",global_index)
+
+    ! number of point received by rank 
+    nb_points_ranks(:)=0
+    DO rp1=1,nbpt_rp1
+      nb_points_ranks(global_ranks(rp1)) = nb_points_ranks(global_ranks(rp1))+1
+    ENDDO
+
+    nb_points_ranks(mpi_rank) = 0 ! not me
+
+    
+    nb_recv_neighbours=0
+    nb_points_to_recv=0
+    DO rank=0, mpi_size-1
+      IF (nb_points_ranks(rank)/=0) THEN
+        nb_recv_neighbours=nb_recv_neighbours+1
+        nb_points_to_recv = nb_points_to_recv + nb_points_ranks(rank)
+      ENDIF
+    ENDDO
+    
+    ALLOCATE(recv_neighbours(nb_recv_neighbours))
+    ALLOCATE(list_rank_to_recv(0:mpi_size-1))
+    k=1
+    DO rank=0, mpi_size-1
+      IF (nb_points_ranks(rank)/=0) THEN
+        recv_neighbours(k)%rank = rank
+        recv_neighbours(k)%n = 0
+        ALLOCATE(recv_neighbours(k)%index(nb_points_ranks(rank)))
+        ALLOCATE(recv_neighbours(k)%ibuffer(nb_points_ranks(rank)))
+        ALLOCATE(recv_neighbours(k)%rbuffer(nb_points_ranks(rank)))
+        list_rank_to_recv(rank) = k
+        k=k+1
+      ENDIF
+    ENDDO    
+    
+    DO rp1=1,nbpt_rp1
+      rank = global_ranks(rp1)
+      IF (rank /= mpi_rank) THEN
+        k = list_rank_to_recv(rank) 
+        recv_neighbours(k)%n=recv_neighbours(k)%n+1
+        recv_neighbours(k)%index(recv_neighbours(k)%n) = rp1
+      ENDIF
+    ENDDO
+
+
+    nb_points_ranks(:) = 0
+    DO k=1,size(recv_neighbours)
+      nb_points_ranks(recv_neighbours(k)%rank) = 1
+    ENDDO
+    
+    CALL MPI_ALLREDUCE(MPI_IN_PLACE , nb_points_ranks, mpi_size, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
+
+    nb_send_neighbours = nb_points_ranks(mpi_rank)
+    ALLOCATE(send_neighbours(nb_send_neighbours))
+
+
+    ALLOCATE(send_request(size(recv_neighbours)))
+    ALLOCATE(send_status(MPI_STATUS_SIZE,size(recv_neighbours)))
+    ALLOCATE(recv_request(size(send_neighbours)))
+    ALLOCATE(recv_status(MPI_STATUS_SIZE, size(send_neighbours)))
+    
+
+    DO k=1,size(recv_neighbours)
+      CALL MPI_ISEND(recv_neighbours(k)%n, 1, MPI_INT_ORCH, recv_neighbours(k)%rank, 0, MPI_COMM_ORCH, send_request(k) , ierr )
+    ENDDO
+
+    DO k=1,size(send_neighbours)
+      CALL MPI_IRecv(send_neighbours(k)%n, 1, MPI_INT_ORCH, MPI_ANY_SOURCE, 0, MPI_COMM_ORCH, recv_request(k) , ierr )
+    ENDDO
+
+    CALL MPI_WAITALL(size(recv_neighbours),send_request, send_status, ierr)
+    CALL MPI_WAITALL(size(send_neighbours),recv_request, recv_status, ierr)
+    
+    DO k=1,size(send_neighbours)
+      send_neighbours(k)%rank = recv_status(MPI_SOURCE, k)
+      ALLOCATE(send_neighbours(k)%index(send_neighbours(k)%n))
+      ALLOCATE(send_neighbours(k)%ibuffer(send_neighbours(k)%n))
+      ALLOCATE(send_neighbours(k)%rbuffer(send_neighbours(k)%n))
+    ENDDO
+
+    DO k=1,size(recv_neighbours)
+      DO ij=1, recv_neighbours(k)%n
+        ind = recv_neighbours(k)%index(ij)
+        recv_neighbours(k)%ibuffer(ij) = global_index(ind)
+      ENDDO
+    ENDDO
+
+    DO k=1,size(recv_neighbours)
+      CALL MPI_ISEND(recv_neighbours(k)%ibuffer, recv_neighbours(k)%n, MPI_INT_ORCH, recv_neighbours(k)%rank, 1, MPI_COMM_ORCH, send_request(k) , ierr )
+    ENDDO
+
+    DO k=1,size(send_neighbours)
+      CALL MPI_IRECV(send_neighbours(k)%index, send_neighbours(k)%n, MPI_INT_ORCH, send_neighbours(k)%rank, 1, MPI_COMM_ORCH, recv_request(k) , ierr )
+    ENDDO
+
+    CALL MPI_WAITALL(size(recv_neighbours),send_request, send_status, ierr)
+    CALL MPI_WAITALL(size(send_neighbours),recv_request, recv_status, ierr)
+
+    global_ranks(:)=1000000000
+    DO r=1,nbpt_r
+      global_ranks(r_to_rp1(r)) = mpi_rank 
+    ENDDO
+    
+    CALL update_halo_real(global_ranks)
+
+  END SUBROUTINE routing_para_initialize
+
+  SUBROUTINE update_halo_integer(field)
+  USE mod_orchidee_para, ONLY :  MPI_COMM_ORCH, MPI_INT_ORCH
+  IMPLICIT NONE
+    INCLUDE 'mpif.h'
+    INTEGER, INTENT(INOUT) :: field(nbpt_rp1)
+    INTEGER :: send_request(size(send_neighbours))
+    INTEGER :: send_status(MPI_STATUS_SIZE,size(send_neighbours))
+    INTEGER :: recv_request(size(recv_neighbours))
+    INTEGER :: recv_status(MPI_STATUS_SIZE, size(recv_neighbours))
+    INTEGER :: ij, k,ind
+    INTEGER :: ierr
+   
+    DO k=1,size(send_neighbours)
+      DO ij=1, send_neighbours(k)%n
+        ind = send_neighbours(k)%index(ij)
+        send_neighbours(k)%ibuffer(ij) = field(ind)
+      ENDDO
+    ENDDO
+
+    DO ij=1,size(send_neighbours)
+      CALL MPI_ISEND(send_neighbours(ij)%ibuffer, send_neighbours(ij)%n, MPI_INT_ORCH, send_neighbours(ij)%rank, 0, MPI_COMM_ORCH, send_request(ij) , ierr )
+    ENDDO
+
+    DO ij=1,size(recv_neighbours)
+      CALL MPI_IRECV(recv_neighbours(ij)%ibuffer, recv_neighbours(ij)%n, MPI_INT_ORCH, recv_neighbours(ij)%rank, 0, MPI_COMM_ORCH, recv_request(ij) , ierr )
+    ENDDO
+
+    CALL MPI_WAITALL(size(send_neighbours),send_request, send_status, ierr)
+    CALL MPI_WAITALL(size(recv_neighbours),recv_request, recv_status, ierr)
+
+    DO k=1,size(recv_neighbours)
+      DO ij=1, recv_neighbours(k)%n
+        ind = recv_neighbours(k)%index(ij)
+        field(ind) = recv_neighbours(k)%ibuffer(ij) 
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE update_halo_integer
+
+  SUBROUTINE update_halo_real(field)
+  USE mod_orchidee_para, ONLY :  MPI_COMM_ORCH, MPI_REAL_ORCH
+  IMPLICIT NONE
+    INCLUDE 'mpif.h'
+    REAL, INTENT(INOUT) :: field(nbpt_rp1)
+    INTEGER :: send_request(size(send_neighbours))
+    INTEGER :: send_status(MPI_STATUS_SIZE,size(send_neighbours))
+    INTEGER :: recv_request(size(recv_neighbours))
+    INTEGER :: recv_status(MPI_STATUS_SIZE, size(recv_neighbours))
+    INTEGER :: ij, k,ind
+    INTEGER :: ierr
+
+    DO k=1,size(send_neighbours)
+      DO ij=1, send_neighbours(k)%n
+        ind = send_neighbours(k)%index(ij)
+        send_neighbours(k)%rbuffer(ij) = field(ind)
+      ENDDO
+    ENDDO
+   
+    DO ij=1,size(send_neighbours)
+      CALL MPI_ISEND(send_neighbours(ij)%rbuffer, send_neighbours(ij)%n, MPI_REAL_ORCH, send_neighbours(ij)%rank, 0, MPI_COMM_ORCH, send_request(ij) , ierr )
+    ENDDO
+
+    DO ij=1,size(recv_neighbours)
+      CALL MPI_IRECV(recv_neighbours(ij)%rbuffer, recv_neighbours(ij)%n, MPI_REAL_ORCH, recv_neighbours(ij)%rank, 0, MPI_COMM_ORCH, recv_request(ij) , ierr )
+    ENDDO
+
+    CALL MPI_WAITALL(size(send_neighbours),send_request, send_status, ierr)
+    CALL MPI_WAITALL(size(recv_neighbours),recv_request, recv_status, ierr)
+
+    DO k=1,size(recv_neighbours)
+      DO ij=1, recv_neighbours(k)%n
+        ind = recv_neighbours(k)%index(ij)
+        field(ind) = recv_neighbours(k)%rbuffer(ij) 
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE update_halo_real
+
+END MODULE routing_native_para
Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90	(revision 8221)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_wrapper.f90	(revision 8227)
@@ -35,4 +35,5 @@
   USE routing_highres
   USE routing_simple
+  USE routing_native_mod
   USE constantes_soil
 
@@ -76,4 +77,6 @@
     ELSEIF(routing_method=='simple') THEN  
        CALL routing_simple_xios_initialize
+    ELSEIF(routing_method=='native') THEN
+       CALL routing_native_xios_initialize
     ENDIF
 
@@ -179,4 +182,20 @@
                                           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 IF(routing_method== 'native') THEN
+
+       CALL routing_native_initialize(    kjit,        nbpt,           index,                 &
+                                          rest_id,     hist_id,        hist2_id,   lalo,      &
+                                          neighbours,  resolution,     contfrac,   stempdiag, &
+                                          returnflow,  reinfiltration, irrigation, riverflow, &
+                                          coastalflow, flood_frac,     flood_res , irrigated_next)
 
        riverflow(:) = zero
@@ -284,5 +303,14 @@
             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)
+            rest_id, hist_id, hist2_id) 
+
+    ELSE IF(routing_method=='native') THEN
+
+       CALL routing_native_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, &
+            soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) 
     ENDIF
 
@@ -331,4 +359,8 @@
        CALL routing_simple_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
 
+    ELSE IF(routing_method=='native') THEN
+
+       CALL routing_native_finalize( kjit, nbpt, rest_id, flood_frac, flood_res )
+
     ENDIF
 
@@ -363,5 +395,8 @@
 
        CALL routing_simple_clear
-
+    
+    ELSE IF(routing_method=='native') THEN
+
+       CALL routing_native_clear
     ENDIF
 
