!!======================================================================

  SUBROUTINE river(xe3)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE river  ***
!&E
!&E ** Purpose : calculation of river discharges (volume in m3/s)
!&E
!&E ** Description : discharges saved inside phiriver array
!&E
!&E ** Called by : main.F90, step.F90, Agrif_User.F90
!&E
!&E ** Channel : 29, 49
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : nbfl,phiriver,iriv,jriv,kriv,nomriviere
!&E                         l_flestat,debmoy,tmfle,amtfle,debit
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !                         Original code
!&E       !  2008-06  (V. Garnier) introduction of phiriver to allow the
!&E       !           prescription of the heat flux relative to river discharges
!&E       !  2009-12  (S. Michel, HOCER) check consistency between river and simulation dates
!&E       !  2010-03  (V. Garnier) correction of phiriver to allow several river
!&E       !                        discharges at the same location
!&E       !  2010-10  (B. Thouvenin) Distribution of the river discharge throughout a thickness
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2011-12  (B. Thouvenin) Correction : Distribution of the river discharge throughout a thickness
!&E       !  2012-03  (S. Petton) Introduction of the csv format
!&E       !  2011-04  (S. Petton) if river_thickness=-1., the amount of river discharge is distributed throughout the water column
!&E       !  2014-07  (L. Pineau) Reading of the hours in *.csv in case of CD-OCO hourly river discharges
!&E       !  2020-01  (F. Grasso) bottom friction z0 varying with the river flow
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d, ONLY : nbouc,nvalmax,h0,h0_g,tchrono,t,indchrono, &
                         iscreenlog,ierrorlog,iwarnlog,l_testcase,z0b
   USE comfluxdo, ONLY : phiriver
   USE comfluxdo, ONLY : iriv,jriv,kriv,nomriviere,l_flestat,nbfl,debmoy,river_thickness
   USE comfluxdo, ONLY : tmfle,amtfle,file_river,temp_riv,l_tempriv
   USE comfluxdo, ONLY : debit 

   USE toolgeom,  ONLY : f_zw
   USE toolgeom,  ONLY : f_dzu,f_zu
   USE comfluxdo, ONLY : phitot_glob 
   USE toolmpi,   ONLY : MASTER,IPROC,MPI_COMM_MARS,ADD_MPI_RSH,BCAST_MPI
   USE mpi
 
   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(IN) :: xe3

!! * Local declarations
   LOGICAL                                    :: ex,l_posit
   CHARACTER(LEN=5)                           :: comment
   CHARACTER(LEN=lchain), DIMENSION(:), ALLOCATABLE :: nom,nomtemp
   CHARACTER(LEN=lchain)                      :: rec
   INTEGER                                    :: i,ii,j,jj,k,kk,eof,ir,IERR_MPI
   INTEGER                                    :: l,flag
   REAL(KIND=rsh)                             :: val1,val2,tool_latlon2i,tool_latlon2j,river_thick
   REAL(KIND=rlg)                             :: tool_datosec,tint1,tint2,dt2,dt1,g,phi,t1
   REAL(kind=rsh)                             :: d3,dz,immer_top,immer_bot,immer_kriv,fdeb
   REAL(kind=rsh)                             :: debinst,phitot,errdeb
   INTEGER :: new_z0

!!----------------------------------------------------------------------
!! * Executable part

   !PRINT*, 'enter RIVER routine, nbouc=',nbouc

! ****  Initialization  ****
!---------------------------
   phitot=0.0

   IF (nbouc == 0 .AND. .NOT.l_testcase) THEN
! save into simu.log
! ------------------

     IF (MASTER) THEN
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) '****************** RIVER.F90 *********************'
     WRITE(iscreenlog,*) '**************************************************'
     WRITE(iscreenlog,*) ' '
     WRITE(iscreenlog,*) 'RIVER FILE :',TRIM(file_river)
     ENDIF

