Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90	(revision 12821)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/domain.F90	(revision 12822)
@@ -44,5 +44,7 @@
    USE lbclnk         ! ocean lateral boundary condition (or mpp link)
    USE lib_mpp        ! distributed memory computing library
-
+!!an45
+!   USE usrdef_nam, ONLY : ln_45machin
+   !
    IMPLICIT NONE
    PRIVATE
@@ -50,5 +52,5 @@
    PUBLIC   dom_init     ! called by nemogcm.F90
    PUBLIC   domain_cfg   ! called by nemogcm.F90
-
+#  include "do_loop_substitute.h90"
    !!-------------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -81,4 +83,6 @@
       INTEGER , DIMENSION(jpi,jpj) ::   ik_top , ik_bot       ! top and bottom ocean level
       REAL(wp), DIMENSION(jpi,jpj) ::   z1_hu_0, z1_hv_0
+      REAL(wp)::   zcoeff                                     ! local real
+
       !!----------------------------------------------------------------------
       !
@@ -154,4 +158,13 @@
          hf_0(:,:) = hf_0(:,:) + e3f_0(:,:,jk) * ssfmask(:,:)     ! CAUTION : only valid in SWE, not with bathymetry
       END DO
+      !
+!!anhf hf_0 = mean(ht_0*tmask) so hf = mimj( ht0 + ssht)
+! ne pas combiner avec an45 tout de suite
+!      DO_2D_10_10
+!         hf_0(ji,jj) = 0.25_wp * (   ht_0(ji,jj+1) * tmask(ji,jj+1,1) + ht_0(ji+1,jj+1) * tmask(ji+1,jj+1,1)   &
+!            &                      + ht_0(ji,jj  ) * tmask(ji,jj  ,1) + ht_0(ji+1,jj  ) * tmask(ji+1,jj  ,1)   )
+!      END_2D
+!      CALL lbc_lnk( 'domain', hf_0, 'F', 1. )      ! Lateral boundary conditions
+!!anhf
       !                                 ! Inverse of reference ocean thickness
       r1_ht_0(:,:) =  ssmask(:,:) / ( ht_0(:,:) + 1._wp -  ssmask(:,:) )
@@ -159,5 +172,22 @@
       r1_hv_0(:,:) = ssvmask(:,:) / ( hv_0(:,:) + 1._wp - ssvmask(:,:) )
       r1_hf_0(:,:) = ssfmask(:,:) / ( hf_0(:,:) + 1._wp - ssfmask(:,:) )
