  MODULE wavediag

   !!======================================================================
   !!                   ***  MODULE wavediag  ***
   !!
   !! ** Description :
   !!
   !! ** Reference :
   !!
   !! ** History :  !  2011-04  (P. Lazure)  Original code
   !!
   !!======================================================================

   !! * Modules used
   USE parameters, ONLY : lchain,rsh,rlg,riolg,riosh,kmax,imin,imax,jmin,jmax, &
                          limin,limax,ljmin,ljmax

   IMPLICIT NONE

#ifdef key_MPI_2D
   INCLUDE 'mpif.h'
#endif

   PRIVATE

   !! * Accessibility

   PUBLIC wave_diag_qeuler_alloc,   &       ! called by alloc.F90
          wave_diag_qeuler_dealloc, &       ! called by dealloc.F90
          wave_diag_qeuler_init,   &       ! called by main.F90
          wave_diag_qeuler_output          ! called by dyn3dxztlssolve.F90, dyn3dyztlssolve.F90

   !! * Shared module variables

   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_advx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_advy
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_corx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_cory
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_vortexx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_vortexy
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_smb1x
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_smb1y
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_grdpx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_grdpy
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_difhx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_difhy
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_difvx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_difvy
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_streamx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_streamy
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_nudtkedz
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_nudpsidz
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_tke
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_psi
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_pk
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_peps
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_prodtke
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_prodpsi
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_dissiptke
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_dissippsi
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_vortx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_vorty
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_vortz
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_dudt
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:,:),PUBLIC :: diag_wf_dvdt
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:),PUBLIC :: diag_wf_xe
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:),PUBLIC :: diag_wf_botfricx
   REAL(KIND=rsh),ALLOCATABLE,DIMENSION(:,:),PUBLIC :: diag_wf_botfricy

   CHARACTER(LEN=lchain) :: name_outdiag_wf_advx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_advy
   CHARACTER(LEN=lchain) :: name_outdiag_wf_corx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_cory
   CHARACTER(LEN=lchain) :: name_outdiag_wf_vortexx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_vortexy
   CHARACTER(LEN=lchain) :: name_outdiag_wf_smb1x
   CHARACTER(LEN=lchain) :: name_outdiag_wf_smb1y
   CHARACTER(LEN=lchain) :: name_outdiag_wf_grdpx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_grdpy
   CHARACTER(LEN=lchain) :: name_outdiag_wf_difhx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_difhy
   CHARACTER(LEN=lchain) :: name_outdiag_wf_difvx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_difvy
   CHARACTER(LEN=lchain) :: name_outdiag_wf_streamx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_streamy
   CHARACTER(LEN=lchain) :: name_outdiag_wf_nudtkedz
   CHARACTER(LEN=lchain) :: name_outdiag_wf_nudpsidz
   CHARACTER(LEN=lchain) :: name_outdiag_wf_tke
   CHARACTER(LEN=lchain) :: name_outdiag_wf_psi
   CHARACTER(LEN=lchain) :: name_outdiag_wf_peps
   CHARACTER(LEN=lchain) :: name_outdiag_wf_pk
   CHARACTER(LEN=lchain) :: name_outdiag_wf_prodtke
   CHARACTER(LEN=lchain) :: name_outdiag_wf_prodpsi
   CHARACTER(LEN=lchain) :: name_outdiag_wf_dissiptke
   CHARACTER(LEN=lchain) :: name_outdiag_wf_dissippsi
   CHARACTER(LEN=lchain) :: name_outdiag_wf_vortx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_vorty
   CHARACTER(LEN=lchain) :: name_outdiag_wf_vortz
   CHARACTER(LEN=lchain) :: name_outdiag_wf_dudt
   CHARACTER(LEN=lchain) :: name_outdiag_wf_dvdt
   CHARACTER(LEN=lchain) :: name_outdiag_wf_xe
   CHARACTER(LEN=lchain) :: name_outdiag_wf_botfricx
   CHARACTER(LEN=lchain) :: name_outdiag_wf_botfricy

   
   LOGICAL,ALLOCATABLE,DIMENSION(:,:),PUBLIC          :: mask_diagx,mask_diagy

   !! * Private variables

   REAL(KIND=rlg)  :: tdebdiagdyn,tfindiagdyn,dtsordiagdyn,tnextout,tcumul,tprecmoy,dtmoy

   LOGICAL,ALLOCATABLE,DIMENSION(:,:)  :: mask_h0_diag,mask_hx_diag,mask_hy_diag

   LOGICAL     :: l_moy_save
   INTEGER     :: nbmoy
   CHARACTER(LEN=lchain) :: filediag
   INTEGER      :: iextmindiag,iextmaxdiag,stepidiag,jextmindiag,jextmaxdiag,stepjdiag,kextmindiag,kextmaxdiag,stepkdiag
   INTEGER      :: igextmindiag,igextmaxdiag,jgextmindiag,jgextmaxdiag
   LOGICAL      :: l_turb,l_botstress,l_xe,l_late

   CHARACTER(LEN=lchain) :: name_outdiag_h0,name_outdiag_xe
   CHARACTER(LEN=lchain) :: name_outdiag_hx,name_outdiag_hy

  CONTAINS

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

  SUBROUTINE wave_diag_qeuler_alloc

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE wave_diag_qeuler_alloc  ***
   !&E
   !&E ** Purpose :
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2009-09 (S. Michel, HOCER) correct i dimension of the arrays
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Local declarations

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

   ALLOCATE(diag_wf_advx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_advy(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_corx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_cory(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_vortexx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_vortexy(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_smb1x(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_smb1y(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_grdpx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_grdpy(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_difhx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_difhy(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_difvx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_difvy(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_streamx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_streamy(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_nudtkedz(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_nudpsidz(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_tke(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_psi(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_peps(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_pk(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_prodtke(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_prodpsi(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_dissiptke(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_dissippsi(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_vortx(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_vorty(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_vortz(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_dudt(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_dvdt(1:kmax,limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_xe(limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_botfricx(limin:limax,ljmin:ljmax) )
   ALLOCATE(diag_wf_botfricy(limin:limax,ljmin:ljmax) )




   ALLOCATE(mask_h0_diag(limin:limax,ljmin:ljmax) )
   ALLOCATE(mask_hx_diag(limin:limax,ljmin:ljmax) )
   ALLOCATE(mask_hy_diag(limin:limax,ljmin:ljmax) )
   ALLOCATE(mask_diagx(limin:limax,ljmin:ljmax) )
   ALLOCATE(mask_diagy(limin:limax,ljmin:ljmax) )

   diag_wf_advx(:,:,:)=0.0_rlg
   diag_wf_advy(:,:,:)=0.0_rlg
   diag_wf_corx(:,:,:)=0.0_rlg
   diag_wf_cory(:,:,:)=0.0_rlg
   diag_wf_vortexx(:,:,:)=0.0_rlg
   diag_wf_vortexy(:,:,:)=0.0_rlg
   diag_wf_smb1x(:,:,:)=0.0_rlg
   diag_wf_smb1y(:,:,:)=0.0_rlg
   diag_wf_grdpx(:,:,:)=0.0_rlg
   diag_wf_grdpy(:,:,:)=0.0_rlg
   diag_wf_difhx(:,:,:)=0.0_rlg
   diag_wf_difhy(:,:,:)=0.0_rlg
   diag_wf_difvx(:,:,:)=0.0_rlg
   diag_wf_difvy(:,:,:)=0.0_rlg
   diag_wf_streamx(:,:,:)=0.0_rlg
   diag_wf_streamy(:,:,:)=0.0_rlg
   diag_wf_nudtkedz(:,:,:)=0.0_rlg
   diag_wf_nudpsidz(:,:,:)=0.0_rlg
   diag_wf_tke(:,:,:)=0.0_rlg
   diag_wf_psi(:,:,:)=0.0_rlg
   diag_wf_pk(:,:,:)=0.0_rlg
   diag_wf_peps(:,:,:)=0.0_rlg
   diag_wf_prodtke(:,:,:)=0.0_rlg
   diag_wf_prodpsi(:,:,:)=0.0_rlg
   diag_wf_dissiptke(:,:,:)=0.0_rlg
   diag_wf_dissippsi(:,:,:)=0.0_rlg
   diag_wf_vortx(:,:,:)=0.0_rlg
   diag_wf_vorty(:,:,:)=0.0_rlg
   diag_wf_vortz(:,:,:)=0.0_rlg
   diag_wf_dudt(:,:,:)=0.0_rlg
   diag_wf_dvdt(:,:,:)=0.0_rlg
   diag_wf_xe(:,:)=0.0_rlg
   diag_wf_botfricx(:,:)=0.0_rlg
   diag_wf_botfricy(:,:)=0.0_rlg


  END SUBROUTINE wave_diag_qeuler_alloc

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

  SUBROUTINE wave_diag_qeuler_dealloc

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE wave_diag_qeuler_dealloc  ***
   !&E
   !&E ** Purpose :
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments

   !! * Local declarations

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

   DEALLOCATE(diag_wf_advx)
   DEALLOCATE(diag_wf_advy)
   DEALLOCATE(diag_wf_corx)
   DEALLOCATE(diag_wf_cory)
   DEALLOCATE(diag_wf_vortexx)
   DEALLOCATE(diag_wf_vortexy)
   DEALLOCATE(diag_wf_smb1x)
   DEALLOCATE(diag_wf_smb1y)
   DEALLOCATE(diag_wf_grdpx)
   DEALLOCATE(diag_wf_grdpy)
   DEALLOCATE(diag_wf_difhx)
   DEALLOCATE(diag_wf_difhy)
   DEALLOCATE(diag_wf_difvx)
   DEALLOCATE(diag_wf_difvy)
   DEALLOCATE(diag_wf_streamx)
   DEALLOCATE(diag_wf_streamy)
   DEALLOCATE(diag_wf_nudtkedz)
   DEALLOCATE(diag_wf_nudpsidz)
   DEALLOCATE(diag_wf_tke)
   DEALLOCATE(diag_wf_psi)
   DEALLOCATE(diag_wf_pk)
   DEALLOCATE(diag_wf_peps)
   DEALLOCATE(diag_wf_prodtke)
   DEALLOCATE(diag_wf_prodpsi)
   DEALLOCATE(diag_wf_dissiptke)
   DEALLOCATE(diag_wf_dissippsi)
   DEALLOCATE(diag_wf_vortx)
   DEALLOCATE(diag_wf_vorty)
   DEALLOCATE(diag_wf_vortz)
   DEALLOCATE(diag_wf_dudt)
   DEALLOCATE(diag_wf_dvdt)
   DEALLOCATE(diag_wf_xe)
   DEALLOCATE(diag_wf_botfricx)
   DEALLOCATE(diag_wf_botfricy)


   DEALLOCATE(mask_h0_diag,mask_hx_diag,mask_hy_diag,mask_diagx,mask_diagy)

  END SUBROUTINE wave_diag_qeuler_dealloc

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

  SUBROUTINE wave_diag_qeuler_init

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE wave_diag_qeuler_init  ***
   !&E
   !&E ** Purpose : Initialize diagnostics of dynamical terms
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :  !  2011-04  (P. Lazure)  Original code
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE comvars2d,   ONLY : lon2d_g,lat2d_g,lon2d,lat2d,ierrorlog
   USE parameters,  ONLY : rsh,rlg,riosh,kmax,imin,imax,limin,limax,ljmin,ljmax,  &
                           rg_valmanq_io,ljminp1,riolg,valmanq,lchain
   USE ionc4,       ONLY : ionc4_createfile,ionc4_gatt,ionc4_createvar,           &
                           ionc4_write,ionc4_write_sig,ionc4_write_sigw,          &
                           ionc4_write_xy,ionc4_write_time,ionc4_vatt_chain
   USE outnetcdf,   ONLY : out_netcdf_conv,out_zaxis_w
   USE output,      ONLY : l_out_nc4par,l_out_pack,l_out_zaxis_w
   USE comvars2d,   ONLY : hx,hy,h0,h0fond,riog_valid_min_h0,riog_valid_max_h0,   &
                           name_out_xe,name_out_h0
   USE comvars2d,   ONLY : ig,id,dx_u,dy_v,jb,jh,htx,hty,tdeb3d
   USE comvars3d,   ONLY : sig
   USE_MPI toolmpi, ONLY : ADD_TAG,IPROC,NPROCS,ADD_ALL_MPI_INT,MPI_COMM_MARS

#if defined key_siggen || defined key_gencoord
   USE comsiggen,   ONLY : b_sig,theta_sig,hc_sig,Csu_sig,Csw_sig
#endif

   !! * Local declarations
   CHARACTER(LEN=lchain) :: chain
   REAL(KIND=riosh),DIMENSION(limin:limax,ljmin:ljmax) :: tecrit
   REAL(KIND=riosh),DIMENSION(0:kmax+1)            :: sigecrit
   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=riolg),DIMENSION(imin:imax)                    :: lon
   REAL(KIND=riosh)                            :: riog_valid_min_diag,riog_valid_max_diag
   REAL(KIND=rlg)        :: tool_datosec
   INTEGER :: i,j,IERR_MPI
   INTEGER :: grp_world,grp_zoom,comm_zoom,ncpu,ip
   INTEGER, DIMENSION(:), ALLOCATABLE  :: numrank,listrank
   REAL(KIND=riosh)                                              :: vecrit

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

   filediag=TRIM('diag_wave_qeuler.nc')

   tdebdiagdyn=tool_datosec('01/01/1900 00:00:00')
   tdebdiagdyn=MAX(tdebdiagdyn, tdeb3d)
   tfindiagdyn=tool_datosec('01/01/2870 00:00:00')
   dtsordiagdyn=3.*3600.0_rlg
   dtsordiagdyn=30.0_rlg
!   dtsordiagdyn=0.00001_rlg
   l_moy_save=.TRUE.
   l_moy_save=.FALSE.

   l_late=.TRUE.
   l_xe=.TRUE.
   l_turb=.TRUE.
   l_botstress=.TRUE.

   igextmindiag=imin
   igextmaxdiag=imax
   iextmindiag=imin
   iextmaxdiag=imax
   stepidiag=1
   jextmindiag=jmin
   jextmaxdiag=jmax
   jgextmindiag=jmin
   jgextmaxdiag=jmax
   stepjdiag=1
#if defined key_MPI_2D
   ! If we do not use the parallel writing, the spatial extension
   ! must be bounded by limin and limax
   IF (.NOT.l_out_nc4par) THEN
     DO WHILE ( iextmindiag < limin )
       iextmindiag=iextmindiag+stepidiag
     END DO
     iextmaxdiag=MIN(limax,iextmaxdiag)
     DO WHILE ( jextmindiag < ljmin )
       jextmindiag=jextmindiag+stepjdiag
     END DO
     jextmaxdiag=MIN(ljmax,jextmaxdiag)
   END IF
#endif
   kextmindiag=1  ! should not be changed
   kextmaxdiag=kmax ! id
   stepkdiag=1 ! id

   riog_valid_min_diag=-10000.0_riosh
   riog_valid_max_diag=10000.0_riosh

   name_outdiag_h0='h0'
   name_outdiag_hx='hx'
   name_outdiag_hy='hy'
   name_outdiag_xe='xe'


   name_outdiag_wf_advx='ADVX'
   name_outdiag_wf_advy='ADVY'
   name_outdiag_wf_corx='CORX'
   name_outdiag_wf_cory='CORY'
   name_outdiag_wf_vortexx='VORTEXX'
   name_outdiag_wf_vortexy='VORTEXY'
   name_outdiag_wf_smb1x='SMB1X'
   name_outdiag_wf_smb1y='SMB1Y'
   name_outdiag_wf_grdpx='GRDPX'
   name_outdiag_wf_grdpy='GRDPY'
   name_outdiag_wf_difhx='DIFHX'
   name_outdiag_wf_difhy='DIFHY'
   name_outdiag_wf_difvx='DIFVX'
   name_outdiag_wf_difvy='DIFVY'
   name_outdiag_wf_streamx='STREAMX'
   name_outdiag_wf_streamy='STREAMY'
   name_outdiag_wf_nudtkedz='NudTKEdz'
   name_outdiag_wf_nudpsidz='NudPSIdz'
   name_outdiag_wf_tke='TKE'
   name_outdiag_wf_psi='PSI'
   name_outdiag_wf_pk='PK'
   name_outdiag_wf_peps='PEPS'
   name_outdiag_wf_prodtke='ProdTKE'
   name_outdiag_wf_prodpsi='ProdPSI'
   name_outdiag_wf_dissiptke='DissipTke'
   name_outdiag_wf_dissippsi='DissipPsi'
   name_outdiag_wf_vortx='VORTX'
   name_outdiag_wf_vorty='VORTY'
   name_outdiag_wf_vortz='VORTZ'
   name_outdiag_wf_dudt='DuDt'
   name_outdiag_wf_dvdt='DvDt'
   name_outdiag_wf_botfricx='BotFricx'
   name_outdiag_wf_botfricy='BotFricy'
!
   tnextout=tdebdiagdyn+dtsordiagdyn
!
! Verifications
! -------------

   IF ( jextmindiag > jextmaxdiag .OR. iextmindiag > iextmaxdiag ) THEN
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) ' Routine WAVE_DIAG_INIT in WAVEDIAG.F90'
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) 'Check your domain : because'
     WRITE(ierrorlog,*) ' jextmindiag > jextmaxdiag '
     WRITE(ierrorlog,*) '           or'
     WRITE(ierrorlog,*) ' iextmindiag > iextmaxdiag '
     WRITE(ierrorlog,*) ' '
     WRITE(ierrorlog,*) ' The simulation is stopped '
     WRITE(ierrorlog,*) ' '
     CALL_MPI MPI_FINALIZE(IERR_MPI)
     STOP
   END IF


! Masks
! -----

   mask_diagx(:,:)=.FALSE.
   DO j=jextmindiag,jextmaxdiag,stepjdiag
     DO i=iextmindiag,iextmaxdiag,stepidiag
       IF ( i >= limin .AND. i <= limax .AND. j >= ljmin .AND. j <= ljmax) THEN
         mask_diagx(i,j)=.TRUE.
       END IF
     END DO
   END DO

   mask_diagy(:,:)=.FALSE.
   DO j=jextmindiag,jextmaxdiag,stepjdiag
     DO i=iextmindiag,iextmaxdiag,stepidiag
       IF ( i >= limin .AND. i <= limax .AND. j >= ljmin .AND. j <= ljmax) THEN
         mask_diagy(i,j)=.TRUE.
       END IF
     END DO
   END DO



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

   mask_h0_diag(:,:)=.TRUE.
   DO j=ljmin,ljmax
     DO i=MAX(limin,ig(j)),MIN(limax,id(j))
       IF ( h0_local(i,j)> h0fond ) mask_h0_diag(i,j)=.FALSE.
     END DO
   END DO

! pascal
! mouais pas clair c t affaire
   mask_hx_diag(:,:)=.TRUE.
   DO j=ljminp1,ljmax
     DO i=MAX(limin,ig(j)-1),MIN(limax,id(j))
       IF ( hx_local(i,j) /= -valmanq -h0fond ) mask_hx_diag(i,j)=.FALSE.
     END DO
   END DO
   mask_hx_diag(:,:)=.FALSE.

   mask_hy_diag(:,:)=.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_local(i,j) /= -valmanq -h0fond) ) mask_hy_diag(i,j)=.FALSE.
     END DO
   END DO


   comm_zoom=0
#ifdef key_MPI_2D
     comm_zoom=MPI_COMM_MARS
     ! 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 ( jextmindiag <= jextmaxdiag .AND. iextmindiag <= iextmaxdiag ) 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 ( jextmindiag <= jextmaxdiag .AND. iextmindiag <= iextmaxdiag ) listrank(1)=iproc
     ! all cpus get the number of active cpus plus sleeping cpus (set to 0)
     CALL MPI_ALLGATHER(listrank(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


   ! to work with acive cpus over a small area
   IF ( jextmindiag <= jextmaxdiag .AND. iextmindiag <= iextmaxdiag ) THEN


   sigecrit(:)=sig(:)
   IF ( l_out_nc4par ) THEN
     ! Create a global file because of the parallel writing
     CALL ionc4_createfile(filediag,lon2d_g,lat2d_g,igextmindiag,igextmaxdiag,stepidiag, &
                           jgextmindiag,jgextmaxdiag,stepjdiag, &
                           kextmindiag,kextmaxdiag,stepkdiag,   &
                           kmax,0,l_out_nc4par = .TRUE.,comm_active=comm_zoom, &
                           zaxis_var=sigecrit(1:kmax))
   ELSE
     ! Create a local file bounded by limin,limax,ljmin,ljmax
     CALL_MPI ADD_TAG(filediag,LEN(filediag))
     CALL ionc4_createfile(filediag,lon2d,lat2d,iextmindiag,iextmaxdiag,stepidiag, &
                           jextmindiag,jextmaxdiag,stepjdiag, &
                           kextmindiag,kextmaxdiag,stepkdiag, &
                           kmax,0,zaxis_var=sigecrit(1:kmax))
   END IF

   ! save attribute of the spatial extension of the domain
   ! -----------------------------------------------------
   CALL ionc4_gatt(filediag,'global_imin',iextmindiag)
   CALL ionc4_gatt(filediag,'global_imax',iextmaxdiag)
   CALL ionc4_gatt(filediag,'global_jmin',jextmindiag)
   CALL ionc4_gatt(filediag,'global_jmax',jextmaxdiag)

   ! Add the vertical axis
   CALL out_zaxis_w(filediag,kextmindiag,kextmaxdiag,stepkdiag,comm_zoom=comm_zoom)

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

#if defined key_siggen || defined key_gencoord

   sigecrit(:)=Csu_sig(:)
   CALL ionc4_createvar(filediag,'Csu_sig',' ','C(s) law at mid layer',dims='z', &
                        standard_name='ocean_s_coordinate_function_at_midlayer', &
                        valid_min=-1.0_riosh,valid_max=10000.0_riosh,fill_value= rg_valmanq_io)
   CALL ionc4_write_sig(filediag,'Csu_sig',sigecrit,kmax)
   CALL ionc4_vatt_chain(filediag,'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))')

   sigecrit(:)=Csw_sig(:)
   CALL ionc4_createvar(filediag,'Csw_sig',' ','C(s) law at interface',dims='z', &
                        standard_name='ocean_s_coordinate_function_at_interface', &
                        valid_min= -1.0_riosh,valid_max= 10000.0_riosh,fill_value= rg_valmanq_io)
   !                        l_auxcoordinate=.false.,location_in_cell='w')
   CALL ionc4_write_sigw(filediag,'Csw_sig',sigecrit,kmax)
   CALL ionc4_vatt_chain(filediag,'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(filediag,'hc','m', 'depth over which the resolution is increased', &
                        standard_name='ocean_s_coordinate_thick_of_surface_resolution',dims='xy', &
                        valid_min= 0.0_riosh,valid_max= 1000.0_riosh,fill_value= rg_valmanq_io,&
                        l_out_nc4par=l_out_nc4par,l_pack=l_out_pack )
   tecrit(limin:limax,ljmin:ljmax)=hc_sig(limin:limax,ljmin:ljmax)
   CALL ionc4_write_xy(filediag,'hc',tecrit,limin,limax,ljmin,ljmax,fill_value=rg_valmanq_io,l_out_nc4par=l_out_nc4par)

   CALL ionc4_createvar(filediag,'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(filediag,'b',vecrit)

   CALL ionc4_createvar(filediag,'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(filediag,'theta',vecrit)

#endif

   ! Add attributes to follow OCO netCDF convention
   CALL out_netcdf_conv(filediag,1,comm_zoom,igextmindiag,igextmaxdiag,jgextmindiag,jgextmaxdiag)

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


! Save the bathymetry
! -------------------

! h0
! --

#if ! defined key_agrif
   ! avec agrif on ne peut compiler car The type of the actual argument differs from the type of the dummy argument.   [ARRAY0]
   !  &,  Agrif_tabvars(606) % var% array0 ,  Agrif_tabvars(605) &
   ! Pas de souci si on utilise riog_valid_min_diag a la place de riog_valid_min_h0 !!!!!!!!!!
   ! Cela m exaspere, et c est un pb recurant, voir mail du 3/3/2010 a Romain Ramel (PREVI_OMPMPI)
   ! Laurent peux-tu faire qq chose STP ou donner une explication rationnelle !!!!

   chain='model_sea_floor_depth_below_geoid'
   CALL ionc4_createvar(filediag,name_outdiag_h0,units='m',&
                        long_name='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 = l_out_pack,&
                        l_out_nc4par = l_out_nc4par)
   tecrit(limin:limax,ljmin:ljmax)=h0_local(limin:limax,ljmin:ljmax)
   CALL ionc4_write_xy(filediag,TRIM(name_outdiag_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(filediag,name_outdiag_hx,units='m',&
                        long_name = '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 = l_out_pack,&
                        l_out_nc4par = l_out_nc4par)
   tecrit(:,:)=hx_local(:,:)
   CALL ionc4_write_xy(filediag,TRIM(name_outdiag_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(filediag,name_outdiag_hy,units='m',&
                        long_name = '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 = l_out_pack,&
                        l_out_nc4par = l_out_nc4par)
   tecrit(:,:)=hy_local(:,:)
   CALL ionc4_write_xy(filediag,TRIM(name_outdiag_hy),tecrit,limin,limax,ljmin,ljmax,&
                       fill_value = -REAL(rg_valmanq_io+h0fond,riosh),l_out_nc4par = l_out_nc4par)

#endif

   IF (l_late) THEN

     CALL ionc4_createvar(filediag,name_outdiag_wf_advx,units= 'm s-2',&
                          long_name = 'x_advection_component',&
                          standard_name = 'x_advection_terms',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_advy,units= 'm s-2',&
                          long_name = 'y_advection_component',&
                          standard_name = 'y_advection_terms',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_corx,units= 'm s-2',&
                          long_name = 'x_coriolis_force_component',&
                          standard_name = 'x_coriolis',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_cory,units= 'm s-2',&
                          long_name = 'y_coriolis_force_component',&
                          standard_name = 'y_coriolis',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_grdpx,units= 'm s-2',&
                          long_name = 'x_pressure_gradient_component',&
                          standard_name = 'x_pressure_gradient',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_grdpy,units= 'm s-2',&
                          long_name = 'y_pressure_gradient_component',&
                          standard_name = 'y_pressure_gradient',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_difhx,units= 'm s-2',&
                          long_name = 'x_horizontal_momentum_diffusion_component',&
                          standard_name = 'x_mtm_hdiffusion',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_difhy,units= 'm s-2',&
                          long_name = 'y_horizontal_momentum_diffusion_component',&
                          standard_name = 'y_mtm_hdiffusion',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_difvx,units= 'm s-2',&
                          long_name = 'x_vertical_momentum_diffusion_component',&
                          standard_name = 'x_mtm_vdiffusion',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_difvy,units= 'm s-2',&
                          long_name = 'y_vertical_momentum_diffusion_component',&
                          standard_name = 'y_mtm_vdiffusion',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_dudt,units= 'm s-2',&
                          long_name = 'x_momentum_trends',&
                          standard_name = 'x_mtm_trends',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_dvdt,units= 'm s-2',&
                          long_name = 'y_momentum_trends',&
                          standard_name = 'y_mtm_trends',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

   END IF

   IF (l_late) THEN
     CALL ionc4_createvar(filediag,name_outdiag_wf_vortexx,units= 'm s-2',&
                          long_name = 'x_vortex_force_component',&
                          standard_name = 'x_vortex_force',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_vortexy,units= 'm s-2',&
                          long_name = 'y_vortex_force_component',&
                          standard_name = 'y_vortex_force',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_smb1x,units= 'm s-2',&
                          long_name = 'x_bernouilli_head_component',&
                          standard_name = 'x_bernouillihead',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_smb1y,units= 'm s-2',&
                          long_name = 'y_bernouilli_head_component',&
                          standard_name = 'y_bernouillihead',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_streamx,units= 'm s-2',&
                          long_name = 'x_streaming_component',&
                          standard_name = 'x_streaming',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_streamy,units= 'm s-2',&
                          long_name = 'y_streaming_component',&
                          standard_name = 'y_streaming',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_botfricx,units= 'm s-2',&
                          long_name = 'x_bottom_friction_component',&
                          standard_name = 'x_bott_fric',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'xyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_botfricy,units= 'm s-2',&
                          long_name = 'y_bottom_friction_component',&
                          standard_name = 'y_bott_fric',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'xyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_vortx,units= 's-1',&
                          long_name = 'x_vorticity_component',&
                          standard_name = 'x_vorticity',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_vorty,units= 's-1',&
                          long_name = 'y_vorticity_component',&
                          standard_name = 'y_vorticity',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_vortz,units= 's-1',&
                          long_name = 'z_vorticity_component',&
                          standard_name = 'z_vorticity',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)
   END IF

   IF (l_late) THEN
     CALL ionc4_createvar(filediag,name_outdiag_wf_tke,units= 'm2.s-2',&
                          long_name = 'turbulent_kinetic_energy',&
                          standard_name = 'TKE',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_psi,units= 'm2.s-3',&
                          long_name = 'turbulent_kinetic_energy_second_moment',&
                          standard_name = 'PSI',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_nudtkedz,units= 'm2.s-3',&
                          long_name = 'turbulent_kinetic_energy_vertical_diffusion',&
                          standard_name = 'tke_vdif',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_nudpsidz,units= 'm2.s-3',&
                          long_name = 'second_turbulent_moment_vertical_diffusion',&
                          standard_name = 'psi_vdif',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_prodtke,units= 'm2.s-3',&
                          long_name = 'turbulent_kinetic_production_term',&
                          standard_name = 'tke_prod',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_prodpsi,units= 'm2.s-3',&
                          long_name = 'turbulent_kinetic_second_moment_production_term',&
                          standard_name = 'psi_prod',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)

     CALL ionc4_createvar(filediag,name_outdiag_wf_dissiptke,units= 'm2.s-3',&
                          long_name = 'turbulent_kinetic_dissipation_term',&
                          standard_name = 'tke_dissip',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)
                          
     CALL ionc4_createvar(filediag,name_outdiag_wf_dissippsi,units= 'm2.s-3',&
                          long_name = 'turbulent_kinetic_second_moment_dissipation_term',&
                          standard_name = 'psi_dissip',&
                          fill_value = rg_valmanq_io,&
                          valid_min = riog_valid_min_diag,&
                          valid_max = riog_valid_max_diag,&
                          dims= 'zxyt',&
                          l_pack = l_out_pack,&
                          l_out_nc4par = l_out_nc4par)
   END IF

   IF (l_xe) &
     CALL ionc4_createvar(filediag,name_outdiag_xe,units= 'm',&
                        long_name = 'sea surface height',&
                        standard_name = 'sea surface height',&
                        fill_value = rg_valmanq_io,&
                        valid_min = riog_valid_min_diag,&
                        valid_max = riog_valid_max_diag,&
                        dims= 'xyt',&
                        l_pack = l_out_pack,&
                        l_out_nc4par = l_out_nc4par)
   END IF  ! ends test on iextmin...

  END SUBROUTINE wave_diag_qeuler_init

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

  SUBROUTINE wave_diag_qeuler_output

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE outputwave_diag_qeuler  ***
   !&E
   !&E ** Purpose : Save diagnostics of dynamical terms
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls :
   !&E
   !&E ** Used ij-arrays :
   !&E
   !&E ** Modified variables :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :  !  2011-04  (P. Lazure)  Original code
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Modules used
   USE parameters,  ONLY : rsh,rlg,riosh,kmax,imin,imax,limin,limax,ljmin,ljmax,  &
                           rg_valmanq_io,ljminp1,riolg,valmanq,lchain
   USE outtools,    ONLY : out_write_xytvar,out_write_zxytvar
   USE ionc4,       ONLY : ionc4_write_time,ionc4_sync
   USE output,      ONLY : l_out_nc4par
   USE comvars2d,   ONLY : hx,hy,h0,h0fond,riog_valid_min_h0,riog_valid_max_h0
   USE comvars2d,   ONLY : ig,id,dx_u,dy_v,jb,jh,htx,hty
   USE comvars2d,   ONLY : t,h0fond
   USE comvars3d,   ONLY : sig
   USE comvarp2d,   ONLY : u,v,ssh
   USE_MPI toolmpi, ONLY : MASTER

   !! * Local declarations
   REAL(KIND=riosh),DIMENSION(imin:imax,jmin:jmax) :: tecrit
   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=riolg),DIMENSION(imin:imax)                    :: lon
   REAL(KIND=riolg)                            :: riog_valid_min_diag,riog_valid_max_diag
   INTEGER :: i,j
   CHARACTER (LEN=19) :: tool_sectodat

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

   PRINT_DBG*, 'enter WAVE_DIAG_QEULER_OUTPUT routine'

   IF ( t < tdebdiagdyn .OR. t > tfindiagdyn )  RETURN

!print*, ' prochaine ecriture:',tool_sectodat(tnextout)

   IF (l_moy_save) THEN
     nbmoy=nbmoy+1

     IF (nbmoy ==1) THEN
       tprecmoy=t
       tcumul=0
     ELSE
       dtmoy=t-tprecmoy
       tcumul=tcumul+dtmoy

!      diagx_wf_1_save(:,:,:)=diagx_wf_1_save(:,:,:)+diagx_wf_1(:,:,:)*dtmoy
!      diagx_wf_2_save(:,:,:)=diagx_wf_2_save(:,:,:)+diagx_wf_2(:,:,:)*dtmoy
!      diagx_wf_3_save(:,:,:)=diagx_wf_3_save(:,:,:)+diagx_wf_3(:,:,:)*dtmoy
!      diagx_wf_5_save(:,:,:)=diagx_wf_5_save(:,:,:)+diagx_wf_5(:,:,:)*dtmoy
!      diagx_wf_6_save(:,:,:)=diagx_wf_6_save(:,:,:)+diagx_wf_6(:,:,:)*dtmoy
!      diagx_wf_7_save(:,:,:)=diagx_wf_7_save(:,:,:)+diagx_wf_7(:,:,:)*dtmoy
!      diagx_wf_10_save(:,:,:)=diagx_wf_10_save(:,:,:)+diagx_wf_10(:,:,:)*dtmoy
!      diagx_wf_12_save(:,:,:)=diagx_wf_12_save(:,:,:)+diagx_wf_12(:,:,:)*dtmoy
!      diagx_wf_13_save(:,:,:)=diagx_wf_13_save(:,:,:)+diagx_wf_13(:,:,:)*dtmoy
!      diagy_wf_1_save(:,:,:)=diagy_wf_1_save(:,:,:)+diagy_wf_1(:,:,:)*dtmoy
!      diagy_wf_2_save(:,:,:)=diagy_wf_2_save(:,:,:)+diagy_wf_2(:,:,:)*dtmoy
!      diagy_wf_6_save(:,:,:)=diagy_wf_6_save(:,:,:)+diagy_wf_6(:,:,:)*dtmoy
!      diagy_wf_4_save(:,:,:)=diagy_wf_4_save(:,:,:)+diagy_wf_4(:,:,:)*dtmoy
!      diagy_wf_7_save(:,:,:)=diagy_wf_7_save(:,:,:)+diagy_wf_7(:,:,:)*dtmoy
!      diagy_wf_10_save(:,:,:)=diagy_wf_10_save(:,:,:)+diagy_wf_10(:,:,:)*dtmoy
!      diagy_wf_12_save(:,:,:)=diagy_wf_12_save(:,:,:)+diagy_wf_12(:,:,:)*dtmoy
!      diagy_wf_13_save(:,:,:)=diagy_wf_13_save(:,:,:)+diagy_wf_13(:,:,:)*dtmoy
!      diagy_dudt_save(:,:,:)=diagy_dudt_save(:,:,:)+diagy_dudt*dtmoy
!      diagy_vel_save(:,:,:)=diagy_vel_save(:,:,:)+diagy_vel*dtmoy
!      diag_boz_save(:,:,:)=diag_boz_save(:,:,:)+diag_boz*dtmoy
!      diag_xe_save(:,:)=diag_xe_save(:,:)+ssh(limin:limax,ljmin:ljmax)*dtmoy
     END IF

   ELSE

     nbmoy=1
     tcumul=1._rlg
!    diagx_wf_1_save=diagx_wf_1
!    diagx_wf_2_save=diagx_wf_2
!    diagx_wf_3_save=diagx_wf_3
!    diagx_wf_5_save=diagx_wf_5
!    diagx_wf_6_save=diagx_wf_6
!    diagx_wf_7_save=diagx_wf_7
!    diagx_wf_10_save=diagx_wf_10
!    diagx_wf_12_save=diagx_wf_12
!    diagx_wf_13_save=diagx_wf_13
!    diagy_wf_1_save=diagy_wf_1
!    diagy_wf_2_save=diagy_wf_2
!    diagy_wf_6_save=diagy_wf_6
!    diagy_wf_4_save=diagy_wf_4
!    diagy_wf_7_save=diagy_wf_7
!    diagy_wf_10_save=diagy_wf_10
!    diagy_wf_12_save=diagy_wf_12
!    diagy_wf_13_save=diagy_wf_13
!    diagy_dudt_save=diagy_dudt
!    diagy_vel_save=diagy_vel
!    diag_boz_save=diag_boz
!    diag_xe_save(:,:)=ssh(limin:limax,ljmin:ljmax)
   END IF


   IF ( t >= tnextout) THEN

     IF_MPI (MASTER) THEN
     print*,' WAVE_DIAG_QEULER_OUTPUT : on va ecrire dans fichier :',TRIM(filediag)
     print*,' date:',tool_sectodat(t)
     print*, 'date sortie prevue',tool_sectodat(tnextout)
     print*, 'longueur moyenne (en h)',tcumul/3600.
     ENDIF_MPI

     tnextout=tnextout+dtsordiagdyn

!    diagx_wf_1_save=diagx_wf_1_save/tcumul
!    diagx_wf_2_save=diagx_wf_2_save/tcumul
!    diagx_wf_3_save=diagx_wf_3_save/tcumul
!    diagx_wf_5_save=diagx_wf_5_save/tcumul
!    diagx_wf_6_save=diagx_wf_6_save/tcumul
!    diagx_wf_7_save=diagx_wf_7_save/tcumul
!    diagx_wf_10_save=diagx_wf_10_save/tcumul
!    diagx_wf_12_save=diagx_wf_12_save/tcumul
!    diagx_wf_13_save=diagx_wf_13_save/tcumul
!    diagy_wf_1_save=diagy_wf_1_save/tcumul
!    diagy_wf_2_save=diagy_wf_2_save/tcumul
!    diagy_wf_6_save=diagy_wf_6_save/tcumul
!    diagy_wf_4_save=diagy_wf_4_save/tcumul
!    diagy_wf_7_save=diagy_wf_7_save/tcumul
!    diagy_wf_10_save=diagy_wf_10_save/tcumul
!    diagy_wf_12_save=diagy_wf_12_save/tcumul
!    diagy_wf_13_save=diagy_wf_13_save/tcumul
!    diagy_dudt_save=diagy_dudt_save/tcumul
!    diagy_vel_save=diagy_vel_save/tcumul
!    diag_boz_save=diag_boz_save/tcumul
!    diag_xe_save=diag_xe_save/tcumul

     CALL ionc4_write_time(filediag,0,t)

     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_advx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_advx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_advy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_advy(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_corx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_corx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_cory,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_cory(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_grdpx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_grdpx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_grdpy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_grdpy(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_difhx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_difhx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_difhy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_difhy(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_difvx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_difvx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_difvy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_difvy(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_dudt,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_dudt(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_dvdt,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_dvdt(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_vortx,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_vortx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_vorty,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_vorty(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_vortz,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_vortz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_xe) CALL out_write_xytvar(filediag,name_outdiag_xe,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_xe(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
!
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_vortexx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_vortexx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_vortexy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_vortexy(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_smb1x,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_smb1x(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_smb1y,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_smb1y(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_xytvar(filediag,name_outdiag_wf_botfricx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_botfricx(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_xytvar(filediag,name_outdiag_wf_botfricy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_botfricy(limin:limax,ljmin:ljmax),0,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_streamx,rg_valmanq_io,mask_hx_diag, &
                                       diag_wf_streamx(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_streamy,rg_valmanq_io,mask_hy_diag, &
                                       diag_wf_streamy(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
!
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_tke,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_tke(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_psi,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_psi(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_nudtkedz,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_nudtkedz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_nudpsidz,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_nudpsidz(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_prodtke,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_prodtke(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_prodpsi,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_prodpsi(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_dissiptke,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_prodtke(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)
     IF (l_late) CALL out_write_zxytvar(filediag,name_outdiag_wf_dissippsi,rg_valmanq_io,mask_h0_diag, &
                                       diag_wf_prodpsi(:,limin:limax,ljmin:ljmax),0,1,l_out_nc4par = l_out_nc4par)


     ! To write the data on the disk and not loose data in case of run crash
     CALL ionc4_sync(filediag)

!    nbmoy=0
!    diagx_wf_1_save=0.0_rsh
!    diagx_wf_2_save=0.0_rsh
!    diagx_wf_3_save=0.0_rsh
!    diagx_wf_5_save=0.0_rsh
!    diagx_wf_6_save=0.0_rsh
!    diagx_wf_7_save=0.0_rsh
!    diagx_wf_10_save=0.0_rsh
!    diagx_wf_12_save=0.0_rsh
!    diagx_wf_13_save=0.0_rsh
!    diagy_wf_1_save=0.0_rsh
!    diagy_wf_2_save=0.0_rsh
!    diagy_wf_4_save=0.0_rsh
!    diagy_wf_6_save=0.0_rsh
!    diagy_wf_7_save=0.0_rsh
!    diagy_wf_10_save=0.0_rsh
!    diagy_wf_12_save=0.0_rsh
!    diagy_wf_13_save=0.0_rsh
!    diagy_dudt_save=0.0_rsh
!    diagy_vel_save=0.0_rsh
!    diag_boz_save=0.0_rsh
!    diag_xe_save=0.0_rsh

   END IF

   PRINT_DBG*, 'exit WAVE_DIAG_QEULER_OUTPUT routine'

  END SUBROUTINE wave_diag_qeuler_output

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

  END MODULE wavediag
