Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90	(revision 7476)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_simple.f90	(revision 7481)
@@ -287,4 +287,6 @@
   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)
@@ -308,4 +310,7 @@
   !$OMP THREADPRIVATE(split_routing)
 
+  REAL(r_std), SAVE                                          :: max_lake_reservoir           !! Maximum limit of water in lake_reservoir [kg/m2]
+  !$OMP THREADPRIVATE(max_lake_reservoir)
+
 
   INTEGER(i_std), PARAMETER :: nb_stations=14
@@ -405,4 +410,5 @@
     REAL(r_std) :: area_mpi(nbp_mpi)
     INTEGER :: basins_count_mpi
+    INTEGER :: nb_coast_points
     INTEGER :: ierr
     LOGICAL :: file_exists
@@ -415,5 +421,5 @@
     CALL gather_omp(contfrac,contfrac_mpi)
     CALL gather_omp(area, area_mpi)
-
+ 
     IF (is_omp_root) THEN
 
@@ -445,5 +451,5 @@
        coast_mask=.FALSE.
        WHERE( contfrac_mpi(:)< 1.-1.e-5) coast_mask(:)=.TRUE.            ! create mask for coastal cells on native grid
-
+       
        INQUIRE(FILE="routing_start.nc", EXIST=file_exists)
 
@@ -485,4 +491,9 @@
        frac_routing_r(:)=0
        WHERE( .NOT. coast_mask(:)) mask_native(:)=1
+
+       nb_coast_points=SUM(1-mask_native)
+       CALL reduce_sum_mpi(nb_coast_points, total_coast_points)
+       CALL bcast_mpi(total_coast_points)
+       
        CALL xios_send_field("mask_native_lake",mask_native)              ! send full land point to XIOS (native grid)
        CALL xios_recv_field("frac_routing_lake_r",frac_routing_r)        ! receive fraction of intersected cell by full land, on routing grid 
@@ -734,5 +745,5 @@
     
     USE xios
-    USE grid, ONLY : area
+    USE grid, ONLY : area 
     IMPLICIT NONE
     INCLUDE "mpif.h"
@@ -779,4 +790,6 @@
     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)                                  :: lake_reservoir_mpi(nbp_mpi) ! cell area         
 
     ! from input model -> routing_grid
@@ -848,5 +861,5 @@
     REAL(r_std) :: basins_riverflow_mpi(0:basins_count)
     REAL(r_std) :: basins_riverflow(0:basins_count)
-
+    REAL(r_std) :: lake_overflow,sum_lake_overflow, total_lake_overflow
     INTEGER :: ierr
 
@@ -873,4 +886,6 @@
     CALL gather_omp(runoff_omp,runoff)
     CALL gather_omp(drainage_omp, drainage)
+    CALL gather_omp(area, area_mpi)
+    CALL gather_omp(lake_reservoir, lake_reservoir_mpi)
 
     IF (is_omp_root) THEN
@@ -1128,5 +1143,21 @@
             " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%"
 