! check the number of rivers
! --------------------------

     nbfl=0
     eof=0

     INQUIRE(file=file_river,exist=ex)
     IF (ex) THEN
       OPEN(49,file = file_river,form='formatted')
       comment='debut'
       DO WHILE (comment /= '*****')
         READ(49,'(a)',iostat=eof) comment
       END DO
       DO WHILE (eof==0)
         READ(49,'(a)',iostat=eof)
         READ(49,*,iostat=eof)
         READ(49,*,iostat=eof)
         READ(49,*,iostat=eof)
         READ(49,*,iostat=eof)
         READ(49,*,iostat=eof)
         READ(49,'(a)',iostat=eof)
         READ(49,*,iostat=eof)
         READ(49,'(a)',iostat=eof)
         READ(49,*,iostat=eof)
         IF (eof==0) nbfl=nbfl+1
         READ(49,*,iostat=eof)
       END DO
       CLOSE(49)

! save into simu.log
! ------------------

       IF (MASTER) THEN
       WRITE(iscreenlog,*) ' '
       WRITE(iscreenlog,*) 'number of river discharges = ',nbfl
       ENDIF

! dynamic allocation
! ------------------

       ALLOCATE(nomriviere(nbfl))
       ALLOCATE(iriv(nbfl),jriv(nbfl),kriv(nbfl),debit(nvalmax,nbfl))
       ALLOCATE(l_flestat(nbfl),debmoy(nbfl),temp_riv(nvalmax,nbfl))
       ALLOCATE(nom(nbfl),nomtemp(nbfl),l_tempriv(nbfl),tmfle(nbfl),amtfle(nbfl))
       ALLOCATE(phiriver(nbfl,kmax))
       ALLOCATE(river_thickness(nbfl))

! read file for initialization
! ----------------------------

       OPEN(49,file = file_river,form='formatted')
       comment='debut'
       DO WHILE (comment /= '*****')
         READ(49,'(a)',iostat=eof) comment
       END DO

       DO k=1,nbfl   ! initialization of each river parametrization

         READ(49,'(a)') nomriviere(k)

         READ(49,*) l_posit
         IF (.NOT.l_posit) THEN

           READ(49,*) g,phi,kriv(k)
           iriv(k)=NINT(tool_latlon2i(g,phi))   ! cval cannot work for key_grid_rotated
           jriv(k)=NINT(tool_latlon2j(g,phi))
           READ(49,*) river_thickness(k)
         ELSE
           READ(49,*) iriv(k),jriv(k),kriv(k)
           READ(49,*) river_thickness(k)
         END IF
         ir=0
         DO WHILE (h0_g(iriv(k),jriv(k))==-valmanq)  ! check the closest wet cell
           ir=ir+1
           IF (MASTER) THEN
           IF (ir==1) WRITE(iwarnlog,*) ' '
           IF (ir==1) WRITE(iwarnlog,*) 'WARNING: The river ',TRIM(nomriviere(k)), ' is shifted '  
           ENDIF
           DO j=MAX(jmin,jriv(k)-ir),MIN(jmax,jriv(k)+ir)
             DO i=MAX(imin,iriv(k)-ir),MIN(imax,iriv(k)+ir)
               IF (h0_g(i,j) /= -valmanq) THEN
                 iriv(k)=i
                 jriv(k)=j
               END IF
             END DO
           END DO
         END DO
         IF (MASTER) THEN
         IF (ir/=0) WRITE(iwarnlog,*) 'into the cell (i,j) : ',iriv(k),jriv(k)
         ENDIF
         IF (kriv(k) > kmax) THEN
           WRITE(ierrorlog,*) ' '
           WRITE(ierrorlog,*) 'ERROR: In RIVER.F90'
           WRITE(ierrorlog,*) 'In river.dat input file, you have chosen kriv > kmax'
           WRITE(ierrorlog,*) 'look at river number',k,',kriv =',kriv(k)
           WRITE(ierrorlog,*) 'Please correct that, the river discharge may not be taken into account'
           WRITE(ierrorlog,*) 'The simulation is stopped'
           CALL MPI_FINALIZE(IERR_MPI)
           STOP
         ELSE IF (kriv(k) == 0) THEN
           WRITE(ierrorlog,*) ' '
           WRITE(ierrorlog,*) 'ERROR: In RIVER.F90'
           WRITE(ierrorlog,*) 'In river.dat input file, you have chosen kriv = 0'
           WRITE(ierrorlog,*) 'look at river number',k,',kriv =',kriv(k)
           WRITE(ierrorlog,*) 'Please correct that, the river discharge must in the water column : 1 <= kriv <=kmax'
           WRITE(ierrorlog,*) 'The simulation is stopped'
           CALL MPI_FINALIZE(IERR_MPI)
           STOP
         END IF

         READ(49,*) l_flestat(k)
         READ(49,*) debmoy(k)
         READ(49,'(a)') nom(k)
         READ(49,*) l_tempriv(k)
         READ(49,'(a)') nomtemp(k)
         READ(49,*) tmfle(k),amtfle(k)
         READ(49,*,END=11,err=11)

       END DO   ! ends loop on rivers
       CLOSE(49)
 11    CONTINUE


