MODULE dynkeg
   !!======================================================================
   !!                       ***  MODULE  dynkeg  ***
   !! Ocean dynamics:  kinetic energy gradient trend
   !!======================================================================
   
   !!----------------------------------------------------------------------
   !!   dyn_keg      : update the momentum trend with the horizontal tke
   !!----------------------------------------------------------------------
   !! * Modules used
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE in_out_manager  ! I/O manager
   USE trddyn_oce      ! dynamics trends diagnostics

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC dyn_keg                ! routine called by step.F90
   
   !! * Substitutions
#  include "vectopt_loop_substitute.h90"
   !!---------------------------------------------------------------------------------
   !!   OPA 9.0 , LODYC-IPSL  (2003)
   !!---------------------------------------------------------------------------------

CONTAINS

   SUBROUTINE dyn_keg( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_keg  ***
      !!
      !! ** Purpose :   Compute the now momentum trend due to the horizontal
      !!      gradient of the horizontal kinetic energy and add it to the 
      !!      general momentum trend.
      !!
      !! ** Method  :   Compute the now horizontal kinetic energy:
      !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
      !!      Take its horizontal gradient and add it to the general momentum
      !!      trend (ua,va).
      !!         ua = ua - 1/e1u di[ zhke ]
      !!         va = va - 1/e2v dj[ zhke ]
      !!
      !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
      !!             - Save the trends in (utrd,vtrd) ('key_trddyn')
      !!
      !! History :
      !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
      !!   7.0  !  97-05  (G. Madec)  Split dynber into dynkeg and dynhpg
      !!   9.0  !  02-07  (G. Madec)  F90: Free form and module
      !!----------------------------------------------------------------------
      !! * modules used     
      USE oce, ONLY :   zhke => ta      ! use ta as 3D workspace   

      !! * Arguments
      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index

      !! * Local declarations
      INTEGER  ::   ji, jj, jk          ! dummy loop indices
      REAL(wp) ::   zua, zva, zu, zv    ! temporary scalars
#if defined key_trddyn_new
      REAL(wp) ::   zuu, zvv            ! temporary scalars
      REAL(wp), DIMENSION(jpi,jpj) ::   &
         zvke, zuke                     ! temporary workspace
#endif
      !!----------------------------------------------------------------------

      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
         IF(lwp) WRITE(numout,*) '~~~~~~~'
      ENDIF
      
      !                                                ! ===============
      DO jk = 1, jpkm1                                 ! Horizontal slab
         !                                             ! ===============
         ! Horizontal kinetic energy at T-point
         DO jj = 2, jpj
            DO ji = fs_2, jpi   ! vector opt.
               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
                            + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
                            + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
               zhke(ji,jj,jk) = zv + zu
#if defined key_trddyn_new
               zvke(ji,jj) = zv
               zuke(ji,jj) = zu
#endif
            END DO  
         END DO  
         
         ! Horizontal gradient of Horizontal kinetic energy
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               ! gradient of kinetic energy
               zua = -( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
               zva = -( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
               ! add to the general momentum trends
               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
               va(ji,jj,jk) = va(ji,jj,jk) + zva
#if defined key_trddyn
               ! add to the general momentum trends
               utrd(ji,jj,jk,2) = zua
               vtrd(ji,jj,jk,2) = zva
#endif
#if defined key_trddyn_new
               zuu = -( zuke(ji+1,jj  ) - zuke(ji,jj) ) / e1u(ji,jj)
               zvv = -( zvke(ji  ,jj+1) - zvke(ji,jj) ) / e2v(ji,jj)
               utrd(ji,jj,jk,2) = zua - zuu
               vtrd(ji,jj,jk,3) = zva - zvv
               utrd(ji,jj,jk,3) = zuu
               vtrd(ji,jj,jk,2) = zvv
#endif
            END DO 
         END DO
         !                                             ! ===============
      END DO                                           !   End of slab
      !                                                ! ===============

      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
         zua = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
         zva = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
         WRITE(numout,*) ' keg  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
         u_ctl = zua   ;   v_ctl = zva
      ENDIF

   END SUBROUTINE dyn_keg

   !!======================================================================
END MODULE dynkeg