Index: branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90
===================================================================
--- branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90	(revision 8320)
+++ branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90	(revision 8365)
@@ -24,9 +24,9 @@
    
   !! PARAMETERS
-  REAL(r_std), SAVE                         :: fast_tcst = 3.0             !! Property of the fast reservoir (day/m)
+  REAL(r_std), SAVE                         :: fast_tcst = 9e-3             !! Property of the fast reservoir (day/km)
   !$OMP THREADPRIVATE(fast_tcst)
-  REAL(r_std), SAVE                         :: slow_tcst = 25.0            !! Property of the slow reservoir (day/m)
+  REAL(r_std), SAVE                         :: slow_tcst = 1.2e-2           !! Property of the slow reservoir (day/km)
   !$OMP THREADPRIVATE(slow_tcst)
-  REAL(r_std), SAVE                         :: stream_tcst = 0.24          !! Property of the stream reservoir (day/m)
+  REAL(r_std), SAVE                         :: stream_tcst = 3e-5           !! Property of the stream reservoir (day/km)
   !$OMP THREADPRIVATE(stream_tcst)
   
@@ -41,5 +41,5 @@
   !$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)
+  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)
 
@@ -47,5 +47,5 @@
   !$OMP THREADPRIVATE(stream_reservoir_r) 
   
-  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:)         !! Water amount in the stream reservoir (kg) - (local routing grid)
   !$OMP THREADPRIVATE(riverflow_mean) 
 
@@ -53,5 +53,5 @@
   !$OMP THREADPRIVATE(coastalflow_mean) 
 
-  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid)
+  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:)        !! Water amount in the stream reservoir (kg) - (local routing grid)
   !$OMP THREADPRIVATE(lakeinflow_mean) 
 
@@ -61,4 +61,6 @@
   REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:)                !! Topographic index of the retention time (m) index - (local routing grid)
   !$OMP THREADPRIVATE(topoind_r)
+  REAL(r_std),SAVE             :: topoind_factor              !! conversion factor for topographic index (merit : m, standard : km), topo index must be in m
+  !$OMP THREADPRIVATE(topoind_factor)
   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)
@@ -69,5 +71,5 @@
   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) 
+  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) 
@@ -85,5 +87,5 @@
   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)
+  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)
@@ -253,8 +255,8 @@
     INTEGER ,INTENT(IN)                     :: kjit           
     INTEGER ,INTENT(IN)                     :: rest_id           
-    INTEGER, INTENT(IN)                     :: nbpt_             !! nb points orchidee grid
+    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(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)
@@ -311,4 +313,7 @@
     LOGICAL     :: coastline_glo(nbp_glo)
     LOGICAL     :: coastline2D_glo(iim_g*jjm_g)
+    LOGICAL     :: mask1d(nbpt)
+    LOGICAL     :: mask1d_glo(nbp_glo)
+    LOGICAL     :: mask2d_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
@@ -316,19 +321,26 @@
     
     is_periodic = global
+    mask1d=.TRUE. ;
+    CALL gather(mask1d,mask1d_glo)
     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
-      
+      contfrac2D_glo(:)=0
+      mask2d_glo=.FALSE.
+      DO i=1,nbp_glo
+        contfrac2D_glo(index_g(i)) = contfrac_glo(i)
+        mask2d_glo(index_g(i))=mask1d_glo(i)
+      ENDDO
+
       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
+          IF (grid_type==unstructured) THEN
+            IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN
+              coastline2D_glo(ij)=.TRUE.
+            ENDIF
+          ELSE
             DO k=-1,1
               DO m=-1,1
@@ -358,6 +370,7 @@
                   CYCLE
                 ENDIF
-                next_ij =  (next_j-1)*iim_g+next_i    
-                IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE.
+                next_ij =  (next_j-1)*iim_g+next_i  
+                IF (.NOT. mask2d_glo(next_ij)) coastline2D_glo(ij)=.TRUE.
+    !            IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE.
               ENDDO
             ENDDO
@@ -376,9 +389,8 @@
   END SUBROUTINE compute_coastline
 
-  
-     
-  SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r)
+ 
+ SUBROUTINE new_routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r)
   USE xios
-  USE grid, ONLY : global
+!  USE grid, ONLY : global
   IMPLICIT NONE
     INCLUDE "mpif.h"
@@ -391,4 +403,5 @@
     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)                :: mask(nbp_mpi)
     REAL(r_std)                :: contfrac_r(ni,nj)
     REAL(r_std)                :: mask_coast(nbp_mpi)
@@ -398,13 +411,20 @@
     REAL(r_std)                :: state_r(ni,nj)
     REAL(r_std)                :: state_rp1(0:ni+1,0:nj+1)
