  MODULE submassbalance

#ifdef key_substance

   !!======================================================================
   !!                   ***  MODULE  submassbalance  ***
   !! Ocean dynamics Bio Sed:  Mass budget in one or seveal zones 
   !!                          and fluxes threw on or several borders
   !! 
   !!  ** History :
   !!        !  2010-07  (B. Thouvenin)  based on sed_massbalance
   !!        !  2011-04  (F. Cayocca, B. Thouvenin) introduction of the fraction of wet cell (fwet)
   !!        !  2015-02  (B. Thouvenin) correction to take into account dissolved substances in sediment
   !!        !  2015-12  (B. Thouvenin, E. Schulz) corrections for complex zones, tdeb_budget, one file_border for each border 
   !!        !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
   !!        !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed 
   !!  
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters,ONLY : rsh,rlg,lchain,imin,imax,jmin,jmax,kmax,limin,limax,ljmin,ljmax,  &
                         liminm1,limaxp2,ljminm1,ljmaxp2,kmax,liminm2,ljminm2,nb_var
   USE comvars2d, ONLY : iscreenlog,fwet
   USE comsubstance,   ONLY : l_bilan,nb_border,name_bilfil,dtout_budget,nv_adv,tdeb_budget,tdebsubs,nvp,nv_fix
   USE output,       ONLY : fileresu,l_filebydate

   IMPLICIT NONE

   PRIVATE

   !! * Accessibility
   PUBLIC sub_budget_readdomain              ! routine called by main
   PUBLIC sub_budget_flxcum                  ! routine called by step3dxy,step3dyx
!   PUBLIC parsub_budget_flxcum                  ! routine called by parsub_eqa (elimine pour nouveau schema - budgets computed with dissolved)
   PUBLIC sub_budget_out                     ! routine called by step

   !! * Shared module variables
   REAL(KIND=rlg),PUBLIC                        :: t_out_budget           ! time output budgets
   INTEGER,DIMENSION(:,:,:),ALLOCATABLE,PUBLIC  :: mask_bud_in
   INTEGER,PUBLIC                               :: nb_zone_budget

   REAL(KIND=rlg),DIMENSION(:,:),ALLOCATABLE    :: bil,stok_sed,stok_wat,stok_wat_fix
   REAL(KIND=rlg),DIMENSION(:,:),ALLOCATABLE    :: flx_obc_cum,flx_obczon_cum
   REAL(KIND=rlg),DIMENSION(:,:),ALLOCATABLE    :: flx_rej_cum,flx_atm_cum,flx_ws_cum
   REAL(KIND=rlg),DIMENSION(:,:),ALLOCATABLE    :: stokw_total,stoks_total,stokwfix_total
   REAL(KIND=rlg),DIMENSION(:,:),ALLOCATABLE    :: flxobc_total,flxobczon_total
   REAL(KIND=rlg),DIMENSION(:,:),ALLOCATABLE    :: flxrej_total,flxatm_total,flxws_total
   INTEGER,DIMENSION(:,:,:),ALLOCATABLE         :: mask_bord_north,mask_bord_east
   INTEGER,DIMENSION(:),ALLOCATABLE             :: iborder_bud
   CHARACTER(LEN=lchain),DIMENSION(:),ALLOCATABLE :: name_bord ,file_budget_out                   
   CHARACTER(LEN=lchain),DIMENSION(:),ALLOCATABLE :: file_border_out                   
   INTEGER                                      :: nb_border_open

#ifdef key_contaminant
   REAL(kind=rlg),DIMENSION(:,:,:),ALLOCATABLE,PUBLIC    :: flx_reac_zone
   REAL(kind=rlg),DIMENSION(:,:,:),ALLOCATABLE,PUBLIC    :: flx_reac_total
   INTEGER, DIMENSION(:)          ,ALLOCATABLE ,PUBLIC   :: nre_act
   INTEGER, DIMENSION(:,:),ALLOCATABLE,PUBLIC            :: ire_act
#endif
   !! * Private variables

  CONTAINS

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

  SUBROUTINE sub_budget_readdomain(xe3)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_budget_readdomain  ***
   !&E
   !&E ** Purpose : read "budget.dat" file describing zones and border 
   !&E              where w want to compute budgets and fluxes
   !&E
   !&E ** Description : usefull for dissolved, biologic, sedimentologic... variables
   !&E
   !&E ** Called by : main.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2010-07 (B. Thouvenin) Original code
   !&E       !  2015-12 (E. Schulz, B. Thouvenin) correction for complexes budgets areas and borders
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,      ONLY :  iscreenlog,ierrorlog,iwarnlog,h0,surf
   USE parameters,ONLY : valmanq,num_testcase
   USE_MPI toolmpi,  ONLY : MASTER,ADD_MPI_RSH
#ifdef key_contaminant
   USE contareaddat,  ONLY : nreacmax,nb_var_re,irkm_var_re
#endif
#if defined key_castest
   USE parameters,ONLY : num_testcase
#endif
   !! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in)    :: xe3

   !! * Local declarations
   LOGICAL                                   :: ex,l_border_close
   INTEGER                                   :: ibord_r,iborder,eof,ind_white,ib,ism1,isp1,         &
                                                        IERR_MPI,nb_seg,is,izon_bud,                &
                                                        ig_bord,id_bord,jb_bord,jh_bord,i,j,ifil
   INTEGER, DIMENSION(:),ALLOCATABLE         :: i1,i2,i3
#ifdef key_contaminant
   INTEGER                                   :: ire,ivr,iv,nreactmax
#endif
   CHARACTER(LEN=lchain)                     :: nombord_r
   CHARACTER(LEN=5)                          :: comment
   CHARACTER(LEN=1),DIMENSION(:),ALLOCATABLE :: rose
   INTEGER,DIMENSION(imin:imax,jmin:jmax)    :: mask_bud_E_tmp,mask_bud_W_tmp
   INTEGER,DIMENSION(imin:imax,jmin:jmax)    :: mask_bud_S_tmp,mask_bud_N_tmp
   INTEGER,DIMENSION(:,:,:),ALLOCATABLE      :: mask_bud_in_tmp
   INTEGER,DIMENSION(:,:,:),ALLOCATABLE      :: mask_bord_north_tmp,mask_bord_east_tmp
   REAL(kind=rsh),DIMENSION(:),ALLOCATABLE      :: surf_tot,surftotale,vol_tot,voltotal

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

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

   IF_MPI (MASTER) THEN
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) '**************** BUDGETS_READ_LIMITS.F90 **********************'
   WRITE(iscreenlog,*) '**************************************************'
   WRITE(iscreenlog,*) ' '
   WRITE(iscreenlog,*) 'file defining budget zones and borders for estimate fluxes :'
   WRITE(iscreenlog,*) TRIM(name_bilfil)
   ENDIF_MPI

   IF (TRIM(name_bilfil) =="" .or. TRIM(name_bilfil) == "all_domain") THEN
      ! file "name_bilfil" is not reading
     ifil=0
     IF_MPI (MASTER) THEN
     WRITE(iscreenlog,*)'only one budget zone is taking into account (all the domain)'
     WRITE(iscreenlog,*)'and any fluxes threw open borders '
     ENDIF_MPI

     nb_zone_budget=1
     nb_border_open=0
     
#if defined key_castest
     IF(num_testcase==33) THEN
     ! test 2 zones for cas test tidalflat2DV
       nb_zone_budget=2
       nb_border=2
     ENDIF
