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

  SUBROUTINE output_diag(tect,x,nb,fileresu1)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE output_diag  ***
   !&E
   !&E ** Purpose :
   !&E
   !&E ** Description : Sortie des variables diagnostique biologique
   !&E
   !&E ** Called by : output_mng.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  09-2007 C. Penard (creation de la routine)
   !&E       !  09-2008 V. Garnier (generalized sigma)
   !&E       !  2008-09 (G. Galibert / ALTRAN) attributes of output files according to NetCDF OCO v1.1 file format
   !&E       !  2009-03 (V. Garnier) attributes according to NetCDF OCO v1.2 file format
   !&E                               saving of generalized sigma parameters
   !&E       !  2011-01 (R. Ramel Alyotech) Ionc library in fortran 90/95, netcdf4 and MPI
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

#ifdef key_biolo

   !! Modules used
   USE parameters,   ONLY : rsh,rlg,riosh,kmax,imin,imax,limin,limax,ljmin,ljmax,  &
                            rg_valmanq_io,ljminp1,riolg,lchain,valmanq
   USE ionc4,        ONLY : ionc4_createfile,ionc4_write_lon,ionc4_gatt,ionc4_createvar,&
                            ionc4_write_sig,ionc4_vatt_standard_name,ionc4_vatt_fill,   &
                            ionc4_write_xy,ionc4_vatt_axis,ionc4_vatt_positive,ionc4_write_time,&
                            ionc4_write_xyt,ionc4_write_t,ionc4_write_zxyt,ionc4_sync,ionc4_close,&
                            ionc4_write_sigw,ionc4_vatt_chain,ionc4_write
   USE outnetcdf,    ONLY : out_netcdf_conv,out_zaxis_w
   USE output,       ONLY : file_output,l_filebydate,stepi,stepj,stepk,   &
                            iextmin,iextmax,jextmin,jextmax,kextmin,kextmax,l_out_zaxis_w,       &
                            igextmin,igextmax,jgextmin,jgextmax,l_out_nc4par,tsort,tfinsor
   USE comvars2d,    ONLY : nbouc,lon2d_g,lat2d_g,lon2d,lat2d,date_ref,h0fond,date_startshift,  &
                            h0,hx,hy,name_out_hx,name_out_hy,name_out_h0,name_out_xe,           &
                            l_modele2d,l_initfromfile,cdate,ig,id,jb,jh,out_dir,suffix,iscreen, &
                            riog_valid_min_h0,riog_valid_max_h0,riog_valid_min_xe,riog_valid_max_xe
   USE comvars3d,    ONLY : sig
   USE cometdyn,     ONLY : tosx,tosy
#if defined key_siggen || defined key_gencoord
   USE comsiggen,    ONLY : hc_sig,Csu_sig,Csw_sig,b_sig,theta_sig
#endif
#ifdef key_substance
   USE parameters,   ONLY : nb_var
   USE comsubstance, ONLY : irk_fil,name_var,long_name_var,standard_name_var, &
                            valid_min_var,valid_max_var,unit_var,nv_state,l_out_log
#endif

   USE bioloinit,    ONLY : irk_diag,ndiag_1d,unit_vardiag,idimv_r,   &
                     ndiag_2d,ndiag_3d_wat,ndiag_tot,name_vardiag,    &
                     diag_1d,diag_2d,diag_3d_wat,l_diag_sortie,       &
                     long_name_vardiag, standard_name_vardiag,        &
                     valid_min_vardiag, valid_max_vardiag
#ifdef key_MPI_2D
   USE toolmpi, ONLY: ADD_TAG,IPROC,NPROCS,ADD_ALL_MPI_INT,MPI_COMM_MARS
#endif
   USE obccombine,ONLY : l_obc_tide

   IMPLICIT NONE

#ifdef key_MPI_2D
   INCLUDE 'mpif.h'
