Ignore:
Timestamp:
2024-01-08T15:20:57+01:00 (12 months ago)
Author:
yann.meurdesoif
Message:

Update of routing native

  • More robust algorithm for correcting routing grid to adjust to coast line model, in global and LAM configuration
  • More robust water conservation
  • Reservoir constant are now given in km/day, independently of the routing grid topoindex unit. So a addtionnal parameter is added to the run.def to specify the type of routing file : routing_file_type = standard => 50km legacy routing file routing_file_type = merit => 2km Merit routing file

YM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90

    r8300 r8365  
    2424    
    2525  !! PARAMETERS 
    26   REAL(r_std), SAVE                         :: fast_tcst = 3.0             !! Property of the fast reservoir (day/m) 
     26  REAL(r_std), SAVE                         :: fast_tcst = 9e-3             !! Property of the fast reservoir (day/km) 
    2727  !$OMP THREADPRIVATE(fast_tcst) 
    28   REAL(r_std), SAVE                         :: slow_tcst = 25.0            !! Property of the slow reservoir (day/m) 
     28  REAL(r_std), SAVE                         :: slow_tcst = 1.2e-2           !! Property of the slow reservoir (day/km) 
    2929  !$OMP THREADPRIVATE(slow_tcst) 
    30   REAL(r_std), SAVE                         :: stream_tcst = 0.24          !! Property of the stream reservoir (day/m) 
     30  REAL(r_std), SAVE                         :: stream_tcst = 3e-5           !! Property of the stream reservoir (day/km) 
    3131  !$OMP THREADPRIVATE(stream_tcst) 
    3232   
    … …  
    4141  !$OMP THREADPRIVATE(fast_reservoir_r) 
    4242 
    43   REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:)         !! Water amount in the slow reservoir (kg) - (local routing grid) 
     43  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:)          !! Water amount in the slow reservoir (kg) - (local routing grid) 
    4444  !$OMP THREADPRIVATE(slow_reservoir_r) 
    4545 
    … …  
    4747  !$OMP THREADPRIVATE(stream_reservoir_r)  
    4848   
    49   REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid) 
     49  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:)         !! Water amount in the stream reservoir (kg) - (local routing grid) 
    5050  !$OMP THREADPRIVATE(riverflow_mean)  
    5151 
    … …  
    5353  !$OMP THREADPRIVATE(coastalflow_mean)  
    5454 
    55   REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:)       !! Water amount in the stream reservoir (kg) - (local routing grid) 
     55  REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:)        !! Water amount in the stream reservoir (kg) - (local routing grid) 
    5656  !$OMP THREADPRIVATE(lakeinflow_mean)  
    5757 
    … …  
    6161  REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:)                !! Topographic index of the retention time (m) index - (local routing grid) 
    6262  !$OMP THREADPRIVATE(topoind_r) 
     63  REAL(r_std),SAVE             :: topoind_factor              !! conversion factor for topographic index (merit : m, standard : km), topo index must be in m 
     64  !$OMP THREADPRIVATE(topoind_factor) 
    6365  INTEGER,SAVE,ALLOCATABLE     :: route_flow_rp1(:)           !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo)    
    6466  !$OMP THREADPRIVATE(route_flow_rp1) 
    … …  
    6971  LOGICAL,SAVE,ALLOCATABLE     :: is_riverflow_r(:)           !! is river flow point - (local routing grid)  
    7072  !$OMP THREADPRIVATE(is_riverflow_r) 
    71   LOGICAL,SAVE,ALLOCATABLE     :: is_coastline(:)           !! is coastline point - (local native grid)  
     73  LOGICAL,SAVE,ALLOCATABLE     :: is_coastline(:)             !! is coastline point - (local native grid)  
    7274  !$OMP THREADPRIVATE(is_coastline) 
    7375  LOGICAL,SAVE,ALLOCATABLE     :: is_streamflow_r(:)          !! is stream flow point - (local routing grid)  
    … …  
    8587  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) 
    8688  !$OMP THREADPRIVATE(routing_weight) 
    87   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) 
     89  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) 
    8890  !$OMP THREADPRIVATE(routing_weight_in) 
    8991  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) 
    … …  
    253255    INTEGER ,INTENT(IN)                     :: kjit            
    254256    INTEGER ,INTENT(IN)                     :: rest_id            
    255     INTEGER, INTENT(IN)                     :: nbpt_             !! nb points orchidee grid 
     257    INTEGER, INTENT(IN)                     :: nbpt_                 !! nb points orchidee grid 
    256258    REAL(r_std), INTENT(IN)                 :: dt_routing                        
    257     REAL(r_std), INTENT(IN)                 :: contfrac(nbpt_)   !! fraction of land 
    258     INTEGER, INTENT(OUT)                    :: nbpt_r_             !! nb points routing native grid 
     259    REAL(r_std), INTENT(IN)                 :: contfrac(nbpt_)       !! fraction of land 
     260    INTEGER, INTENT(OUT)                    :: nbpt_r_               !! nb points routing native grid 
    259261    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) 
    260262    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) 
    … …  
    311313    LOGICAL     :: coastline_glo(nbp_glo) 
    312314    LOGICAL     :: coastline2D_glo(iim_g*jjm_g) 
     315    LOGICAL     :: mask1d(nbpt) 
     316    LOGICAL     :: mask1d_glo(nbp_glo) 
     317    LOGICAL     :: mask2d_glo(iim_g*jjm_g) 
    313318    INTEGER(i_std) :: i,j,ij,next_i,next_j,next_ij,m,k 
    314319    REAL(r_std), PARAMETER :: epsilon=1e-5 
    … …  
    316321     
    317322    is_periodic = global 
     323    mask1d=.TRUE. ; 
     324    CALL gather(mask1d,mask1d_glo) 
    318325    CALL gather(contfrac,contfrac_glo) 
    319     contfrac2D_glo(:)=0 
    320     DO i=1,nbp_glo 
    321       contfrac2D_glo(index_g(i)) = contfrac_glo(i) 
    322     ENDDO 
     326     
    323327     
    324328    IF (is_mpi_root .AND. is_omp_root) THEN 
    325        
     329      contfrac2D_glo(:)=0 
     330      mask2d_glo=.FALSE. 
     331      DO i=1,nbp_glo 
     332        contfrac2D_glo(index_g(i)) = contfrac_glo(i) 
     333        mask2d_glo(index_g(i))=mask1d_glo(i) 
     334      ENDDO 
     335 
    326336      coastline2D_glo(:)=.FALSE. 
    327337      DO j=1,jjm_g 
    328338        DO i=1,iim_g 
    329339          ij=(j-1)*iim_g+i 
    330           IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN 
    331             coastline2D_glo(ij)=.TRUE. 
    332           ELSE IF (contfrac2D_glo(ij)>=1-epsilon .AND. grid_type/=unstructured ) THEN 
     340          IF (grid_type==unstructured) THEN 
     341            IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN 
     342              coastline2D_glo(ij)=.TRUE. 
     343            ENDIF 
     344          ELSE 
    333345            DO k=-1,1 
    334346              DO m=-1,1 
    … …  
    358370                  CYCLE 
    359371                ENDIF 
    360                 next_ij =  (next_j-1)*iim_g+next_i     
    361                 IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE. 
     372                next_ij =  (next_j-1)*iim_g+next_i   
     373                IF (.NOT. mask2d_glo(next_ij)) coastline2D_glo(ij)=.TRUE. 
     374    !            IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE. 
    362375              ENDDO 
    363376            ENDDO 
    … …  
    376389  END SUBROUTINE compute_coastline 
    377390 
    378    
    379       
    380   SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) 
     391  
     392 SUBROUTINE new_routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) 
    381393  USE xios 
    382   USE grid, ONLY : global 
     394!  USE grid, ONLY : global 
    383395  IMPLICIT NONE 
    384396    INCLUDE "mpif.h" 
    … …  
    391403    REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj)        ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean 
    392404     
     405    REAL(r_std)                :: mask(nbp_mpi) 
    393406    REAL(r_std)                :: contfrac_r(ni,nj) 
    394407    REAL(r_std)                :: mask_coast(nbp_mpi) 
    … …  
    398411    REAL(r_std)                :: state_r(ni,nj) 
    399412    REAL(r_std)                :: state_rp1(0:ni+1,0:nj+1) 
    400      
     413    REAL(r_std)                :: done_rp1(0:ni+1,0:nj+1) 
     414    REAL(r_std)                :: trip_out_r(0:ni+1,0:nj+1) 
     415 
    401416    INTEGER, PARAMETER ::       is_ter=0 
    402417    INTEGER, PARAMETER ::       is_coast=1 
    403418    INTEGER, PARAMETER ::       is_oce=2 
     419    INTEGER, PARAMETER ::       riverflow=99 
     420    INTEGER, PARAMETER ::       coastalflow=98 
     421    INTEGER, PARAMETER ::       lakeinflow=97 
     422    INTEGER, PARAMETER ::       norouting=100 
     423 
    404424    INTEGER            :: updated 
    405425    INTEGER            :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr 
    406426    INTEGER            :: ibegin, jbegin, ni_glo, nj_glo 
    407427    REAL(r_std)        :: lonvalue_1d(ni), latvalue_1d(nj) 
    408     REAL(r_std),PARAMETER :: epsilon=1e-3 
     428    REAL(r_std),PARAMETER :: epsilon=1e-5 
    409429  
    410430    TYPE(xios_duration) :: ts  
    … …  
    419439    CHARACTER(LEN=20)  :: calendar_type 
    420440    LOGICAL :: true=.TRUE. 
    421  
     441    LOGICAL :: global=.FALSE. 
     442    LOGICAL :: ok 
     443 
     444      CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo)    ! get routing domain dimension 
     445      CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d)    ! get routing domain dimension 
     446       
     447      ! manage non periodic boundary in case of LAM       
     448      IF (.NOT. global) THEN 
     449        IF (ibegin==0) THEN 
     450          i=1 ; 
     451          DO j=1,nj 
     452            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 
     453            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 
     454          ENDDO 
     455        ENDIF 
     456 
     457        IF (ibegin+ni==ni_glo) THEN 
     458          i=ni ; 
     459          DO j=1,nj 
     460            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 
     461            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 
     462          ENDDO 
     463        ENDIF 
     464         
     465        IF (jbegin==0) THEN 
     466          j=1 ; 
     467          DO i=1,ni 
     468            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 
     469            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 
     470          ENDDO 
     471        ENDIF 
     472 
     473        IF (jbegin+nj==nj_glo) THEN 
     474          j=nj ; 
     475          DO i=1,ni 
     476            IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 
     477            IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 
     478          ENDDO 
     479        ENDIF 
     480      ENDIF 
     481 
     482      mask = 1 
     483      CALL xios_send_field("routing_contfrac", mask)  
     484      CALL xios_recv_field("routing_contfrac_r", contfrac_r) 
     485 
     486      CALL xios_send_field("trip_ext_r",trip_extended_r) 
     487      CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1) 
     488 
     489      mask_coast=0 
     490      WHERE(coastline) mask_coast=1 
     491      CALL xios_send_field("mask_coastline",mask_coast) 
     492      CALL xios_recv_field("frac_coastline_r",frac_coast_r) 
     493       
     494      DO j=1,nj 
     495        DO i=1,ni 
     496          IF (frac_coast_r(i,j)>epsilon) THEN 
     497            state_r(i,j)=is_coast 
     498          ELSE  
     499            IF (contfrac_r(i,j)> epsilon) THEN  
     500             state_r(i,j)=is_ter 
     501            ELSE 
     502              state_r(i,j)=is_oce 
     503            ENDIF 
     504          ENDIF 
     505        ENDDO 
     506      ENDDO 
     507 
     508 
     509      DO j=1,nj 
     510        DO i=1,ni 
     511          IF (trip_r(i,j)>=100) THEN 
     512            trip_r(i,j)=norouting 
     513          ENDIF 
     514        ENDDO 
     515      ENDDO           
     516 
     517      CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH) 
     518      CALL xios_orchidee_change_context("orchidee_init_routing") 
     519      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 
     520      CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d)                           ! get routing domain dimension 
     521      CALL xios_close_context_definition() 
     522      CALL xios_update_calendar(1) 
     523 
     524      CALL xios_send_field("trip_r_init",trip_r) 
     525       
     526      CALL xios_send_field("trip_r",trip_r) 
     527      CALL xios_recv_field("trip_rp1",trip_rp1) 
     528      
     529      CALL xios_send_field("state_r",state_r) 
     530      CALL xios_recv_field("state_rp1",state_rp1) 
     531       
     532!      DO j=1,nj 
     533!        DO i=1,ni 
     534!          IF (trip_rp1(i,j)<=8) THEN 
     535!            CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj) 
     536!            IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow 
     537!          ENDIF 
     538!        ENDDO 
     539!      ENDDO           
     540 
     541      state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 
     542      trip_r(1:ni,1:nj)  = trip_rp1(1:ni,1:nj) 
     543 
     544      updated=1 
     545      it=2 
     546      DO WHILE(updated>0) 
     547        updated=0 
     548 
     549        CALL xios_update_calendar(it) 
     550         
     551        CALL xios_send_field("trip_r", trip_r) 
     552        CALL xios_recv_field("trip_rp1",trip_rp1) 
     553      
     554        CALL xios_send_field("state_r",state_r) 
     555        CALL xios_recv_field("state_rp1",state_rp1) 
     556 
     557        state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 
     558        trip_r(1:ni,1:nj)  = trip_rp1(1:ni,1:nj) 
     559         
     560        ! first => riverflow and coastalflow routing to ocean must be routed to coast => to long  
     561        DO j=1,nj 
     562          DO i=1,ni 
     563            CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj, ok) 
     564            IF (ok) THEN 
     565              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 
     566                trip_r(i,j) = trip_rp1(nexti,nextj) 
     567                updated=updated + 1 
     568              ENDIF 
     569            ENDIF 
     570          ENDDO 
     571        ENDDO 
     572 
     573        DO j=1,nj 
     574          DO i=1,ni 
     575            IF (state_rp1(i,j)==is_oce) THEN 
     576              IF (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow) THEN 
     577                trip_r(i,j) = norouting 
     578                updated=updated+1 
     579              ENDIF  
     580            ENDIF 
     581          ENDDO 
     582        ENDDO      
     583 
     584        CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 
     585        it = it +1 
     586      ENDDO 
     587 
     588      ! second prolongate riverflow and coastalflow to coast if they are in land  
     589      updated=1 
     590      DO WHILE(updated>0) 
     591        updated=0 
     592 
     593        CALL xios_update_calendar(it) 
     594         
     595        CALL xios_send_field("trip_r", trip_r) 
     596        CALL xios_recv_field("trip_rp1",trip_rp1) 
     597      
     598        CALL xios_send_field("state_r",state_r) 
     599        CALL xios_recv_field("state_rp1",state_rp1) 
     600 
     601        state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 
     602        trip_r(1:ni,1:nj)  = trip_rp1(1:ni,1:nj) 
     603         
     604        DO j=1,nj 
     605          DO i=1,ni 
     606            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 
     607              CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj, ok) 
     608              IF (ok) THEN 
     609                IF (state_rp1(nexti,nextj)==is_ter .OR. state_rp1(nexti,nextj)==is_coast) THEN 
     610                  trip_r(i,j) = trip_extended_rp1(i,j) 
     611                  topoind_r(i,j) = 1e-10  ! flow instantanesly for now but better to take the topind of the incomming flux 
     612                  updated=updated+1 
     613                ENDIF 
     614              ENDIF 
     615            ENDIF 
     616          ENDDO 
     617        ENDDO 
     618 
     619        DO j=1,nj 
     620          DO i=1,ni 
     621            IF (state_rp1(i,j)==is_ter .OR. state_rp1(i,j)==is_coast) THEN 
     622              DO k=-1,1 
     623                DO m=-1,1 
     624                  IF (k==0 .AND. m==0) CYCLE 
     625                  prevj = j+k ; previ = i+m 
     626                  CALL next_trip(trip_extended_rp1, ni, nj, previ, prevj, nexti, nextj, ok) 
     627                  IF (ok) THEN 
     628                    IF (nexti==i .AND. nextj==j) THEN 
     629                      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 
     630                        IF (trip_r(i,j)/=riverflow) THEN 
     631                          trip_r(i,j)=trip_rp1(previ,prevj) 
     632                        ELSE  
     633                          trip_r(i,j)=riverflow 
     634                        ENDIF 
     635                        updated=updated+1 
     636                      ENDIF 
     637                    ENDIF 
     638                  ENDIF 
     639                ENDDO 
     640              ENDDO 
     641            ENDIF 
     642          ENDDO 
     643        ENDDO 
     644         
     645        CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 
     646        it = it +1 
     647      ENDDO 
     648 
     649      CALL xios_update_calendar(it) 
     650         
     651      CALL xios_send_field("trip_r", trip_r) 
     652      CALL xios_recv_field("trip_rp1",trip_rp1) 
     653      
     654      CALL xios_send_field("state_r",state_r) 
     655      CALL xios_recv_field("state_rp1",state_rp1)       
     656      DO j=1,nj 
     657        DO i=1,ni 
     658          IF (state_rp1(i,j)==is_coast) THEN 
     659            IF (trip_r(i,j)==norouting) trip_r(i,j) = coastalflow 
     660          ENDIF 
     661        ENDDO 
     662      ENDDO 
     663       
     664!      it = it +1 
     665!      CALL xios_update_calendar(it) 
     666!         
     667!      CALL xios_send_field("trip_r", trip_r) 
     668!      CALL xios_recv_field("trip_rp1",trip_rp1) 
     669!      
     670!      CALL xios_send_field("state_r",state_r) 
     671!      CALL xios_recv_field("state_rp1",state_rp1) 
     672! 
     673!      DO j=1,nj 
     674!        DO i=1,ni 
     675!          IF (state_rp1(i,j)==is_oce) THEN 
     676!            trip_r(i,j) = norouting 
     677!          ENDIF 
     678!        ENDDO 
     679!      ENDDO 
     680       
     681      it = it +1 
     682      CALL xios_update_calendar(it) 
     683         
     684      CALL xios_send_field("trip_r", trip_r) 
     685      CALL xios_recv_field("trip_rp1",trip_rp1) 
     686      
     687      CALL xios_send_field("state_r",state_r) 
     688      CALL xios_recv_field("state_rp1",state_rp1) 
     689! check 
     690      updated=0 
     691      DO j=1,nj 
     692        DO i=1,ni 
     693          IF (state_rp1(i,j)==is_ter) THEN 
     694            IF (trip_r(i,j) == coastalflow .OR. trip_r(i,j) == riverflow .OR. trip_r(i,j) == norouting) THEN 
     695              updated=updated+1 
     696              trip_r(i,j) = lakeinflow 
     697            ENDIF 
     698          ELSE IF (state_rp1(i,j)==is_coast) THEN 
     699            IF (trip_r(i,j) == norouting) THEN 
     700              updated=updated+1 
     701              trip_r(i,j) = coastalflow 
     702            ENDIF 
     703          ENDIF 
     704        ENDDO 
     705      ENDDO 
     706 
     707 
     708      CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 
     709      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" 
     710      CALL MPI_BARRIER(MPI_COMM_ORCH, ierr) 
     711 
     712      updated=0 
     713      DO j=1,nj 
     714        DO i=1,ni 
     715          CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok) 
     716          IF (ok) THEN 
     717            IF (trip_rp1(nexti,nextj)==norouting) THEN 
     718              updated=updated+1 
     719              PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to nothing => ",ibegin+nexti,jbegin+nextj, & 
     720                     " (trip(i,j)=",trip_rp1(nexti,nextj),")" 
     721            ENDIF 
     722          ENDIF 
     723        ENDDO 
     724      ENDDO 
     725       
     726      CALL MPI_BARRIER(MPI_COMM_ORCH, ierr) 
     727      CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 
     728      IF (updated/=0) THEN 
     729        IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are not correctly routed (routed to nothing)" 
     730        CALL xios_context_finalize 
     731        STOP 
     732      ENDIF 
     733 
     734 
     735      updated=0 
     736      DO j=1,nj 
     737        DO i=1,ni 
     738          CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok) 
     739          IF (ok) THEN 
     740            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 
     741              updated=updated+1 
     742              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, & 
     743                     " (trip(i,j)=",trip_rp1(nexti,nextj),")" 
     744            ENDIF 
     745          ENDIF 
     746        ENDDO 
     747      ENDDO 
     748       
     749      CALL MPI_BARRIER(MPI_COMM_ORCH, ierr) 
     750      CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 
     751      IF (updated/=0) THEN 
     752        IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are routed to coastalflow/riverflow/lakeinflow points that are not on land grid" 
     753        CALL xios_context_finalize 
     754        STOP 
     755      ENDIF 
     756 
     757 
     758 
     759      it = it +1 
     760      CALL xios_update_calendar(it) 
     761         
     762      CALL xios_send_field("trip_r", trip_r) 
     763      CALL xios_recv_field("trip_rp1",trip_rp1) 
     764      
     765      CALL xios_send_field("state_r",state_r) 
     766      CALL xios_recv_field("state_rp1",state_rp1) 
     767      state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 
     768      trip_r(1:ni,1:nj)  = trip_rp1(1:ni,1:nj) 
     769 
     770      CALL xios_send_field("trip_r_final",trip_r) 
     771       
     772      CALL xios_context_finalize 
     773      CALL xios_orchidee_change_context("orchidee") 
     774     
     775  CONTAINS 
     776 
     777    SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj, ok) 
     778    IMPLICIT NONE 
     779      REAL, INTENT(IN)     :: trip(0:ni+1,0:nj+1) 
     780      INTEGER, INTENT(IN)  :: ni 
     781      INTEGER, INTENT(IN)  :: nj 
     782      INTEGER, INTENT(IN)  :: i 
     783      INTEGER, INTENT(IN)  :: j 
     784      INTEGER, INTENT(OUT) :: nexti 
     785      INTEGER, INTENT(OUT) :: nextj 
     786      LOGICAL, INTENT(OUT) :: ok 
     787 
     788      ok=.TRUE.  
     789       
     790      SELECT CASE(NINT(trip(i,j))) 
     791        CASE (1) 
     792           nextj=j-1 ; nexti=i    ! north 
     793        CASE (2) 
     794           nextj=j-1 ; nexti=i+1  ! north-east 
     795        CASE (3) 
     796           nextj=j ;   nexti=i+1  ! east 
     797        CASE (4) 
     798           nextj=j+1 ;  nexti=i+1 ! south-east 
     799        CASE (5) 
     800           nextj=j+1 ;  nexti=i   ! south 
     801        CASE(6) 
     802           nextj=j+1 ;  nexti=i-1 ! south-west 
     803        CASE (7) 
     804           nextj=j ;   nexti=i-1  ! west 
     805        CASE (8) 
     806           nextj=j-1 ;  nexti=i-1 ! north-west 
     807        CASE DEFAULT 
     808          ok=.FALSE. 
     809       END SELECT 
     810        
     811     END SUBROUTINE next_trip 
     812      
     813  END SUBROUTINE new_routing_flow_correct_riverflow 
     814 
     815 
     816 
     817      
     818  SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) 
     819  USE xios 
     820!  USE grid, ONLY : global 
     821  IMPLICIT NONE 
     822    INCLUDE "mpif.h" 
     823    INTEGER, INTENT(IN)        :: ni                      ! INPUT : size i (longitude) of the local routing domain 
     824    INTEGER, INTENT(IN)        :: nj                      ! INPUT : size i (latitude)  of the local routing domain 
     825    REAL(r_std), INTENT(IN)    :: contfrac(nbp_mpi)       ! INPUT : continental fraction (unittless) 
     826    LOGICAL, INTENT(IN)        :: coastline(nbp_mpi)      ! INPUT/OUTPUT : coastline mask 
     827    REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj)           ! INPUT/OUTPUT : diection of flow, which will be modified by the routine 
     828    REAL(r_std), INTENT(INOUT)    :: trip_extended_r(ni,nj)  ! INPUT  : direction of flow extended to ocean points 
     829    REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj)        ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean 
     830     
     831    REAL(r_std)                :: contfrac_r(ni,nj) 
     832    REAL(r_std)                :: mask_coast(nbp_mpi) 
     833    REAL(r_std)                :: frac_coast_r(ni,nj) 
     834    REAL(r_std)                :: trip_extended_rp1(0:ni+1,0:nj+1) 
     835    REAL(r_std)                :: trip_rp1(0:ni+1,0:nj+1) 
     836    REAL(r_std)                :: state_r(ni,nj) 
     837    REAL(r_std)                :: state_rp1(0:ni+1,0:nj+1) 
     838     
     839    INTEGER, PARAMETER ::       is_ter=0 
     840    INTEGER, PARAMETER ::       is_coast=1 
     841    INTEGER, PARAMETER ::       is_oce=2 
     842    INTEGER            :: updated 
     843    INTEGER            :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr 
     844    INTEGER            :: ibegin, jbegin, ni_glo, nj_glo 
     845    REAL(r_std)        :: lonvalue_1d(ni), latvalue_1d(nj) 
     846    REAL(r_std),PARAMETER :: epsilon=1e-5 
     847  
     848    TYPE(xios_duration) :: ts  
     849    TYPE(xios_domain) :: domain_hdl 
     850    TYPE(xios_domaingroup) :: domain_def_hdl 
     851    TYPE(xios_field) :: field_hdl 
     852    TYPE(xios_fieldgroup) :: field_def_hdl 
     853    TYPE(xios_file) :: file_hdl 
     854    TYPE(xios_filegroup) :: file_def_hdl 
     855    TYPE(xios_expand_domain) :: domain_expand_hdl 
     856    TYPE(xios_date) :: start_date, time_origin 
     857    CHARACTER(LEN=20)  :: calendar_type 
     858    LOGICAL :: true=.TRUE. 
     859    LOGICAL :: global=.FALSE. 
    422860 
    423861      CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo)    ! get routing domain dimension 
    … …  
    493931        ENDDO 
    494932      ENDDO           
    495  
    496 ! create a new XIOS context to solve iterative process of correction 
    497 ! so we can use filter to update halo 
    498 !      CALL xios_get_calendar_type(calendar_type) 
    499 !      CALL xios_get_start_date(start_date) 
    500 !      CALL xios_get_time_origin(time_origin) 
    501 ! 
    502 !      CALL xios_context_initialize("orchidee_routing_correction",MPI_COMM_ORCH) 
    503 !      CALL xios_orchidee_change_context("orchidee_routing_correction") 
    504 !     ts.second=1 
    505 !     calendar_type="gregorian"  
    506 !     CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts ) 
    507  
    508 !     CALL xios_get_handle("domain_definition",domain_def_hdl) 
    509 !     CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain")  
    510 !     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 
    511 !     CALL xios_set_attr(domain_hdl, lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d)                           ! set routing domain dimension 
    512 !     CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain_expand")  
    513 !     CALL xios_set_attr(domain_hdl, domain_ref="routing_domain") 
    514 !     CALL xios_add_child(domain_hdl,domain_expand_hdl) 
    515 !     CALL xios_set_attr(domain_expand_hdl, type="edge", i_periodic=.TRUE., j_periodic=.TRUE.) 
    516 !   
    517 !     CALL xios_get_handle("field_definition",field_def_hdl) 
    518 !     CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_init") 
    519 !     CALL xios_set_attr(field_hdl,domain_ref="routing_domain") 
    520 !     CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_final") 
    521 !     CALL xios_set_attr(field_hdl,domain_ref="routing_domain") 
    522  
    523 !     CALL xios_add_child(field_def_hdl,field_hdl,"trip_r") 
    524 !     CALL xios_set_attr(field_hdl,domain_ref="routing_domain") 
    525 !     CALL xios_add_child(field_def_hdl,field_hdl,"trip_rp1") 
    526 !     CALL xios_set_attr(field_hdl, field_ref="trip_r", domain_ref="routing_domain_expand", read_access=.TRUE.) 
    527  
    528 !     CALL xios_add_child(field_def_hdl,field_hdl,"state_r") 
    529 !     CALL xios_set_attr(field_hdl,domain_ref="routing_domain") 
    530 !     CALL xios_add_child(field_def_hdl,field_hdl,"state_rp1") 
    531 !     CALL xios_set_attr(field_hdl, field_ref="state_r", domain_ref="routing_domain_expand", read_access=.TRUE.) 
    532 !    
    533 !     CALL xios_get_handle("file_definition",file_def_hdl) 
    534 !     CALL xios_add_child(file_def_hdl,file_hdl,"routing_correction") 
    535 !     CALL xios_set_attr(file_hdl,type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.) 
    536 !     CALL xios_add_child(file_hdl, field_hdl) 
    537 !     CALL xios_set_attr(field_hdl, field_ref="trip_r", operation="instant") 
    538 !     CALL xios_add_child(file_hdl, field_hdl) 
    539 !     CALL xios_set_attr(field_hdl, field_ref="state_r", operation="instant") 
    540 !     CALL xios_add_child(file_hdl, field_hdl) 
    541 !     CALL xios_set_attr(field_hdl, field_ref="trip_r_init", operation="once") 
    542 !     CALL xios_add_child(file_hdl, field_hdl) 
    543 !     CALL xios_set_attr(field_hdl, field_ref="state_r_init", operation="once") 
    544  
    545933 
    546934      CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH) 
    … …  
    8311219    LOGICAL :: file_exists 
    8321220    REAL(r_std) :: epsilon = 1e-5 
    833  
     1221    CHARACTER(LEN=255) :: routing_file_type 
    8341222!_ ================================================================================================================================ 
    8351223    ! 
    … …  
    8791267    !Config Units = [days] 
    8801268 
     1269    CALL getin_p('SLOW_TCST', slow_tcst) 
     1270 
     1271    routing_file_type = "standard" 
     1272    CALL getin_p("routing_file_type", routing_file_type)  
     1273 
     1274    IF (TRIM(routing_file_type)=="standard") THEN 
     1275      topoind_factor=1000 
     1276    ELSE IF (TRIM(routing_file_type)=="merit") THEN 
     1277      topoind_factor=1 
     1278    ELSE 
     1279      CALL ipslerr_p(3,'routing_flow_init_local', & 
     1280        'Error in call routing_flow_init_local','getin routing_file_type : bad value, must be <standard> or <merit>','') 
     1281    ENDIF 
     1282 
    8811283    split_routing=1 
    8821284    CALL getin_p("SPLIT_ROUTING",split_routing) 
    … …  
    8921294 
    8931295    IF (is_omp_root) THEN 
    894  
    895        WHERE(contfrac_mpi<=epsilon) contfrac_mpi=0 
    896        WHERE(contfrac_mpi>=epsilon) contfrac_mpi=1 
    897         
    898        diag=0 
    899        WHERE(coastline_mpi) diag=1 
    900        CALL xios_send_field("is_coastline", diag) 
    901         
     1296       WHERE(contfrac_mpi <= epsilon) contfrac_mpi=0 
     1297       WHERE(contfrac_mpi >= 1-epsilon) contfrac_mpi=1 
     1298 
    9021299       CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj)    ! get routing domain dimension 
    9031300 
    … …  
    9331330 
    9341331       ALLOCATE(trip_extended_r(ni,nj)) 
     1332 
     1333      ! correction on coastline in case of LAM, the cells near the frontier are considered as coastline  
     1334      ! => interpolate model grid to routing grid and reinterpolate to model grid. Fractionnal cells will be considered as coastline too 
     1335    !  diag=1 
     1336    !  CALL xios_send_field("one", diag) 
     1337    !  CALL xios_recv_field("tmp1_coastline", diag_r) 
     1338    !  CALL xios_send_field("tmp2_coastline", diag_r) 
     1339    !  CALL xios_recv_field("lam_coastline", diag) 
     1340    !  WHERE(diag<1-epsilon) is_coastline=.TRUE. 
     1341 
     1342       diag=0 
     1343       WHERE(is_coastline) diag=1 
     1344       CALL xios_send_field("is_coastline", diag) 
    9351345 
    9361346       is_lakeinflow_r(:) = .FALSE. 
    … …  
    9601370       CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file 
    9611371 
    962        CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r)  
     1372!       CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r)  
     1373       CALL new_routing_flow_correct_riverflow(ni, nj, contfrac_mpi, is_coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r)  
    9631374    
    9641375       CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj))        ! send to xios trip array to update for halo 
    … …  
    10411452       CALL xios_send_field("routing_mask_r", diag_r) 
    10421453       CALL xios_recv_field("frac_routing", diag) 
     1454       WHERE (diag < epsilon) diag=0 
    10431455       routing_weight_in (:) = 0 
    10441456       unrouted_weight (:) = 0 
    … …  
    10491461      ! compute the number of coast cells to distribute lakeinflow onto 
    10501462       diag=0 
    1051        WHERE (coastline_mpi) diag=1 
     1463       WHERE (is_coastline) diag=1 
    10521464       nb_coast_points=SUM(diag) 
    10531465       CALL reduce_sum_mpi(nb_coast_points, total_coast_points) 
    … …  
    10671479 
    10681480       diag(:) = 0 
    1069        WHERE (coastline_mpi) diag=1 
     1481       WHERE (is_coastline) diag=1 
    10701482       CALL xios_send_field("mask_coastal",diag) 
    10711483       CALL xios_recv_field("frac_coastal_r",diag_r) 
    1072        WHERE (diag_r > 100.) ! missing_value 
     1484       WHERE (diag_r > 100) ! missing_value 
    10731485         diag_r=0 
    10741486       ELSE WHERE (diag_r < epsilon)  
    … …  
    10801492 
    10811493       diag(:) = 0 
    1082        WHERE (.NOT. coastline_mpi) diag=1 
     1494       WHERE (.NOT. is_coastline) diag=1 
    10831495       CALL xios_send_field("mask_lake",diag) 
    10841496       CALL xios_recv_field("frac_lake_r",diag_r) 
    1085        WHERE (diag_r > 100.) ! missing_value 
     1497       WHERE (diag_r > 100) ! missing_value 
    10861498         diag_r=0. 
    10871499       ELSEWHERE (diag_r < epsilon)  
    … …  
    11051517             weight_coast_to_lake_r(ij) = 1./frac_lake_r(ij) 
    11061518             weight_coast_to_coast_r(ij) = 0 
    1107            ELSE IF (frac_coast_r(ij)>0 ) THEN 
    1108              weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij)) 
     1519           ELSE IF (frac_coast_r(ij)>0) THEN 
     1520             weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij)) 
    11091521             weight_coast_to_lake_r(ij) = 0 
     1522           ELSE 
    11101523           ENDIF 
    11111524         ELSE IF (is_lakeinflow_r(ij)) THEN 
    … …  
    11171530             weight_lake_to_lake_r(ij) = 0 
    11181531           ELSE IF (frac_lake_r(ij)>0 ) THEN 
    1119              weight_lake_to_lake_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij)) 
     1532             weight_lake_to_lake_r(ij) = 1./(frac_lake_r(ij)) 
    11201533             weight_lake_to_coast_r(ij) = 0 
    11211534           ENDIF 
    … …  
    14321845      WHERE(.NOT. routing_mask_r) drainage_r = 0 
    14331846 
    1434       CALL MPI_ALLREDUCE(sum(runoff*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
     1847      CALL MPI_ALLREDUCE(sum(runoff*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
    14351848      CALL MPI_ALLREDUCE(sum(runoff_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
    14361849      IF (sum_water_after/=0) THEN 
    … …  
    14411854      ENDIF 
    14421855 
    1443       CALL MPI_ALLREDUCE(sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
     1856      CALL MPI_ALLREDUCE(sum(drainage*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
    14441857      CALL MPI_ALLREDUCE(sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
    14451858      IF (sum_water_after/=0) THEN 
    … …  
    14511864           
    14521865 
    1453       CALL MPI_ALLREDUCE(sum(runoff*routing_weight)+sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
    1454       CALL MPI_ALLREDUCE(sum(runoff_r) +sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
     1866      CALL MPI_ALLREDUCE(sum((runoff+drainage)*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
     1867      CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 
    14551868 
    14561869      IF (is_mpi_root) PRINT *,"Loose water by interpolation ;  before : ", sum_water_before," ; after : ",sum_water_after,  & 
    … …  
    14701883          IF ( routing_mask_r(ig) ) THEN 
    14711884 
    1472             flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),& 
     1885            flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),& 
    14731886                   & fast_reservoir_r(ig)-min_sechiba) 
    14741887            fast_flow_r(ig) = MAX(flow, zero) 
    14751888          ! 
    1476             flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),& 
     1889            flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),& 
    14771890                   & slow_reservoir_r(ig)-min_sechiba) 
    14781891            slow_flow_r(ig) = MAX(flow, zero) 
    14791892          ! 
    1480             flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)/1000.)*stream_tcst*one_day/(dt_routing/split_routing)),& 
     1893            flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*stream_tcst*one_day/(dt_routing/split_routing)),& 
    14811894                   & stream_reservoir_r(ig)-min_sechiba) 
    14821895            stream_flow_r(ig) = MAX(flow, zero) 
    … …  
    16412054      ENDIF    
    16422055        
    1643       coastalflow=coastalflow+flow_coast 
     2056      coastalflow=coastalflow+flow_coast+(runoff+drainage)*unrouted_weight 
    16442057      lakeinflow=lakeinflow+flow_lake 
    16452058 
    16462059 
    1647       WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in 
    1648       WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in 
     2060!      WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in 
     2061!      WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in 
    16492062 
    16502063      sum_water_after = sum(coastalflow(:) + riverflow(:) + lakeinflow(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:)) 
Note: See TracChangeset for help on using the changeset viewer.