  !*******************************************************************

  SUBROUTINE diagnostic(u,v,uz,vz,xe,sal,temp,        &
#ifdef key_substance
                                         cv_wat,      &
#endif
                                         ect,boz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE diagnostic  ***
   !&E
   !&E ** Purpose : This routine is dedicated to users : they can
   !&E              implement any diagnoses from prognostic variables
   !&E              sal, temp, ect , boz, u, v, xe, uz, vz...
   !&E
   !&E ** Description : comments on variables to be used by the user
   !&E              to implement his own diagnosis.
   !&E
   !&E    available variables in each cell
   !&E
   !&E       ---- V(i,j) --- phi(i,j)
   !&E       ^                  |
   !&E       |                  |
   !&E       |                  |
   !&E       |     XE(i,j)     U(i,j)
   !&E       |                  |
   !&E       |                  |
   !&E       <---dx_v(i,j-1)--->|
   !&E
   !&E    For a cell of indexes (i,j) in the grid,
   !&E    the meridional velocity V(i,j) or VZ(k,i,j),
   !&E    the zonal velocity U(i,j) or UZ(k,i,j),
   !&E    the sea surface heigth XE(i,j)
   !&E    or the tracers TEMP(k,i,j), SAL(k,i,j) or cv_wat(k,i,j)
   !&E    are not estimated at the same location.
   !&E
   !&E
   !&E       |<-- dx_v(i,j)--->|            |--------|--------|
   !&E       |                 |            |        |        |
   !&E       |                 |            |        |        |
   !&E       |<----dx(i,j)---->|            |     dy(i,j)   dy_u(i,j)
   !&E       |                 |            |        |        |
   !&E       |                 |            |        |        |
   !&E       |-----------------|            |--------|--------|
   !&E
   !&E
   !&E    In the middle of the cell,
   !&E    the zonal step (in meter) or the width of the cell is dx(i,j)
   !&E    the meridional step (in meter) or the height of the cell is dy(i,j).
   !&E    At the U location,
   !&E    the zonal step (in meter) is dx_u(i,j)
   !&E    the meridional step (in meter) is dy_u(i,j).
   !&E    At the V location,
   !&E    the zonal step (in meter) is dx_v(i,j)
   !&E    the meridional step (in meter) is dy_v(i,j).
   !&E
   !&E
   !&E    surf(i,j) is the surface of the cell (i,j)
   !&E    hex(i,j) is the flowing thickness at location U(i,j)
   !&E    hey(i,j) is the flowing thickness at location V(i,j)
   !&E    h0(i,j) +xe(i,j) is the total water depth at the center of the cell
   !&E
   !&E                                                                          ^
   !&E                                                                          |
   !&E  o--- wz(k,i,j) ----o         o--- wz(k,i,j) ----o    o--- wz(k,i,j)---dsigw(k)
   !&E                     |                            ^                       |
   !&E                     |                            |                       |
   !&E       sal(k,i,j)  uz(k,i,j)                    dsigu(k)     temp(k,i,j)    o
   !&E                     |                            |                       ^
   !&E                     |                            |                       |
   !&E  o-- wz(k-1,i,j) ---o          -- wz(k-1,i,j) ---o    o-- wz(k-1,i,j)--dsigw(k-1)
   !&E                     |                            ^                       |
   !&E                     |                            |                       |
   !&E      sal(k-1,i,j) uz(k-1,i,j)               dsigu(k-1)      sal(k-1,i,j)   o
   !&E                     |                            |                       |
   !&E                     |                            |                       |
   !&E
   !&E
   !&E
   !&E    sal(k,i,j) is the salinity in cell(i,j) (middle of the cell)
   !&E    temp(k,i,j) is the remperature in cell(i,j) (middle of the cell)
   !&E    cv_wat(l,k,i,j) is the tracer concentration l in cell(i,j) (middle of the cell)
   !&E    wz(k,i,j) is the (sigma) vertical velocity at the interface between layers
   !&E              k and k+1 (middle of the cell)
   !&E
   !&E    0.5*dt is the time step (dt is the time step between two consecutive estimates
   !&E             or column)
   !&E
   !&E ** Called by : main, step, Agrif_User
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !                                 Original code
   !&E       !  2007-09  (V. Garnier)          homogeneneity of dx dy...
   !&E       !  2008-05  (V. Garnier T. Odaka) Introduction of MPI
   !&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
   !&E       !  2009-03  (V. Garnier) MPI2D
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------

#include "toolcpp.h"

   !! * Used modules
   USE parameters
   USE comvars2d,  ONLY : dt,nbouc,t,h0,surf,dx_f,dy_f,cdate,dx
   USE comvars2d,  ONLY : frot_lin,tdeb
   USE comvars3d,  ONLY : dsigu,sig,nzinit
   USE comdiag,    ONLY : file_diag,tdebdiag,tfindiag,t0
   USE_MPI toolmpi,ONLY : MASTER,MAX_ALL_MPI_RSH,MIN_ALL_MPI_RSH,ADD_ALL_MPI_RSH
   USE toolgeom,   ONLY : f_dzu,f_zu

   IMPLICIT NONE

  !! * Arguments

   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in)             :: u,v
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in)             :: xe
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in)        :: uz,vz
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in)        :: sal,temp
   REAL(KIND=rsh), DIMENSION(kmax,limin:limax,ljmin:ljmax), INTENT(in)                :: ect
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in)        :: boz
#ifdef key_substance
   REAL(KIND=rsh), DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: cv_wat