! save into simu.log and initialize the discharge amount
! ------------------------------------------------------

       DO k = 1,nbfl

         IF (MASTER) THEN
         WRITE(iscreenlog,*) ' '
         WRITE(iscreenlog,*) 'RIVER NAME        : ', TRIM(nomriviere(k))
         WRITE(iscreenlog,*) 'LOCATION i,j,k    : ',iriv(k),jriv(k),kriv(k)
         WRITE(iscreenlog,*) 'RIVER DEPTH       : ',river_thickness(k)
         ENDIF

! initialize the river discharge amount
! -------------------------------------

         IF (l_flestat(k)) THEN
           IF (MASTER) THEN
           WRITE(iscreenlog,*) 'MEAN FLOW         : ',debmoy(k) 
           ENDIF
         ELSE
           IF (MASTER) THEN
           WRITE(iscreenlog,*) 'FLOW FILE         : ',TRIM(nom(k))
           ENDIF
           OPEN(29,file = nom(k),form='formatted')
           i=1

! check the first available date in river discharge file
! ------------------------------------------------------
           IF (index(nom(k),'.dat') /= 0  .OR. index(nom(k),'.txt') /= 0 ) THEN

             READ(29,*)
             READ(29,*)
             READ(29,*)
             READ(29,*)
             READ(29,'(a)',END=999) rec
             kk=INDEX(rec,'/')-2
             READ(rec(kk+19:),*) val1
             tint1 = tool_datosec(rec(kk:kk+18))

           ELSE IF (index(nom(k),'.csv') /= 0) THEN

             READ(29,*)
             READ(29,'(a)',END=999) rec

             l=len_trim(rec)
             DO WHILE (l==0)
                READ(29,'(a)',END=999) rec
                l=len_trim(rec)
             ENDDO
             IF(rec(l:l)<'0'.and.rec(l:l) >'9') l=l-1
             READ(rec(l:l),*)flag
             DO WHILE( flag==9 .OR. flag==4 ) ! pour les mauvaises lignes
                READ(29,'(a)',END=999) rec
                l=len_trim(rec)
                DO WHILE (l.eq.0)
                   WRITE(*,*) 'ligne blanche'
                   READ(29,'(a)',END=999) rec
                   l=len_trim(rec)
                ENDDO
                IF (rec(l:l)<'0'.and.rec(l:l) >'9') l=l-1
                READ(rec(l:l),*)flag
             ENDDO

             kk=index(rec,'/')-4 ! qd on a une donnee valide lecture de la date
             IF (rec(kk+10:kk+10)==',') THEN
                 tint1 = tool_datosec(''//rec(kk+8:kk+9)//'/'//rec(kk+5:kk+6)//'/'//rec(kk:kk+3)//' 12:00:00')
             ELSE
                 tint1 = tool_datosec(''//rec(kk+8:kk+9)//'/'//rec(kk+5:kk+6)//'/'//rec(kk:kk+3)//' '//rec(kk+11:kk+18))
             ENDIF

             kk=index(rec(1:l-6),',',.true.)
             READ(rec(kk+1:l-6),*)val1

           ELSE

             WRITE(ierrorlog,*) ' '
             WRITE(ierrorlog,*) 'ERROR: In RIVER.F90'
             WRITE(ierrorlog,*) 'In river.dat input file, river file names must end with :'
             WRITE(ierrorlog,*) '     - .txt or .dat for MARS format'
             WRITE(ierrorlog,*) '     - .csv for CDOCO format'
             WRITE(ierrorlog,*) ' look at river ',trim(nom(k))
             WRITE(ierrorlog,*) 'The simulation is stopped'
             CALL MPI_FINALIZE(IERR_MPI)
             STOP

           END IF

! If the simulation starts before the first data record,
! the river flow is set to the first available data.
           t1=tchrono(i)
           DO WHILE (t1 < tint1.AND.i<=nvalmax)
             debit(i,k)=val1
             i=i+1
             IF (i<=nvalmax) t1=tchrono(i)
           END DO
           IF (i > 1) THEN
             IF (MASTER) THEN
             WRITE(iwarnlog,*) ' '
             WRITE(iwarnlog,*) 'WARNING: The flow data for the river ',TRIM(nomriviere(k))
             WRITE(iwarnlog,*) 'is not available from the beginning of the simulation.'
             WRITE(iwarnlog,*) 'The first available value from the file : ',TRIM(nom(k))
             WRITE(iwarnlog,*) 'corresponds to the date : ',rec(kk:kk+18)
             WRITE(iwarnlog,*) 'It is equal to : ',val1,' m3/s'
             WRITE(iwarnlog,*) 'and will be used at the beginning of the simulation.'
             ENDIF
           END IF

! time interpolation of river discharge data
           IF (i < nvalmax) THEN
             IF (index(nom(k),'.dat') /= 0  .OR. index(nom(k),'.txt') /= 0) THEN
               eof=0
               DO WHILE (eof==0)
                  READ(29,'(a)',END=999) rec
                  kk=INDEX(rec,'/')-2
                  READ(rec(kk+19:),*) val2
                  tint2 = tool_datosec(rec(kk:kk+18))
                  DO WHILE ((tchrono(i)-tint1)*(tchrono(i)-tint2) <= 0.0_rlg)
                    dt1=(tchrono(i)-tint1)/(tint2-tint1)
                    dt2=(tint2-tchrono(i))/(tint2-tint1)
                    debit(i,k)=val1*dt2+val2*dt1
                    i=i+1
                    IF (i > nvalmax) GOTO 999
                  END DO 
                  tint1=tint2
                  val1=val2
               END DO
             ELSE
               eof=0
               DO WHILE (eof==0)
                  READ(29,'(a)',END=999) rec
                  l=len_trim(rec)
                  DO WHILE (l==0)
                    READ(29,'(a)',END=999) rec
                    l=len_trim(rec)
                 ENDDO
                 IF (rec(l:l)<'0'.AND.rec(l:l) >'9') THEN
                    l=l-1
                 ENDIF
                 READ(rec(l:l),*)flag
                 DO WHILE((flag.EQ.9).or.(flag.EQ.4))
                    READ(29,'(a)',END=999) rec
                    l=len_trim(rec)
                    DO WHILE (l==0)
                       READ(29,'(a)',END=999) rec
                       l=len_trim(rec)
                    ENDDO
                    IF(rec(l:l)<'0'.AND.rec(l:l) >'9') THEN
                       l=l-1
                    ENDIF
                    READ(rec(l:l),*)flag
                 ENDDO
! lecture de la date
                 kk=index(rec,'/')-4
                 IF (rec(kk+10:kk+10)==',') THEN
                    tint2 = tool_datosec(''//rec(kk+8:kk+9)//'/'//rec(kk+5:kk+6)//'/'//rec(kk:kk+3)//' 12:00:00')
                 ELSE
                    tint2 = tool_datosec(''//rec(kk+8:kk+9)//'/'//rec(kk+5:kk+6)//'/'//rec(kk:kk+3)//' '//rec(kk+11:kk+18))
                 ENDIF
! lecture du debit
                 kk=index(rec(1:l-6),',',.true.)
                 READ(rec(kk+1:l-6),*)val2

                 DO WHILE ((tchrono(i)-tint1)*(tchrono(i)-tint2) <= 0.0_rlg)
                   dt1=(tchrono(i)-tint1)/(tint2-tint1)
                   dt2=(tint2-tchrono(i))/(tint2-tint1)
                   debit(i,k)=val1*dt2+val2*dt1
                   i=i+1
                   IF (i > nvalmax) GOTO 999
                 END DO 
                 tint1=tint2
                 val1=val2
               END DO
             END IF
           END IF

 999       CONTINUE

! If the simulation ends after the last data record,
! set the river flow to the last available data.
           IF (i < nvalmax) THEN
             IF (MASTER) THEN
             WRITE(iwarnlog,*) ' '
             WRITE(iwarnlog,*) 'WARNING: The flow data for the river ',TRIM(nomriviere(k))
             WRITE(iwarnlog,*) 'is not available to the end of the simulation.'
             WRITE(iwarnlog,*) 'The last available value from the file : ',TRIM(nom(k))
             WRITE(iwarnlog,*) 'corresponds to the date : ',rec(kk:kk+18)
             WRITE(iwarnlog,*) 'It is equal to : ',val2,' m3/s'
             WRITE(iwarnlog,*) 'and will be used until the end of the simulation.'
             ENDIF
             debit(i:nvalmax,k)=val2
           END IF

           CLOSE(29)

         END IF   ! ends test on l_flestat

! initialize the river temperature
! --------------------------------

         IF (.NOT.l_tempriv(k)) THEN
           IF (MASTER) THEN
           WRITE(iscreenlog,*) 'MEAN/STD TEMP     : ',tmfle(k),amtfle(k)
           ENDIF
         ELSE

! read file of river temperature
! ------------------------------

           IF (MASTER) THEN
           WRITE(iscreenlog,*) 'TEMPERATURE FILE  : ',TRIM(nomtemp(k))
           ENDIF
           OPEN(29,file = nomtemp(k),form='formatted')
           i=1

! check the first available date in river discharge file
! ------------------------------------------------------

           READ(29,*)
           READ(29,*)
           READ(29,*)
           READ(29,*)
           READ(29,'(a)',err=998,END=998) rec
           kk=INDEX(rec,'/')-2
           READ(rec(kk+19:),*) val1
           tint1 = tool_datosec(rec(kk:kk+18))

           t1=tchrono(i)

! If the simulation starts before the first data record,
! the river temperature is set to the first available data.
           DO WHILE (t1 < tint1 .AND. i<=nvalmax)
             temp_riv(i,k)=val1
             i=i+1
             IF (i<=nvalmax) t1=tchrono(i)
           END DO
           IF (i > 1) THEN
             IF (MASTER) THEN
             WRITE(iwarnlog,*) ' '
             WRITE(iwarnlog,*) 'WARNING: The emperature data for the river ',TRIM(nomriviere(k))
             WRITE(iwarnlog,*) 'is not available from the beginning of the simulation.'
             WRITE(iwarnlog,*) 'The first available value from the file : ',TRIM(nomtemp(k))
             WRITE(iwarnlog,*) 'corresponds to the date : ',rec(kk:kk+18)
             WRITE(iwarnlog,*) 'It is equal to : ',val1,' deg.C'
             WRITE(iwarnlog,*) 'and will be used at the beginning of the simulation.'
             ENDIF
           END IF

! time interpolation of river temperature data
           IF (i < nvalmax) THEN
             eof=0
             DO WHILE (eof==0)
               READ(29,'(a)',END=998) rec
               kk=INDEX(rec,'/')-2
               READ(rec(kk+19:),*) val2
               tint2 = tool_datosec(rec(kk:kk+18))
               DO WHILE ((tchrono(i)-tint1)*(tchrono(i)-tint2) <= 0.0_rlg)
                 dt1=(tchrono(i)-tint1)/(tint2-tint1)
                 dt2=(tint2-tchrono(i))/(tint2-tint1)
                 temp_riv(i,k)=val1*dt2+val2*dt1
                 i=i+1
                 IF (i > nvalmax) GOTO 998
               END DO
               tint1=tint2
               val1=val2
             END DO
           END IF

 998       CONTINUE

! If the simulation ends after the last data record,
! set the river temperature to the last available data.
           IF (i < nvalmax) THEN
             IF (MASTER) THEN
             WRITE(iwarnlog,*) ' '
             WRITE(iwarnlog,*) 'WARNING: The temperature data for the river ',TRIM(nomriviere(k))
             WRITE(iwarnlog,*) 'is not available to the end of the simulation.'
             WRITE(iwarnlog,*) 'The last available value from the file : ',TRIM(nomtemp(k))
             WRITE(iwarnlog,*) 'corresponds to the date : ',rec(kk:kk+18)
             WRITE(iwarnlog,*) 'It is equal to : ',val2,' deg.C'
             WRITE(iwarnlog,*) 'and will be used until the end of the simulation.'
             temp_riv(i:nvalmax,k)=val2
             ENDIF
           END IF

           CLOSE(29)

         END IF
       END DO   ! ends loop on rivers

       DEALLOCATE(nom,nomtemp)

     ELSE

       IF (MASTER) THEN
       WRITE(iwarnlog,*) ''
       WRITE(iwarnlog,*) 'WARNING : RIVER routine'
       WRITE(iwarnlog,*) 'The ',TRIM(file_river),' file does not exist.'
       WRITE(iwarnlog,*) 'No river discharge is taking into account'
       ENDIF

     END IF   ! ends test on file availability

   END IF   ! ends test on nbouc

! **** End of initialization part ****
! ------------------------------------





! ****  Current progress  ****
! ----------------------------


   DO k = 1,nbfl
     ii=iriv(k)
     jj=jriv(k)
     phiriver(k,:)=0.0_rsh
     debinst=0.0_rsh
     IF (jj >= ljmin .AND. jj <= ljmax .AND. ii >= limin .AND. ii <= limax) THEN
       d3=h0(ii,jj)+xe3(ii,jj)
       kk=kriv(k)
! river thickness estimated downwards from the interface of the layer kriv
! if river_thickness=-1., the amount of river discharge is distributed throughout the water column
! if discharge at the bottom, kriv is set to 1 by the user in river.dat,
! none river thickness is taken into account
! if river_thickness=0., none river thickness is taken into account
! verification of the initialisation of river_thickness
       IF (river_thickness(k) == 0.0_rsh) THEN
         river_thick = f_dzu(h0(ii,jj),xe3(ii,jj),kk,ii,jj)
       ELSE IF (river_thickness(k) > f_zw(h0(ii,jj),xe3(ii,jj),kk,ii,jj) + h0(ii,jj)) THEN
         river_thick = h0(ii,jj) +f_zw(h0(ii,jj),xe3(ii,jj),kk,ii,jj)  ! the depth f_zw is negative
       ELSE
         river_thick = river_thickness(k)
       END IF
       IF (l_flestat(k)) THEN
         debinst=debmoy(k)


       ELSE
         dt1=(t-tchrono(indchrono))/(tchrono(indchrono+1)-tchrono(indchrono))
         debinst=dt1*debit(indchrono+1,k)+(1.d0-dt1)*debit(indchrono,k)
       END IF

! wetting
! if river_thickness=-1., the amount of river discharge is distributed throughout the water column
       IF (d3 == 0.0_rsh .OR. river_thick==-1.0_rsh) THEN
         phiriver(k,:)=debinst/kmax
       ELSE  
! immersion from the sea surface. The water thickness is xe - z (z being negative)
         immer_kriv= - (f_zw(h0(ii,jj),xe3(ii,jj),kk  ,ii,jj)-xe3(ii,jj))
         immer_bot = - (f_zw(h0(ii,jj),xe3(ii,jj),kk-1,ii,jj)-xe3(ii,jj))
         immer_top=immer_kriv
!  test on kk==kriv(k) to make sure you enter the loop in case of rounding errors in estimates of river_thick,immer_bot,immer_kriv
         DO WHILE ( ((immer_bot - immer_kriv) <= river_thick .AND. kk/=0) .OR. (kk==kriv(k)) )
           immer_bot = - (f_zw(h0(ii,jj),xe3(ii,jj),kk-1,ii,jj)-xe3(ii,jj))
           dz = f_dzu(h0(ii,jj),xe3(ii,jj),kk,ii,jj)
           IF ((immer_bot -immer_kriv) <= river_thick) THEN
             fdeb=dz/river_thick
           ELSE
             fdeb=(river_thick-immer_top+immer_kriv)/river_thick
           END IF
           phiriver(k,kk)=debinst*fdeb
           immer_top = immer_bot
           kk = kk-1
         END DO
       ENDIF
!verification
       phitot=SUM(phiriver(k,:))
       errdeb=debinst*1.0e-5
       IF (phitot < debinst-errdeb .OR. phitot > debinst+errdeb) THEN
         write(ierrorlog,*)'PROBLEM during the evaluation of river fluxes'
         write(ierrorlog,*)'in subroutine river.F90'
         write(ierrorlog,*)'river number',k
         write(ierrorlog,*)' reading total input flux=',debinst,' river_thickness=',river_thickness
         write(ierrorlog,*)' different from computing total input flux=',phitot
         write(ierrorlog,*)'flux into each level k', phiriver(k,:)
         write(ierrorlog,*)
         write(ierrorlog,*)'it could be happen if no water in the cell and d3<0 '
         write(ierrorlog,*)' here d3=', d3
         write(ierrorlog,*)'if d3<0 choose the key_dyn_wetdry_fct'
         CALL MPI_FINALIZE(IERR_MPI)
         STOP
       ENDIF
       IF (l_tempriv(k)) &
         tmfle(k)=dt1*temp_riv(indchrono+1,k)+(1.d0-dt1)*temp_riv(indchrono,k)
     ENDIF
   END DO
!Modif MPI pour partager le debit avec tous les procs
      phitot_glob=0.0_rsh
      phitot=SUM(phiriver(1,:)) ! 1=LOIRE
!$OMP barrier
!$OMP master
      CALL ADD_MPI_RSH(phitot,phitot_glob)
!CALL MPI_BCAST(phitot_glob,1,MPI_REAL,0,MPI_COMM_MARS,ierr_mpi)
      CALL BCAST_MPI(phitot_glob)
!$OMP end master
!$OMP barrier
!$OMP flush(phitot_glob)
!FIN MODIF MPI

!update Z0 with phitot in the estuary vFG (11/06/2020)
!write(111,*)'phitot',IPROC,phitot_glob
   new_z0=1
   IF (new_z0==1) THEN

   DO j=50,90
     IF (ljmin <= j .AND. j <= ljmax) THEN
!!! Mouth section: 1e-3
       DO i=130,200
         IF (limin <= i .AND. i <= limax) THEN
           z0b(i,j) = 1e-3-((j-68)/12.0_rsh*4e-4)
         ENDIF
       END DO
!!! Transition: 1e-3 -> 1e-4
       DO i=201,210
         IF (limin <= i .AND. i <= limax) THEN
           z0b(i,j) = 0.001-(i-201)*1e-4
         ENDIF
       END DO
!!! Mid section: 1e-6 pour Q=100 -> 1e-5 pour Q=1000 -> 1e-4 pour Q=2000 -> 2e-4 pour Q=3000
       DO i=211,265
         IF (limin <= i .AND. i <= limax) THEN
           z0b(i,j) = 1e-6+min((phitot_glob/100.0_rsh-1.0_rsh)*1e-6,9e-6)+max(phitot_glob/100.0_rsh-10.0_rsh,0.0_rsh)*1e-5
         ENDIF
       END DO
!!! Embanked section 1: 1e-5 pour Q=0 -> 7e-3 pour Q=1000 -> 3.6e-2 pour Q=2000
       DO i=266,280
         IF (limin <= i .AND. i <= limax) THEN
           z0b(i,j) = 1e-5+phitot_glob/100.0_rsh*3e-4+max(phitot_glob/50.0_rsh-10.0_rsh,0.0_rsh)*3e-4+max(phitot_glob/100.0_rsh-10.0_rsh,0.0_rsh)*2e-3
         ENDIF
       END DO
!!! Embanked section 2: 1e-5 pour Q=0 -> 7e-3 pour Q=1000 -> 4e-2 pour Q=2000
       DO i=281,360
         IF (limin <= i .AND. i <= limax) THEN
           z0b(i,j) = 1e-5+phitot_glob/100.0_rsh*3e-4+max(phitot_glob/150.0_rsh-10.0_rsh,0.0_rsh)*1e-4
         ENDIF
       END DO
!!! Up-river section: 3e-3
       DO i=361,imax
         IF (limin <= i .AND. i <= limax) THEN
           z0b(i,j) = 3e-3
         ENDIF
       END DO
     ENDIF
   END DO

   ENDIF
   !PRINT*, 'exit RIVER routine'

  END SUBROUTINE river