#endif

   !! * 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,kk,ifois
   INTEGER               :: grp_world,grp_zoom,comm_zoom,ncpu,ip,ierr_mpi
   INTEGER, DIMENSION(:), ALLOCATABLE  :: numrank,listrank
   INTEGER, DIMENSION(1) :: listrank_nb
   CHARACTER(LEN=lchain) :: file_diag_ecosys_wat
   CHARACTER(LEN=lchain) :: chain
   CHARACTER(LEN=5) :: vardims
   CHARACTER (len=19)    :: tool_sectodat,cdateout
   LOGICAL,         DIMENSION(limin:limax,ljmin:ljmax)      :: mask_h0,mask_hx,mask_hy
   REAL(KIND=rsh)        :: xx,yy
   REAL(KIND=riosh)      :: vecrit
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)      :: tecrit
   REAL(KIND=rlg)        :: decrit
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax)      :: h0_local,hx_local,hy_local
   REAL(KIND=riosh),DIMENSION(kmax,limin:limax,ljmin:ljmax) :: tzecrit
   REAL(KIND=riosh),DIMENSION(0:kmax+1)                     :: sigecrit
   REAL(KIND=riolg),DIMENSION(imin:imax)                    :: lon

   !DATA ifois /-1/
   !!----------------------------------------------------------------------
   !! * Executable part

   PRINT_DBG*, 'enter OUTPUT_DIAG routine'
  ! Masks
   ! -----

   mask_h0(:,:)=.TRUE.
   DO j=ljminp1,ljmax
     DO i=MAX(limin,ig(j)),MIN(limax,id(j))
       ! For heights smaller than h0fond, the cell is considered as wet
       IF ( h0(i,j)+x(i,j) > h0fond ) mask_h0(i,j)=.FALSE.
     END DO
   END DO

   mask_hx(:,:)=.TRUE.
   DO j=ljminp1,ljmax
     DO i=MAX(limin,ig(j)-1),MIN(limax,id(j))
       IF ( hx(i,j) /= -valmanq -h0fond ) mask_hx(i,j)=.FALSE.
     END DO
   END DO

   mask_hy(:,:)=.TRUE.
   DO j=ljmin,ljmax
     DO i=MAX(limin,ig(j)),MIN(limax,id(j))
       IF ( (j >= jb(i)-1 .AND. j <= jh(i)) .AND. (hy(i,j) /= -valmanq -h0fond) ) mask_hy(i,j)=.FALSE.
     END DO
   END DO

   comm_zoom=0
#ifdef key_MPI_2D
   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_MPI ( 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
#endif

   h0_local=h0(limin:limax,ljmin:ljmax)
   hx_local=hx(limin:limax,ljmin:ljmax)
   hy_local=hy(limin:limax,ljmin:ljmax)

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

     !ifois=ifois+1

! modif juin2014 : definition of file_diag_ecosys_wat at each tsort because of using OMP
       IF (l_filebydate(nb) ) THEN
   !     file_diag_ecosys_wat=TRIM(fileresu1)//'_'//TRIM(suffix)//'_diag'
   !     file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//'_'//cdate(7:10)//cdate(4:5)//cdate(1:2)//'T'//cdate(12:13)
   !     file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//cdate(15:16)//cdate(18:19)//'Z'
   !     file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//'.nc'
         cdateout=tool_sectodat(tect)
         file_diag_ecosys_wat=TRIM(out_dir)//TRIM(fileresu1)//'_diag'//TRIM(suffix)
         file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//'_'//cdateout(7:10)//cdateout(4:5)//cdateout(1:2)//'T'
         file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//cdateout(12:13)//cdateout(15:16)//cdateout(18:19)//'Z'
         file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//'.nc'
         IF_MPI (.NOT. l_out_nc4par) CALL ADD_TAG(file_diag_ecosys_wat,LEN(file_diag_ecosys_wat))

       ELSE
        file_diag_ecosys_wat=TRIM(fileresu1)
        kk=INDEX(file_diag_ecosys_wat,'.nc')
        file_diag_ecosys_wat=file_diag_ecosys_wat(1:kk-1)
        file_diag_ecosys_wat=TRIM(file_diag_ecosys_wat)//'_diag.nc'
       END IF


     !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

!---------------------------------
! creation du fichier de resultats
!---------------------------------

       WRITE(*,*) 'Creation du fichier result - Diag-bio'


!      DO i=imin,imax
!        lon(i)=gwac+(i-1)*dgac  !+360.0
!      END DO
       sigecrit(:)=sig(:)
       IF (l_out_nc4par) THEN
       ! To allow the different cpus to write at the same time into a global
       ! file, set l_out_nc4par= .TRUE.
       ! The dimension of the file must be global (igextmin)
         CALL ionc4_createfile(file_diag_ecosys_wat,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), &
                                     kextmin(nb),kextmax(nb),stepk(nb),   &
                                     kmax,0,l_out_nc4par= .TRUE.,comm_active=comm_zoom, &
                                     zaxis_var=sigecrit(1:kmax))
       ELSE
         ! If none l_out_nc4par argument, l_out_nc4par is considered as false
         ! Each cpu create its own output file inside the local domain (iextmin)
         CALL ionc4_createfile(file_diag_ecosys_wat,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), &
                       kextmin(nb),kextmax(nb),stepk(nb),kmax,0,zaxis_var=sigecrit(1:kmax))
       ENDIF

