Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8177)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 8178)
@@ -25,30 +25,29 @@
    !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain
-   USE dynadv, ONLY: ln_dynadv_vec
-   USE zdf_oce         ! ocean vertical physics
-   USE zdfdrg          ! ocean vertical physics: top/bottom friction
-   USE ldftra          ! lateral physics: eddy diffusivity coef.
-   USE ldfdyn          ! lateral physics: eddy viscosity   coef.
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbc_ice         ! Surface boundary condition: ice fields
-   USE icb_oce         ! Icebergs
-   USE icbdia          ! Iceberg budgets
-   USE sbcssr          ! restoring term toward SST/SSS climatology
-   USE phycst          ! physical constants
-   USE zdfmxl          ! mixed layer
-   USE dianam          ! build name of file (routine)
-!   USE zdfddm          ! vertical  physics: double diffusion
-   USE diahth          ! thermocline diagnostics
-   USE wet_dry         ! wetting and drying
-   USE sbcwave         ! wave parameters
+   USE oce            ! ocean dynamics and tracers 
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE dianam         ! build name of file (routine)
+   USE diahth         ! thermocline diagnostics
+   USE dynadv   , ONLY: ln_dynadv_vec
+   USE icb_oce        ! Icebergs
+   USE icbdia         ! Iceberg budgets
+   USE ldftra         ! lateral physics: eddy diffusivity coef.
+   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE sbcssr         ! restoring term toward SST/SSS climatology
+   USE sbcwave        ! wave parameters
+   USE wet_dry        ! wetting and drying
+   USE zdf_oce        ! ocean vertical physics
+   USE zdfdrg         ! ocean vertical physics: top/bottom friction
+   USE zdfmxl         ! mixed layer
    !
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager  ! I/O manager
-   USE diatmb          ! Top,middle,bottom output
-   USE dia25h          ! 25h Mean output
-   USE iom
-   USE ioipsl
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager ! I/O manager
+   USE diatmb         ! Top,middle,bottom output
+   USE dia25h         ! 25h Mean output
+   USE iom            ! 
+   USE ioipsl         ! 
 
 #if defined key_lim2
@@ -194,6 +193,6 @@
       ENDIF
          
-      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
-      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
+      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
+      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
       IF ( iom_use("sbu") ) THEN
          DO jj = 1, jpj
@@ -206,6 +205,6 @@
       ENDIF
       
-      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
-      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
+      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
+      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
       IF ( iom_use("sbv") ) THEN
          DO jj = 1, jpj
@@ -229,7 +228,7 @@
       ENDIF
 
-      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
-      CALL iom_put( "avs" , avs                        )    ! S vert. eddy diff. coef.
-      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
+      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
+      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
+      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
 
       IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
@@ -247,7 +246,7 @@
          END DO
          CALL lbc_lnk( z2d, 'T', 1. )
-         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
+         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
          z2d(:,:) = SQRT( z2d(:,:) )
-         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
+         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
       ENDIF
          
@@ -262,5 +261,5 @@
             END DO
          END DO
-         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
+         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
       ENDIF
 
@@ -274,23 +273,22 @@
             END DO
          END DO
-         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
+         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
       ENDIF
       !
       IF ( iom_use("eken") ) THEN