-
+   
+   !! 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,nbp_mpi
+         lake_overflow = MAX(0., lake_reservoir_mpi(ig) - max_lake_reservoir*area_mpi(ig))
+         lake_reservoir_mpi(ig) = lake_reservoir_mpi(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_mpi(sum_lake_overflow,total_lake_overflow)
+      CALL bcast_mpi(total_lake_overflow)
+
+      WHERE(coast_mask) coastalflow = coastalflow + total_lake_overflow/total_coast_points
 
     ENDIF ! is_omp_root
@@ -1135,4 +1166,5 @@
     CALL scatter_omp(coastalflow,coastalflow_omp)
     CALL scatter_omp(lakeinflow,lakeinflow_omp)
+    CALL scatter_omp(lake_reservoir_mpi,lake_reservoir)
     CALL scatter_omp(fast_diag_mpi,fast_diag)
     CALL scatter_omp(slow_diag_mpi,slow_diag)
@@ -1166,4 +1198,5 @@
     pond_diag(:) =  zero
     irrigation(:) = zero
+
 
 
@@ -1733,4 +1766,65 @@
     ENDIF
 
+    !
+    ! Write restart variables
+    !
+    tmp_day(1) = time_counter
+    IF (is_root_prc) CALL restput (rest_id, 'routingcounter', 1, 1, 1, kjit, tmp_day)
+
+    CALL restput_p (rest_id, 'routingarea', nbp_glo, nbasmax, 1, kjit, routing_area, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'routetogrid', nbp_glo, nbasmax, 1, kjit, REAL(route_togrid,r_std), 'scatter', &
+         nbp_glo, index_g)
+    CALL restput_p (rest_id, 'routetobasin', nbp_glo, nbasmax, 1, kjit, REAL(route_tobasin,r_std), 'scatter', &
+         nbp_glo, index_g)
+    CALL restput_p (rest_id, 'basinid', nbp_glo, nbasmax, 1, kjit, REAL(global_basinid,r_std), 'scatter', &
+         nbp_glo, index_g)
+    CALL restput_p (rest_id, 'topoindex', nbp_glo, nbasmax, 1, kjit, topo_resid, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'fastres', nbp_glo, nbasmax, 1, kjit, fast_reservoir, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'slowres', nbp_glo, nbasmax, 1, kjit, slow_reservoir, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'streamres', nbp_glo, nbasmax, 1, kjit, stream_reservoir, 'scatter',nbp_glo,index_g)
+    CALL restput_p (rest_id, 'floodres', nbp_glo, nbasmax, 1, kjit, flood_reservoir, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'floodh', nbp_glo, 1, 1, kjit, flood_height, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'flood_frac_bas', nbp_glo, nbasmax, 1, kjit, flood_frac_bas, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'pond_frac', nbp_glo, 1, 1, kjit, pond_frac, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'flood_frac', nbp_glo, 1, 1, kjit, flood_frac, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'flood_res', nbp_glo, 1, 1, kjit, flood_res, 'scatter', nbp_glo, index_g)
+
+    CALL restput_p (rest_id, 'lakeres', nbp_glo, 1, 1, kjit, lake_reservoir, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'pondres', nbp_glo, 1, 1, kjit, pond_reservoir, 'scatter',  nbp_glo, index_g)
+
+    CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'returnflow', nbp_glo, 1, 1, kjit, returnflow_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'reinfiltration', nbp_glo, 1, 1, kjit, reinfiltration_mean, 'scatter',  nbp_glo, index_g)
+    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, 'hydrographs', nbp_glo, 1, 1, kjit, hydrographs, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'slowflow_diag', nbp_glo, 1, 1, kjit, slowflow_diag, 'scatter',  nbp_glo, index_g)
+    !
+    ! Keep track of the accumulated variables
+    !
+    CALL restput_p (rest_id, 'floodout_route', nbp_glo, 1, 1, kjit, floodout_mean, 'scatter',  nbp_glo, index_g)
+    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)
+    CALL restput_p (rest_id, 'transpot_route', nbp_glo, 1, 1, kjit, transpot_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'precip_route', nbp_glo, 1, 1, kjit, precip_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'humrel_route', nbp_glo, 1, 1, kjit, humrel_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'totnobio_route', nbp_glo, 1, 1, kjit, totnobio_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'k_litt_route', nbp_glo, 1, 1, kjit, k_litt_mean, 'scatter',  nbp_glo, index_g)
+    CALL restput_p (rest_id, 'vegtot_route', nbp_glo, 1, 1, kjit, vegtot_mean, 'scatter',  nbp_glo, index_g)
+
+    IF ( do_irrigation ) THEN
+       CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter',  nbp_glo, index_g)
+       CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter',  nbp_glo, index_g)
+    ENDIF
+
+    IF ( do_floodplains ) THEN
+       CALL restput_p (rest_id, 'floodplains', nbp_glo, 1, 1, kjit, floodplains, 'scatter',  nbp_glo, index_g)
+    ENDIF
+    IF ( doswamps ) THEN
+       CALL restput_p (rest_id, 'swamp', nbp_glo, 1, 1, kjit, swamp, 'scatter',  nbp_glo, index_g)
+    ENDIF
+
+
+
   END SUBROUTINE routing_simple_finalize
 
@@ -1950,4 +2044,12 @@
     CALL getin_p("POND_CRI", pondcri)
 
+    !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)
 
     ! In order to simplify the time cascade check that dt_routing