!     CALL ionc4_write_lon(file_diag_ecosys_wat,lon(iextmin(nb):iextmax(nb)),iextmin(nb),iextmax(nb))

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

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

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

!     CALL ionc4_torigine(file_diag_ecosys_wat,date_ref)
   ! Add the vertical axis at w location
       IF (l_out_zaxis_w(nb)) CALL out_zaxis_w(file_diag_ecosys_wat,kextmin(nb),kextmax(nb),stepk(nb),comm_zoom=comm_zoom)


       CALL ionc4_createvar(file_diag_ecosys_wat,'SIG',' ','sigma variable',&
#if defined key_siggen || defined key_gencoord
                          standard_name= 'ocean_s_variable',&
#else
                          standard_name= 'ocean_sigma_variable',&
#endif                          
                          l_auxcoordinate=.false.,&
                          dims='z',valid_min= -1.0_riosh,valid_max= 0.0_riosh,fill_value= rg_valmanq_io)
       CALL ionc4_write_sig(file_diag_ecosys_wat,'SIG',sigecrit,kmax)

#if defined key_siggen || defined key_gencoord
     
       sigecrit(:)=Csu_sig(:)
       CALL ionc4_createvar(file_diag_ecosys_wat,'Csu_sig',' ','C(s) law at mid layer',dims='z', &
                           standard_name='ocean_s_coordinate_function_at_midlayer', &
                           valid_min=-2.0_riosh,valid_max=2.0_riosh,fill_value= rg_valmanq_io, &
                           l_auxcoordinate=.false.)
       CALL ionc4_write_sig(file_diag_ecosys_wat,'Csu_sig',sigecrit,kmax)
       CALL ionc4_vatt_chain(file_diag_ecosys_wat,'Csu_sig','formula_definition', &
          'C(k) = (1-b)*sinh(a*s(k))/sinh(a) + b*[tanh(a*(s(k)+0.5))-tanh(0.5*a)]/(2*tanh(0.5*a))')

       IF (l_out_zaxis_w(nb)) THEN
         sigecrit(:)=Csw_sig(:)
         CALL ionc4_createvar(file_diag_ecosys_wat,'Csw_sig',' ','C(s) law at interface',dims='z',location_in_cell='w', &
                           standard_name='ocean_s_coordinate_function_at_interface', &
                           valid_min= -2.0_riosh,valid_max= 2.0_riosh,fill_value= rg_valmanq_io, &
                           l_auxcoordinate=.false.)
         CALL ionc4_write_sigw(file_diag_ecosys_wat,'Csw_sig',sigecrit,kmax)
         CALL ionc4_vatt_chain(file_diag_ecosys_wat,'Csw_sig','formula_definition', &
          'C_w(k) = (1-b)*sinh(a*s_w(k))/sinh(a) + b*[tanh(a*(s_w(k)+0.5))-tanh(0.5*a)]/(2*tanh(0.5*a))')
       END IF

       !sigecrit(:)=Csw_sig(:)
       !CALL ionc4_createvar(file_diag_ecosys_wat,'Csw_sig',' ','C(s) law at interface',dims='z', &
       !                      standard_name='ocean_s_coordinate_function_at_interface', &
       !                      valid_min= -2.0_riosh,valid_max= 2.0_riosh,fill_value= rg_valmanq_io, &
       !                      l_auxcoordinate=.false.,location_in_cell='w')
       !CALL ionc4_write_sigw(file_diag_ecosys_wat,'Csw_sig',sigecrit,kmax)
       !CALL ionc4_vatt_chain(file_diag_ecosys_wat,'Csw_sig','formula_definition', &
       !     'C_w(k) = (1-b)*sinh(a*s_w(k))/sinh(a) + b*[tanh(a*(s_w(k)+0.5))-tanh(0.5*a)]/(2*tanh(0.5*a))')
     
       CALL ionc4_createvar(file_diag_ecosys_wat,'hc','m', 'depth over which the resolution is increased', &
                           standard_name='ocean_s_coordinate_thick_of_surface_resolution',dims='xy', &
                           valid_min= -20.0_riosh,valid_max= 500.0_riosh,fill_value= rg_valmanq_io,&
                           l_out_nc4par=l_out_nc4par,l_pack=.false. )
       tecrit(limin:limax,ljmin:ljmax)=hc_sig(limin:limax,ljmin:ljmax)
       CALL ionc4_write_xy(file_diag_ecosys_wat,'hc',tecrit,limin,limax,ljmin,ljmax,fill_value=rg_valmanq_io, &
                                              l_out_nc4par=l_out_nc4par)
     
       CALL ionc4_createvar(file_diag_ecosys_wat,'b',' ','bottom control parameter',  &
                           standard_name='ocean_s_coordinate_bottom_control', &
                           valid_min=0.0_riosh,valid_max=1.0_riosh,fill_value=rg_valmanq_io)
       vecrit=b_sig
       CALL ionc4_write(file_diag_ecosys_wat,'b',vecrit)
     
       CALL ionc4_createvar(file_diag_ecosys_wat,'theta',' ','surface control parameter', &
                           standard_name ='ocean_s_coordinate_surface_control', &
                           valid_min=0.0_riosh,valid_max=20.0_riosh,fill_value=rg_valmanq_io)
       vecrit=theta_sig
       CALL ionc4_write(file_diag_ecosys_wat,'theta',vecrit)
