Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5327)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5328)
@@ -7,7 +7,7 @@
    !!            7.0  !  1997-05  (G. Madec)  Split dynber into dynkeg and dynhpg
    !!  NEMO      1.0  !  2002-07  (G. Madec)  F90: Free form and module
-   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option
+   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option 
    !!----------------------------------------------------------------------
-
+   
    !!----------------------------------------------------------------------
    !!   dyn_keg      : update the momentum trend with the horizontal tke
@@ -29,15 +29,15 @@
 
    PUBLIC   dyn_keg    ! routine called by step module
-
+   
    INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2  = 0   !: 2nd order centered scheme (standard scheme)
    INTEGER, PARAMETER, PUBLIC  ::   nkeg_HW  = 1   !: Hollingsworth et al., QJRMS, 1983
    !
    REAL(wp) ::   r1_48 = 1._wp / 48._wp   !: =1/(4*2*6)
-
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.6 , NEMO Consortium (2015)
-   !! $Id$
+   !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -49,9 +49,9 @@
       !!
       !! ** Purpose :   Compute the now momentum trend due to the horizontal
-      !!      gradient of the horizontal kinetic energy and add it to the
+      !!      gradient of the horizontal kinetic energy and add it to the 
       !!      general momentum trend.
       !!
-      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that
-      !!      conserve kinetic energy. Compute the now horizontal kinetic energy
+      !! ** Method  : * kscheme = nkeg_C2 : 2nd order centered scheme that 
+      !!      conserve kinetic energy. Compute the now horizontal kinetic energy 
       !!         zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
       !!              * kscheme = nkeg_HW : Hollingsworth correction following
@@ -59,5 +59,5 @@
       !!         zhke = 1/6 [ mi-1(  2 * un^2 + ((un(j+1)+un(j-1))/2)^2  )
       !!                    + mj-1(  2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2  ) ]
-      !!
+      !!      
       !!      Take its horizontal gradient and add it to the general momentum
       !!      trend (ua,va).
@@ -72,15 +72,15 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
-      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme
+      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zu, zv       ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:)   :: zhke
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
       !
-      CALL wrk_alloc( jpi,jpj,   zhke )
+      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
       !
       IF( kt == nit000 ) THEN
@@ -92,16 +92,14 @@
       IF( l_trddyn ) THEN           ! Save ua and va trends
          CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
-         ztrdu(:,:,:) = ua(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:)
+         ztrdu(:,:,:) = ua(:,:,:) 
+         ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
-
-      zhke(:,:) = 0._wp
-
-      DO jk = 1, jpkm1
-
-         SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
-         !
-         CASE ( nkeg_C2 )                          !--  Standard scheme  --!
-         
+      
+      zhke(:,:,jpk) = 0._wp
+      
+      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
+      !
+      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
+         DO jk = 1, jpkm1
             DO jj = 2, jpj
                DO ji = fs_2, jpi   ! vector opt.
@@ -110,11 +108,12 @@
                   zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
                      &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
-                  zhke(ji,jj) = 0.25_wp * ( zv + zu )
-               END DO
+                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+               END DO  
             END DO
+         END DO
          !
-         CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
-
-            DO jj = 2, jpjm1
+      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1       
                DO ji = fs_2, jpim1   ! vector opt.
                   zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
@@ -127,18 +126,19 @@
                      &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
                      &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
-                  zhke(ji,jj) = r1_48 * ( zv + zu )
-               END DO
+                  zhke(ji,jj,jk) = r1_48 * ( zv + zu )
+               END DO  
             END DO
-            CALL lbc_lnk( zhke, 'T', 1. )
-            !
-         END SELECT
+         END DO
+         CALL lbc_lnk( zhke, 'T', 1. )
          !
+      END SELECT
+      !
+      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ) - zhke(ji,jj) ) / e1u(ji,jj)
-               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1) - zhke(ji,jj) ) / e2v(ji,jj)
-            END DO
+               ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+               va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+            END DO 
          END DO
-
       END DO
       !
@@ -153,5 +153,5 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi,jpj,   zhke )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