-         z3d(:,:,jk) = 0._wp                               !      kinetic energy 
+         z3d(:,:,jk) = 0._wp 
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zztmp  = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-                  z3d(ji,jj,jk) = 0.25_wp * zztmp * (                                    &
-                     &            un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
-                     &          + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
-                     &          + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
-                     &          + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
+                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
+                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
+                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
+                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
                END DO
             END DO
          END DO
          CALL lbc_lnk( z3d, 'T', 1. )
-         CALL iom_put( "eken", z3d )           
+         CALL iom_put( "eken", z3d )                 ! kinetic energy
       ENDIF
       !
@@ -304,10 +302,10 @@
             z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
-         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
-         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum
+         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
+         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
       ENDIF
       
       IF( iom_use("u_heattr") ) THEN
-         z2d(:,:) = 0.e0 
+         z2d(:,:) = 0._wp 
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -318,5 +316,5 @@
          END DO
          CALL lbc_lnk( z2d, 'U', -1. )
-         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
+         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
       ENDIF
 
@@ -331,5 +329,5 @@
          END DO
          CALL lbc_lnk( z2d, 'U', -1. )
-         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
+         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
       ENDIF
 
@@ -340,5 +338,5 @@
             z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
          END DO
-         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
+         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
       ENDIF
       
@@ -353,9 +351,9 @@
          END DO
          CALL lbc_lnk( z2d, 'V', -1. )
-         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
+         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
       ENDIF
 
       IF( iom_use("v_salttr") ) THEN
-         z2d(:,:) = 0.e0 
+         z2d(:,:) = 0._wp 
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
@@ -366,22 +364,19 @@
          END DO
          CALL lbc_lnk( z2d, 'V', -1. )
-         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
-      ENDIF
-
-      ! Vertical integral of temperature
+         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
+      ENDIF
+
       IF( iom_use("tosmint") ) THEN
-         z2d(:,:)=0._wp
+         z2d(:,:) = 0._wp
          DO jk = 1, jpkm1
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
                END DO
             END DO
          END DO
          CALL lbc_lnk( z2d, 'T', -1. )
-         CALL iom_put( "tosmint", z2d ) 
-      ENDIF
-
-      ! Vertical integral of salinity
+         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
+      ENDIF
       IF( iom_use("somint") ) THEN
          z2d(:,:)=0._wp
@@ -389,21 +384,18 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
                END DO
             END DO
          END DO
          CALL lbc_lnk( z2d, 'T', -1. )
-         CALL iom_put( "somint", z2d ) 
-      ENDIF
-
-      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2)
-      !
-
-      IF (ln_diatmb) THEN      ! If we want tmb values 
-         CALL dia_tmb 
-      ENDIF 
-      IF (ln_dia25h) THEN
-         CALL dia_25h( kt )
-      ENDIF 
+         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
+      ENDIF
+
+      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
+      !
+
+      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values 
+          
+      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
 
       IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
@@ -729,5 +721,5 @@
          CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
             &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
-         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
+         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
             &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
 
@@ -856,5 +848,5 @@
       CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
       CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
-      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
+      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
       IF( ln_zdfddm ) THEN
          CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 8177)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90	(revision 8178)
@@ -37,11 +37,11 @@
 
    !                      ! 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)
+   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, PARAMETER, PUBLIC ::   np_lap    = 10   ,   np_blp    = 20  !: iso-level operator
+   INTEGER, PARAMETER, PUBLIC ::   np_lap_i  = 11                       !: iso-neutral or geopotential operator
 
-   INTEGER ::   nldf   ! type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
+   INTEGER, PUBLIC ::   nldf   !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals)
 
    !! * Substitutions
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 8177)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90	(revision 8178)
@@ -37,4 +37,6 @@
    PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akzu, akzv   !: vertical component of rotated lateral viscosity
+   
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso) 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      -
@@ -53,6 +55,6 @@
       !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
-         &      zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
+      ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
+         &      akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
          !
       IF( dyn_ldf_iso_alloc /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
@@ -99,8 +101,7 @@
       !!
       !! ** Action :
-      !!        Update (ua,va) arrays with the before geopotential biharmonic
-      !!      mixing trend.
-      !!        Update (avmu,avmv) to accompt for the diagonal vertical component
-      !!      of the rotated operator in dynzdf module
+      !!       -(ua,va) updated with the before geopotential harmonic mixing trend
+      !!       -(akzu,akzv) to accompt for the diagonal vertical component
+      !!                    of the rotated operator in dynzdf module
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
@@ -144,12 +145,6 @@
          CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
          CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
- 
-!!bug
-         IF( kt == nit000 ) then
-            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
-               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
-         endif
-!!end
-      ENDIF
+         !
+       ENDIF
 
       !                                                ! ===============
@@ -365,7 +360,7 @@
                            + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
-               ! update avmu (add isopycnal vertical coefficient to avmu)
-               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
-               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
+               ! vertical mixing coefficient (akzu)
+               ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               akzu(ji,jj,jk) = ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0
             END DO
          END DO
@@ -391,7 +386,7 @@
                   &        + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
                   &                    +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
-               ! update avmv (add isopycnal vertical coefficient to avmv)
-               ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
-               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0
+               ! vertical mixing coefficient (akzv)
+               ! Note: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0
+               akzv(ji,jj,jk) = ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0
             END DO
          END DO
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 8177)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 8178)
@@ -6,9 +6,9 @@
    !! History :  1.0  !  2005-11  (G. Madec)  Original code
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
-   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dyn_zdf       : Update the momentum trend with the vertical diffusion
+   !!            4.0  !  2017-06  (G. Madec) remove the explicit time-stepping option + avm at t-point
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   dyn_zdf       : compute the after velocity through implicit calculation of vertical mixing
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
@@ -18,5 +18,7 @@
    USE zdf_oce        ! ocean vertical physics variables
    USE zdfdrg         ! vertical physics: top/bottom drag coef.