#endif

  !! * Local declarations

   INTEGER,PARAMETER :: isortie=4
   INTEGER           :: i,j,k
   REAL(KIND=rsh)    :: ke,ape,surface,ens,smin,smax,tool_latlon2j
   REAL(KIND=rsh)    :: curl,d,z,roa,dz,dv,ro,tool_latlon2i
   REAL(KIND=rsh)    :: omega,periode,zeta0
   REAL(KIND=rlg)                                            :: dtm,rms
   REAL(KIND=rsh), DIMENSION(kmax) :: utot
   COMPLEX(KIND=rsh) :: ccosh,csinh,gamma,ub,B,A,Ubar,gammaH,gammaz,iomegat,ic=(0,1)


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

   dtm=0.5_rlg*dt

   IF (nbouc == 0) THEN
     t0=t
     ! file_diag=trim(root)//trim(file_diag)//trim(suffix)//'.txt'
     IF_MPI (MASTER) THEN
     PRINT*,'CREATE DIAGNOSTICS OUTPUT FILE : ',TRIM(file_diag)
     OPEN(55,file=file_diag,form='formatted')
     ENDIF_MPI
   END IF

   IF (t >= tdebdiag .AND. t <= tfindiag ) THEN


   !! test cases part
     IF (l_testcase) THEN

       IF (num_testcase==1) THEN

         ! kinetic energy
         roa=1025.0_rsh+0.78_rsh*(34.85_rsh-33.75_rsh)
         ke=0.0_rsh
         DO j=ljminp1,ljmax
           DO i=liminp1,limax
             IF (h0(i,j) > -valmanq) THEN
               d=h0(i,j)+xe(i,j)
               DO k=1,kmax
                 dz=f_dzu(h0(i,j),xe(i,j),k,i,j)
                 ke=ke+0.5_rsh*roa*surf(i,j)*dz*(uz(k,i,j)**2+vz(k,i,j)**2)
               END DO
             END IF
           END DO
         END DO
         CALL_MPI ADD_ALL_MPI_RSH(ke)

         ! available potential energy
         ape=0.0_rsh
         DO j=ljminp1,ljmax
           DO i=liminp1,limax
             IF (h0(i,j) > -valmanq) THEN
               d=h0(i,j)+xe(i,j)
               DO k=1,kmax
                 ! immersion from the sea surface. The water thickness is xe - f_zu (f_zu being negative)
                 z=-(f_zu(h0(i,j),xe(i,j),k,i,j)-xe(i,j))
                 dz=f_dzu(h0(i,j),xe(i,j),k,i,j)
                 dv=dz*surf(i,j)
                 ro=1025.0_rsh+0.78_rsh*(sal(k,i,j)-33.75_rsh)
                 ape=ape+grav*dv*z*(ro-roa)
               END DO
             END IF
           END DO
         END DO
         CALL_MPI ADD_ALL_MPI_RSH(ape)

         ! enstrophy
         ens=0.0_rsh
         DO j=ljminp1,ljmaxm1
           DO i=liminp1,limaxm1
             IF (h0(i,j) > -valmanq) THEN
               d=h0(i,j)+xe(i,j)
               DO k=1,kmax
                 dz=f_dzu(h0(i,j),xe(i,j),k,i,j)
                 dv=dz*surf(i,j)
                 curl=(vz(k,i+1,j)-vz(k,i,j))/dx_f(i,j)-(uz(k,i,j+1)-uz(k,i,j))/dy_f(i,j)
                 ens=ens+dv*curl*curl
               END DO
             END IF
           END DO
         END DO
         CALL_MPI ADD_ALL_MPI_RSH(ens)

         ! surface of the fresh eddy
         surface=0.0_rsh
         DO j=ljminp1,ljmax
           DO i=liminp1,limax
              IF (sal(kmax,i,j) < 34.839_rsh) surface=surface+surf(i,j)
           END DO
         END DO
         CALL_MPI ADD_ALL_MPI_RSH(surface)

         ! Minimum and maximum of salinity
         smin=40.0_rsh
         smax=0.0_rsh
         DO j=ljminp1,ljmax
           DO i=liminp1,limax
             DO k=1,kmax
               smin=MIN(smin,sal(k,i,j))
               smax=MAX(smax,sal(k,i,j))
             END DO
           END DO
         END DO
         CALL_MPI MAX_ALL_MPI_RSH(smax)
         CALL_MPI MIN_ALL_MPI_RSH(smin)

         IF_MPI (MASTER) THEN
         WRITE(55,155) (t-t0)/86400.0_rlg,cdate,ke,ape,ens,surface,smin,smax
         ENDIF_MPI

       ELSE IF(num_testcase==9) THEN