#endif



!ecriture de la bathymetrie
!--------------------------
       chain='model_sea_floor_depth_below_geoid'
       CALL ionc4_createvar(file_diag_ecosys_wat,name_out_h0,units='m',&
                           long_name='model bathymetry relative to the mean level',&
                           standard_name = TRIM(chain),&
                           fill_value = -rg_valmanq_io,&
                           valid_min = riog_valid_min_h0,&
                           valid_max = riog_valid_max_h0,&
                           dims = 'xy',&
                           l_pack = .false.,&
                           l_out_nc4par = l_out_nc4par) 
       tecrit(:,:)=h0_local(:,:)
       CALL ionc4_write_xy(file_diag_ecosys_wat,TRIM(name_out_h0),tecrit,limin,limax,ljmin,ljmax,&
                          fill_value=-rg_valmanq_io,l_out_nc4par = l_out_nc4par )

! hx
!---
       chain='model_sea_floor_depth_below_geoid_at_u_location'
       CALL ionc4_createvar(file_diag_ecosys_wat,name_out_hx,units='m',&
                           long_name = 'model bathymetry at u-location relative to the mean level',&
                           standard_name = TRIM(chain),&
                           fill_value = -REAL(rg_valmanq_io+h0fond,riosh),&
                           valid_min = riog_valid_min_h0,&
                           valid_max = riog_valid_max_h0,&
                           dims = 'xy',&
                           l_pack = .false.,&
                           l_out_nc4par = l_out_nc4par)
       tecrit(:,:)=hx_local(:,:)
       CALL ionc4_write_xy(file_diag_ecosys_wat,TRIM(name_out_hx),tecrit,limin,limax,ljmin,ljmax ,&
                         fill_value = -REAL(rg_valmanq_io+h0fond,riosh),l_out_nc4par = l_out_nc4par)

! hy
!---
       chain='model_sea_floor_depth_below_geoid_at_v_location'
       CALL ionc4_createvar(file_diag_ecosys_wat,name_out_hy,units='m',&
                           long_name = 'model bathymetry at v-location relative to the mean level',&
                           standard_name = TRIM(chain),&
                           fill_value = -REAL(rg_valmanq_io+h0fond,riosh),&
                           valid_min = riog_valid_min_h0,&
                           valid_max = riog_valid_max_h0,&
                           dims = 'xy',&
                           l_pack = .false.,&
                           l_out_nc4par = l_out_nc4par)
       tecrit(:,:)=hy_local(:,:)
       CALL ionc4_write_xy(file_diag_ecosys_wat,TRIM(name_out_hy),tecrit,limin,limax,ljmin,ljmax,&
                         fill_value = -REAL(rg_valmanq_io+h0fond,riosh),l_out_nc4par = l_out_nc4par)

! creation des variables qui varient au cours du temps
!--------------------------------------------------------

! save sea surface height
!--------------------------------------------
      CALL ionc4_createvar(file_diag_ecosys_wat,name_out_xe,units= 'm',&
                         long_name = 'sea surface height',&
                         standard_name = 'sea_surface_height_above_geoid',&
                         fill_value = rg_valmanq_io,&
                         valid_min = riog_valid_min_xe,&
                         valid_max = riog_valid_max_xe,&
                         dims= 'xyt',&
                         l_pack = .false.,&
                         l_out_nc4par = l_out_nc4par)