-   USE dynadv   , ONLY: ln_dynadv_vec ! dynamics: advection form
+   USE dynadv    ,ONLY: ln_dynadv_vec    ! dynamics: advection form
+   USE dynldf    ,ONLY: nldf, np_lap_i   ! dynamics: type of lateral mixing 
+   USE dynldf_iso,ONLY: akzu, akzv       ! dynamics: vertical component of rotated lateral mixing 
    USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
    USE trd_oce        ! trends: ocean variables
@@ -26,5 +28,4 @@
    USE lib_mpp        ! MPP library
    USE prtctl         ! Print control
-   USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
 
@@ -49,47 +50,5 @@
       !!                  ***  ROUTINE dyn_zdf  ***
       !!
-      !! ** Purpose :   compute the vertical ocean dynamics physics.
-      !!---------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf')
-      !
-      !                             !* set time step
-      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
-      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
-      ENDIF
-
-      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
-         ztrdu(:,:,:) = ua(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:)
-      ENDIF
-      !                             !* compute lateral mixing trend and add it to the general trend
-      !
-      CALL dyn_zdf_imp( kt, r2dt )
-      !
-      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
-         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
-         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
-      ENDIF
-      !                                          ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
-         !
-      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf')
-      !
-   END SUBROUTINE dyn_zdf
-
-
-   SUBROUTINE dyn_zdf_imp( kt, p2dt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_zdf_imp  ***
-      !!                   
-      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
+      !! ** Purpose :   compute the trend due to the vert. momentum diffusion
       !!              together with the Leap-Frog time stepping using an 
       !!              implicit scheme.
@@ -100,24 +59,24 @@
       !!               - update the after velocity with the implicit vertical mixing.
       !!      This requires to solver the following system: 
-      !!         ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ]
+      !!         ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ]
       !!      with the following surface/top/bottom boundary condition:
       !!      surface: wind stress input (averaged over kt-1/2 & kt+1/2)
       !!      top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90)
       !!
-      !! ** Action :   (ua,va) after velocity 
+      !! ** Action :   (ua,va)   after velocity 
       !!---------------------------------------------------------------------
       INTEGER , INTENT(in) ::  kt     ! ocean time-step index
-      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step
-      !
-      INTEGER  ::   ji, jj, jk    ! dummy loop indices
-      INTEGER  ::   ikbu, ikbv    ! local integers
-      REAL(wp) ::   zzwi, ze3ua   ! local scalars
-      REAL(wp) ::   zzws, ze3va   !   -      -
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwd, zws
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
-      !
-      IF( kt == nit000 ) THEN
+      !
+      INTEGER  ::   ji, jj, jk         ! dummy loop indices
+      INTEGER  ::   iku, ikv           ! local integers
+      REAL(wp) ::   zzwi, ze3ua, zdt   ! local scalars
+      REAL(wp) ::   zzws, ze3va        !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::  zwi, zwd, zws   ! 3D workspace 
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdu, ztrdv   !  -      -
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf')
+      !
+      IF( kt == nit000 ) THEN       !* initialization
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
@@ -128,54 +87,36 @@
          ENDIF
       ENDIF