#ifdef key_toseelater
         periode=44712.0_rlg
         omega=2.0_rsh*pi/periode
         zeta0=0.001_rsh
         Ubar=ic*omega*zeta0/h0(3,2)*SUM(dx(2:isortie,2))
         gamma=(1-ic)*sqrt(omega/(2.0_rsh*nzinit))
         gammaH=gamma*h0(3,2)
         B=frot_lin*ccosh(gammaH)-gamma*nzinit*csinh(gammaH)*(ic*frot_lin/(omega*h0(3,2))-1.0_rsh)
         A=-Ubar*frot_lin/B
         ub=-Ubar/(h0(3,2)*omega)*(frot_lin*h0(3,2)*omega*ccosh(gammaH) &
                                        -ic*frot_lin*gamma*nzinit*csinh(gammaH))/B
         DO k=1,kmax
           z=h0(3,2)*(-1+REAL(k-1)/REAL(kmax-1))
           gammaz=gamma*z
           iomegat=-ic*omega*(t-tdeb)
           utot(k)=REAL((A*ccosh(gammaz)+ic*frot_lin/(omega*h0(3,2))*(ub+Ubar)+Ubar)*EXP(iomegat),rsh)
         END DO
         rms=SUM(ABS(utot(:)-uz(:,isortie,2)))/SUM(ABS(utot(:)))
         WRITE(55,155) (t-t0)/86400.0_rlg,cdate,rms,u(isortie,2),REAL(Ubar*EXP(iomegat),rsh),uz(:,isortie,2),utot(:)
#endif
       END IF   ! ends test on num_testcase

   !! general part
     ELSE     ! not a test case

           ! Here the user can implement any diagnoses, the variable xe,u,uz,sal... are global
           ! In order to save time, the user should comment all the lines  ALLOCATE... CALL correct_rsh...
           ! DEALLOCATE...  treating variables he does not use




#ifdef key_MPI_2D
   ! a ajouter pour sortir un fichier CALL ADD_TAG(filename,LEN(filename))
#endif


     END IF   ! ends test on test case

   END IF   ! ends test on period


  !! write diagnoses of integrated variables
!  IF (t >= tfindiag) THEN
!  END IF

  155 FORMAT(e15.7,1x,a19,1000(e13.5,1x),i4)

  END SUBROUTINE diagnostic

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

  COMPLEX FUNCTION csinh(zz)
   USE parameters
   complex(kind=rsh),INTENT(in) :: zz
   complex(kind=rsh)            :: zzm
   zzm=-zz
   csinh=(EXP(zz)-EXP(zzm))/2.0_rsh
  END FUNCTION csinh

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

  COMPLEX FUNCTION ccosh(zz)
   USE parameters
   complex(kind=rsh),INTENT(in) :: zz
   complex(kind=rsh)            :: zzm
   zzm=-zz
   ccosh=(EXP(zz)+EXP(zzm))/2.0_rsh
  END FUNCTION ccosh