#endif

   ELSE    
     ! first reading   
     eof=0
     INQUIRE(file=name_bilfil,exist=ex)
     IF (ex) THEN
       ifil=1
       OPEN(49,file = name_bilfil,form='formatted')
       comment='debut'
       DO WHILE (comment /= '*****')
         READ(49,'(a)',iostat=eof) comment
       END DO
       nb_zone_budget=0
       nb_border_open=0
       DO iborder=1,nb_border
         READ(49,*,iostat=eof)ibord_r
         IF (eof==0) THEN
           READ(49,*)
           READ(49,*)l_border_close
           IF (l_border_close) THEN
             nb_zone_budget=nb_zone_budget+1
           ELSE
              nb_border_open=nb_border_open+1
           END IF
           READ(49,*)nb_seg
           DO is=1,nb_seg
              READ(49,*)
           END DO
           READ(49,*)
         ELSE
           WRITE(ierrorlog,*) 'ERROR : the total number of border for the '
           WRITE(ierrorlog,*) 'estimation of fluxes and budgets is =',nb_border
           WRITE(ierrorlog,*) '(read in parasubs.txt)'
           WRITE(ierrorlog,*) 'It is not the good number in the data file =',TRIM(name_bilfil)
           CALL_MPI MPI_FINALIZE(IERR_MPI)
           STOP
         END IF
       END DO
       CLOSE(49)

     ELSE
       WRITE(ierrorlog,*) 'ERROR : the data file for budget does not exist'
       WRITE(ierrorlog,*) 'you have given l_bilan=.true. in parasubs.txt'
       WRITE(ierrorlog,*) 'and the name for the data file =',TRIM(name_bilfil)
       WRITE(ierrorlog,*) 'build or rename your data file'
       WRITE(ierrorlog,*) 'or set l_bilan=.false.' 
       WRITE(ierrorlog,*) 'if you don t want the calculation of budget'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF     ! fin du test sur l existence du fichier
     
   END IF     ! fin du test sur un fichier ou bilan sur tout le domaine


    ! allocation and initialize masks, stock and fluxes
   ALLOCATE(name_bord(nb_border))
   ALLOCATE(mask_bord_north(nb_border,liminm1:limaxp2,ljminm1:ljmaxp2))
   ALLOCATE(mask_bord_east(nb_border,liminm1:limaxp2,ljminm1:ljmaxp2))
   ALLOCATE(mask_bord_north_tmp(nb_border,imin:imax,jmin:jmax))
   ALLOCATE(mask_bord_east_tmp(nb_border,imin:imax,jmin:jmax))
   ALLOCATE(flx_obc_cum(nb_border_open,-1:nv_adv))
   ALLOCATE(flxobc_total(nb_border_open,-1:nv_adv))
   ALLOCATE(file_border_out(nb_border_open))
   IF (nb_zone_budget > 0) THEN
     ALLOCATE(file_budget_out(nb_zone_budget))
     ALLOCATE(iborder_bud(nb_zone_budget))
     ALLOCATE(mask_bud_in(nb_zone_budget,liminm1:limaxp2,ljminm1:ljmaxp2))
     ALLOCATE(mask_bud_in_tmp(nb_zone_budget,imin:imax,jmin:jmax))
     ALLOCATE(flx_rej_cum(nb_zone_budget,1:nv_adv))
     ALLOCATE(flx_atm_cum(nb_zone_budget,1:nv_adv))
     ALLOCATE(flx_obczon_cum(nb_zone_budget,1:nv_adv))
     ALLOCATE(flxrej_total(nb_zone_budget,1:nv_adv))
     ALLOCATE(flxatm_total(nb_zone_budget,1:nv_adv))
     ALLOCATE(flxobczon_total(nb_zone_budget,1:nv_adv))
#ifdef key_sedim_MUSTANG
     ALLOCATE(stok_sed(nb_zone_budget,1:nv_adv))
     ALLOCATE(flx_ws_cum(nb_zone_budget,1:nv_adv))
     ALLOCATE(stoks_total(nb_zone_budget,1:nv_adv))
     ALLOCATE(flxws_total(nb_zone_budget,1:nv_adv))
#endif
     ALLOCATE(stok_wat(nb_zone_budget,1:nv_adv))
     IF(nv_fix > 0 ) THEN
       ALLOCATE(stok_wat_fix(nb_zone_budget,1:nv_fix))
       ALLOCATE(stokwfix_total(nb_zone_budget,1:nv_fix))
     ENDIF
     ALLOCATE(bil(nb_zone_budget,1:nv_adv))
     ALLOCATE(stokw_total(nb_zone_budget,1:nv_adv))
     ALLOCATE(surf_tot(nb_zone_budget))
     ALLOCATE(surftotale(nb_zone_budget))
     ALLOCATE(vol_tot(nb_zone_budget))
     ALLOCATE(voltotal(nb_zone_budget))
   END IF

   !initialization
   IF (nb_zone_budget > 0) THEN
     mask_bud_in(:,:,:)=0
     mask_bud_in_tmp(:,:,:)=0
     flx_rej_cum(:,:)=0.0_rlg
     flx_atm_cum(:,:)=0.0_rlg
#ifdef key_sedim_MUSTANG
      flx_ws_cum(:,:)=0.0_rlg
#endif
      flx_obczon_cum(:,:)=0.0_rlg
   ENDIF
   mask_bord_north_tmp(:,:,:)=0
   mask_bord_east_tmp(:,:,:)=0
   mask_bord_north_tmp(:,:,:)=0
   mask_bord_east_tmp(:,:,:)=0
   flx_obc_cum(:,:)=0.0_rlg
   dtout_budget=dtout_budget*3600.0_rlg
   WRITE(*,*)'dtout_budget',dtout_budget
   t_out_budget=tdeb_budget
    IF(tdeb_budget < tdebsubs)t_out_budget=tdebsubs
    
   IF (ifil==0) THEN
#if defined key_castest
     IF(num_testcase==33) THEN
       ! test zone profonde
       iborder=1
       name_bord(iborder)="zone profonde"
       izon_bud=1
       iborder_bud(izon_bud)=1
       mask_bud_in_tmp(1,imin+5:99,jmin+2:jmax-1)=1
       mask_bord_north_tmp(1,imin+5:99,jmax-1)=-1
       mask_bord_north_tmp(1,imin+5:99,jmin+1)=1
       mask_bord_east_tmp(1,imin+4,jmin+2:jmax-1)=1
       mask_bord_east_tmp(1,99,jmin+2:jmax-1)=-1

     ! test zone peu profonde
       iborder=2
       name_bord(iborder)="zone peu profonde"
       izon_bud=2
       iborder_bud(izon_bud)=2
       mask_bud_in_tmp(2,150:imax-1,jmin+2:jmax-1)=1
       mask_bord_north_tmp(2,150:imax-1,jmax-1)=-1
       mask_bord_north_tmp(2,150:imax-1,jmin+1)=1
       mask_bord_east_tmp(2,149,jmin+2:jmax-1)=1
       mask_bord_east_tmp(2,imax-1,jmin+2:jmax-1)=-1
      ELSE
#endif
       iborder=1
       name_bord(iborder)="all_the_domain"
       izon_bud=1
       iborder_bud(izon_bud)=1
       mask_bud_in_tmp(1,imin+2:imax-1,jmin+2:jmax-1)=1
       mask_bord_north_tmp(1,imin+2:imax-1,jmax-1)=-1
       mask_bord_north_tmp(1,imin+2:imax-1,jmin+1)=1
       mask_bord_east_tmp(1,imin+1,jmin+2:jmax-1)=1
       mask_bord_east_tmp(1,imax-1,jmin+2:jmax-1)=-1
#if defined key_castest
      ENDIF