!     CALL ionc4_vatt_axis(file_diag_ecosys_wat,TRIM(name_out_xe),'Z')
!     CALL ionc4_vatt_positive(file_diag_ecosys_wat,TRIM(name_out_xe),'up')
!--------------------------------------------


       DO l = 1,ndiag_tot
! test si on ecrit cette variable diagnostique ou pas
         IF (l_diag_sortie(l)) THEN
! Variable diagnostique 1d
           IF (idimv_r(l)==1) THEN
             vardims = 't'
! Variable diagnostique 2d	   
           ELSEIF  (idimv_r(l)==2) THEN
             vardims = 'xyt'
! Variable diagnostique 3d		   
           ELSEIF (idimv_r(l)==3) THEN
             IF (l_modele2d) THEN
               vardims = 'xyt'
             ELSE
               vardims = 'zxyt'
             ENDIF
           ENDIF

           unit_vardiag(l)=TRIM(unit_vardiag(l))
           CALL ionc4_createvar(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
                              units = TRIM(unit_vardiag(l)),&
                              long_name = TRIM(long_name_vardiag(l)),&
                              standard_name = TRIM(standard_name_vardiag(l)),&
                              fill_value = rg_valmanq_io,&
                              valid_min = REAL(valid_min_vardiag(l),riosh),&
                              valid_max = REAL(valid_max_vardiag(l),riosh),&
                              dims = TRIM(vardims),&
                              l_pack = .FALSE., &
                              l_out_nc4par = l_out_nc4par)
         END IF
       END DO


!------------------------------------
! sauvegarde de la condition initiale
!------------------------------------

!     IF ( (.NOT.l_initfromfile .AND. date_startshift == 0.0_rlg) .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) .OR. &
!         (l_initfromfile .AND. nbouc== 0) ) THEN

!       CALL ionc4_write_time(file_diag_ecosys_wat,0,tect)

! sauvegarde des xe
!       tecrit(:,:)=x(limin:limax,ljmin:ljmax)
!       WHERE(h0_local(:,:)+x(limin:limax,ljmin:ljmax) <= h0fond) tecrit(:,:)=rg_valmanq_io
!       CALL ionc4_write_xyt(file_diag_ecosys_wat,TRIM(name_out_xe),tecrit,limin,limax,ljmin,ljmax,0,&
!                            fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)


!       DO l = 1,ndiag_tot
! test si on ecrit cette variable diagnostique ou pas
!         IF (l_diag_sortie(l)) THEN
! Variable diagnostique 1d
!           IF (idimv_r(l)==1) THEN
!            vecrit = diag_1d(irk_diag(l))
!            CALL ionc4_write_t(file_diag_ecosys_wat,TRIM(name_vardiag(l)),vecrit,0, fill_value = rg_valmanq_io)

! Variable diagnostique 2d
!           ELSEIF (idimv_r(l)==2) THEN
!             tecrit(:,:) = diag_2d(irk_diag(l),limin:limax,ljmin:ljmax)  
!             IF (l_modele2d) THEN
!               WHERE(h0_local(:,:)+x(limin:limax,ljmin:ljmax) <= h0fond) tecrit(:,:) = rg_valmanq_io
!             ELSE
!               WHERE(h0_local(:,:)==-rg_valmanq_io) tecrit(:,:) = rg_valmanq_io
!             ENDIF
!             CALL ionc4_write_xyt(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
!                                  tecrit,limin,limax,ljmin,ljmax,0,&
!                                  fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)

! Variable diagnostique 3d
!           ELSEIF (idimv_r(l)==3) THEN
!             IF (l_modele2d) THEN
!               tecrit(:,:) = diag_3d_wat(irk_diag(l),1,limin:limax,ljmin:ljmax)
!               WHERE(h0_local(:,:)+x(limin:limax,ljmin:ljmax) <= h0fond) tecrit(:,:)=rg_valmanq_io
!               CALL ionc4_write_xyt(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
!                                    tecrit,limin,limax,ljmin,ljmax,0,&
!                                    fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)
!             ELSE
!               tzecrit(:,:,:) = diag_3d_wat(irk_diag(l),:,limin:limax,ljmin:ljmax)
!               DO j=ljmin,ljmax
!                 DO i=limin,limax
!                   IF (h0_local(i,j)==-rg_valmanq_io) tzecrit(:,i,j)=rg_valmanq_io
!                 END DO
!               END DO
!               CALL ionc4_write_zxyt(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
!                                     tzecrit,limin,limax,ljmin,ljmax,1,kmax,0,&
!                                     fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)
!             END IF
!           END IF
!         END IF
!       END DO

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