-      !
-      !              !==  Time step dynamics  ==!
-      !
-      IF( ln_dynadv_vec .OR. ln_linssh ) THEN      ! applied on velocity
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
-            va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
-         END DO
-      ELSE                                         ! applied on thickness weighted velocity
+      !                             !* set time step
+      IF( neuler == 0 .AND. kt == nit000     ) THEN   ;   r2dt =      rdt   ! = rdt (restart with Euler time stepping)
+      ELSEIF(               kt <= nit000 + 1 ) THEN   ;   r2dt = 2. * rdt   ! = 2 rdt (leapfrog)
+      ENDIF
+
+      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
+         ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 
+         ztrdu(:,:,:) = ua(:,:,:)
+         ztrdv(:,:,:) = va(:,:,:)
+      ENDIF
+      !
+      !              !==  RHS: Leap-Frog time stepping on all trends but the vertical mixing  ==!   (put in ua,va)
+      !
+      !                    ! time stepping except vertical diffusion
+      IF( ln_dynadv_vec .OR. ln_linssh ) THEN   ! applied on velocity
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)
+            va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)
+         END DO
+      ELSE                                      ! applied on thickness weighted velocity
          DO jk = 1, jpkm1
             ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
-               &          + p2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
+               &          + r2dt * e3u_n(:,:,jk) * ua(:,:,jk)  ) / e3u_a(:,:,jk) * umask(:,:,jk)
             va(:,:,jk) = (         e3v_b(:,:,jk) * vb(:,:,jk)  &
-               &          + p2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
-         END DO
-      ENDIF
-      !
-      !              !==  Apply semi-implicit bottom friction  ==!
-      !
-      ! Only needed for semi-implicit bottom friction setup. The explicit
-      ! bottom friction has been included in "u(v)a" which act as the R.H.S
-      ! column vector of the tri-diagonal matrix equation
-      !
-      IF( ln_drgimp ) THEN
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               ikbu = mbku(ji,jj)       ! ocean bottom level at u- and v-points 
-               ikbv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
-               avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1)
-               avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1)
-            END DO
-         END DO
-         IF ( ln_isfcav ) THEN
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  ikbu = miku(ji,jj)       ! ocean top level at u- and v-points 
-                  ikbv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
-                  ! top Cd is masked (=0 outside cavities) no need of test on mik>=2
-                  avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu)
-                  avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv)
-               END DO
-            END DO
-         END IF
-      ENDIF
-      !
-      ! With split-explicit free surface, barotropic stress is treated explicitly
-      ! Update velocities at the bottom.
-      ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 
-      !            not lead to the effective stress seen over the whole barotropic loop. 
-      ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
+               &          + r2dt * e3v_n(:,:,jk) * va(:,:,jk)  ) / e3v_a(:,:,jk) * vmask(:,:,jk)
+         END DO
+      ENDIF
+      !                    ! add top/bottom friction 
+      !     With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
+      !     J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 
+      !                not lead to the effective stress seen over the whole barotropic loop. 
+      !     G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a
       IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
          DO jk = 1, jpkm1        ! remove barotropic velocities
@@ -185,10 +126,10 @@
          DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ikbu = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
-               ikbv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
-               ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-               va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
+               iku = mbku(ji,jj)         ! ocean bottom level at u- and v-points 
+               ikv = mbkv(ji,jj)         ! (deepest ocean u- and v-points)
+               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
+               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
+               ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
+               va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
             END DO
          END DO
@@ -196,10 +137,10 @@
             DO jj = 2, jpjm1        
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ikbu = miku(ji,jj)         ! top ocean level at u- and v-points 
-                  ikbv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
-                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu)
-                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv)
-                  ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
-                  va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
+                  iku = miku(ji,jj)         ! top ocean level at u- and v-points 
+                  ikv = mikv(ji,jj)         ! (first wet ocean u- and v-points)
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)
+                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)
+                  ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
+                  va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
                END DO
             END DO
@@ -209,21 +150,37 @@
       !              !==  Vertical diffusion on u  ==!
       !