-      
+      !
+!!an45 Ligne de cote a 45deg : e1e2t *= ( mi(umask) + mj(vmask) ) /2
+!!                             idem pour e1e2f
+!      DO_2D_10_10
+!      zcoeff = 0.25_wp * (   umask(ji,jj+1,1) + umask(ji+1,jj+1,1)   &
+!         &                 + vmask(ji,jj  ,1) + vmask(ji+1,jj  ,1)   )
+!         IF ( zcoeff /= 0._wp )   THEN
+!               e1e2t(ji,jj) =    e1e2t(ji,jj) * zcoeff
+!            r1_e1e2t(ji,jj) = r1_e1e2t(ji,jj) / zcoeff
+!         ENDIF
+!      END_2D
+!      WRITE(numout,*) '   an45 half T cell e1e2t '
+!      zcoeff = 0.25_wp * (   umask(ji,jj+1,1) + umask(ji+1,jj+1,1)   &
+!         &                 + vmask(ji,jj  ,1) + vmask(ji+1,jj  ,1)   )
+!         IF ( zcoeff /= 0._wp )   THEN
+!               e1e2t(ji,jj) =    e1e2t(ji,jj) * zcoeff
+!            r1_e1e2t(ji,jj) = r1_e1e2t(ji,jj) / zcoeff
+!!an45
       !           !==  time varying part of coordinate system  ==!
       !
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 12822)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90	(revision 12822)
@@ -0,0 +1,153 @@
+MODULE dynadv
+   !!==============================================================================
+   !!                       ***  MODULE  dynadv  ***
+   !! Ocean active tracers:  advection scheme control
+   !!==============================================================================
+   !! History :  1.0  !  2006-11  (G. Madec)  Original code
+   !!            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 
+   !!            4.0  !  2017-07  (G. Madec)  add a linear dynamics option
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dyn_adv      : compute the momentum advection trend 
+   !!   dyn_adv_init : control the different options of advection scheme
+   !!----------------------------------------------------------------------
+   USE dom_oce         ! ocean space and time domain
+   USE dynadv_cen2     ! centred flux form advection      (dyn_adv_cen2 routine)
+   USE dynadv_ubs      ! UBS flux form advection          (dyn_adv_ubs  routine)
+   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
+   USE timing          ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dyn_adv       ! routine called by step module
+   PUBLIC dyn_adv_init  ! routine called by opa  module
+ 
+   !                                   !!* namdyn_adv namelist *
+   LOGICAL, PUBLIC ::   ln_dynadv_OFF   !: linear dynamics (no momentum advection)
+   LOGICAL, PUBLIC ::   ln_dynadv_vec   !: vector form
+   INTEGER, PUBLIC ::      nn_dynkeg       !: scheme of grad(KE): =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
+   
+   INTEGER, PUBLIC ::   n_dynadv   !: choice of the formulation and scheme for momentum advection
+   !                               !  associated indices:
+   INTEGER, PUBLIC, PARAMETER ::   np_LIN_dyn = 0   ! no advection: linear dynamics
+   INTEGER, PUBLIC, PARAMETER ::   np_VEC_c2  = 1   ! vector form : 2nd order centered scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_FLX_c2  = 2   ! flux   form : 2nd order centered scheme
+   INTEGER, PUBLIC, PARAMETER ::   np_FLX_ubs = 3   ! flux   form : 3rd order Upstream Biased Scheme
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_adv  ***
+      !!                
+      !! ** Purpose :   compute the ocean momentum advection trend.
+      !!
+      !! ** Method  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the advection term following n_dynadv
+      !!
+      !!      NB: in flux form advection (ln_dynadv_cen2 or ln_dynadv_ubs=T) 
+      !!      a metric term is add to the coriolis term while in vector form 
+      !!      it is the relative vorticity which is added to coriolis term
+      !!      (see dynvor module).
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'dyn_adv' )
+      !
+      SELECT CASE( n_dynadv )    !==  compute advection trend and add it to general trend  ==!
+      CASE( np_VEC_c2  )     
+         CALL dyn_keg     ( kt, nn_dynkeg,      Kmm, puu, pvv, Krhs )    ! vector form : horizontal gradient of kinetic energy
+         CALL dyn_zad     ( kt,                 Kmm, puu, pvv, Krhs )    ! vector form : vertical advection
+      CASE( np_FLX_c2  ) 
+         CALL dyn_adv_cen2( kt,                 Kmm, puu, pvv, Krhs )    ! 2nd order centered scheme
+      CASE( np_FLX_ubs )   
+         CALL dyn_adv_ubs ( kt,            Kbb, Kmm, puu, pvv, Krhs )    ! 3rd order UBS      scheme (UP3)
+      END SELECT
+      !
+      IF( ln_timing )   CALL timing_stop( 'dyn_adv' )
+      !
+   END SUBROUTINE dyn_adv
+
+   
+   SUBROUTINE dyn_adv_init
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_adv_init  ***
+      !!                
+      !! ** Purpose :   Control the consistency between namelist options for 
+      !!              momentum advection formulation & scheme and set n_dynadv
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio, ios   ! Local integer
+      !
+      NAMELIST/namdyn_adv/ ln_dynadv_OFF, ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2, ln_dynadv_ubs
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dyn_adv_init : choice/control of the momentum advection scheme'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
+      READ  ( numnam_ref, namdyn_adv, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_adv in reference namelist' )
+      READ  ( numnam_cfg, namdyn_adv, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_adv in configuration namelist' )
+      IF(lwm) WRITE ( numond, namdyn_adv )
+
+      IF(lwp) THEN                    ! Namelist print
+         WRITE(numout,*) '   Namelist namdyn_adv : chose a advection formulation & scheme for momentum'
+         WRITE(numout,*) '      linear dynamics : no momentum advection          ln_dynadv_OFF  = ', ln_dynadv_OFF
+         WRITE(numout,*) '      Vector form: 2nd order centered scheme           ln_dynadv_vec  = ', ln_dynadv_vec
+!!an45
+         WRITE(numout,*) '         with Hollingsworth scheme (=1) or not (=0,2)      nn_dynkeg  = ', nn_dynkeg
+!!an45
+         WRITE(numout,*) '      flux form: 2nd order centred scheme              ln_dynadv_cen2 = ', ln_dynadv_cen2
+         WRITE(numout,*) '                 3rd order UBS scheme                  ln_dynadv_ubs  = ', ln_dynadv_ubs
+      ENDIF
+
+      ioptio = 0                      ! parameter control and set n_dynadv
+      IF( ln_dynadv_OFF  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_LIN_dyn   ;   ENDIF
+      IF( ln_dynadv_vec  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_VEC_c2    ;   ENDIF
+      IF( ln_dynadv_cen2 ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_FLX_c2    ;   ENDIF
+      IF( ln_dynadv_ubs  ) THEN   ;   ioptio = ioptio + 1   ;   n_dynadv = np_FLX_ubs   ;   ENDIF
+
+      IF( ioptio /= 1 )   CALL ctl_stop( 'choose ONE and only ONE advection scheme' )
+!!an45
+      IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_C2_wpo .AND. nn_dynkeg /= nkeg_HW )   & 
+         &   CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
+!!an45
+
+      IF(lwp) THEN                    ! Print the choice
+         WRITE(numout,*)
+         SELECT CASE( n_dynadv )
+         CASE( np_LIN_dyn )   ;   WRITE(numout,*) '   ==>>>   linear dynamics : no momentum advection used'
+         CASE( np_VEC_c2  )   ;   WRITE(numout,*) '   ==>>>   vector form : keg + zad + vor is used' 
+            IF( nn_dynkeg == nkeg_C2     )   WRITE(numout,*) '              with Centered standard keg scheme'
+!!an45
+            IF( nn_dynkeg == nkeg_C2_wpo )   WRITE(numout,*) '              with Centered standard keg scheme (wet point only)'
+!!an45
+            IF( nn_dynkeg == nkeg_HW     )   WRITE(numout,*) '              with Hollingsworth keg scheme'
+         CASE( np_FLX_c2  )   ;   WRITE(numout,*) '   ==>>>   flux form   : 2nd order scheme is used'
+         CASE( np_FLX_ubs )   ;   WRITE(numout,*) '   ==>>>   flux form   : UBS       scheme is used'
+         END SELECT
+      ENDIF
+      !
+   END SUBROUTINE dyn_adv_init
+
+  !!======================================================================
+END MODULE dynadv
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 12822)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90	(revision 12822)
@@ -0,0 +1,159 @@
+MODULE dynkeg
+   !!======================================================================
+   !!                       ***  MODULE  dynkeg  ***
+   !! Ocean dynamics:  kinetic energy gradient trend
+   !!======================================================================
+   !! 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 
+   !!----------------------------------------------------------------------
+   
+   !!----------------------------------------------------------------------
+   !!   dyn_keg      : update the momentum trend with the horizontal tke
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   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
+   USE timing          ! Timing
+   USE bdy_oce         ! ocean open boundary conditions
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   dyn_keg    ! routine called by step module
+   
+   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2     = 0   !: 2nd order centered scheme (standard scheme)
+!!an45
+   INTEGER, PARAMETER, PUBLIC  ::   nkeg_C2_wpo = 2   !: 2nd order centered scheme (wet point only : coastline at 45 degrees)
+   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 "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: dynkeg.F90 12377 2020-02-12 14:39:06Z acc $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
+      !!----------------------------------------------------------------------
+      !!                  ***  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  : * 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 + ((u(j+1)+u(j-1))/2)^2  )
+      !!                    + mj-1(  2 * vn^2 + ((v(i+1)+v(i-1))/2)^2  ) ]
+      !!      
+      !!      Take its horizontal gradient and add it to the general momentum
+      !!      trend.
+      !!         u(rhs) = u(rhs) - 1/e1u di[ zhke ]
+      !!         v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
+      !!
+      !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) 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 )  ::  kscheme          ! =0/1/2   type of KEG scheme 
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs        ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
+      !
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::   zu, zv                   ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('dyn_keg')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
+         IF(lwp) WRITE(numout,*) '~~~~~~~'
+      ENDIF
+
+      IF( l_trddyn ) THEN           ! Save the input trends
+         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
+         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
+      ENDIF
+      
+      zhke(:,:,jpk) = 0._wp
+
+      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
+!!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2)
+      CASE ( nkeg_C2_wpo )                          !--  Standard scheme  --!
+         DO_3D_01_01( 1, jpkm1 )
+            zu =  (   puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
+               &    + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)   ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) )
+            zv =  (   pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   &
+               &    + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm)   ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) )
+            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+         END_3D
+!!an45         
+      !
+      CASE ( nkeg_C2 )                          !--  Standard scheme  --!
+         DO_3D_01_01( 1, jpkm1 )
+            zu =    puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)   &
+               &  + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm)
+            zv =    pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)   &
+               &  + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm)
+            zhke(ji,jj,jk) = 0.25_wp * ( zv + zu )
+         END_3D
+      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
+         DO_3D_00_00( 1, jpkm1 )
+            zu = 8._wp * ( puu(ji-1,jj  ,jk,Kmm) * puu(ji-1,jj  ,jk,Kmm)    &
+               &         + puu(ji  ,jj  ,jk,Kmm) * puu(ji  ,jj  ,jk,Kmm) )  &
+               &   +     ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) )   &
+               &   +     ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) ) * ( puu(ji  ,jj-1,jk,Kmm) + puu(ji  ,jj+1,jk,Kmm) )
+               !
+            zv = 8._wp * ( pvv(ji  ,jj-1,jk,Kmm) * pvv(ji  ,jj-1,jk,Kmm)    &
+               &         + pvv(ji  ,jj  ,jk,Kmm) * pvv(ji  ,jj  ,jk,Kmm) )  &
+               &  +      ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) )   &
+               &  +      ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) ) * ( pvv(ji-1,jj  ,jk,Kmm) + pvv(ji+1,jj  ,jk,Kmm) )
+            zhke(ji,jj,jk) = r1_48 * ( zv + zu )
+         END_3D
+         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
+         !
+      END SELECT 
+      !
+      DO_3D_00_00( 1, jpkm1 )
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj  ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji  ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)
+      END_3D
+      !
+      IF( l_trddyn ) THEN                 ! save the Kinetic Energy trends for diagnostic
+         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt, Kmm )
+         DEALLOCATE( ztrdu , ztrdv )
+      ENDIF
+      !
+      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg  - Ua: ', mask1=umask,   &
+         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      !
+      IF( ln_timing )   CALL timing_stop('dyn_keg')
+      !
+   END SUBROUTINE dyn_keg
+
+   !!======================================================================
+END MODULE dynkeg
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90	(revision 12822)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90	(revision 12822)
@@ -0,0 +1,238 @@
+MODULE dynldf_lap_blp
+   !!======================================================================
+   !!                   ***  MODULE  dynldf_lap_blp  ***
+   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian)
+   !!======================================================================
+   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
+   !!           4.0  ! 2020-04  (A. Nasser, G. Madec)  Add symmetric mixing tensor 
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator
+   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
+   USE ldfslp         ! iso-neutral slopes 
+   USE zdf_oce        ! ocean vertical physics
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
+   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
+!!anSYM  
+   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_rot  = 1         ! div-rot   laplacian 
+   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_sym  = 2         ! symmetric laplacian (Griffies&Hallberg 2000)
+   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian)
+   
+   INTEGER, PUBLIC, PARAMETER ::   ln_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist)
+!!anSYM
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE dyn_ldf_lap  ***
+      !!                       
+      !! ** Purpose :   Compute the before horizontal momentum diffusive 
+      !!      trend and add it to the general trend of momentum equation.
+      !!
+      !! ** Method  :   The Laplacian operator apply on horizontal velocity is 
+      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
+      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
+      !!
+      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
+      !!
+      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ 
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
+      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zsign        ! local scalars
+      REAL(wp) ::   zua, zva     ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms
+      !!----------------------------------------------------------------------
+      !
+!!anSYM TO BE ADDED : reading of laplacian operator (ln_dynldf_lap_typ -> to be written nn_) shall be added in dyn_ldf_init 
+!!                 as the writing   
+!!                 and an integer as np_dynldf_lap for instance taken as argument by dyn_ldf_lap call in dyn_ldf
+      IF( kt == nit000 .AND. lwp ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '                  ln_dynldf_lap_typ = ', ln_dynldf_lap_typ
+         SELECT CASE( ln_dynldf_lap_typ )             ! print the choice of operator
+         CASE( np_dynldf_lap_rot )   ;   WRITE(numout,*) '   ==>>>   div-rot   laplacian'
+         CASE( np_dynldf_lap_sym )   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (covariant form)'
+         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (simple form)'
+         END SELECT
+      ENDIF
+      !
+      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
+      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
+      ENDIF
+      !
+      SELECT CASE( ln_dynldf_lap_typ )  
+         !              
+         CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==!
+            !   
+            DO jk = 1, jpkm1                                 ! Horizontal slab
+               !                                             
+               DO_2D_01_01
+               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
+!!gm open question here : e3f  at before or now ?    probably now... 
+!!gm note that ahmf has already been multiplied by fmask
+                  zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
+                     &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+                     &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
+                  !                                      ! ahm * div        (computed from 2 to jpi/jpj)
+!!gm note that ahmt has already been multiplied by tmask
+                  zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         &
+                     &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
+                     &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
+               END_2D
+               !
+               DO_2D_00_00
+                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                             &
+                     &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
+                     &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                       )
+                     !
+                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                             &
+                     &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
+                     &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                       )
+               END_2D
+               !
+            END DO                                           !   End of slab
+            !
+         CASE ( np_dynldf_lap_sym )       !==  Symmetric form  ==!   (Griffies&Hallberg 2000)
+            !
+            DO jk = 1, jpkm1                                 ! Horizontal slab
+               !
+               DO_2D_01_01
+                  !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
+                  zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              &
+                     &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             &
+                     &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   &
+                     &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             &
+                     &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   ) 
+                  !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
+                  zten(ji,jj)    = ahmt(ji,jj,jk)                                                       &
+                     &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         &
+                     &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   &
+                     &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         &
+                     &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )   
+               END_2D
+               !
+               DO_2D_00_00
+                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               &
+                     &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       &
+                     &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &                                                    
+                     &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
+                     &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )   
+                  !
+                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               &
+                     &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
+                     &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     &
+                     &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       &
+                     &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) )
+                   !
+               END_2D
+               !
+            END DO                                           !   End of slab
+            !
+         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (naive way)
+            !
+            DO jk = 1, jpkm1                                 ! Horizontal slab
+               !
+               DO_2D_01_01
+                  !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
+                  zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                           &
+                     &     * (   r1_e2f(ji-1,jj-1) * ( pu(ji-1,jj  ,jk) - pu(ji-1,jj-1,jk)  )   &
+                     &         + r1_e1f(ji-1,jj-1) * ( pv(ji  ,jj-1,jk) - pv(ji-1,jj-1,jk)  )   ) 
+                  !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
+                  zten(ji,jj)    = ahmt(ji,jj,jk)                                       &
+                     &     * (   r1_e1t(ji,jj) * ( pu(ji,jj,jk) - pu(ji-1,jj  ,jk)  )   &
+                     &         - r1_e2t(ji,jj) * ( pv(ji,jj,jk) - pv(ji  ,jj-1,jk)  )   )   
+               END_2D
+               !
+               DO_2D_00_00
+                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
+                     &    * (   zten(ji+1,jj  ) * e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)              &
+                     &        - zten(ji  ,jj  ) * e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm)              &
+                     &        + zshe(ji  ,jj  ) * e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                  &
+                     &        - zshe(ji  ,jj-1) * e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)                  )    
+                  !
+                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
+                     &    * (   zshe(ji  ,jj  ) * e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                  &
+                     &        - zshe(ji-1,jj  ) * e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)                  &
+                     &        - zten(ji  ,jj+1) * e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)              &
+                     &        + zten(ji  ,jj  ) * e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm)              )
+                   !
+               END_2D
+               !
+            END DO                                           !   End of slab
+            !  
+         CASE DEFAULT                                     ! error
+            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for ln_dynldf_lap_typ'  )
+         END SELECT
+         !
+      !
+   END SUBROUTINE dyn_ldf_lap
+
+
+   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dyn_ldf_blp  ***
+      !!                    
+      !! ** Purpose :   Compute the before lateral momentum viscous trend 
+      !!              and add it to the general trend of momentum equation.
+      !!
+      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
+      !!      operator applied to before field (forward in time).
+      !!      It is computed by two successive calls to dyn_ldf_lap routine
+      !!
+      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
+      !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
+      !!----------------------------------------------------------------------
+      !
+      IF( kt == nit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+      !
+      zulap(:,:,:) = 0._wp
+      zvlap(:,:,:) = 0._wp
+      !
+      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
+      !
+      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
+      !
+      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
+      !
+   END SUBROUTINE dyn_ldf_blp
+
+   !!======================================================================
+END MODULE dynldf_lap_blp
Index: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90
===================================================================
--- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90	(revision 12822)
+++ NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/ldfdyn.F90	(revision 12822)
@@ -0,0 +1,525 @@
+MODULE ldfdyn
+   !!======================================================================
+   !!                       ***  MODULE  ldfdyn  ***
+   !! Ocean physics:  lateral viscosity coefficient 
+   !!=====================================================================
+   !! History :  OPA  ! 1997-07  (G. Madec)  multi dimensional coefficients
+   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
+   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
+   !!                 !                                  add velocity dependent coefficient and optional read in file
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   ldf_dyn_init  : initialization, namelist read, and parameters control
+   !!   ldf_dyn       : update lateral eddy viscosity coefficients at each time step 
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers   
+   USE dom_oce         ! ocean space and time domain 
+   USE phycst          ! physical constants
+   USE ldfslp          ! lateral diffusion: slopes of mixing orientation
+   USE ldfc1d_c2d      ! lateral diffusion: 1D and 2D cases
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O module for ehanced bottom friction file
+   USE timing          ! Timing
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ldf_dyn_init   ! called by nemogcm.F90
+   PUBLIC   ldf_dyn        ! called by step.F90
+
+   !                                    !!* Namelist namdyn_ldf : lateral mixing on momentum *
+   LOGICAL , PUBLIC ::   ln_dynldf_OFF   !: No operator (i.e. no explicit diffusion)
+   LOGICAL , PUBLIC ::   ln_dynldf_lap   !: laplacian operator
+   LOGICAL , PUBLIC ::   ln_dynldf_blp   !: bilaplacian operator
+   LOGICAL , PUBLIC ::   ln_dynldf_lev   !: iso-level direction
+   LOGICAL , PUBLIC ::   ln_dynldf_hor   !: horizontal (geopotential) direction
+!  LOGICAL , PUBLIC ::   ln_dynldf_iso   !: iso-neutral direction                        (see ldfslp)
+   INTEGER , PUBLIC ::   nn_ahm_ijk_t    !: choice of time & space variations of the lateral eddy viscosity coef.
+   !                                        !  time invariant coefficients:  aht = 1/2  Ud*Ld   (lap case) 
+   !                                           !                             bht = 1/12 Ud*Ld^3 (blp case)
+   REAL(wp), PUBLIC ::   rn_Uv                 !: lateral viscous velocity  [m/s]
+   REAL(wp), PUBLIC ::   rn_Lv                 !: lateral viscous length    [m]
+   !                                        ! Smagorinsky viscosity  (nn_ahm_ijk_t = 32) 
+   REAL(wp), PUBLIC ::   rn_csmc               !: Smagorinsky constant of proportionality 
+   REAL(wp), PUBLIC ::   rn_minfac             !: Multiplicative factor of theorectical minimum Smagorinsky viscosity
+   REAL(wp), PUBLIC ::   rn_maxfac             !: Multiplicative factor of theorectical maximum Smagorinsky viscosity
+   !                                        ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T)
+   REAL(wp), PUBLIC ::   rn_ahm_b              !: lateral laplacian background eddy viscosity  [m2/s]
+
+   !                                    !!* Parameter to control the type of lateral viscous operator
+   INTEGER, PARAMETER, PUBLIC ::   np_ERROR  =-10                       !: error in setting the operator
+   INTEGER, PARAMETER, PUBLIC ::   np_no_ldf = 00                       !: without operator (i.e. no lateral viscous trend)
+   !                          !!      laplacian     !    bilaplacian    !
+   INTEGER, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator
+   !
+   INTEGER           , PUBLIC ::   nldf_dyn         !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
+   LOGICAL           , PUBLIC ::   l_ldfdyn_time    !: flag for time variation of the lateral eddy viscosity coef.
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ahmt, ahmf   !: eddy viscosity coef. at T- and F-points [m2/s or m4/s]
+   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dtensq       !: horizontal tension squared         (Smagorinsky only)
+   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dshesq       !: horizontal shearing strain squared (Smagorinsky only)
+   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   esqt, esqf   !: Square of the local gridscale (e1e2/(e1+e2))**2           
+
+   REAL(wp) ::   r1_2    = 0.5_wp            ! =1/2
+   REAL(wp) ::   r1_4    = 0.25_wp           ! =1/4
+   REAL(wp) ::   r1_8    = 0.125_wp          ! =1/8
+   REAL(wp) ::   r1_12   = 1._wp / 12._wp    ! =1/12
+   REAL(wp) ::   r1_288  = 1._wp / 288._wp   ! =1/( 12^2 * 2 )
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id$ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ldf_dyn_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_dyn_init  ***
+      !!                   
+      !! ** Purpose :   set the horizontal ocean dynamics physics
+      !!
+      !! ** Method  :   the eddy viscosity coef. specification depends on:
+      !!              - the operator:
+      !!             ln_dynldf_lap = T     laplacian operator
+      !!             ln_dynldf_blp = T   bilaplacian operator
+      !!              - the parameter nn_ahm_ijk_t:
+      !!    nn_ahm_ijk_t  =  0 => = constant
+      !!                  = 10 => = F(z) :     = constant with a reduction of 1/4 with depth 
+      !!                  =-20 => = F(i,j)     = shape read in 'eddy_viscosity.nc' file
+      !!                  = 20    = F(i,j)     = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case)
+      !!                  =-30 => = F(i,j,k)   = shape read in 'eddy_viscosity.nc'  file
+      !!                  = 30    = F(i,j,k)   = 2D (case 20) + decrease with depth (case 10)
+      !!                  = 31    = F(i,j,k,t) = F(local velocity) (  |u|e  /12   laplacian operator
+      !!                                                           or |u|e^3/12 bilaplacian operator )
+      !!                  = 32    = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)  
+      !!                                                           (   L^2|D|      laplacian operator
+      !!                                                           or  L^4|D|/8  bilaplacian operator )
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk                     ! dummy loop indices
+      INTEGER  ::   ioptio, ierr, inum, ios, inn   ! local integer
+      REAL(wp) ::   zah0, zah_max, zUfac           ! local scalar
+      REAL(wp) ::   zsum                           ! local scalar
+      CHARACTER(len=5) ::   cl_Units               ! units (m2/s or m4/s)
+      !!
+      NAMELIST/namdyn_ldf/ ln_dynldf_OFF, ln_dynldf_lap, ln_dynldf_blp,   &   ! type of operator
+         &                 ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso,   &   ! acting direction of the operator
+         &                 nn_ahm_ijk_t , rn_Uv    , rn_Lv,   rn_ahm_b,   &   ! lateral eddy coefficient
+         &                 rn_csmc      , rn_minfac    , rn_maxfac            ! Smagorinsky settings
+      !!----------------------------------------------------------------------
+      !
+      READ  ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist' )
+
+      READ  ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist' )
+      IF(lwm) WRITE ( numond, namdyn_ldf )
+
+      IF(lwp) THEN                      ! Parameter print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ldf_dyn : lateral momentum physics'
+         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) '   Namelist namdyn_ldf : set lateral mixing parameters'
+         !
+         WRITE(numout,*) '      type :'
+         WRITE(numout,*) '         no explicit diffusion                ln_dynldf_OFF = ', ln_dynldf_OFF
+         WRITE(numout,*) '         laplacian operator                   ln_dynldf_lap = ', ln_dynldf_lap
+         WRITE(numout,*) '         bilaplacian operator                 ln_dynldf_blp = ', ln_dynldf_blp
+         !
+         WRITE(numout,*) '      direction of action :'
+         WRITE(numout,*) '         iso-level                            ln_dynldf_lev = ', ln_dynldf_lev
+         WRITE(numout,*) '         horizontal (geopotential)            ln_dynldf_hor = ', ln_dynldf_hor
+         WRITE(numout,*) '         iso-neutral                          ln_dynldf_iso = ', ln_dynldf_iso
+         !
+         WRITE(numout,*) '      coefficients :'
+         WRITE(numout,*) '         type of time-space variation         nn_ahm_ijk_t  = ', nn_ahm_ijk_t
+         WRITE(numout,*) '         lateral viscous velocity  (if cst)      rn_Uv      = ', rn_Uv, ' m/s'
+         WRITE(numout,*) '         lateral viscous length    (if cst)      rn_Lv      = ', rn_Lv, ' m'
+         WRITE(numout,*) '         background viscosity (iso-lap case)     rn_ahm_b   = ', rn_ahm_b, ' m2/s'
+         !
+         WRITE(numout,*) '      Smagorinsky settings (nn_ahm_ijk_t  = 32) :'
+         WRITE(numout,*) '         Smagorinsky coefficient              rn_csmc       = ', rn_csmc
+         WRITE(numout,*) '         factor multiplier for eddy visc.'
+         WRITE(numout,*) '            lower limit (default 1.0)         rn_minfac    = ', rn_minfac
+         WRITE(numout,*) '            upper limit (default 1.0)         rn_maxfac    = ', rn_maxfac
+      ENDIF
+
+      !
+      !           !==  type of lateral operator used  ==!   (set nldf_dyn)
+      !           !=====================================!
+      !
+      nldf_dyn = np_ERROR
+      ioptio = 0
+      IF( ln_dynldf_OFF ) THEN   ;   nldf_dyn = np_no_ldf   ;   ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_dynldf_lap ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
+      IF( ln_dynldf_blp ) THEN   ;                              ioptio = ioptio + 1   ;   ENDIF
+      IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' )
+      !
+      IF(.NOT.ln_dynldf_OFF ) THEN     !==  direction ==>> type of operator  ==!
+         ioptio = 0
+         IF( ln_dynldf_lev )   ioptio = ioptio + 1
+         IF( ln_dynldf_hor )   ioptio = ioptio + 1
+         IF( ln_dynldf_iso )   ioptio = ioptio + 1
+         IF( ioptio /= 1   )   CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' )
+         !
+         !                             ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals
+         ierr = 0
+         IF( ln_dynldf_lap ) THEN         ! laplacian operator
+            IF( ln_zco ) THEN                   ! z-coordinate
+               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
+            ENDIF
+            IF( ln_zps ) THEN                   ! z-coordinate with partial step
+               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level              (no rotation)
+               IF ( ln_dynldf_hor )   nldf_dyn = np_lap     ! iso-level              (no rotation)
+               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
+            ENDIF
+            IF( ln_sco ) THEN                   ! s-coordinate
+               IF ( ln_dynldf_lev )   nldf_dyn = np_lap     ! iso-level = horizontal (no rotation)
+               IF ( ln_dynldf_hor )   nldf_dyn = np_lap_i   ! horizontal             (   rotation)
+               IF ( ln_dynldf_iso )   nldf_dyn = np_lap_i   ! iso-neutral            (   rotation)
+            ENDIF
+         ENDIF
+         !
+         IF( ln_dynldf_blp ) THEN         ! bilaplacian operator
+            IF( ln_zco ) THEN                   ! z-coordinate
+               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
+               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level = horizontal (no rotation)
+               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
+            ENDIF
+            IF( ln_zps ) THEN                   ! z-coordinate with partial step
+               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
+               IF( ln_dynldf_hor )   nldf_dyn = np_blp   ! iso-level              (no rotation)
+               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
+            ENDIF
+            IF( ln_sco ) THEN                   ! s-coordinate
+               IF( ln_dynldf_lev )   nldf_dyn = np_blp   ! iso-level              (no rotation)
+               IF( ln_dynldf_hor )   ierr = 2            ! horizontal             (   rotation)
+               IF( ln_dynldf_iso )   ierr = 2            ! iso-neutral            (   rotation)
+            ENDIF
+         ENDIF
+         !
+         IF( ierr == 2 )   CALL ctl_stop( 'rotated bi-laplacian operator does not exist' )
+         !
+         IF( nldf_dyn == np_lap_i )   l_ldfslp = .TRUE.  ! rotation require the computation of the slopes
+         !
+      ENDIF
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         SELECT CASE( nldf_dyn )
+         CASE( np_no_ldf )   ;   WRITE(numout,*) '   ==>>>   NO lateral viscosity'
+         CASE( np_lap    )   ;   WRITE(numout,*) '   ==>>>   iso-level laplacian operator'
+         CASE( np_lap_i  )   ;   WRITE(numout,*) '   ==>>>   rotated laplacian operator with iso-level background'
+         CASE( np_blp    )   ;   WRITE(numout,*) '   ==>>>   iso-level bi-laplacian operator'
+         END SELECT
+         WRITE(numout,*)
+      ENDIF
+      
+      !
+      !           !==  Space/time variation of eddy coefficients  ==!
+      !           !=================================================!
+      !
+      l_ldfdyn_time = .FALSE.                ! no time variation except in case defined below
+      !
+      IF( ln_dynldf_OFF ) THEN
+         IF(lwp) WRITE(numout,*) '   ==>>>   No viscous operator selected. ahmt and ahmf are not allocated'
+         RETURN
+         !
+      ELSE                                   !==  a lateral diffusion operator is used  ==!
+         !
+         !                                         ! allocate the ahm arrays
+         ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays')
+         !
+         ahmt(:,:,:) = 0._wp                       ! init to 0 needed 
+         ahmf(:,:,:) = 0._wp
+         !
+         !                                         ! value of lap/blp eddy mixing coef.
+         IF(     ln_dynldf_lap ) THEN   ;   zUfac = r1_2 *rn_Uv   ;   inn = 1   ;   cl_Units = ' m2/s'   !   laplacian
+         ELSEIF( ln_dynldf_blp ) THEN   ;   zUfac = r1_12*rn_Uv   ;   inn = 3   ;   cl_Units = ' m4/s'   ! bilaplacian
+         ENDIF
+         zah0    = zUfac *    rn_Lv**inn              ! mixing coefficient
+         zah_max = zUfac * (ra*rad)**inn              ! maximum reachable coefficient (value at the Equator)
+         !
+         SELECT CASE(  nn_ahm_ijk_t  )             !* Specification of space-time variations of ahmt, ahmf
+         !
+         CASE(   0  )      !==  constant  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity. = constant = ', zah0, cl_Units
+            ahmt(:,:,1:jpkm1) = zah0
+            ahmf(:,:,1:jpkm1) = zah0
+            !
+         CASE(  10  )      !==  fixed profile  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( depth )'
+            IF(lwp) WRITE(numout,*) '           surface viscous coef. = constant = ', zah0, cl_Units
+            ahmt(:,:,1) = zah0                        ! constant surface value
+            ahmf(:,:,1) = zah0
+            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )
+            !
+         CASE ( -20 )      !== fixed horizontal shape read in file  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j) read in eddy_viscosity.nc file'
+            CALL iom_open( 'eddy_viscosity_2D.nc', inum )
+            CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) )
+            CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) )
+            CALL iom_close( inum )
+            DO jk = 2, jpkm1
+               ahmt(:,:,jk) = ahmt(:,:,1)
+               ahmf(:,:,jk) = ahmf(:,:,1)
+            END DO
+            !
+         CASE(  20  )      !== fixed horizontal shape  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)'
+            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Lv = Max(e1,e2)'
+            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
+            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
+            !
+         CASE( -30  )      !== fixed 3D shape read in file  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file'
+            CALL iom_open( 'eddy_viscosity_3D.nc', inum )
+            CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt )
+            CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf )
+            CALL iom_close( inum )
+            !
+         CASE(  30  )       !==  fixed 3D shape  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth )'
+            IF(lwp) WRITE(numout,*) '           using a fixed viscous velocity = ', rn_Uv  ,' m/s   and   Ld = Max(e1,e2)'
+            IF(lwp) WRITE(numout,*) '           maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, '  for e1=1°)'
+            CALL ldf_c2d( 'DYN', zUfac      , inn        , ahmt, ahmf )         ! surface value proportional to scale factor^inn
+            CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )  ! reduction with depth
+            !
+         CASE(  31  )       !==  time varying 3D field  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
+            IF(lwp) WRITE(numout,*) '           proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)'
+            !
+            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
+            !
+         CASE(  32  )       !==  time varying 3D field  ==!
+            IF(lwp) WRITE(numout,*) '   ==>>>   eddy viscosity = F( latitude, longitude, depth , time )'
+            IF(lwp) WRITE(numout,*) '           proportional to the local deformation rate and gridscale (Smagorinsky)'
+            !
+            l_ldfdyn_time = .TRUE.     ! will be calculated by call to ldf_dyn routine in step.F90
+            !
+            !                          ! allocate arrays used in ldf_dyn. 
+            ALLOCATE( dtensq(jpi,jpj,jpk) , dshesq(jpi,jpj,jpk) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr )
+            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays')
+            !
+            DO_2D_11_11
+               esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
+               esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
+            END_2D
+            !
+         CASE DEFAULT
+            CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm')
+         END SELECT
+         !
+         IF( .NOT.l_ldfdyn_time ) THEN             !* No time variation 
+            IF(     ln_dynldf_lap ) THEN                 !   laplacian operator (mask only)
+               ahmt(:,:,1:jpkm1) =       ahmt(:,:,1:jpkm1)   * tmask(:,:,1:jpkm1)
+               WRITE(numout,*) ' ahmt tmask '
+!! mask tension at the coast (equivalent of ghostpoints at T)
+!              DO jk = 1, jpk
+!                 DO jj = 1, jpjm1
+!                    DO ji = 1, jpim1      ! NO vector opt.
+!                       ! si sum(fmask)==3 = mouillé (on touche pas)
+!                       ! si sum = 2 alors on met a 0 zsum = fmask + fmask + fmask,.. et si zsum < 2 -> 0 sinon = 1
+!                       zsum =   fmask(ji,jj  ,jk) + fmask(ji+1,jj  ,jk)   &
+!                          &   + fmask(ji,jj+1,jk) + fmask(ji+1,jj+1,jk)
+!                       IF ( zsum < 2._wp )   ahmt(ji,jj,jk) = 0
+!                       !
+!                       !ahmt(ji,jj,jk) = ahmt(ji,jj,jk) * fmask(ji,jj  ,jk) * fmask(ji+1,jj  ,jk)   &
+!                       !   &                            * fmask(ji,jj+1,jk) * fmask(ji+1,jj+1,jk)
+!                    END DO
+!                 END DO
+!              END DO
+!              ahmt(jpi,:,1:jpkm1) =  0._wp
+!              ahmt(:,jpj,1:jpkm1) =  0._wp
+!              WRITE(numout,*) '   an45 ahmt x0'
+
+               ahmf(:,:,1:jpkm1) =       ahmf(:,:,1:jpkm1)   * fmask(:,:,1:jpkm1)
+               WRITE(numout,*) ' ahmf fmask '
+!!an apply no slip at the coast (ssfmask = 1 within the domain and at the coast contrary to fmask in free slip)
+!               DO jk = 1, jpkm1
+!                  ahmf(:,:,jk) =    ahmf(:,:,jk) * ( 2._wp * ssfmask(:,:) - fmask(:,:,jk) )
+!               END DO
+!               WRITE(numout,*) '   an45 ahmf x2'
+
+            ELSEIF( ln_dynldf_blp ) THEN                 ! bilaplacian operator (square root + mask)
+               ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1)
+               ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1)
+            ENDIF
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE ldf_dyn_init
+
+
+   SUBROUTINE ldf_dyn( kt, Kbb )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_dyn  ***
+      !! 
+      !! ** Purpose :   update at kt the momentum lateral mixing coeff. (ahmt and ahmf)
+      !!
+      !! ** Method  :   time varying eddy viscosity coefficients:
+      !!
+      !!    nn_ahm_ijk_t = 31    ahmt, ahmf = F(i,j,k,t) = F(local velocity) 
+      !!                         ( |u|e /12  or  |u|e^3/12 for laplacian or bilaplacian operator )
+      !!
+      !!    nn_ahm_ijk_t = 32    ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky)  
+      !!                         ( L^2|D|    or  L^4|D|/8  for laplacian or bilaplacian operator )
+      !!
+      !! ** note    :    in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian
+      !! ** action  :    ahmt, ahmf   updated at each time step
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! time step index
+      INTEGER, INTENT(in) ::   Kbb  ! ocean time level indices
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zemax   ! local scalar (option 31)
+      REAL(wp) ::   zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb     ! local scalar (option 32)
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('ldf_dyn')
+      !
+      SELECT CASE(  nn_ahm_ijk_t  )       !== Eddy vicosity coefficients ==!
+      !
+      CASE(  31  )       !==  time varying 3D field  ==!   = F( local velocity )
+         !
+         IF( ln_dynldf_lap   ) THEN        ! laplacian operator : |u| e /12 = |u/144| e
+            DO jk = 1, jpkm1
+               DO_2D_00_00
+                  zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
+                  zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
+                  zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
+                  ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
+               END_2D
+               DO_2D_10_10
+                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, Kbb) * uu(ji  ,jj+1,jk, Kbb) + vv(ji+1,jj  ,jk, Kbb) * vv(ji+1,jj  ,jk, Kbb)
+                  zu2pv2_ij    = uu(ji  ,jj  ,jk, Kbb) * uu(ji  ,jj  ,jk, Kbb) + vv(ji  ,jj  ,jk, Kbb) * vv(ji  ,jj  ,jk, Kbb)
+                  zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
+                  ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk)      ! 288= 12*12 * 2
+               END_2D
+            END DO
+         ELSEIF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e
+            DO jk = 1, jpkm1
+               DO_2D_00_00
+                  zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
+                  zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
+                  zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
+                  ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax  ) * zemax * tmask(ji,jj,jk)
+               END_2D
+               DO_2D_10_10
+                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, Kbb) * uu(ji  ,jj+1,jk, Kbb) + vv(ji+1,jj  ,jk, Kbb) * vv(ji+1,jj  ,jk, Kbb)
+                  zu2pv2_ij    = uu(ji  ,jj  ,jk, Kbb) * uu(ji  ,jj  ,jk, Kbb) + vv(ji  ,jj  ,jk, Kbb) * vv(ji  ,jj  ,jk, Kbb)
+                  zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
+                  ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax  ) * zemax * fmask(ji,jj,jk)
+               END_2D
+            END DO
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1.,  ahmf, 'F', 1. )
+         !
+         !
+      CASE(  32  )       !==  time varying 3D field  ==!   = F( local deformation rate and gridscale ) (Smagorinsky)
+         !
+         IF( ln_dynldf_lap .OR. ln_dynldf_blp  ) THEN        ! laplacian operator : (C_smag/pi)^2 L^2 |D|
+            !
+            zcmsmag   = (rn_csmc/rpi)**2                                            ! (C_smag/pi)^2
+            zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 12._wp * 12._wp * zcmsmag ) ! lower limit stability factor scaling
+            zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rn_Dt )               ! upper limit stability factor scaling
+            IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo          ! provide |U|L^3/12 lower limit instead 
+            !                                                                       ! of |U|L^3/16 in blp case
+            DO jk = 1, jpkm1
+               !
+               DO_2D_00_00
+                  zdb =    ( uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) -  uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) )  &
+                       &                      * r1_e1t(ji,jj) * e2t(ji,jj)                           &
+                       & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) -  vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) )  &
+                       &                      * r1_e2t(ji,jj) * e1t(ji,jj)
+                  dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)
+               END_2D
+               !
+               DO_2D_10_10
+                  zdb =   (  uu(ji,jj+1,jk,Kbb) * r1_e1u(ji,jj+1) -  uu(ji,jj,jk,Kbb) * r1_e1u(ji,jj) )  &
+                       &                        * r1_e2f(ji,jj)   * e1f(ji,jj)                       &
+                       & + ( vv(ji+1,jj,jk,Kbb) * r1_e2v(ji+1,jj) -  vv(ji,jj,jk,Kbb) * r1_e2v(ji,jj) )  &
+                       &                        * r1_e1f(ji,jj)   * e2f(ji,jj)
+                  dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)
+               END_2D
+               !
+            END DO
+            !
+            CALL lbc_lnk_multi( 'ldfdyn', dtensq, 'T', 1. )  ! lbc_lnk on dshesq not needed
+            !
+            DO jk = 1, jpkm1
+              !
+               DO_2D_00_00
+                  !
+                  zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
+                  zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
+                  !
+                  zdelta         = zcmsmag * esqt(ji,jj)                                        ! L^2 * (C_smag/pi)^2
+                  ahmt(ji,jj,jk) = zdelta * SQRT(          dtensq(ji  ,jj,jk) +                         &
+                     &                            r1_4 * ( dshesq(ji  ,jj,jk) + dshesq(ji  ,jj-1,jk) +  &
+                     &                                     dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) )
+                  ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
+                  ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
+                  !
+               END_2D
+               !
+               DO_2D_10_10
+                  !
+                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, kbb) * uu(ji  ,jj+1,jk, kbb) + vv(ji+1,jj  ,jk, kbb) * vv(ji+1,jj  ,jk, kbb)
+                  zu2pv2_ij    = uu(ji  ,jj  ,jk, kbb) * uu(ji  ,jj  ,jk, kbb) + vv(ji  ,jj  ,jk, kbb) * vv(ji  ,jj  ,jk, kbb)
+                  !
+                  zdelta         = zcmsmag * esqf(ji,jj)                                        ! L^2 * (C_smag/pi)^2
+                  ahmf(ji,jj,jk) = zdelta * SQRT(          dshesq(ji  ,jj,jk) +                         &
+                     &                            r1_4 * ( dtensq(ji  ,jj,jk) + dtensq(ji  ,jj+1,jk) +  &
+                     &                                     dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) )
+                  ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac  * |U|L/2
+                  ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk),                                    zdelta * zstabf_up )   ! Impose upper limit == maxfac  * L^2/(4*2dt)
+                  !
+               END_2D
+               !
+            END DO
+            !
+         ENDIF
+         !
+         IF( ln_dynldf_blp ) THEN      ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8)
+            !                          !                      = sqrt( A_lap_smag L^2/8 )
+            !                          ! stability limits already applied to laplacian values
+            !                          ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4
+            DO jk = 1, jpkm1
+               DO_2D_00_00
+                  ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
+               END_2D
+               DO_2D_10_10
+                  ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
+               END_2D
+            END DO
+            !
+         ENDIF
+         !
+         CALL lbc_lnk_multi( 'ldfdyn', ahmt, 'T', 1. , ahmf, 'F', 1. )
+         !
+      END SELECT
+      !
+      CALL iom_put( "ahmt_2d", ahmt(:,:,1) )   ! surface u-eddy diffusivity coeff.
+      CALL iom_put( "ahmf_2d", ahmf(:,:,1) )   ! surface v-eddy diffusivity coeff.
+      CALL iom_put( "ahmt_3d", ahmt(:,:,:) )   ! 3D      u-eddy diffusivity coeff.
+      CALL iom_put( "ahmf_3d", ahmf(:,:,:) )   ! 3D      v-eddy diffusivity coeff.
+      !
+      IF( ln_timing )   CALL timing_stop('ldf_dyn')
+      !
+   END SUBROUTINE ldf_dyn
+
+   !!======================================================================
+END MODULE ldfdyn
