MODULE diamke
   !!==============================================================================
   !!                       ***  MODULE  diacfl  ***
   !! Output Kinetic Energy diagnostics to ocean output
   !!==============================================================================
   !! History :  1.0  !  2010-03  (E. Blockley)  Original code
   !!----------------------------------------------------------------------
#if defined key_diamke
   !!----------------------------------------------------------------------
   !!   dia_mke        : 
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and active tracers
   USE dom_oce         ! ocean space and time domain
   USE lib_mpp         ! distribued memory computing
   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
   USE in_out_manager  ! I/O manager

   IMPLICIT NONE
   PRIVATE

   LOGICAL, PUBLIC, PARAMETER   ::   lk_diamke = .TRUE.    ! KE diag flag

   INTEGER  :: numke                                       ! outfile unit
   CHARACTER(LEN=50) :: clname="ke_diagnostics.ascii"      ! ascii filename

   PUBLIC   dia_mke       ! routine called by step.F90

   !! * Substitutions
#  include "domzgr_substitute.h90"
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
   !! $Id$
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------


CONTAINS


   SUBROUTINE dia_mke ( )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dia_mke  ***
      !!
      !! ** Purpose :  Compute the total Kinetic Energy at each depth level,  
      !!               output to ascii file 'ke_diagnostics.ascii' and output
      !!               total KE to numout
      !!----------------------------------------------------------------------

      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmke        ! Kinetic Energy array
      REAL(wp), DIMENSION(jpk)         :: zmke_h      ! Total horizontal Kinetic Energy array

      INTEGER  :: ji, jj, jk                          ! dummy loop indices


      zmke(:,:,:) = 0.0
      zmke_h(:)   = 0.0

      IF( lwp ) THEN
         WRITE(numout,*)
         WRITE(numout,*) 'dia_mke     : Outputting Kinetic Energy diagnostics to '//TRIM(clname)
         WRITE(numout,*) '~~~~~~~~~~~~'
         WRITE(numout,*)

         ! create output ascii file
         CALL ctl_opn( numke, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, 1 )
         WRITE(numke,*) 'Level      Depth   KE per unit depth'
         WRITE(numke,*) '            (m)      (tera J / m)   '
         WRITE(numke,*) '************************************'
      ENDIF

      ! calculate Kinetic Energy at T points
      DO jk = 1, jpk
         DO jj = 2, jpj
            DO ji = 2, jpi

               zmke(ji,jj,jk) = 0.5 * rhop(ji,jj,jk) *                        &  ! 1/2 * density
                           &    e1t(ji,jj) * e2t(ji,jj) * e3t(ji,jj,jk) *     &  ! T point volume
                           &    ( 0.25 * (un(ji-1,jj,jk) + un(ji,jj,jk))**2   &  ! U^2 at T point
                           &    + 0.25 * (vn(ji-1,jj,jk) + vn(ji,jj,jk))**2 )    ! V^2 at T point

            END DO
         END DO         
      END DO
      CALL lbc_lnk( zmke(:,:,:), 'T', 1. )

      ! sum over each level and write KE per unit depth profiles to ascii file
      DO jk = 1, jpk-1
         zmke_h(jk) = SUM( zmke(:,:,jk) * tmask(:,:,jk) )
         IF( lk_mpp )  CALL  mpp_sum( zmke_h(jk) )
         IF( lwp ) WRITE(numke,FMT='(3x,i3,4x,f7.2,10x,f12.4)') jk, fsdept(1,1,jk+1), zmke_h(jk) * 1.0e-12_wp / e3t(1,1,jk)
      END DO

      IF( lwp ) THEN
         ! write total out to ascii file
         WRITE(numout,FMT='(12x,a31,f12.7,a23)') 'Total Kinetic Energy at nitend ', SUM( zmke_h(:) ) * 1.0e-15_wp, ' * 1e15 J (peta Joules)' 
         WRITE(numout,*)

         ! and ocean output file
         WRITE(numke,*)
         WRITE(numke,*) 'Total Kinetic Energy at nitend :' 
         WRITE(numke,FMT='(5x,f12.7,a23)') SUM( zmke_h(:) ) * 1.0e-15_wp,' * 1e15 J (peta Joules)'
         CLOSE(numke)
      ENDIF

   END SUBROUTINE dia_mke

#else
   !!----------------------------------------------------------------------
   !!   Default option :                                       Dummy Module
   !!----------------------------------------------------------------------
   LOGICAL, PUBLIC, PARAMETER   ::   lk_diamke = .FALSE.   ! KE diag flag
CONTAINS
   SUBROUTINE dia_mke( )           ! Empty routine
      WRITE(*,*) 'dia_mke: : You should not have seen this print! error?'
   END SUBROUTINE dia_mke
#endif

   !!======================================================================
END MODULE diamke