-      ! Matrix and second member construction
-      ! bottom boundary condition: both zwi and zws must be masked as avmu can take
-      ! non zero value at the ocean bottom depending on the bottom friction used.
-      !
-      DO jk = 1, jpkm1        ! Matrix
-         DO jj = 2, jpjm1 
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
-               zzwi = - p2dt * avmu(ji,jj,jk  ) / ( ze3ua * e3uw_n(ji,jj,jk  ) )
-               zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) )
-               zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk  )
-               zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1        ! Surface boundary conditions
+      !                    !* Matrix construction
+      zdt = r2dt * 0.5
+      IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu)
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
+                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) + akzu(ji,jj,jk  ) )   &
+                     &         / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) )   &
+                     &         / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                  zwi(ji,jj,jk) = zzwi
+                  zws(ji,jj,jk) = zzws
+                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+               END DO
+            END DO
+         END DO
+      ELSE                          ! standard case
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1 
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk)   ! after scale factor at T-point
+                  zzwi = - zdt * ( avm(ji+1,jj,jk  ) + avm(ji,jj,jk  ) ) / ( ze3ua * e3uw_n(ji,jj,jk  ) ) * wumask(ji,jj,jk  )
+                  zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)
+                  zwi(ji,jj,jk) = zzwi
+                  zws(ji,jj,jk) = zzws
+                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jj = 2, jpjm1     !* Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
@@ -231,5 +188,31 @@
          END DO
       END DO
-
+      !
+      !              !==  Apply semi-implicit bottom friction  ==!
+      !
+      !     Only needed for semi-implicit bottom friction setup. The explicit
+      !     bottom friction has been included in "u(v)a" which act as the R.H.S
+      !     column vector of the tri-diagonal matrix equation
+      !
+      IF ( ln_drgimp ) THEN      ! implicit bottom friction
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               iku = mbku(ji,jj)       ! ocean bottom level at u- and v-points
+               ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
+               zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
+            END DO
+         END DO
+         IF ( ln_isfcav ) THEN   ! top friction (always implicit)
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  !!gm   top Cd is masked (=0 outside cavities) no need of test on mik>=2  ==>> it has been suppressed
+                  iku = miku(ji,jj)       ! ocean top level at u- and v-points 
+                  ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)   ! after scale factor at T-point
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
+               END DO
+            END DO
+         END IF
+      ENDIF
+      !
       ! Matrix inversion starting from the first level
       !-----------------------------------------------------------------------
@@ -258,5 +241,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ze3ua =  ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 
-            ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
+            ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   &
                &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
          END DO
@@ -285,21 +268,37 @@
       !              !==  Vertical diffusion on v  ==!
       !
-      ! Matrix and second member construction
-      ! bottom boundary condition: both zwi and zws must be masked as avmv can take
-      ! non zero value at the ocean bottom depending on the bottom friction used
-      !
-      DO jk = 1, jpkm1        ! Matrix
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
-               zzwi = - p2dt * avmv (ji,jj,jk  ) / ( ze3va * e3vw_n(ji,jj,jk  ) )
-               zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) )
-               zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
-               zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
-               zwd(ji,jj,jk) = 1._wp - zzwi - zzws
-            END DO
-         END DO
-      END DO
-      DO jj = 2, jpjm1        ! Surface boundary conditions
+      !                       !* Matrix construction
+      zdt = r2dt * 0.5
+      IF( nldf == np_lap_i ) THEN   ! rotated lateral mixing: add its vertical mixing (akzu)
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1   
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
+                  zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) + akzv(ji,jj,jk  ) )   &
+                     &         / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                  zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) )   &
+                     &         / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
+                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
+                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+               END DO
+            END DO
+         END DO
+      ELSE                          ! standard case
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1   
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk)   ! after scale factor at T-point
+                  zzwi = - zdt * ( avm(ji,jj+1,jk  )+ avm(ji,jj,jk  ) ) / ( ze3va * e3vw_n(ji,jj,jk  ) ) * wvmask(ji,jj,jk  )
+                  zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)
+                  zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk  )
+                  zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
+                  zwd(ji,jj,jk) = 1._wp - zzwi - zzws
+               END DO
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jj = 2, jpjm1        !* Surface boundary conditions
          DO ji = fs_2, fs_jpim1   ! vector opt.
             zwi(ji,jj,1) = 0._wp
@@ -307,4 +306,28 @@
          END DO
       END DO
+      !              !==  Apply semi-implicit top/bottom friction  ==!
+      !
+      !     Only needed for semi-implicit bottom friction setup. The explicit
+      !     bottom friction has been included in "u(v)a" which act as the R.H.S
+      !     column vector of the tri-diagonal matrix equation
+      !
+      IF( ln_drgimp ) THEN
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ikv = mbkv(ji,jj)       ! (deepest ocean u- and v-points)
+               ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
+               zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va           
+            END DO
+         END DO
+         IF ( ln_isfcav ) THEN
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  ikv = mikv(ji,jj)       ! (first wet ocean u- and v-points)
+                  ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)   ! after scale factor at T-point
+                  zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va
+               END DO
+            END DO
+         ENDIF
+      ENDIF
 
       ! Matrix inversion
