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


      SUBROUTINE sed_output(nb,tect,x,fileresu1)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sed_output  ***
!&E
!&E ** Purpose :
!&E
!&E ** Description :
!&E
!&E ** Called by : output.F90
!&E
!&E ** External calls :
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2010-07 (V. Garnier) Original code extracted from sedim
!&E       !  01-2011 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
!!   ------------
   USE parameters,  ONLY : lchain,imin,imax,jmin,jmax,ksdmin,ksdmax,rsh,rlg,riosh,  &
                           limin,limax,ljmin,ljmax,liminp1,limaxm1,ljminp1,ljmaxm1, &
                           liminm1,ljminm1,limaxp1,ljmaxp1,limaxp2,ljmaxp2,         &
                           rg_valmanq_io,valmanq,kmax
   USE ionc4,       ONLY : ionc4_createfile,ionc4_gatt,ionc4_createvar,ionc4_write_sig,&
                           ionc4_vatt_missing,ionc4_write_time,ionc4_write_xyt,&
                           ionc4_write_zxyt,ionc4_write_xy,ionc4_sync,ionc4_close
   USE comsubstance, ONLY : nv_adv,name_var,irk_fil,long_name_var,nvpc,unit_modif_mudbio_N2dw,nvp
   USE comvars2d,    ONLY : lon2d_g,lat2d_g,lon2d,lat2d,date_ref,h0,ig,id,jb,jh,iscreenlog,&
                            h0fond,out_dir,cdate,suffix,nbouc,name_out_h0,name_out_xe,     &
                            riog_valid_min_h0,riog_valid_max_h0,name_out_temp,name_out_sal,&
                            name_out_hx,name_out_hy,hx,hy
   USE toolmpi,  ONLY : ADD_TAG,IPROC,NPROCS,ADD_ALL_MPI_INT,MPI_COMM_MARS,MASTER
   USE outtools,     ONLY : out_write_xytvar,out_write_zxytvar,out_grid
   USE outnetcdf,    ONLY : out_netcdf_conv
   USE output,       ONLY : iextmin,iextmax,jextmin,jextmax,     &
                            igextmin,igextmax,jgextmin,jgextmax, &
                            stepi,stepj,stepk,l_filebydate,tsort,&
                            tfinsor,l_out_nc4par,l_out_pack
   USE comMUSTANG,    ONLY : dzs,cv_sed,tenfon,tenfonc,tenfonw,dzsmaxuni,l_dzsmaxuni,   &
                            ksmi,ksma,l_morphocoupl,filoutsed,ksmauni,hseduni,         &
                            name_out_eptot,name_out_nbniv,name_out_dzs,name_out_tenfon,&
                            name_out_tenfonc,name_out_tenfonw,l_outsed_saltemp,        &
                            riog_valid_min_eptot,riog_valid_max_eptot,riog_valid_min_nbniv, &
                            riog_valid_max_nbniv,riog_valid_min_dzs,riog_valid_max_dzs,     &
                            riog_valid_min_tenfon,riog_valid_max_tenfon,ep_nivsed_outp1,    &
                            nk_nivsed_out,ep_nivsed_out,nivsed_out,choice_nivsed_out,epmax_nivsed_out

   USE wave,         ONLY : riog_valid_min_wave_hs,riog_valid_max_wave_hs, &
                            riog_valid_min_wave_tp,riog_valid_max_wave_tp, &
                            name_out_wave_hs,name_out_wave_tp
   USE wave,         ONLY : name_wave_hs,name_wave_tp,wave_tp,wave_hs


   IMPLICIT NONE


   INCLUDE 'mpif.h'


!! * Arguments
   INTEGER       ,       INTENT(in) :: nb
   REAL(KIND=rlg),       INTENT(in) :: tect
   REAL(KIND=rsh),DIMENSION(     limin:limax,ljmin:ljmax),INTENT(in)     :: x
   CHARACTER(LEN=lchain),INTENT(in) :: fileresu1

!! * Local declarations
!!   ------------------
   INTEGER               :: i,j,k,l,lnam,kniv,finkniv,nk_nivsed_outlu
   INTEGER               :: grp_world,grp_zoom,comm_zoom,ncpu,ip,ierr_mpi,nv_use,iv0
   INTEGER, DIMENSION(:), ALLOCATABLE  :: numrank,listrank
   INTEGER, DIMENSION(1) :: listrank_nb
   CHARACTER(len=4)       :: nomc
   CHARACTER(len=lchain)      :: chaine
   CHARACTER(len=lchain)  :: name_out_cvsed
   CHARACTER(LEN=lchain) :: chain
   LOGICAL,         DIMENSION(limin:limax,ljmin:ljmax)    :: mask_h0
   REAL(KIND=rsh),  DIMENSION(limin:limax,ljmin:ljmax)    :: var2D
   REAL(KIND=rsh),  DIMENSION(limin:limax,ljmin:ljmax)    :: epsedtot
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)    :: h0_local,hx_local,hy_local
   REAL(kind=riosh)::        riog_valid_min_cvsed,riog_valid_max_cvsed
   REAL(KIND=riosh),DIMENSION(:,:,:), ALLOCATABLE         :: tecrit3D
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)    :: tecrit2D
   REAL(KIND=rsh),DIMENSION(:), ALLOCATABLE               :: niv
   REAL(KIND=rsh)                                         :: eptmp,eptmp1,cstmp,epintbottom,dzs_estim,unitmudbinv