-    
+    REAL(r_std)                :: done_rp1(0:ni+1,0:nj+1)
+    REAL(r_std)                :: trip_out_r(0:ni+1,0:nj+1)
+
     INTEGER, PARAMETER ::       is_ter=0
     INTEGER, PARAMETER ::       is_coast=1
     INTEGER, PARAMETER ::       is_oce=2
+    INTEGER, PARAMETER ::       riverflow=99
+    INTEGER, PARAMETER ::       coastalflow=98
+    INTEGER, PARAMETER ::       lakeinflow=97
+    INTEGER, PARAMETER ::       norouting=100
+
     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
+    REAL(r_std),PARAMETER :: epsilon=1e-5
  
     TYPE(xios_duration) :: ts 
@@ -419,5 +439,423 @@
     CHARACTER(LEN=20)  :: calendar_type
     LOGICAL :: true=.TRUE.
-
+    LOGICAL :: global=.FALSE.
+    LOGICAL :: ok
+
+      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) = coastalflow
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
+          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) = coastalflow
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
+          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) = coastalflow
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
+          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) = coastalflow
+            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow
+          ENDDO
+        ENDIF
+      ENDIF
+
+      mask = 1
+      CALL xios_send_field("routing_contfrac", mask) 
+      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 (frac_coast_r(i,j)>epsilon) THEN
+            state_r(i,j)=is_coast
+          ELSE 
+            IF (contfrac_r(i,j)> epsilon) THEN 
+             state_r(i,j)=is_ter
+            ELSE
+              state_r(i,j)=is_oce
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO
+
+
+      DO j=1,nj
+        DO i=1,ni
+          IF (trip_r(i,j)>=100) THEN
+            trip_r(i,j)=norouting
+          ENDIF
+        ENDDO
+      ENDDO          
+
+      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)
+        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)
+        
+        ! first => riverflow and coastalflow routing to ocean must be routed to coast => to long 
+        DO j=1,nj
+          DO i=1,ni
+            CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj, ok)
+            IF (ok) THEN
+              IF (state_rp1(nexti,nextj)==is_oce .AND. (trip_rp1(nexti,nextj)==riverflow .OR. trip_rp1(nexti,nextj)==coastalflow .OR. trip_rp1(nexti,nextj)==lakeinflow)) THEN
+                trip_r(i,j) = trip_rp1(nexti,nextj)
+                updated=updated + 1
+              ENDIF
+            ENDIF
+          ENDDO
+        ENDDO
+
+        DO j=1,nj
+          DO i=1,ni
+            IF (state_rp1(i,j)==is_oce) THEN
+              IF (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow) THEN
+                trip_r(i,j) = norouting
+                updated=updated+1
+              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
+
+      ! second prolongate riverflow and coastalflow to coast if they are in land 
+      updated=1
+      DO WHILE(updated>0)
+        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=1,nj
+          DO i=1,ni
+            IF (state_rp1(i,j)==is_ter .AND. (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow)) THEN
+              CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj, ok)
+              IF (ok) THEN
+                IF (state_rp1(nexti,nextj)==is_ter .OR. state_rp1(nexti,nextj)==is_coast) THEN
+                  trip_r(i,j) = trip_extended_rp1(i,j)
+                  topoind_r(i,j) = 1e-10  ! flow instantanesly for now but better to take the topind of the incomming flux
+                  updated=updated+1
+                ENDIF
+              ENDIF
+            ENDIF
+          ENDDO
+        ENDDO
+
+        DO j=1,nj
+          DO i=1,ni
+            IF (state_rp1(i,j)==is_ter .OR. state_rp1(i,j)==is_coast) THEN
+              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_extended_rp1, ni, nj, previ, prevj, nexti, nextj, ok)
+                  IF (ok) THEN
+                    IF (nexti==i .AND. nextj==j) THEN
+                      IF (state_rp1(previ,prevj)==is_ter .AND. (trip_rp1(previ,prevj)==riverflow .OR. trip_rp1(previ,prevj)==coastalflow .OR. trip_rp1(previ,prevj)==lakeinflow)) THEN
+                        IF (trip_r(i,j)/=riverflow) THEN
+                          trip_r(i,j)=trip_rp1(previ,prevj)
+                        ELSE 
+                          trip_r(i,j)=riverflow
+                        ENDIF
+                        updated=updated+1
+                      ENDIF
+                    ENDIF
+                  ENDIF
+                ENDDO
+              ENDDO
+            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_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)      
+      DO j=1,nj
+        DO i=1,ni
+          IF (state_rp1(i,j)==is_coast) THEN
+            IF (trip_r(i,j)==norouting) trip_r(i,j) = coastalflow
+          ENDIF
+        ENDDO
+      ENDDO
+      
+!      it = it +1
+!      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)
+!
+!      DO j=1,nj
+!        DO i=1,ni
+!          IF (state_rp1(i,j)==is_oce) THEN
+!            trip_r(i,j) = norouting
+!          ENDIF
+!        ENDDO
+!      ENDDO
+      
+      it = it +1
+      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)
+! check
+      updated=0
+      DO j=1,nj
+        DO i=1,ni
+          IF (state_rp1(i,j)==is_ter) THEN
+            IF (trip_r(i,j) == coastalflow .OR. trip_r(i,j) == riverflow .OR. trip_r(i,j) == norouting) THEN
+              updated=updated+1
+              trip_r(i,j) = lakeinflow
+            ENDIF
+          ELSE IF (state_rp1(i,j)==is_coast) THEN
+            IF (trip_r(i,j) == norouting) THEN
+              updated=updated+1
+              trip_r(i,j) = coastalflow
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO
+
+
+      CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
+      IF (updated/=0 .AND. is_mpi_root) PRINT*,"WARNING : ",updated," riverflow or coastal flow are not on the coast, ie in middle of land ==> converted into lakeinflow"
+      CALL MPI_BARRIER(MPI_COMM_ORCH, ierr)
+
+      updated=0
+      DO j=1,nj
+        DO i=1,ni
+          CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok)
+          IF (ok) THEN
+            IF (trip_rp1(nexti,nextj)==norouting) THEN
+              updated=updated+1
+              PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to nothing => ",ibegin+nexti,jbegin+nextj, &
+                     " (trip(i,j)=",trip_rp1(nexti,nextj),")"
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO
+      
+      CALL MPI_BARRIER(MPI_COMM_ORCH, ierr)
+      CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
+      IF (updated/=0) THEN
+        IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are not correctly routed (routed to nothing)"
+        CALL xios_context_finalize
+        STOP
+      ENDIF
+
+
+      updated=0
+      DO j=1,nj
+        DO i=1,ni
+          CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok)
+          IF (ok) THEN
+            IF (state_rp1(nexti,nextj)==is_oce .AND. (trip_rp1(nexti,nextj)==riverflow .OR. trip_rp1(nexti,nextj)==coastalflow .OR. trip_rp1(nexti,nextj)==lakeinflow )) THEN
+              updated=updated+1
+              PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to coastalflow/riverflow/lakeinflow  point that is not on land grid => ",ibegin+nexti,jbegin+nextj, &
+                     " (trip(i,j)=",trip_rp1(nexti,nextj),")"
+            ENDIF
+          ENDIF
+        ENDDO
+      ENDDO
+      
+      CALL MPI_BARRIER(MPI_COMM_ORCH, ierr)
+      CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr)
+      IF (updated/=0) THEN
+        IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are routed to coastalflow/riverflow/lakeinflow points that are not on land grid"
+        CALL xios_context_finalize
+        STOP
+      ENDIF
+
+
+
+      it = it +1
+      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)
+
+      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, ok)
+    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
+      LOGICAL, INTENT(OUT) :: ok
+
+      ok=.TRUE. 
+      
+      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
+          ok=.FALSE.
+       END SELECT
+       
+     END SUBROUTINE next_trip
+     
+  END SUBROUTINE new_routing_flow_correct_riverflow
+
+
+
+     
+  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-5
+ 
+    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.
+    LOGICAL :: global=.FALSE.
 
       CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo)    ! get routing domain dimension