@@ -334,5 +357,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.          
             ze3va =  ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 
-            va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
+            va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   &
                &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
          END DO
@@ -359,7 +382,17 @@
       END DO
       !
-      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp')
-      !
-   END SUBROUTINE dyn_zdf_imp
+      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics
+         ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)
+         ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)
+         CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )
+         DEALLOCATE( ztrdu, ztrdv ) 
+      ENDIF
+      !                                          ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
+         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+         !
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf')
+      !
+   END SUBROUTINE dyn_zdf
 
    !!==============================================================================
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 8177)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 8178)
@@ -6,5 +6,5 @@
    !! history :  1.0  !  2002-06  (G. Madec)  Original code
    !!            3.2  !  2009-07  (G. Madec)  addition of avm
-   !!            4.0  !  2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
+   !!            4.0  !  2017-05  (G. Madec)  avm and drag coef. defined at t-point
    !!----------------------------------------------------------------------
    USE par_oce        ! ocean parameters
@@ -44,12 +44,9 @@
 
 
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm, avt, avs  !: vertical mixing coefficient (w-point) [m2/s]
-!!gm
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avmu , avmv    !: vertical viscosity coef at uw- & vw-pts  [m2/s]
-!!gm
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k, avt_k   !: avm, avt computed by turbulent closure alone
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile
-   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   avm, avt, avs  !: vertical mixing coefficients (w-point) [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avm_k , avt_k  !: Kz computed by turbulent closure alone [m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy          [m2/s2]
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:)     ::   avmb , avtb    !: background profile of avm and avt      [m2/s]
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:)   ::   avtb_2d        !: horizontal shape of background Kz profile [-]
 
    !!----------------------------------------------------------------------
@@ -65,8 +62,7 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( avm  (jpi,jpj,jpk) , avt  (jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,   &
-         &      avm_k(jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   & 
-         &      avmb(jpk) , avtb(jpk) ,  avtb_2d(jpi,jpj) ,                    &
-         &      avmu(jpi,jpj,jpk), avmv(jpi,jpj,jpk)    , STAT = zdf_oce_alloc )
+      ALLOCATE( avm (jpi,jpj,jpk) , avm_k(jpi,jpj,jpk) , avs(jpi,jpj,jpk) ,   &
+         &      avt (jpi,jpj,jpk) , avt_k(jpi,jpj,jpk) , en (jpi,jpj,jpk) ,   & 
+         &      avmb(jpk)         , avtb(jpk)          , avtb_2d(jpi,jpj) , STAT = zdf_oce_alloc )
          !
       IF( zdf_oce_alloc /= 0 )   CALL ctl_warn('zdf_oce_alloc: failed to allocate arrays')
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8177)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8178)
@@ -29,5 +29,6 @@
    !!            3.7  !  2014-10  (G. Madec)  LDF simplication 
    !!             -   !  2014-12  (G. Madec) remove KPP scheme
-   !!             -   !  2015-11  (J. Chanut) free surface simplification
+   !!             -   !  2015-11  (J. Chanut) free surface simplification (remove filtered free surface)
+   !!            4.0  !  2017-05  (G. Madec)  introduction of the vertical physics manager (zdfphy)
    !!----------------------------------------------------------------------
 
@@ -36,7 +37,4 @@
    !!----------------------------------------------------------------------
    USE step_oce         ! time stepping definition modules
-!!gm to be removed when removing avmu, avmv
-   USE zdf_oce        ! ocean vertical physics variables
-!!gm
    !
    USE iom              ! xIOs server
@@ -48,5 +46,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -131,18 +129,4 @@
                          CALL zdf_phy( kstp )         ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
 
-
-!!gm  ===>>>   TO BE REMOVED   (require changes in zdf_oce, dynzdf(_imp,_exp), dynldf_iso, diawri)
-      DO jk = 1, jpkm1            !* vertical eddy viscosity at wu- and wv-points
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
-               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
-!!gm end
-
-
       !  LATERAL  PHYSICS
       !