!       END IF

     ENDIF

     !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

! On enregistre le temps
     CALL ionc4_write_time(file_diag_ecosys_wat,0,tect)

! sauvegarde des xe
       tecrit(:,:)=x(limin:limax,ljmin:ljmax)
       WHERE(h0_local(:,:)+x(limin:limax,ljmin:ljmax) <= h0fond) tecrit(:,:)=rg_valmanq_io
       CALL ionc4_write_xyt(file_diag_ecosys_wat,TRIM(name_out_xe),tecrit,limin,limax,ljmin,ljmax,0,&
                          fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)

! On enregistre les variables diagnostiques

       DO l = 1,ndiag_tot
! test si on ecrit cette variable diagnostique ou pas
         IF (l_diag_sortie(l)) THEN
! Variable diagnostique 1d
           IF (idimv_r(l)==1) THEN
            vecrit = diag_1d(irk_diag(l))
            CALL ionc4_write_t(file_diag_ecosys_wat,TRIM(name_vardiag(l)),vecrit,0, fill_value = rg_valmanq_io)

! Variable diagnostique 2d
           ELSEIF (idimv_r(l)==2) THEN
             tecrit(:,:) = diag_2d(irk_diag(l),limin:limax,ljmin:ljmax)  
             IF (l_modele2d) THEN
               WHERE(h0_local(:,:)+x(limin:limax,ljmin:ljmax) <= h0fond) tecrit(:,:) = rg_valmanq_io
             ELSE
               WHERE(h0_local(:,:)==-rg_valmanq_io) tecrit(:,:) = rg_valmanq_io
             ENDIF
             CALL ionc4_write_xyt(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
                                  tecrit,limin,limax,ljmin,ljmax,0,&
                                  fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)

! Variable diagnostique 3d
           ELSEIF (idimv_r(l)==3) THEN
             IF (l_modele2d) THEN
               tecrit(:,:) = diag_3d_wat(irk_diag(l),1,limin:limax,ljmin:ljmax)
               WHERE(h0_local(:,:)+x(limin:limax,ljmin:ljmax) <= h0fond) tecrit(:,:)=rg_valmanq_io
               CALL ionc4_write_xyt(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
                                    tecrit,limin,limax,ljmin,ljmax,0,&
                                    fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)
             ELSE
               tzecrit(:,:,:) = diag_3d_wat(irk_diag(l),:,limin:limax,ljmin:ljmax)
               DO j=ljmin,ljmax
                 DO i=limin,limax
                   IF (h0_local(i,j)==-rg_valmanq_io) tzecrit(:,i,j)=rg_valmanq_io
                 END DO
               END DO
               CALL ionc4_write_zxyt(file_diag_ecosys_wat,TRIM(name_vardiag(l)), &
                                     tzecrit,limin,limax,ljmin,ljmax,1,kmax,0,&
                                     fill_value = rg_valmanq_io, l_out_nc4par = l_out_nc4par)
             END IF
           END IF
         END IF
       END DO

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

     END IF

   ELSE IF ( nbouc == 0 .OR. (l_filebydate(nb) .AND. tect >= tsort(nb)) ) THEN
     WRITE(iscreen,*) 'Output file ',TRIM(file_diag_ecosys_wat),' is not created'
     WRITE(iscreen,*) 'Have a look at line iextmin,iextmax,jextmin,jextmax'
     WRITE(iscreen,*) 'in file :',TRIM(file_output)
   END IF

   PRINT_DBG*, 'exit OUTPUT_DIAG routine'

#else
   !! * Modules
   USE parameters,   ONLY : rsh,rlg,limin,limax,ljmin,ljmax
   USE comvars2d,    ONLY : iscreen
   !! * Arguments
   INTEGER       ,INTENT(in) :: nb
   REAL(KIND=rlg),INTENT(in) :: tect
   REAL(KIND=rsh),DIMENSION(limin:limax,ljmin:ljmax),INTENT(in) :: x
   WRITE(iscreen,*) 'OUTPUT_DIAG subroutine should not be called '
   WRITE(iscreen,*) 'if the CPP key "key_biolo" is not activated.'
   WRITE(iscreen,*) 'Simulation stopped.'
   STOP

#endif

  END SUBROUTINE output_diag

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