Index: /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5253)
+++ /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 5254)
@@ -803,4 +803,5 @@
 !-----------------------------------------------------------------------
    ln_dynadv_vec = .true.  !  vector form (T) or flux form (F)
+   nn_dynkeg     = 0       !  scheme for grad(KE): =0  C2  ; =1  Hollingsworth correction
    ln_dynadv_cen2= .false. !  flux form - 2nd order centered scheme
    ln_dynadv_ubs = .false. !  flux form - 3rd order UBS      scheme
Index: /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
===================================================================
--- /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5253)
+++ /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90	(revision 5254)
@@ -5,5 +5,6 @@
    !!==============================================================================
    !! History :  1.0  !  2006-11  (G. Madec)  Original code
-   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!            3.3  !  2010-10  (C. Ethe, G. Madec)  reorganisation of initialisation phase
+   !!            3.6  !  2015-05  (N. Ducousso, G. Madec)  add Hollingsworth scheme as an option 
    !!----------------------------------------------------------------------
 
@@ -17,4 +18,5 @@
    USE dynkeg          ! kinetic energy gradient          (dyn_keg      routine)
    USE dynzad          ! vertical advection               (dyn_zad      routine)
+   !
    USE in_out_manager  ! I/O manager
    USE lib_mpp         ! MPP library
@@ -25,7 +27,9 @@
 
    PUBLIC dyn_adv       ! routine called by step module
-   PUBLIC dyn_adv_init  ! routine called by opa module
+   PUBLIC dyn_adv_init  ! routine called by opa  module
  
+   !                                    !* namdyn_adv namelist *
    LOGICAL, PUBLIC ::   ln_dynadv_vec   !: vector form flag
+   INTEGER, PUBLIC ::   nn_dynkeg       !: scheme of kinetic energy gradient: =0 C2 ; =1 Hollingsworth
    LOGICAL, PUBLIC ::   ln_dynadv_cen2  !: flux form - 2nd order centered scheme flag
    LOGICAL, PUBLIC ::   ln_dynadv_ubs   !: flux form - 3rd order UBS scheme flag
@@ -38,5 +42,5 @@
 #  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)
@@ -63,16 +67,16 @@
       SELECT CASE ( nadv )                  ! compute advection trend and add it to general trend
       CASE ( 0 )     
-                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
-                      CALL dyn_zad     ( kt )    ! vector form : vertical advection
+                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad     ( kt )               ! vector form : vertical advection
       CASE ( 1 )     
-                      CALL dyn_keg     ( kt )    ! vector form : horizontal gradient of kinetic energy
-                      CALL dyn_zad_zts ( kt )    ! vector form : vertical advection with sub-timestepping
+                      CALL dyn_keg     ( kt, nn_dynkeg )    ! vector form : horizontal gradient of kinetic energy
+                      CALL dyn_zad_zts ( kt )               ! vector form : vertical advection with sub-timestepping
       CASE ( 2 ) 
-                      CALL dyn_adv_cen2( kt )    ! 2nd order centered scheme
+                      CALL dyn_adv_cen2( kt )               ! 2nd order centered scheme
       CASE ( 3 )   
-                      CALL dyn_adv_ubs ( kt )    ! 3rd order UBS      scheme
+                      CALL dyn_adv_ubs ( kt )               ! 3rd order UBS      scheme
       !
-      CASE (-1 )                                 ! esopa: test all possibility with control print
-                      CALL dyn_keg     ( kt )
+      CASE (-1 )                                            ! esopa: test all possibility with control print
+                      CALL dyn_keg     ( kt, nn_dynkeg )
                       CALL dyn_zad     ( kt )
                       CALL dyn_adv_cen2( kt )
@@ -92,10 +96,9 @@
       !!              momentum advection formulation & scheme and set nadv
       !!----------------------------------------------------------------------
-      INTEGER ::   ioptio
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      !!
-      NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
+      INTEGER ::   ioptio, ios   ! Local integer
+      !
+      NAMELIST/namdyn_adv/ ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts
       !!----------------------------------------------------------------------
-
+      !
       REWIND( numnam_ref )              ! Namelist namdyn_adv in reference namelist : Momentum advection scheme
       READ  ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901)
@@ -112,8 +115,9 @@
          WRITE(numout,*) '~~~~~~~~~~~'
          WRITE(numout,*) '       Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
-         WRITE(numout,*) '          Vector/flux form (T/F)             ln_dynadv_vec  = ', ln_dynadv_vec
-         WRITE(numout,*) '          2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
-         WRITE(numout,*) '          3rd order UBS advection scheme     ln_dynadv_ubs  = ', ln_dynadv_ubs
-         WRITE(numout,*) '      Sub timestepping of vertical advection ln_dynzad_zts  = ', ln_dynzad_zts
+         WRITE(numout,*) '          Vector/flux form (T/F)                           ln_dynadv_vec  = ', ln_dynadv_vec
+         WRITE(numout,*) '          = 0 standard scheme  ; =1 Hollingsworth scheme   nn_dynkeg      = ', nn_dynkeg
+         WRITE(numout,*) '          2nd order centred advection scheme               ln_dynadv_cen2 = ', ln_dynadv_cen2
+         WRITE(numout,*) '          3rd order UBS advection scheme                   ln_dynadv_ubs  = ', ln_dynadv_ubs
+         WRITE(numout,*) '          Sub timestepping of vertical advection           ln_dynzad_zts  = ', ln_dynzad_zts
       ENDIF
 
@@ -126,7 +130,7 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' )
       IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec )   &
-          CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
-      IF( ln_dynzad_zts .AND. ln_isfcav )   &
-          CALL ctl_stop( 'Sub timestepping of vertical advection does not work with ln_isfcav = .TRUE.' )
+         CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' )
+      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW )   &  
+         CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
 
       !                               ! Set nadv
@@ -139,6 +143,10 @@
       IF(lwp) THEN                    ! Print the choice
          WRITE(numout,*)
-         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used'
+         IF( nadv ==  0 )   WRITE(numout,*) '         vector form : keg + zad + vor is used' 
          IF( nadv ==  1 )   WRITE(numout,*) '         vector form : keg + zad_zts + vor is used'
+         IF( nadv ==  0 .OR. nadv ==  1 ) THEN
+            IF( nn_dynkeg == nkeg_C2  )   WRITE(numout,*) 'with Centered standard keg scheme'
+            IF( nn_dynkeg == nkeg_HW  )   WRITE(numout,*) 'with Hollingsworth keg scheme'
+         ENDIF
          IF( nadv ==  2 )   WRITE(numout,*) '         flux form   : 2nd order scheme is used'
          IF( nadv ==  3 )   WRITE(numout,*) '         flux form   : UBS       scheme is used'
Index: /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5253)
+++ /branches/2015/dev_r5224_CNRS_OPA_Hollingsworth/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 5254)
@@ -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 
    !!----------------------------------------------------------------------
    
@@ -18,4 +19,5 @@
    !
    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
@@ -28,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)
@@ -37,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_keg( kt )
+   SUBROUTINE dyn_keg( kt, kscheme )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_keg  ***
@@ -45,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).
@@ -54,6 +67,10 @@
       !! ** Action : - Update the (ua, va) with the hor. ke gradient trend
       !!             - 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
@@ -63,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)
@@ -104,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_dyn( ztrdu, ztrdv, jpdyn_keg, kt )
-         CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
+         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdu, ztrdv )
       ENDIF
       !
@@ -139,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