@@ -493,54 +931,4 @@
         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)
@@ -831,5 +1219,5 @@
     LOGICAL :: file_exists
     REAL(r_std) :: epsilon = 1e-5
-
+    CHARACTER(LEN=255) :: routing_file_type
 !_ ================================================================================================================================
     !
@@ -879,4 +1267,18 @@
     !Config Units = [days]
 
+    CALL getin_p('SLOW_TCST', slow_tcst)
+
+    routing_file_type = "standard"
+    CALL getin_p("routing_file_type", routing_file_type) 
+
+    IF (TRIM(routing_file_type)=="standard") THEN
+      topoind_factor=1000
+    ELSE IF (TRIM(routing_file_type)=="merit") THEN
+      topoind_factor=1
+    ELSE
+      CALL ipslerr_p(3,'routing_flow_init_local', &
+        'Error in call routing_flow_init_local','getin routing_file_type : bad value, must be <standard> or <merit>','')
+    ENDIF
+
     split_routing=1
     CALL getin_p("SPLIT_ROUTING",split_routing)
@@ -892,12 +1294,7 @@
 
     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)
-       
+       WHERE(contfrac_mpi <= epsilon) contfrac_mpi=0
+       WHERE(contfrac_mpi >= 1-epsilon) contfrac_mpi=1
+
        CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)    ! get routing domain dimension
 