!DATA ifois /-1/
!INTEGER :: ifois=-1
   
!! * Executable part
!!   ---------------

   !PRINT*, 'enter sed_output'

! Masks
! -----

   h0_local(:,:)=h0(limin:limax,ljmin:ljmax)
   mask_h0(:,:)=.TRUE.
   DO j=ljmin,ljmax
     DO i=MAX(limin,ig(j)),MIN(limax,id(j))
!IF ( h0_local(i,j)+x(i,j) > h0fond ) mask_h0(i,j)=.FALSE.
       IF ( h0_local(i,j) .NE. -valmanq .AND. ksma(i,j) > 0 ) mask_h0(i,j)=.FALSE.
     END DO
   END DO

   comm_zoom=0

   comm_zoom=MPI_COMM_MARS
   IF ( (nbouc == 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) THEN
! what a code ! Not a clever code probably but not clear in my head at the moment
! Count the number of processors active over the area to be saved
     ncpu=0
     IF ( jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb) ) ncpu=ncpu+1
     CALL ADD_ALL_MPI_INT(ncpu)
! Save the number of active processors
     ALLOCATE (numrank(ncpu))
     ALLOCATE (listrank(nprocs))
     numrank(:)=0   ! must be initialized if the proc 0 is active on the saved zoom
     listrank(:)=0
! which cpu is active ?
     IF ( jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb) ) listrank(1)=iproc
! all cpus get the number of active cpus plus sleeping cpus (set to 0)
     listrank_nb(1) = listrank(1)
     CALL MPI_ALLGATHER(listrank_nb(1),1,MPI_INTEGER,listrank,1,MPI_INTEGER,MPI_COMM_MARS,ierr_mpi)
! all cpus keep the active cpus only
     IF (ncpu/=nprocs) THEN
       ip=0
       DO i=1,nprocs
         IF ( listrank(i) /= 0 ) THEN
         ip=ip+1
         numrank(ip)=listrank(i)
         ENDIF
       END DO
! Group relative to the MPI_COMM_MARS
       CALL MPI_COMM_GROUP(MPI_COMM_MARS,grp_world,ierr_mpi)
! Create the group of active cpus
       CALL MPI_GROUP_INCL(grp_world,ncpu,numrank,grp_zoom,ierr_mpi)
! Create the group of active cpus
       CALL MPI_COMM_CREATE(MPI_COMM_MARS,grp_zoom,comm_zoom,ierr_mpi)
     ENDIF
   ENDIF

       IF(l_outsed_saltemp) THEN
         iv0=-1
       ELSE
         iv0=1
       ENDIF

       nv_use=nv_adv

    