#endif


   ELSE
     OPEN(49,file = name_bilfil,form='formatted')
     comment='debut'
     DO WHILE (comment /= '*****')
         READ(49,'(a)',iostat=eof) comment
     END DO
       izon_bud=0
     ! FIRST open borders
     DO iborder=1,nb_border_open
       READ(49,*,iostat=eof)ibord_r
       READ(49,'(a)',iostat=eof) nombord_r
       ind_white=INDEX(nombord_r,' ')
       name_bord(iborder)=TRIM(ADJUSTL(ADJUSTR(nombord_r(1:ind_white))))
       READ(49,*)l_border_close
       IF (l_border_close) THEN
         WRITE(ierrorlog,*) 'ERROR in the data file for budget :',TRIM(name_bilfil)
         WRITE(ierrorlog,*) 'you must give first all the open borders '
         WRITE(ierrorlog,*) 'before closed borders '
         CALL_MPI MPI_FINALIZE(IERR_MPI)
         STOP   
       END IF
       READ(49,*)nb_seg
       ALLOCATE(i1(nb_seg),i2(nb_seg),i3(nb_seg))
       ALLOCATE(rose(nb_seg))
       DO is=1,nb_seg
         READ(49,*)i1(is),i2(is),i3(is),rose(is)
       ENDDO
       DO is=1,nb_seg
         ism1=is-1
         isp1=is+1
         IF (rose(is) == 'N') THEN
           ! verification
           IF_MPI (MASTER) THEN
           IF(is > 1) THEN
             IF( rose(ism1) /= 'N' .AND. (i1(is) .NE. i1(ism1) .OR. i3(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of open border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment NORD number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i3 of segment',is,'must be = to i1 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
             ENDIF
           ENDIF
           ENDIF_MPI
           IF (i3(is) >= jmin .and.i3(is) <= jmax) THEN
             ig_bord=MIN0(i1(is),i2(is))
             id_bord=MAX0(i1(is),i2(is))
             DO i=MAX0(imin,ig_bord+1),MIN0(imax,id_bord)
               mask_bord_north_tmp(iborder,i,i3(is))=-1    
             END DO
           END IF
         ELSE IF (rose(is) == 'S') THEN
            ! verification
           IF_MPI (MASTER) THEN
           IF(is > 1) THEN
             IF(rose(ism1) /= 'S' .AND. (i1(is) .NE. i1(ism1) .OR. i3(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of open border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment SOUTH number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i3 of segment',is,'must be = to i1 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
             ENDIF
           ENDIF
           ENDIF_MPI
           IF (i3(is)-1 >= jmin .and.i3(is)-1 <= jmax) THEN
             ig_bord=MIN0(i1(is),i2(is))
             id_bord=MAX0(i1(is),i2(is))
             DO i=MAX0(imin,ig_bord+1),MIN0(imax,id_bord)
               mask_bord_north_tmp(iborder,i,i3(is))=1
             END DO
           END IF
         ELSE IF (rose(is) == 'E') THEN
            ! verification
           IF_MPI (MASTER) THEN
           IF(is > 1) THEN
             IF(rose(ism1) /= 'E' .AND. (i1(is) .NE. i2(ism1) .OR. i2(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of open border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment East number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i2 of segment',is,'must be = to i2 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
             ENDIF
           ENDIF
           ENDIF_MPI
           IF (i1(is) >= imin .and.i1(is) <= imax) THEN
             jb_bord=MIN0(i3(is),i2(is))
             jh_bord=MAX0(i3(is),i2(is))
             DO j=MAX0(jmin,jb_bord+1),MIN0(jmax,jh_bord)
               mask_bord_east_tmp(iborder,i1(is),j)=-1
             END DO
           END IF
         ELSE IF (rose(is) == 'W') THEN
            ! verification
           IF_MPI (MASTER) THEN
           IF(is > 1) THEN
             IF(rose(ism1) /= 'W' .AND. (i1(is) .NE. i2(ism1) .OR. i2(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of open border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment West number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i2 of segment',is,'must be = to i2 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
             ENDIF
           ENDIF
           ENDIF_MPI
           IF (i1(is)-1 >= imin .and.i1(is)-1 <= imax) THEN
             jb_bord=MIN0(i3(is),i2(is))
             jh_bord=MAX0(i3(is),i2(is))
             DO j=MAX0(jmin,jb_bord+1),MIN0(jmax,jh_bord)
               mask_bord_east_tmp(iborder,i1(is),j)=1
             END DO
           END IF
         END IF
       END DO
       DEALLOCATE(i1,i2,i3,rose)
       READ(49,*)

     END DO

     DO iborder=nb_border_open+1,nb_border
       mask_bud_S_tmp(:,:)=0
       mask_bud_N_tmp(:,:)=0
       mask_bud_W_tmp(:,:)=0
       mask_bud_E_tmp(:,:)=0
       READ(49,*,iostat=eof)ibord_r
       READ(49,'(a)',iostat=eof) nombord_r
       ind_white=INDEX(nombord_r,' ')
       name_bord(iborder)=TRIM(ADJUSTL(ADJUSTR(nombord_r(1:ind_white))))
       READ(49,*)l_border_close
       izon_bud=izon_bud+1
       iborder_bud(izon_bud)=iborder
       READ(49,*)nb_seg
       ALLOCATE(i1(nb_seg),i2(nb_seg),i3(nb_seg))
       ALLOCATE(rose(nb_seg))
       DO is=1,nb_seg
         READ(49,*)i1(is),i2(is),i3(is),rose(is)
       ENDDO
       DO is=1,nb_seg
         ism1=is-1
         IF(is==1)ism1=nb_seg
         isp1=is+1
         IF(is==nb_seg)isp1=1
         
         IF (rose(is) == 'S') THEN
            ! verification
           IF_MPI (MASTER) THEN
           IF(rose(ism1) /= 'S' .AND. (i1(is) .NE. i1(ism1) .OR. i3(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of closed border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment SOUTH number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i3 of segment',is,'must be = to i1 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
           ENDIF
           ENDIF_MPI
           IF (i3(is)-1 >= jmin .and.i3(is)-1 <= jmax) THEN
             ig_bord=MIN0(i1(is),i2(is))
             id_bord=MAX0(i1(is),i2(is))
             DO i=MAX0(imin,ig_bord+1),MIN0(imax,id_bord)
               mask_bord_north_tmp(iborder,i,i3(is))=1
               DO j=jmin,i3(is)
                 mask_bud_S_tmp(i,j)=mask_bud_S_tmp(i,j)-1
               END DO
               DO j=i3(is)+1,jmax
                 mask_bud_S_tmp(i,j)=mask_bud_S_tmp(i,j)+1
               END DO
             END DO
           END IF
         ELSE IF (rose(is) == 'N') THEN
           ! verification
           IF_MPI (MASTER) THEN
           IF(rose(ism1) /= 'N' .AND. (i1(is) .NE. i1(ism1) .OR. i3(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of closed border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment NORD number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i3 of segment',is,'must be = to i1 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
           ENDIF
           ENDIF_MPI
           IF (i3(is) >= jmin .and.i3(is) <= jmax) THEN
             ig_bord=MIN0(i1(is),i2(is))
             id_bord=MAX0(i1(is),i2(is))
             DO i=MAX0(imin,ig_bord+1),MIN0(imax,id_bord)
               mask_bord_north_tmp(iborder,i,i3(is))=-1    
               DO j=jmin,i3(is)
                 mask_bud_N_tmp(i,j)=mask_bud_N_tmp(i,j)+1
               END DO
               DO j=i3(is)+1,jmax
                 mask_bud_N_tmp(i,j)=mask_bud_N_tmp(i,j)-1
               END DO
             END DO
           END IF
         ELSE IF (rose(is) == 'E') THEN
            ! verification
           IF_MPI (MASTER) THEN
           IF(rose(ism1) /= 'E' .AND. (i1(is) .NE. i2(ism1) .OR. i2(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of closed border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment East number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i2 of segment',is,'must be = to i2 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
           ENDIF
           ENDIF_MPI
           IF (i1(is) >= imin .and.i1(is) <= imax) THEN
             jb_bord=MIN0(i3(is),i2(is))
             jh_bord=MAX0(i3(is),i2(is))
             DO j=MAX0(jmin,jb_bord+1),MIN0(jmax,jh_bord)
               mask_bord_east_tmp(iborder,i1(is),j)=-1
               DO i=imin,i1(is)
                 mask_bud_E_tmp(i,j)=mask_bud_E_tmp(i,j)+1
               END DO
               DO i=i1(is)+1,imax
                 mask_bud_E_tmp(i,j)=mask_bud_E_tmp(i,j)-1
               END DO
             END DO
           END IF
         ELSE IF (rose(is) == 'W') THEN
            ! verification
           IF_MPI (MASTER) THEN
           IF(rose(ism1) /= 'W' .AND. (i1(is) .NE. i2(ism1) .OR. i2(is) .NE. i3(ism1))) THEN
               write(ierrorlog,*)'Definition of closed border number :',iborder,' name ',TRIM(name_bord(iborder))
               write(ierrorlog,*)'segment West number :',is,'i1,i2,i3 ',i1(is),i2(is),i3(is)
               write(ierrorlog,*)'is not contiguous to the previous one whom i1,i2,i3 are :', i1(ism1),i2(ism1),i3(ism1)
               write(ierrorlog,*)'i1 and i2 of segment',is,'must be = to i2 and i3 respectively of previous segment'
               CALL_MPI MPI_FINALIZE(IERR_MPI)
               STOP
           ENDIF
           ENDIF_MPI
           IF (i1(is)-1 >= imin .and.i1(is)-1 <= imax) THEN
             jb_bord=MIN0(i3(is),i2(is))
             jh_bord=MAX0(i3(is),i2(is))
             DO j=MAX0(jmin,jb_bord+1),MIN0(jmax,jh_bord)
               mask_bord_east_tmp(iborder,i1(is),j)=1
               DO i=imin,i1(is)
                 mask_bud_W_tmp(i,j)=mask_bud_W_tmp(i,j)-1
               END DO
               DO i=i1(is)+1,imax
                 mask_bud_W_tmp(i,j)=mask_bud_W_tmp(i,j)+1
               END DO
             END DO
           END IF
         END IF
       END DO
       DO i=imin,imax
         DO j=jmin,jmax
           IF (mask_bud_W_tmp(i,j)+mask_bud_E_tmp(i,j)==2  .and. &
               mask_bud_N_tmp(i,j)+mask_bud_S_tmp(i,j)==2) THEN
              mask_bud_in_tmp(izon_bud,i,j)=1
           ELSE
             mask_bud_in_tmp(izon_bud,i,j)=0
           END IF
         END DO
       END DO
       DEALLOCATE(i1,i2,i3,rose)

       READ(49,*)

     END DO
   END IF

  !! computing total surface and volume of each budget zone 
   DO is=1,nb_zone_budget
      ib=is+nb_border_open
      surf_tot(is)=0._rsh
      vol_tot(is)=0._rsh
       DO i=limin,limax
         DO j=ljmin,ljmax
           IF (h0(i,j).NE.-valmanq) THEN
             surf_tot(is)=surf_tot(is)+ mask_bud_in_tmp(is,i,j)*surf(i,j)
             vol_tot(is)=vol_tot(is)+ mask_bud_in_tmp(is,i,j)*surf(i,j)*(h0(i,j)+xe3(i,j))
           END IF
         END DO
       END DO
       CALL_MPI ADD_MPI_RSH(surf_tot(is),surftotale(is))
       CALL_MPI ADD_MPI_RSH(vol_tot(is),voltotal(is))
       IF_MPI (MASTER) write(*,*)'zone,',is,TRIM(ADJUSTL(ADJUSTR(name_bord(ib))))
       IF_MPI (MASTER) write(*,*)'total surface zone,',is,limin,limax,ljmin,ljmax,surf_tot(is),surftotale(is)
       IF_MPI (MASTER) write(*,*)'total volume zone',is,limin,limax,ljmin,ljmax,vol_tot(is),voltotal(is)
   END DO
            
#ifdef key_MPI_2D
   !  one keeps only mask into the zone corresponding to each processor
   mask_bord_east(:,liminm1:limaxp2,ljminm1:ljmaxp2) = mask_bord_east_tmp(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   mask_bord_north(:,liminm1:limaxp2,ljminm1:ljmaxp2) = mask_bord_north_tmp(:,liminm1:limaxp2,ljminm1:ljmaxp2)
   IF (nb_zone_budget > 0) THEN
      mask_bud_in(:,liminm1:limaxp2,ljminm1:ljmaxp2) = mask_bud_in_tmp(:,liminm1:limaxp2,ljminm1:ljmaxp2)   
   ENDIF
   ! writing in a file for verification of mask
   IF_MPI (MASTER) THEN

#else
   mask_bord_east(:,:,:) = mask_bord_east_tmp(:,:,:)
   mask_bord_north(:,:,:) = mask_bord_north_tmp(:,:,:)
   IF (nb_zone_budget > 0) THEN
     mask_bud_in(:,:,:) = mask_bud_in_tmp(:,:,:)
     surftotale(:)=surf_tot(:)   
     voltotal(:)=vol_tot(:)   
   ENDIF
#endif
   OPEN(9999,file='ficverif_zone_budget')
   WRITE(9999,*)'Verification of masks for budget'
   DO is=1,nb_zone_budget
     ib=is+nb_border_open
     WRITE(9999,*) 
     WRITE(9999,*)'zone ',is
     WRITE(9999,*)'nom ',TRIM(ADJUSTL(ADJUSTR(name_bord(ib))))
     WRITE(9999,*)'total surface  (m2) ',surftotale(is)
     WRITE(9999,*)'total initial volume  (m3) ',voltotal(is)
     WRITE(9999,*)'mask =1 if cell into the budget zone' 
     WRITE(9999,*)'mask =0 if cell out of the budget zone' 
     WRITE(9999,101)(i,i=imin,imax)
     DO j=jmax,jmin,-1
       WRITE(9999,100)j,(mask_bud_in_tmp(is,i,j),i=imin,imax)
     END DO
     WRITE(9999,*)'mask W and E' 
     WRITE(9999,*)'mask =1 if border cell W' 
     WRITE(9999,*)'mask =-1 if border cell E' 
     WRITE(9999,*)'mask =0 if not a border cell' 
     iborder=iborder_bud(is)
     WRITE(9999,101)(i,i=imin,imax)
     DO j=jmax,jmin,-1
       WRITE(9999,100)j,(mask_bord_east_tmp(iborder,i,j),i=imin,imax)
     END DO
     WRITE(9999,*)'mask S and N' 
     WRITE(9999,*)'mask =1 if border cell S' 
     WRITE(9999,*)'mask =-1 if border cell N' 
     WRITE(9999,*)'mask =0 if not a border cell' 
     WRITE(9999,101)(i,i=imin,imax)
     DO j=jmax,jmin,-1
       WRITE(9999,100)j,(mask_bord_north_tmp(iborder,i,j),i=imin,imax)
     END DO
   END DO
  100 FORMAT(i3,': ',1000i3)
  101 FORMAT(' i=  ',1000i3)
   CLOSE(9999)

   ENDIF_MPI

#ifdef key_contaminant
   ALLOCATE(flx_reac_zone(nb_zone_budget,nreacmax,nb_var))
   flx_reac_zone(:,:,:)=0._rlg
   ALLOCATE(nre_act(nv_adv))
   ALLOCATE(ire_act(nv_adv,nreacmax))
   nre_act(:)=0
   ire_act(:,:)=0
   DO iv=1,nv_adv
    nre_act(iv)=0
    DO ire=1,nreacmax
     DO ivr=1,nb_var_re(ire)
       IF(iv == irkm_var_re(ire,ivr))then
         nre_act(iv)=nre_act(iv)+1
         ire_act(iv,nre_act(iv))=ire
       ENDIF
     ENDDO
    ENDDO
   ENDDO
   nreactmax=maxval(nre_act)
   ALLOCATE(flx_reac_total(nb_zone_budget,nv_adv,nreactmax))
#endif

  END SUBROUTINE sub_budget_readdomain

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

  SUBROUTINE sub_budget_flxcum(cv_wat,xe3)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_budget_flxcum  ***
   !&E
   !&E ** Purpose : flux cumulating threw each border line defined by the user 
   !&E               for dissolved substances
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step3dxy, step3dyx
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2010-07 (B. Thouvenin) based on sed_massbalance 
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2014-12 (B.Thouvenin) modifications to use OPENMP
   !&E       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate) - 
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,   ONLY : kmax,valmanq
   USE comsubstance, ONLY : fx,fy,phicon,phicon_atm,l_subs2D
#ifdef key_sedim_MUSTANG
   USE comsubstance, ONLY : typdiss
   USE comMUSTANG,   ONLY : flx_s2w,flx_w2s,flx_w2s_sum
   USE comvars2d,    ONLY : halfdt
#endif
   USE comvars2d,    ONLY : t,h0,dx_v,dy,dx,dy_u,surf,hm,l_tssubadv_over_smalldepth,h0fond
   USE comvars3d,    ONLY : dsigu

   !! * Arguments
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in)    :: xe3
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_wat

   !! * Local declarations
   INTEGER                                 :: iv,i,j,k,iz,iborder,ib
   INTEGER :: IERR_MPI
   REAL(KIND=rsh)                          :: fxij,fyij,d3,phi

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

!$OMP PARALLEL DEFAULT(NONE) SHARED(ljmin,ljmax,limin,limax,h0,xe3,hm,h0fond       &
!$OMP& ,nvp,nv_adv,fx,fy,dsigu,dy_u,dx_v,nb_border_open,flx_obc_cum,flx_rej_cum    &
#ifdef key_sedim_MUSTANG
!$OMP& ,flx_ws_cum,halfdt,surf,flx_s2w,flx_w2s,flx_w2s_sum,typdiss                 &
#endif
!$OMP& ,flx_obczon_cum,flx_atm_cum,mask_bord_north,mask_bord_east,nb_zone_budget   &
!$OMP& ,iborder_bud,mask_bud_in,phicon,phicon_atm,cv_wat                           &
!$OMP& ,l_tssubadv_over_smalldepth,l_subs2D) &
!$OMP& PRIVATE(i,j,k,ib,d3,iv,fxij,fyij,iborder,iz,phi)  

!$OMP DO REDUCTION(+:flx_obc_cum,flx_rej_cum,flx_obczon_cum    &
#ifdef key_sedim_MUSTANG
!$OMP& ,flx_ws_cum        &
#endif
!$OMP& ,flx_atm_cum)
   DO j=ljmin,ljmax
     DO i=limin,limax
       IF (h0(i,j).NE.-valmanq) THEN
         d3=h0(i,j)+xe3(i,j)
         IF(d3 > hm) THEN

           DO k=1,kmax
             DO iv=1,nv_adv
               IF (l_subs2D(iv)) THEN
                   IF (k==1) THEN
                     phi =SUM(phicon(iv,:,i,j))
                     ! sur toute la hauteur d eau
                     fxij=fx(iv,k,i,j)*dy_u(i,j)
                     fyij=fy(iv,k,i,j)*dx_v(i,j)
                   ELSE
                     fxij=0.0_rsh
                     fyij=0.0_rsh
                     phi=0.0_rsh
                   ENDIF
               ELSE
                   fxij=fx(iv,k,i,j)*dsigu(k)*dy_u(i,j)
                   fyij=fy(iv,k,i,j)*dsigu(k)*dx_v(i,j)
                   phi=REAL(phicon(iv,k,i,j),rlg)
               ENDIF
               DO ib=1,nb_border_open
                 flx_obc_cum(ib,iv)=flx_obc_cum(ib,iv)   &
                  + mask_bord_north(ib,i,j)*fyij  &
                  + mask_bord_east(ib,i,j)*fxij 
               END DO
               DO iz=1,nb_zone_budget
                 iborder=iborder_bud(iz)
                 flx_rej_cum(iz,iv)=flx_rej_cum(iz,iv)   &
                  + mask_bud_in(iz,i,j)*REAL(phi,rlg)
                 flx_obczon_cum(iz,iv)=flx_obczon_cum(iz,iv)   &
                  + mask_bord_north(iborder,i,j)*REAL(fyij,rlg)  &
                  + mask_bord_east(iborder,i,j)*REAL(fxij,rlg) 
               END DO
             END DO     
           END DO
       
           DO iv=1,nv_adv
             DO iz=1,nb_zone_budget
               flx_atm_cum(iz,iv)=flx_atm_cum(iz,iv)  &
                  + mask_bud_in(iz,i,j)*REAL(phicon_atm(iv,i,j),rlg)
#ifdef key_sedim_MUSTANG
              ! appelle avant effdep et apres step3d.. ==>> flx_w2s(iv,i,j) exprime en m/s
              !     mais flx_s2w(iv,i,j) en kg/m2/s
                  flx_ws_cum(iz,iv)=flx_ws_cum(iz,iv)                       &
                   + mask_bud_in(iz,i,j)*REAL(surf(i,j),rlg)*         &
              ! for dissolved substances : flx_w2s explicit and included in flx_s2w (M/m2/s)
              ! for particulate substances : flx_w2s cumule sur ndt_part sous pas de temps : flx_w2s_sum (M/m2)  
              !                              and   flx_s2w (M/m2/s)
                   (halfdt*REAL(flx_s2w(iv,i,j),rlg)-(1-typdiss(iv))*REAL(flx_w2s_sum(iv,i,j),rlg))
#endif
             END DO     
           END DO

         ELSE
           IF (l_tssubadv_over_smalldepth .AND. d3 > h0fond) THEN
              DO iv=1,nv_adv
                IF (l_subs2D(iv)) THEN
                   phi =SUM(phicon(iv,:,i,j))
                   fxij=fx(iv,1,i,j)*dy_u(i,j)
                   fyij=fy(iv,1,i,j)*dx_v(i,j)
                ELSE
                   fxij=SUM(fx(iv,:,i,j)*dsigu(:))*dy_u(i,j)
                   fyij=SUM(fy(iv,:,i,j)*dsigu(:))*dx_v(i,j)
                   phi =SUM(phicon(iv,:,i,j))
                ENDIF
                DO ib=1,nb_border_open
                  flx_obc_cum(ib,iv)=flx_obc_cum(ib,iv)   &
                   + mask_bord_north(ib,i,j)*fyij  &
                   + mask_bord_east(ib,i,j)*fxij 
                END DO
                DO iz=1,nb_zone_budget
                  iborder=iborder_bud(iz)
                  flx_rej_cum(iz,iv)=flx_rej_cum(iz,iv)   &
                   + mask_bud_in(iz,i,j)*REAL(phi,rlg)
                  flx_obczon_cum(iz,iv)=flx_obczon_cum(iz,iv)   &
                   + mask_bord_north(iborder,i,j)*REAL(fyij,rlg)  &
                   + mask_bord_east(iborder,i,j)*REAL(fxij,rlg) 
                  flx_atm_cum(iz,iv)=flx_atm_cum(iz,iv)  &
                   + mask_bud_in(iz,i,j)*REAL(phicon_atm(iv,i,j),rlg)
#ifdef key_sedim_MUSTANG
                   flx_ws_cum(iz,iv)=flx_ws_cum(iz,iv)                       &
                      + mask_bud_in(iz,i,j)*REAL(surf(i,j),rlg)*             &
              ! dissolved substances : flx_w2s explicit and included in flx_s2w
              ! particulate substances  : flx_w2s explicit and included in flx_s2w
                     (halfdt*REAL(flx_s2w(iv,i,j),rlg)-(1-typdiss(iv))*REAL(flx_w2s_sum(iv,i,j),rlg))
#endif
                END DO     
              END DO
           
           END IF
         END IF
       END IF

     END DO
   END DO
!$OMP END DO
!$OMP END PARALLEL
   
  END SUBROUTINE sub_budget_flxcum

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

  SUBROUTINE sub_budget_out(tect,concwat   &
#ifdef key_sedim_MUSTANG
                                        ,concsed    &
#endif
                                                ,xeloc)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_budget_out  ***
   !&E
   !&E ** Purpose : mass budget computation and writing on output files
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : step
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2010-07 (B. Thouvenin) based on sed_massbalance 
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,   ONLY : nb_var,kmax,ljmin,ljmax,valmanq
   USE comsubstance, ONLY : nvp,nvpc,nv_tot,l_subs2d
   USE comvars2d,    ONLY : t,dt,h0,surf,dx_v,dy,hx,hy
   USE comvars3d,    ONLY : dsigu
#ifdef key_sedim_MUSTANG
   USE parameters,   ONLY : ksdmin,ksdmax
   USE comMUSTANG,   ONLY : ksmi,ksma,dzs,poro
#endif
#if defined key_siggen || defined key_gencoord
   USE toolgeom,    ONLY : f_hzu
#endif
   USE_MPI toolmpi,  ONLY : MASTER,ADD_MPI_RLG

   !! * Arguments
   REAL(KIND=rlg),INTENT(IN)                                         :: tect
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(IN) :: concwat
#ifdef key_sedim_MUSTANG
   REAL(KIND=rsh),DIMENSION(-1:nv_tot,ksdmin:ksdmax,limin:limax,ljmin:ljmax),INTENT(IN) :: concsed
#endif
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT(IN)                     :: xeloc

   !! * Local declarations
   INTEGER                                 :: iv,i,j,k,ksmin,ksmax,iz,iborder,ib
   INTEGER                                 :: IERR_MPI
   REAL(KIND=rlg)                          :: voltot
#ifdef key_sedim_MUSTANG
   REAL(KIND=rlg),DIMENSION(ksdmin:ksdmax) :: vol
#endif
#ifdef key_contaminant
   INTEGER                                 :: ire
#endif

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

   t_out_budget=t_out_budget+dtout_budget

!$OMP PARALLEL DEFAULT(NONE) SHARED(ljmin,ljmax,limin,limax,nv_adv,nv_fix,h0,xeloc       &
!$OMP&  ,t_out_budget,stok_wat,stok_wat_fix,dsigu,surf,fwet,mask_bud_in,concwat,l_subs2D &
#ifdef key_sedim_MUSTANG
!$OMP& ,stok_sed,ksmi,ksma,dzs,nvp,concsed,poro                         &
#endif
!$OMP&  ,nb_zone_budget)                         &
!$OMP& PRIVATE(i,j,k,iv,iz                       &
#ifdef key_sedim_MUSTANG
!$OMP& ,ksmin,ksmax,vol                          &
#endif
!$OMP& ,voltot)  

   ! stocks and budgets and sinksources
!$OMP SINGLE
   stok_wat(:,1:nv_adv)=0.0_rlg
   IF(nv_fix>0)stok_wat_fix(:,1:nv_fix)=0.0_rsh
#ifdef key_sedim_MUSTANG
   stok_sed(:,1:nv_adv)=0.0_rlg
#endif
!$OMP END SINGLE

!$OMP DO REDUCTION(+:stok_wat    &
#ifdef key_sedim_MUSTANG
!$OMP& ,stok_sed                &
#endif
!$OMP& ,stok_wat_fix)
   DO j=ljmin,ljmax
     DO i=limin,limax

       IF (h0(i,j).NE.-valmanq) THEN
         DO k=1,kmax
#if defined key_siggen || defined key_gencoord
           voltot=REAL(dsigu(k)*f_hzu(h0(i,j),xeloc(i,j),k,i,j)*surf(i,j)*fwet(i,j),rlg)
#else
           voltot=REAL(dsigu(k)*(h0(i,j)+xeloc(i,j))*surf(i,j)*fwet(i,j),rlg)
#endif

           DO iv=1,nv_adv
             DO iz=1,nb_zone_budget
                 IF (l_subs2D(iv)) THEN
                     IF(k==1)stok_wat(iz,iv)= &
                        stok_wat(iz,iv)+mask_bud_in(iz,i,j)*  &
                             REAL(concwat(iv,k,i,j)*(h0(i,j)+xeloc(i,j))*surf(i,j)*fwet(i,j),rlg)
                 ELSE
                     stok_wat(iz,iv)=stok_wat(iz,iv)+mask_bud_in(iz,i,j)*  &
                                    REAL(concwat(iv,k,i,j)*voltot,rlg)
                 ENDIF
             END DO     
           END DO     
         END DO     
         IF(nv_fix>0) THEN
            DO iv=1,nv_fix
             DO iz=1,nb_zone_budget
               stok_wat_fix(iz,iv)=stok_wat_fix(iz,iv)+mask_bud_in(iz,i,j)*  &
                                  REAL(concwat(nv_adv+iv,1,i,j)*surf(i,j),rlg)
             END DO
            END DO     
         ENDIF

#ifdef key_sedim_MUSTANG
         ksmin=ksmi(i,j)
         ksmax=ksma(i,j)
         DO k=ksmin,ksmax
           vol(k)=REAL(surf(i,j)*dzs(k,i,j),rlg)
           DO iv=1,nvp
             DO iz=1,nb_zone_budget
               stok_sed(iz,iv)=stok_sed(iz,iv)+mask_bud_in(iz,i,j)*  &
                               REAL(concsed(iv,k,i,j),rlg)*vol(k)
             END DO
           END DO
           DO iv=nvp+1,nv_adv
             DO iz=1,nb_zone_budget
               stok_sed(iz,iv)=stok_sed(iz,iv)+mask_bud_in(iz,i,j)*  &
                               REAL(concsed(iv,k,i,j)*poro(k,i,j),rlg)*vol(k)
             END DO
           END DO
         END DO    
#endif
       END IF

     END DO
   END DO
!$OMP END DO
!$OMP END PARALLEL


#ifdef key_MPI_2D
   DO iv=1,nv_adv
     DO iz=1,nb_zone_budget
       CALL_MPI ADD_MPI_RLG(stok_wat(iz,iv),stokw_total(iz,iv))
       CALL_MPI ADD_MPI_RLG(flx_rej_cum(iz,iv),flxrej_total(iz,iv))
       CALL_MPI ADD_MPI_RLG(flx_atm_cum(iz,iv),flxatm_total(iz,iv))
       CALL_MPI ADD_MPI_RLG(flx_obczon_cum(iz,iv),flxobczon_total(iz,iv))
#ifdef key_contaminant
       DO ire=1,nre_act(iv)
          CALL_MPI ADD_MPI_RLG(flx_reac_zone(iz,ire_act(iv,ire),iv),flx_reac_total(iz,iv,ire))
       END DO
#endif
     END DO
     DO ib=1,nb_border_open
       CALL_MPI ADD_MPI_RLG(flx_obc_cum(ib,iv),flxobc_total(ib,iv))
     END DO
   END DO
#ifdef key_sedim_MUSTANG
   DO iv=1,nv_adv
     DO iz=1,nb_zone_budget
        CALL_MPI ADD_MPI_RLG(flx_ws_cum(iz,iv),flxws_total(iz,iv))
        CALL_MPI ADD_MPI_RLG(stok_sed(iz,iv),stoks_total(iz,iv))
     END DO
   END DO
#endif
 IF (nv_fix > 0) THEN
   DO iv=1,nv_fix
     DO iz=1,nb_zone_budget
      CALL_MPI ADD_MPI_RLG(stok_wat_fix(iz,iv),stokwfix_total(iz,iv))
      IF(stokwfix_total(iz,iv) < 1.e-14)stokwfix_total(iz,iv)=0.0_rlg
     ENDDO
   ENDDO
 ENDIF
#endif

   IF_MPI (MASTER) THEN
     CALL sub_budget_output(tect)
   ENDIF_MPI

  END SUBROUTINE sub_budget_out


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

  SUBROUTINE sub_budget_output(tect)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE sub_budget_output  ***
   !&E
   !&E ** Purpose : mass budget writing into output files
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : sub_budget_out 
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2010-07 (B. Thouvenin) based on sed_massbalance 
   !&E       !  2015-11 (E. Schulz,B. Thouvenin) write name area in result file and time from 01/01/1900 
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : iscreenlog,ierrorlog,iwarnlog,out_dir,suffix
   USE comsubstance, ONLY : name_var,irk_fil,iout_budget
#ifdef key_contaminant
   USE contareaddat,  ONLY : name_reac
#endif
   USE_MPI toolmpi,  ONLY : MASTER

   !! * Arguments
   REAL(KIND=rlg),INTENT(IN)                                         :: tect

   !! * Local declarations
   INTEGER                                 :: iz,inumfic,i,iv,lnam,ifois
   INTEGER                                 :: nfoiv,ifoiv,iv1,iv2
#ifdef key_contaminant
   INTEGER                                 :: ire
#endif
   CHARACTER(LEN=2)                        :: name_zone
   CHARACTER(LEN=1)                        :: name_fois
   CHARACTER(LEN=lchain)                   :: fileoutbal,fileresu1              
   CHARACTER(LEN=1000)                     :: comment1,comment2              
   REAL(KIND=rlg)                          :: tec_h

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

          !! one file for each zone of budget (stock, and fluxes
  !! one file for all the borders (fluxes only)

  !IF ((tect-tdebsubs) < dtout_budget ) THEN
  IF (iout_budget == -1 ) THEN
     iout_budget=0
  !! first create output files
     inumfic=589
     DO iz=1,nb_zone_budget
       IF (iz <10) THEN
         WRITE(name_zone,200)iz
       ELSE
         WRITE(name_zone,201)iz
       END IF
 200   FORMAT(i1)
 201   FORMAT(i2)
 300   FORMAT(a)
       nfoiv=nv_adv/10+1
       DO ifoiv=1,nfoiv
         inumfic=inumfic+1
         WRITE(name_fois,200)ifoiv
         fileresu1=fileresu(1)
         IF (l_filebydate(1)) THEN
           fileresu1=TRIM(out_dir)//TRIM(fileresu1)//'_'//TRIM(suffix)//'.nc'
         ENDIF
         lnam=index(fileresu1,'.nc')-1
         file_budget_out(iz)=TRIM(fileresu1(1:lnam))//'_bil'//TRIM(name_zone)//'_'//TRIM(name_fois)//'.csv'
         fileoutbal=file_budget_out(iz)
         IF_MPI (MASTER) THEN
         WRITE(iscreenlog,*) 
         WRITE(iscreenlog,*) 'name of output file for budgets into zone ',iz
         WRITE(iscreenlog,*) TRIM(fileoutbal)
         WRITE(iscreenlog,*)
         ENDIF_MPI
       
         OPEN(inumfic,file=fileoutbal,form='formatted', access='sequential')
         WRITE(inumfic,*)'RESULTS FOR ZONE :;',iz,';',TRIM(ADJUSTL(ADJUSTR(name_bord(nb_border_open+iz))))
         WRITE(inumfic,*)'Bud   = global budget (would must be constant if conservative variable)'
         WRITE(inumfic,*)'StWat   = Stock Water'
         WRITE(inumfic,*)'Flx obc = FLux threw obc zone (>if in)'
         WRITE(inumfic,*)'Flx IN  = FLux river & discharge & atmosphere'
         WRITE(inumfic,*)'Flx atm  = FLux atmosphere only'
#ifdef key_sedim_MUSTANG
         WRITE(inumfic,*)'StSed   = Stock Sediment'
         WRITE(inumfic,*)'Flx WS  = FLux interface water/sediment'
#endif
#ifdef key_contaminant
         WRITE(inumfic,*)'Flre_name_reac   = Reactif Flux for process name_reac '
#endif
         comment1='time;'
         comment2='hours since 1900-01-01 00:00:00;'
        ! DO iv=1,nv_adv 
         iv1=(ifoiv-1)*10+1
         iv2=min(ifoiv*10,nv_adv)
         DO iv=iv1,iv2
           ifois=5
#ifdef key_sedim_MUSTANG
           ifois=ifois+2
#endif
#ifdef key_contaminant
           ifois=ifois+nre_act(iv)
#endif
           DO i=1,ifois
             comment1=TRIM(ADJUSTL(ADJUSTR(comment1)))//TRIM(ADJUSTL(ADJUSTR(name_var(irk_fil(iv))))) //';'
           END DO 
           comment2=TRIM(ADJUSTL(ADJUSTR(comment2)))//'Bud;StWat;'//&
#ifdef key_sedim_MUSTANG
           'StSed; Flx WS;'// &
#endif
           'Flx obc; Flx IN;Flx atm;'
#ifdef key_contaminant
           DO ire=1,nre_act(iv)
             comment2=TRIM(ADJUSTL(ADJUSTR(comment2)))//'FLre_'//TRIM(ADJUSTL(ADJUSTR(name_reac(ire_act(iv,ire)))))//';'
           ENDDO
#endif
         END DO
         DO iv=1,nv_fix
           comment1=TRIM(ADJUSTL(ADJUSTR(comment1)))//TRIM(ADJUSTL(ADJUSTR(name_var(irk_fil(nv_adv+iv))))) //';'
           comment2=TRIM(ADJUSTL(ADJUSTR(comment2)))//'StWat;'
         ENDDO
         WRITE(inumfic,300) TRIM(ADJUSTL(ADJUSTR(comment1)))
         WRITE(inumfic,300) TRIM(ADJUSTL(ADJUSTR(comment2)))
       END DO
     END DO
  
     IF (nb_border_open > 0) THEN
       inumfic=990
       IF (nb_zone_budget==0) THEN
         fileresu1=fileresu(1)
         IF (l_filebydate(1)) THEN
           fileresu1=TRIM(out_dir)//TRIM(fileresu1)//'_'//TRIM(suffix)//'.nc'
         ENDIF
         lnam=index(fileresu1,'.nc')-1
       ENDIF
       write(*,*)'nom du file_border :', lnam,fileresu1
       DO iz=1,nb_border_open
         IF (iz <10) THEN
           WRITE(name_zone,200)iz
         ELSE
           WRITE(name_zone,201)iz
         END IF
         inumfic=inumfic+1
         file_border_out(iz)=TRIM(fileresu1(1:lnam))//'_border'//TRIM(name_zone)//'.csv'
         fileoutbal=file_border_out(iz)
         write(*,*)'nom du file_border :', fileoutbal
         IF_MPI (MASTER) THEN
           WRITE(iscreenlog,*) 
           WRITE(iscreenlog,*) 'name of output file for fluxes threw borders : ',TRIM(fileoutbal)
           WRITE(iscreenlog,*) 
         ENDIF_MPI

         OPEN(inumfic,file=fileoutbal,status='replace',form='formatted', access='sequential')
         WRITE(inumfic,*)'RESULTS of FLUXES THREW BORDERS '
         comment1='time;'
         comment2='hours since 1900-01-01 00:00:00;'
         DO iv=1,nv_adv 
           comment1=TRIM(ADJUSTL(ADJUSTR(comment1)))//TRIM(ADJUSTL(ADJUSTR(name_var(irk_fil(iv))))) //';'
           comment2=TRIM(ADJUSTL(ADJUSTR(comment2)))//TRIM(ADJUSTL(ADJUSTR(name_bord(iz)))) //';' 
         END DO 
         WRITE(inumfic,300) TRIM(ADJUSTL(ADJUSTR(comment1)))
         WRITE(inumfic,300) TRIM(ADJUSTL(ADJUSTR(comment2)))
       END DO ! iz=1,nb_border_open
     END IF 
   END IF
  
  !! general case : computing budget and writing into files
 
   bil(:,:)=0.0_rlg 
   inumfic=589
!  E. Schulz :  write time in hours since tref = 1900-01-01 00:00:00
!               to standardise post-processing and to avoid mistakes)
   tec_h=tect/3600.0_rlg
   nfoiv=nv_adv/10+1

! MPI : the MASTER knows the total computed from each processor
#ifdef key_MPI_2D
   DO iz=1,nb_zone_budget
     bil(iz,1:nv_adv)=stokw_total(iz,1:nv_adv)     &
#ifdef key_sedim_MUSTANG
                    + stoks_total(iz,1:nv_adv)     &
#endif
                    - flxobczon_total(iz,1:nv_adv) &
                    - flxrej_total(iz,1:nv_adv)
#ifdef key_contaminant
     DO iv=1,nv_adv
     DO ire=1,nre_act(iv)
          bil(iz,iv)=bil(iz,iv)-flx_reac_total(iz,iv,ire)
     ENDDO
     ENDDO
#endif
     DO ifoiv=1,nfoiv
       inumfic=inumfic+1
       iv1=(ifoiv-1)*10+1
       iv2=min(ifoiv*10,nv_adv)
!toseelater not possible to keep that we do not write the result ! 
#if ! defined key_gfortran
!!!#if ! defined key_agrif
       IF(nv_fix == 0 .or. iv2 < nv_adv ) THEN
         WRITE(inumfic,100) tec_h,( bil(iz,iv),stokw_total(iz,iv),      &
#ifdef key_sedim_MUSTANG
                        stoks_total(iz,iv),flxws_total(iz,iv),      &
#endif
                        flxobczon_total(iz,iv),flxrej_total(iz,iv), &
                        flxatm_total(iz,iv) &
#ifdef key_contaminant
                       ,(flx_reac_total(iz,iv,ire),ire=1,nre_act(iv))  &
#endif
                        ,iv=iv1,iv2)
       ELSE
          WRITE(inumfic,100) tec_h,(bil(iz,iv),stokw_total(iz,iv),      &
#ifdef key_sedim_MUSTANG
                        stoks_total(iz,iv),flxws_total(iz,iv),      &
#endif
                        flxobczon_total(iz,iv),flxrej_total(iz,iv), &
                        flxatm_total(iz,iv) &
#ifdef key_contaminant
                       ,(flx_reac_total(iz,iv,ire),ire=1,nre_act(iv))  &
#endif
                        ,iv=iv1,iv2),(stokwfix_total(iz,iv),iv=1,nv_fix)
       ENDIF
!!!#endif
#endif

     END DO
   END DO
        
   IF (nb_border_open > 0) THEN
     inumfic=990
     DO iz=1,nb_border_open
       inumfic=inumfic+1
       WRITE(inumfic,100)tec_h,(flxobc_total(iz,iv),iv=1,nv_adv)
     END DO 
   END IF

! not MPI 
#else
   DO iz=1,nb_zone_budget
     bil(iz,1:nv_adv)=stok_wat(iz,1:nv_adv)       &
#ifdef key_sedim_MUSTANG
                    + stok_sed(iz,1:nv_adv)       &
#endif
                    - flx_obczon_cum(iz,1:nv_adv) &
                    - flx_rej_cum(iz,1:nv_adv)
#ifdef key_contaminant
     DO iv=1,nv_adv
     DO ire=1,nre_act(iv)
          bil(iz,iv)=bil(iz,iv)-flx_reac_zone(iz,ire_act(iv,ire),iv)
          flx_reac_total(iz,iv,ire)=flx_reac_zone(iz,ire_act(iv,ire),iv)
     ENDDO
     ENDDO
#endif
     DO ifoiv=1,nfoiv
       inumfic=inumfic+1
       iv1=(ifoiv-1)*10+1
       iv2=min(ifoiv*10,nv_adv)
!toseelater not possible to keep that we do not write the result ! 
#if ! defined key_gfortran
#if ! defined key_agrif
       IF(nv_fix == 0 .or. iv2 < nv_adv ) THEN
         WRITE(inumfic,100) tec_h,((bil(iz,iv),stok_wat(iz,iv),        &
#ifdef key_sedim_MUSTANG
                        stok_sed(iz,iv),flx_ws_cum(iz,iv),         &
#endif
                        flx_obczon_cum(iz,iv),flx_rej_cum(iz,iv),  &
                        flx_atm_cum(iz,iv) &
#ifdef key_contaminant
                       ,(flx_reac_total(iz,iv,ire),ire=1,nre_act(iv))  &
#endif
                       ),iv=iv1,iv2)
       ELSE
          WRITE(inumfic,100) tec_h,((bil(iz,iv),stok_wat(iz,iv),      &
#ifdef key_sedim_MUSTANG
                        stok_sed(iz,iv),flx_ws_cum(iz,iv),      &
#endif
                        flx_obczon_cum(iz,iv),flx_rej_cum(iz,iv), &
                        flx_atm_cum(iz,iv) &
#ifdef key_contaminant
                       ,(flx_reac_total(iz,iv,ire),ire=1,nre_act(iv))  &
#endif
                        ),iv=iv1,iv2),((stok_wat_fix(iz,iv)),iv=1,nv_fix)
       ENDIF
#endif
#endif
     END DO
   END DO
  
   IF (nb_border_open > 0) THEN
     inumfic=990
     DO iz=1,nb_border_open
       inumfic=inumfic+1
       WRITE(inumfic,100)tec_h,(flx_obc_cum(iz,iv),iv=1,nv_adv)
     END DO 
   ENDIF
#endif


100 FORMAT(120(e17.10,1x,';',1x))
  
  END SUBROUTINE sub_budget_output

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

#endif

 END MODULE