@@ -933,4 +1330,17 @@
 
        ALLOCATE(trip_extended_r(ni,nj))
+
+      ! correction on coastline in case of LAM, the cells near the frontier are considered as coastline 
+      ! => interpolate model grid to routing grid and reinterpolate to model grid. Fractionnal cells will be considered as coastline too
+    !  diag=1
+    !  CALL xios_send_field("one", diag)
+    !  CALL xios_recv_field("tmp1_coastline", diag_r)
+    !  CALL xios_send_field("tmp2_coastline", diag_r)
+    !  CALL xios_recv_field("lam_coastline", diag)
+    !  WHERE(diag<1-epsilon) is_coastline=.TRUE.
+
+       diag=0
+       WHERE(is_coastline) diag=1
+       CALL xios_send_field("is_coastline", diag)
 
        is_lakeinflow_r(:) = .FALSE.
@@ -960,5 +1370,6 @@
        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 routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r) 
+       CALL new_routing_flow_correct_riverflow(ni, nj, contfrac_mpi, is_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
@@ -1041,4 +1452,5 @@
        CALL xios_send_field("routing_mask_r", diag_r)
        CALL xios_recv_field("frac_routing", diag)
+       WHERE (diag < epsilon) diag=0
        routing_weight_in (:) = 0
        unrouted_weight (:) = 0
@@ -1049,5 +1461,5 @@
       ! compute the number of coast cells to distribute lakeinflow onto
        diag=0
-       WHERE (coastline_mpi) diag=1
+       WHERE (is_coastline) diag=1
        nb_coast_points=SUM(diag)
        CALL reduce_sum_mpi(nb_coast_points, total_coast_points)
@@ -1067,8 +1479,8 @@
 
        diag(:) = 0
-       WHERE (coastline_mpi) diag=1
+       WHERE (is_coastline) diag=1
        CALL xios_send_field("mask_coastal",diag)
        CALL xios_recv_field("frac_coastal_r",diag_r)
-       WHERE (diag_r > 100.) ! missing_value
+       WHERE (diag_r > 100) ! missing_value
          diag_r=0
        ELSE WHERE (diag_r < epsilon) 
@@ -1080,8 +1492,8 @@
 
        diag(:) = 0
-       WHERE (.NOT. coastline_mpi) diag=1
+       WHERE (.NOT. is_coastline) diag=1
        CALL xios_send_field("mask_lake",diag)
        CALL xios_recv_field("frac_lake_r",diag_r)
-       WHERE (diag_r > 100.) ! missing_value
+       WHERE (diag_r > 100) ! missing_value
          diag_r=0.
        ELSEWHERE (diag_r < epsilon) 
@@ -1105,7 +1517,8 @@
              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))
+           ELSE IF (frac_coast_r(ij)>0) THEN
+             weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij))
              weight_coast_to_lake_r(ij) = 0
+           ELSE
            ENDIF
          ELSE IF (is_lakeinflow_r(ij)) THEN
@@ -1117,5 +1530,5 @@
              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_lake_r(ij) = 1./(frac_lake_r(ij))
              weight_lake_to_coast_r(ij) = 0
            ENDIF
@@ -1432,5 +1845,5 @@
       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*(routing_weight-unrouted_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
@@ -1441,5 +1854,5 @@
       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*(routing_weight-unrouted_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
@@ -1451,6 +1864,6 @@
           
 
-      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)
+      CALL MPI_ALLREDUCE(sum((runoff+drainage)*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)
+      CALL MPI_ALLREDUCE(sum(runoff_r+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,  &
@@ -1470,13 +1883,13 @@
           IF ( routing_mask_r(ig) ) THEN
 
-            flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),&
+            flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/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)),&
+            flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/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)),&
+            flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*stream_tcst*one_day/(dt_routing/split_routing)),&
                    & stream_reservoir_r(ig)-min_sechiba)
             stream_flow_r(ig) = MAX(flow, zero)
@@ -1641,10 +2054,10 @@
       ENDIF   
        
-      coastalflow=coastalflow+flow_coast
+      coastalflow=coastalflow+flow_coast+(runoff+drainage)*unrouted_weight
       lakeinflow=lakeinflow+flow_lake
 
 
-      WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in
-      WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in
+!      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(:))