!IF (ifois==-1 .AND. jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb)) THEN
   IF (nbouc==0 .AND. jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb)) THEN
    IF (MASTER) THEN
           WRITE(iscreenlog,*)
           WRITE(iscreenlog,*)'***************************************************************'
           WRITE(iscreenlog,*)'*****************     SEDIMENT OUTPUT  ************************'
           WRITE(iscreenlog,*)'***************************************************************'
           WRITE(iscreenlog,*)
    ENDIF
    nk_nivsed_outlu=nk_nivsed_out
    IF(choice_nivsed_out == 1 ) THEN
       ALLOCATE(nivsed_out(0:ksdmax+1))
       nivsed_out(0)=0
       DO k=ksdmin,ksdmax
             nivsed_out(k)=ksdmax+1-k
       END DO
       nivsed_out(ksdmax+1)=0
       IF (MASTER)THEN
           WRITE(iscreenlog,*)'results in sediment will be save on all the sediment layers (ksdmax)'
       ENDIF
    ELSE IF(choice_nivsed_out == 2 ) THEN
       IF(nk_nivsed_outlu > ksdmax) THEN
          nk_nivsed_out=ksdmax
       ENDIF
       ALLOCATE(nivsed_out(0:nk_nivsed_out+1))
       nivsed_out(0)=0
       DO k=1,nk_nivsed_out
         nivsed_out(k)=nk_nivsed_out+1-k
       END DO
       nivsed_out(nk_nivsed_out+1)=0       
       IF (MASTER) THEN
          WRITE(iscreenlog,*)'results in sediment will be save only for the first ', nk_nivsed_out, &
                              ' layers from sediment surface'
       ENDIF
    ELSE IF(choice_nivsed_out == 3 ) THEN
       dzs_estim=MIN(dzsmaxuni,hseduni/ksmauni)
       nk_nivsed_out = MIN(ksdmax, INT(epmax_nivsed_out / dzs_estim) +3) 
       ALLOCATE(nivsed_out(0:nk_nivsed_out+1))
       nivsed_out(0)=0
       DO k=1,nk_nivsed_out
         nivsed_out(k)=nk_nivsed_out+1-k
       END DO
       nivsed_out(nk_nivsed_out+1)=0
       IF (MASTER) THEN
           WRITE(iscreenlog,*)'results in sediment will be save from sediment surface till ',  &
                                epmax_nivsed_out*1000,'mmeters (integration below)'
           WRITE(iscreenlog,*)'the last layer saved in the sediment output file will be that  ',  &
                              'the bottom of which exceeds the desired maximum thickness '          
           WRITE(iscreenlog,*)'number of sediment layers in output file :',nk_nivsed_out,'INT(',epmax_nivsed_out,'/',dzs_estim,'+3)'
           IF (.NOT. l_dzsmaxuni) THEN
             WRITE(iscreenlog,*)
             WRITE(iscreenlog,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
             WRITE(iscreenlog,*)'ATTENTION : the maximum number of saved sed. layers has been evaluated from dzsmaxuni=',dzsmaxuni
             WRITE(iscreenlog,*)'            if dzsmax is not uniform, it could be smaller than dzsmaxuni and then there may be points where '
             WRITE(iscreenlog,*)'             this number is not sufficient to describe the maximum thickness'
             WRITE(iscreenlog,*)
             WRITE(iscreenlog,*)'This option is not recommended with l_dzsmaxuni=.False.'
             WRITE(iscreenlog,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
           ENDIF
       ENDIF
     ELSE IF(choice_nivsed_out == 4 ) THEN
       IF(nk_nivsed_outlu > 5 ) THEN
            IF (MASTER) THEN
               WRITE(iscreenlog,*)'ERROR in namelist namsedoutput (parasedim.txt) '
               WRITE(iscreenlog,*)'nk_nivsed_out must be <=5'
               WRITE(iscreenlog,*)'nk_nivsed_out is automatically set to 5 '
               WRITE(iscreenlog,*)'and an latter layer (6th) is added to integrate till the bottom sediment  '
            ENDIF
            nk_nivsed_outlu = 5
       ENDIF        
       nk_nivsed_out =  nk_nivsed_outlu+1   
       ALLOCATE(ep_nivsed_outp1(nk_nivsed_out))
       ep_nivsed_outp1(1:nk_nivsed_outlu)=ep_nivsed_out(1:nk_nivsed_outlu)
       ep_nivsed_outp1(nk_nivsed_out)=10.0_rsh  
       IF (MASTER) THEN
             WRITE(iscreenlog,*)'results in sediment will be save on ',nk_nivsed_out-1, &
                 'integrated layers whom thickness are constant and given by user - first is  sediment surface'
             WRITE(iscreenlog,*)'the ',nk_nivsed_out,'th layer will be an integrated layer till the bottom' 
       ENDIF  
       ALLOCATE(nivsed_out(0:nk_nivsed_out+1))
       nivsed_out(0)=0
       DO k=1,nk_nivsed_out
           nivsed_out(k)=nk_nivsed_out+1-k
       END DO
       nivsed_out(nk_nivsed_out+1)=0
      ENDIF

   ENDIF
 
   IF ( jextmin(nb) <= jextmax(nb) .AND. iextmin(nb) <= iextmax(nb) ) THEN

!ifois=ifois+1

!IF ( (ifois == 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) THEN
      IF ( (nbouc == 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) THEN
!
!       sig2(0)=-0.001
!       do k=ksdmin,ksdmax
!         sig2(k)=(k-1)/(ksdmax-ksdmin)
!       end do
!       sig2(ksdmax+1)=1.0
!
!
!
        IF (l_filebydate(nb)) THEN
          filoutsed=TRIM(out_dir)//TRIM(fileresu1)//'_'//TRIM(suffix)//'_sed'
          filoutsed=TRIM(filoutsed)//'_'//cdate(7:10)//cdate(4:5)//cdate(1:2)//cdate(12:13)//cdate(15:16)//cdate(18:19)
          filoutsed=TRIM(filoutsed)//'.nc'
        ELSE
          lnam=index(fileresu1,'.nc')-1
          filoutsed=TRIM(fileresu1(1:lnam))//'_sed.nc'
          write(*,*)filoutsed
        END IF
        IF( .NOT. l_out_nc4par) CALL ADD_TAG(filoutsed,LEN(filoutsed))
        IF (.NOT. l_out_nc4par) THEN
          CALL ionc4_createfile(filoutsed, &
                           lon2d(iextmin(nb):iextmax(nb):stepi(nb),jextmin(nb):jextmax(nb):stepj(nb)), &
                           lat2d(iextmin(nb):iextmax(nb):stepi(nb),jextmin(nb):jextmax(nb):stepj(nb)), &
                           iextmin(nb),iextmax(nb),stepi(nb),   &
                           jextmin(nb),jextmax(nb),stepj(nb),   &
                           1,nk_nivsed_out,1,nk_nivsed_out,0)
        ELSE
          CALL ionc4_createfile(filoutsed, &
                           lon2d_g(igextmin(nb):igextmax(nb):stepi(nb),jgextmin(nb):jgextmax(nb):stepj(nb)), &
                           lat2d_g(igextmin(nb):igextmax(nb):stepi(nb),jgextmin(nb):jgextmax(nb):stepj(nb)), &
                           igextmin(nb),igextmax(nb),stepi(nb), &
                           jgextmin(nb),jgextmax(nb),stepj(nb), &
                           1,nk_nivsed_out,1,nk_nivsed_out,0,l_out_nc4par = .TRUE.,comm_active=comm_zoom)
        ENDIF

! save attribute of the spatial extension of the domain
! -----------------------------------------------------
     CALL ionc4_gatt(filoutsed,'global_imin',igextmin(nb))
     CALL ionc4_gatt(filoutsed,'global_imax',igextmax(nb))
     CALL ionc4_gatt(filoutsed,'global_jmin',jgextmin(nb))
     CALL ionc4_gatt(filoutsed,'global_jmax',jgextmax(nb))

! Add attributes to follow OCO netCDF convention
     CALL out_netcdf_conv(filoutsed,nb,comm_zoom,igextmin(nb),igextmax(nb),jgextmin(nb),jgextmax(nb))

     IF (l_out_nc4par .AND. comm_zoom/=MPI_COMM_MARS) CALL MPI_COMM_FREE(comm_zoom,ierr_mpi)

! indispensable pour concatenation pas terrible a revoir
     CALL ionc4_createvar(filoutsed,'SIG',' ','sigma variable',&

                          standard_name= 'ocean_sigma_variable',&

                          l_auxcoordinate=.false.,&
                          valid_min= -1.0_riosh,valid_max= 0.0_riosh,fill_value= rg_valmanq_io,&
                          dims= 'z', l_out_nc4par = l_out_nc4par)
     CALL ionc4_write_sig(filoutsed,'SIG',nivsed_out,nk_nivsed_out)

     IF (l_morphocoupl) THEN
        CALL ionc4_createvar(filoutsed,TRIM(name_out_h0),'m',                   &
                             long_name = 'bathymetry relative to the mean level',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_h0,&
                             valid_max = riog_valid_max_h0,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
        CALL ionc4_createvar(filoutsed,TRIM(name_out_hx),'m',                   &
                             long_name = 'bathymetry x relative to the mean level',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_h0,&
                             valid_max = riog_valid_max_h0,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
        CALL ionc4_createvar(filoutsed,TRIM(name_out_hy),'m',                   &
                             long_name = 'bathymetry y relative to the mean level',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_h0,&
                             valid_max = riog_valid_max_h0,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
     ELSE
        CALL ionc4_createvar(filoutsed,TRIM(name_out_h0),'m',                   &
                             long_name = 'bathymetry relative to the mean level',&
                             dims = 'xy',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_h0,&
                             valid_max = riog_valid_max_h0,&
                             l_pack = .FALSE. , &
                             l_out_nc4par = l_out_nc4par)
        CALL ionc4_write_xy(filoutsed,name_out_h0,h0_local,limin,limax,ljmin,ljmax,&
                            fill_value = -rg_valmanq_io, l_out_nc4par = l_out_nc4par)
     END IF

        CALL ionc4_createvar(filoutsed,name_out_eptot,'m',                           &
                             long_name = 'total thickness of sediment',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_eptot,&
                             valid_max = riog_valid_max_eptot,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)

        CALL ionc4_createvar(filoutsed,name_out_nbniv,'quantite',                           &
                             long_name = 'number of sediment layers',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_nbniv,&
                             valid_max = riog_valid_max_nbniv,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)

        CALL ionc4_createvar(filoutsed,name_out_dzs,'m',                           &
                             long_name = 'thickness of sediment layer',&
                             dims = 'zxyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_dzs,&
                             valid_max = riog_valid_max_dzs,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
                    

        DO l = iv0,nv_use ! pour sauver T et S
          IF (l > 0) THEN
              name_out_cvsed=TRIM(name_var(irk_fil(l)))
              chaine=TRIM(ADJUSTL(ADJUSTR(long_name_var(irk_fil(l)))))
          ELSE IF (l==-1) THEN
              name_out_cvsed=TRIM(name_out_temp)
              chaine=name_out_cvsed
          ELSE IF (l==0) THEN
              name_out_cvsed=TRIM(name_out_sal)
              chaine=name_out_cvsed
          ENDIF

          CALL ionc4_createvar(filoutsed,name_out_cvsed,'u.s.i',                           &
                               long_name = chaine,&
                               dims = 'zxyt',&
                               fill_value = -rg_valmanq_io,&
!cval                               valid_min = riog_valid_min_cvsed,&
!cval                               valid_max = riog_valid_max_cvsed,&
                               l_pack = .FALSE. ,&
                               l_out_nc4par = l_out_nc4par)
        ENDDO
        
        CALL ionc4_createvar(filoutsed,name_out_tenfon,'N.m-2',                           &
                             long_name = 'bottom shear stress for erosion ',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_tenfon,&
                             valid_max = riog_valid_max_tenfon,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)


        CALL ionc4_createvar(filoutsed,name_out_tenfonc,'N.m-2',                           &
                             long_name = 'bottom shear stress due to the current',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_tenfon,&
                             valid_max = riog_valid_max_tenfon,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
        
        CALL ionc4_createvar(filoutsed,name_out_tenfonw,'N.m-2',                           &
                             long_name = 'bottom shear stress due to waves',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_tenfon,&
                             valid_max = riog_valid_max_tenfon,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
        
        CALL ionc4_createvar(filoutsed,name_out_wave_tp,'s',                           &
                             long_name = 'peak period',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_wave_tp,&
                             valid_max = riog_valid_max_wave_tp,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)
        
        CALL ionc4_createvar(filoutsed,name_out_wave_hs,'m',                           &
                             long_name = 'significant wave height',&
                             dims = 'xyt',&
                             fill_value = -rg_valmanq_io,&
                             valid_min = riog_valid_min_wave_hs,&
                             valid_max = riog_valid_max_wave_hs,&
                             l_pack = .FALSE. ,&
                             l_out_nc4par = l_out_nc4par)



     END IF

!IF ( (ifois > 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb) .AND. tect <= tfinsor(nb)) ) THEN
!IF ( (nbouc > 0 .AND. .NOT.l_filebydate(nb)) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb) .AND. tect <= tfinsor(nb)) ) THEN
     IF ( .NOT.l_filebydate(nb) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb) .AND. tect <= tfinsor(nb)) ) THEN
     
       CALL ionc4_write_time(filoutsed,0,tect)
       IF (l_morphocoupl) THEN
         DO j = ljmin,ljmax
           DO i = limin,limax
               h0_local(i,j)=h0(i,j)
               hx_local(i,j)=hx(i,j)
               hy_local(i,j)=hy(i,j)
           ENDDO
         ENDDO
         CALL ionc4_write_xyt(filoutsed,name_out_h0,h0_local,limin,limax,ljmin,ljmax,0,&
                              fill_value=-rg_valmanq_io, l_out_nc4par = l_out_nc4par)
         CALL ionc4_write_xyt(filoutsed,name_out_hx,hx_local,limin,limax,ljmin,ljmax,0,&
                              fill_value=-rg_valmanq_io, l_out_nc4par = l_out_nc4par)
         CALL ionc4_write_xyt(filoutsed,name_out_hy,hy_local,limin,limax,ljmin,ljmax,0,&
                              fill_value=-rg_valmanq_io, l_out_nc4par = l_out_nc4par)
       END IF
    
       DO j = ljmin,ljmax
         DO i = limin,limax
           IF(mask_h0(i,j)) THEN
             var2D(i,j)=-rg_valmanq_io
           ELSE 
             var2D(i,j)=ksma(i,j)
           ENDIF
         END DO
       END DO
       CALL out_write_xytvar(filoutsed,name_out_nbniv,-rg_valmanq_io,mask_h0, &
                             var2D(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

       DO j = ljmin,ljmax
         DO i = limin,limax
           DO k = ksma(i,j)+1,ksdmax 
             dzs(k,i,j)=-rg_valmanq_io
           END DO
           var2D(i,j)=0.0_rsh
           epsedtot(i,j)=0.0_rsh
           IF(mask_h0(i,j)) THEN
             var2D(i,j)=-rg_valmanq_io
             epsedtot(i,j)=-rg_valmanq_io
           ELSE
             DO k = ksdmin,ksma(i,j)
               var2D(i,j)=var2D(i,j)+dzs(k,i,j)
               epsedtot(i,j)=epsedtot(i,j)+dzs(k,i,j)
             END DO
           ENDIF
         END DO
       END DO
       CALL out_write_xytvar(filoutsed,name_out_eptot,-rg_valmanq_io,mask_h0, &
                             var2D(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

       ALLOCATE(tecrit3D(nk_nivsed_out,limin:limax,ljmin:ljmax))
       tecrit3D(:,:,:)=0.0_rsh
       DO j = ljmin,ljmax
         DO i = limin,limax
           IF(mask_h0(i,j)) THEN
             tecrit3D(:,i,j)=-rg_valmanq_io
           ELSE 
            IF(choice_nivsed_out == 1) THEN
! all layers are saved - k=1 bottom layer
              IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                  tecrit3D(ksdmin:ksma(i,j),i,j)=dzs(ksdmin:ksma(i,j),i,j)
              ELSE
                tecrit3D(:,i,j)=-rg_valmanq_io
              END IF
            ELSE IF ( choice_nivsed_out == 2) THEN
! nk_nivsed_out layers are saved (from surface) - k=1 most deep layer
                IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                  DO kniv=1,nk_nivsed_out
                    k=ksma(i,j)-nk_nivsed_out+kniv
                    IF(k >= 1) THEN
                       tecrit3D(kniv,i,j)=dzs(k,i,j)
                    ELSE
                       tecrit3D(kniv,i,j)=-rg_valmanq_io
                    ENDIF
                  ENDDO
                ELSE
                  tecrit3D(:,i,j)=-rg_valmanq_io
                END IF
            ELSE IF ( choice_nivsed_out == 3) THEN
! all layers till a max thickness given in parasedim.txt - k=1 : surface layer
! last layer is the layer whoch the bottom level exceeds the maximum thickness
              IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                eptmp=0.0_rsh
                finkniv=0
                tecrit3D(:,i,j)=0.
                kniv=1
                DO k = ksma(i,j),ksmi(i,j),-1
                  eptmp=eptmp+dzs(k,i,j)
                  IF(eptmp < epmax_nivsed_out) THEN
                    tecrit3D(kniv,i,j)=dzs(k,i,j)
                    kniv=kniv+1
                  ELSE
                    finkniv=1
                    IF(kniv <= nk_nivsed_out) THEN
                      tecrit3D(kniv,i,j)=tecrit3D(kniv,i,j)+dzs(k,i,j)
                    ENDIF
                  ENDIF
                END DO
              ELSE
                tecrit3D(:,i,j)=-rg_valmanq_io
              END IF
            ELSE IF ( choice_nivsed_out == 4) THEN
! a fixed number af layers is saved whith constant and fixed thickness given in parasedim.txt - k=1 : surface layer
! the last layer in addition integrates the results till the bottom sediment
              IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                eptmp=0.0_rsh
                kniv=1
                finkniv=0
                tecrit3D(:,i,j)=-rg_valmanq_io
                k=ksma(i,j)
                DO WHILE (k >= ksdmin .AND. finkniv==0)
                  eptmp=eptmp+dzs(k,i,j)
                  IF(eptmp < ep_nivsed_outp1(kniv) .AND. kniv==nk_nivsed_out ) THEN
                    tecrit3D(kniv,i,j)=epsedtot(i,j)-SUM(tecrit3D(1:kniv-1,i,j))
                    finkniv=1
                  ELSE IF (eptmp < ep_nivsed_outp1(kniv) )THEN
                    tecrit3D(kniv,i,j)=ep_nivsed_outp1(kniv)
                    eptmp=eptmp-ep_nivsed_outp1(kniv)
                    kniv=kniv+1  
                    IF(eptmp > 0 ) k=k+1
                  ELSE
                    DO WHILE (eptmp > ep_nivsed_outp1(kniv).AND. finkniv==0 ) 
                      eptmp=eptmp-ep_nivsed_outp1(kniv)
                      IF(kniv < nk_nivsed_out) THEN
                         tecrit3D(kniv,i,j)=ep_nivsed_outp1(kniv)
                         kniv=kniv+1
                      ELSE                                ! kniv=nk_nivsed_out
                         tecrit3D(kniv,i,j)=epsedtot(i,j)-SUM(tecrit3D(1:kniv-1,i,j))
                         finkniv=1
                      ENDIF
                    ENDDO
                  ENDIF
                  k=k-1
                END DO
                IF(k==0 .AND. finkniv==0 ) THEN
                  IF(eptmp > 0.0_rsh) THEN
                    tecrit3D(kniv,i,j)=eptmp
                    finkniv=1
                  ELSE
                    tecrit3D(kniv-1,i,j)=tecrit3D(kniv-1,i,j)+eptmp
                    finkniv=1
                  ENDIF
                ENDIF
              ELSE
                tecrit3D(:,i,j)=-rg_valmanq_io
              END IF
            ENDIF
           ENDIF
         END DO 
       END DO
       IF (l_out_nc4par) THEN
         CALL ionc4_write_zxyt(filoutsed,TRIM(ADJUSTL(ADJUSTR(name_out_dzs))),tecrit3D,    &
                limin,limax,ljmin,ljmax,1,nk_nivsed_out,0,fill_value=-rg_valmanq_io,l_out_nc4par = l_out_nc4par)
       ELSE
         CALL ionc4_write_zxyt(filoutsed,TRIM(ADJUSTL(ADJUSTR(name_out_dzs))),tecrit3D,    &
                limin,limax,ljmin,ljmax,1,nk_nivsed_out,0,fill_value=-rg_valmanq_io)
       ENDIF
 
       DO l = iv0,nv_use
         IF (l > 0) THEN
              name_out_cvsed=TRIM(name_var(irk_fil(l)))
              unitmudbinv=1._rsh/unit_modif_mudbio_N2dw(irk_fil(l))
         ELSE IF (l==-1) THEN
              name_out_cvsed=TRIM(name_out_temp)
              unitmudbinv=1._rsh
         ELSE IF (l==0) THEN
              name_out_cvsed=TRIM(name_out_sal)
              unitmudbinv=1._rsh
         ENDIF
         tecrit3D(:,:,:)=0.0_rsh

         DO j = ljmin,ljmax
           DO i = limin,limax
             IF(mask_h0(i,j)) THEN
               tecrit3D(:,i,j)=-rg_valmanq_io
             ELSE 
              IF(choice_nivsed_out == 1) THEN
! all layers are saved - k=1 bottom layer
                IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                  tecrit3D(ksdmin:ksma(i,j),i,j)=cv_sed(l,ksdmin:ksma(i,j),i,j)*unitmudbinv
                ELSE
                  tecrit3D(:,i,j)=-rg_valmanq_io
                END IF
              ELSE IF (choice_nivsed_out == 2) THEN
! nk_nivsed_out layers are saved (from surface) - k=1 most deep layer
                IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                   DO kniv=1,nk_nivsed_out
                    k=ksma(i,j)-nk_nivsed_out+kniv
                    IF (k >= 1) THEN
                      tecrit3D(kniv,i,j)=cv_sed(l,k,i,j)*unitmudbinv
                    ELSE
                      tecrit3D(kniv,i,j)=-rg_valmanq_io
                    ENDIF
                  ENDDO
                ELSE
                  tecrit3D(:,i,j)=-rg_valmanq_io
                END IF
              ELSE IF (choice_nivsed_out == 3) THEN
! all layers till a max thickness given in parasedim.txt
! last layer is the layer which the bottom level exceeds the maximum thickness
                IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                  eptmp=0.0_rsh
                  finkniv=0
                  tecrit3D(:,i,j)=0.0_rsh
                  kniv=1
                  epintbottom=0.0_rsh
                  DO k = ksma(i,j),ksmi(i,j),-1
                    eptmp=eptmp+dzs(k,i,j)
                    IF(eptmp < epmax_nivsed_out) THEN
                      tecrit3D(kniv,i,j)=cv_sed(l,k,i,j)*unitmudbinv
                      kniv=kniv+1
                    ELSE
                      finkniv=1
                      epintbottom=epintbottom+dzs(k,i,j)
                      IF(kniv <= nk_nivsed_out) THEN
                        tecrit3D(kniv,i,j)=tecrit3D(kniv,i,j)+cv_sed(l,k,i,j)*dzs(k,i,j)*unitmudbinv
                      ENDIF
                    ENDIF
                  END DO
                  IF(finkniv==1) tecrit3D(kniv,i,j)=tecrit3D(kniv,i,j)/epintbottom
                ELSE
                  tecrit3D(:,i,j)=-rg_valmanq_io
                END IF
              ELSE IF (choice_nivsed_out == 4) THEN
! a fixed number af layers is saved whith constant et fixed thickness given in parasedim.txt
! the last layer in addition integrates the results till the bottom sediment
                IF(ksma(i,j) >= 1 .AND. ksma(i,j) /= INT(rg_valmanq_io)) THEN
                  eptmp1=0.0_rsh
                  eptmp=0.0_rsh
                  cstmp=0.0_rsh
                  kniv=1
                  finkniv=0
                  tecrit3D(:,i,j)=-rg_valmanq_io
                  k=ksma(i,j)
                  DO WHILE (k >= ksdmin .and. finkniv ==0)
                    eptmp=eptmp+dzs(k,i,j)
                    IF(eptmp < ep_nivsed_outp1(kniv) .AND. kniv==nk_nivsed_out) THEN
                      tecrit3D(kniv,i,j)=(cstmp+SUM(cv_sed(l,ksdmin:k,i,j)*dzs(ksdmin:k,i,j)*unitmudbinv))  &
                                            /(SUM(dzs(ksdmin:k,i,j))+eptmp-dzs(k,i,j))
                      finkniv=1
                    ELSE
                     IF(eptmp <= ep_nivsed_outp1(kniv)) THEN
                         cstmp=cstmp+dzs(k,i,j)*cv_sed(l,k,i,j)*unitmudbinv 
                         eptmp1=eptmp
                     ENDIF
                     DO WHILE (eptmp > ep_nivsed_outp1(kniv) .and. finkniv ==0 ) 
                        cstmp=cstmp+(ep_nivsed_outp1(kniv)-eptmp1)*cv_sed(l,k,i,j)*unitmudbinv
                        eptmp=eptmp-ep_nivsed_outp1(kniv)
                        IF(kniv < nk_nivsed_out) THEN
                          tecrit3D(kniv,i,j)=cstmp/(ep_nivsed_outp1(kniv))
                          kniv=kniv+1
                          IF(eptmp > ep_nivsed_outp1(kniv)) THEN
                             cstmp=0.0_rsh
                             eptmp1=0.0_rsh
                          ELSE
                             eptmp1=eptmp
                             cstmp=eptmp1*cv_sed(l,k,i,j)*unitmudbinv
                          ENDIF
                        ELSE                                ! kniv=nk_nivsed_out
                          IF(k==1) THEN
                            tecrit3D(kniv,i,j)=cv_sed(l,k,i,j)*unitmudbinv
                          ELSE
                            tecrit3D(kniv,i,j)=(cstmp+SUM(cv_sed(l,ksdmin:k-1,i,j)*dzs(ksdmin:k-1,i,j)*unitmudbinv))  &
                                                 /(eptmp+SUM(dzs(ksdmin:k-1,i,j)))
                          ENDIF
                          finkniv=1
                        ENDIF
                      ENDDO 
                    ENDIF
                    k=k-1
                  END DO
                  IF(k==0 .AND. finkniv==0 .AND. eptmp > 0.0_rsh) THEN
                    tecrit3D(kniv,i,j)=cv_sed(l,1,i,j)*unitmudbinv
                    finkniv=1
                  ENDIF
                ELSE
                  tecrit3D(:,i,j)=-rg_valmanq_io
                END IF
              ENDIF
             ENDIF
           END DO 
         END DO
 
         IF (l_out_nc4par) THEN
           CALL ionc4_write_zxyt(filoutsed,TRIM(ADJUSTL(ADJUSTR(name_out_cvsed))),tecrit3D,    &
                limin,limax,ljmin,ljmax,1,nk_nivsed_out,0,fill_value=-rg_valmanq_io,l_out_nc4par = l_out_nc4par)
         ELSE
           CALL ionc4_write_zxyt(filoutsed,TRIM(ADJUSTL(ADJUSTR(name_out_cvsed))),tecrit3D,    &
                limin,limax,ljmin,ljmax,1,nk_nivsed_out,0,fill_value=-rg_valmanq_io)
         ENDIF
       ENDDO
       DEALLOCATE(tecrit3D)

       CALL out_write_xytvar(filoutsed,name_out_tenfon,-rg_valmanq_io,mask_h0, &
                             tenfon(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)


       CALL out_write_xytvar(filoutsed,name_out_tenfonc,-rg_valmanq_io,mask_h0,&
                          tenfonc(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(filoutsed,name_out_tenfonw,-rg_valmanq_io,mask_h0,&
                          tenfonw(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(filoutsed,name_out_wave_tp,-rg_valmanq_io,mask_h0,&
                          wave_tp(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
       CALL out_write_xytvar(filoutsed,name_out_wave_hs,-rg_valmanq_io,mask_h0,&
                          wave_hs(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)

     END IF

     IF (l_filebydate(nb)) THEN
       CALL ionc4_close(filoutsed)
     ELSE
! To write the data on the disk and not loose data in case of run crash
       CALL ionc4_sync(filoutsed)
     ENDIF

   END IF

   10 FORMAT('C',i2.2)

   !PRINT*, 'exit sed_output'



  END SUBROUTINE sed_output

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