Index: branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 4737)
+++ branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5837)
@@ -4,7 +4,8 @@
    !! Ocean dynamics:  kinetic energy gradient trend
    !!======================================================================
-   !! 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
+   !! History :  1.0  !  1987-09  (P. Andrich, M.-A. Foujols)  Original code
+   !!            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 
    !!----------------------------------------------------------------------
    
@@ -14,7 +15,9 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE trdmod          ! ocean dynamics trends 
-   USE trdmod_oce      ! ocean variables trends
+   USE trd_oce         ! trends: ocean variables
+   USE trddyn          ! trend manager: dynamics
+   !
    USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         ! MPP library
    USE prtctl          ! Print control
@@ -27,8 +30,13 @@
    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.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.6 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -36,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_keg( kt )
+   SUBROUTINE dyn_keg( kt, kscheme )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_keg  ***
@@ -44,6 +52,12 @@
       !!      general momentum trend.
       !!
-      !! ** Method  :   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
+      !!      Arakawa (2001). The now horizontal kinetic energy is given by:
+      !!         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).
@@ -52,8 +66,12 @@
       !!
       !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
-      !!             - save this trends (l_trddyn=T) for post-processing
+      !!             - send this trends to trd_dyn (l_trddyn=T) for post-processing
+      !!
+      !! ** References : Arakawa, A., International Geophysics 2001.
+      !!                 Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      !!
+      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
+      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
+      !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       REAL(wp) ::   zu, zv       ! temporary scalars
@@ -62,40 +80,60 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_keg')
+      IF( nn_timing == 1 )   CALL timing_start('dyn_keg')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zhke )
+      CALL wrk_alloc( jpi,jpj,jpk,   zhke )
       !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend'
+         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
          IF(lwp) WRITE(numout,*) '~~~~~~~'
       ENDIF
 
       IF( l_trddyn ) THEN           ! Save ua and va trends
-         CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         CALL wrk_alloc( jpi,jpj,jpk,   ztrdu, ztrdv )
          ztrdu(:,:,:) = ua(:,:,:) 
          ztrdv(:,:,:) = va(:,:,:) 
       ENDIF
       
-      !                                                ! ===============
-      DO jk = 1, jpkm1                                 ! Horizontal slab
-         !                                             ! ===============
-         DO jj = 2, jpj         ! Horizontal kinetic energy at T-point
-            DO ji = fs_2, jpi   ! vector opt.
-               zu = 0.25 * (  un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-                  &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)  )
-               zv = 0.25 * (  vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-                  &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)  )
-               zhke(ji,jj,jk) = zv + zu
-!!gm simplier coding  ==>>   ~ faster
-!    don't forget to suppress local zu zv scalars
-!               zhke(ji,jj,jk) = 0.25 * (   un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-!                  &                      + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
-!                  &                      + vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-!                  &                      + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   )
-!!gm end <<==
-            END DO  
-         END DO  
-         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
+      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.
+                  zu =    un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
+                     &  + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
+                  zv =    vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
+                     &  + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
+                  zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+               END DO  
+            END DO
+         END DO
+         !
+      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)    &
+                     &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
+                     &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
+                     &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
+                     !
+                  zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
+                     &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
+                     &  +      ( 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,jk) = r1_48 * ( zv + zu )
+               END DO  
+            END DO
+         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  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
@@ -103,34 +141,11 @@
             END DO 
          END DO
-!!gm idea to be tested  ==>>   is it faster on scalar computers ?
-!         DO jj = 2, jpjm1       ! add the gradient of kinetic energy to the general momentum trends
-!            DO ji = fs_2, fs_jpim1   ! vector opt.
-!               ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj  ,jk) * un(ji+1,jj  ,jk)   &
-!                  &                                   + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk)   &
-!                  &                                   + vn(ji+1,jj  ,jk) * vn(ji+1,jj  ,jk)   &
-!                  !
-!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-!                  &                                   - vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)   &
-!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e1u(ji,jj)
-!                  !
-!               va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * (   un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk)   &
-!                  &                                   + un(ji  ,jj+1,jk) * un(ji  ,jj+1,jk)   &
-!                  &                                   + vn(ji  ,jj+1,jk) * vn(ji  ,jj+1,jk)   &
-!                  !
-!                  &                                   - un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)   &
-!                  &                                   - un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)   &
-!                  &                                   - vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)   ) / e2v(ji,jj)
-!            END DO 
-!         END DO
-!!gm en idea            <<==
-         !                                             ! ===============
-      END DO                                           !   End of slab
-      !                                                ! ===============
-
-      IF( l_trddyn ) THEN      ! save the Kinetic Energy trends for diagnostic
+      END DO
+      !
+      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
          ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
          ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
-         CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )
-         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
+         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
       ENDIF
       !
@@ -138,7 +153,7 @@
          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zhke )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zhke )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('dyn_keg')
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_keg')
       !
    END SUBROUTINE dyn_keg
