Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/AMM12/EXP00/namelist_cfg	(revision 8160)
@@ -340,8 +340,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       ! coefficients
    rn_avm0     =   0.1e-6     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/C1D_PAPA/EXP00/namelist_cfg	(revision 8160)
@@ -279,8 +279,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       ! coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg	(revision 8160)
@@ -268,8 +268,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       ! coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg	(revision 8160)
@@ -222,8 +222,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       !  Coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_cfg	(revision 8160)
@@ -259,8 +259,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       !  Coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8160)
@@ -914,8 +914,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       ! coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/ISOMIP/EXP00/namelist_cfg	(revision 8160)
@@ -339,7 +339,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-   !
    !                       ! coefficients
    rn_avm0     =   1.e-3     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/LOCK_EXCHANGE/EXP00/namelist_cfg	(revision 8160)
@@ -201,7 +201,4 @@
    ln_zdfnpc   = .false.      !  Non-Penetrative Convective algorithm
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-   !
    !                       ! coefficients
    rn_avm0     =   1.e-4      !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/OVERFLOW/EXP00/namelist_cfg	(revision 8160)
@@ -198,7 +198,4 @@
    ln_zdfnpc   = .false.      !  Non-Penetrative Convective algorithm
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-   !
    !                       ! coefficients
    rn_avm0     =   1.e-4      !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/CONFIG/TEST_CASES/WAD/EXP00/namelist_cfg	(revision 8160)
@@ -359,8 +359,4 @@
    ln_zdfswm   = .false.      ! surface  wave-induced mixing            (T => ln_wave=ln_sdw=T )
    !
-   !                       ! time-stepping
-   ln_zdfexp   = .false.      ! split-explicit (T) or implicit (F) scheme
-      nn_zdfexp=    3            !  number of sub-timestep for ln_zdfexp=T
-   !
    !                       ! coefficients
    rn_avm0     =   1.2e-4     !  vertical eddy viscosity   [m2/s]       (background Kz if ln_zdfcst=F)
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 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90	(revision 8160)
@@ -6,15 +6,17 @@
    !! 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
-   !!   dyn_zdf_init  : initializations of the vertical diffusion scheme
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
+   USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain variables 
+   USE sbc_oce        ! surface boundary condition: ocean
    USE zdf_oce        ! ocean vertical physics variables
-   USE dynzdf_exp     ! vertical diffusion: explicit (dyn_zdf_exp     routine)
-   USE dynzdf_imp     ! vertical diffusion: implicit (dyn_zdf_imp     routine)
+   USE zdfdrg         ! vertical physics: top/bottom drag coef.
+   USE dynadv   , ONLY: ln_dynadv_vec ! dynamics: advection form
    USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
    USE trd_oce        ! trends: ocean variables
@@ -30,13 +32,12 @@
    PRIVATE
 
-   PUBLIC   dyn_zdf       !  routine called by step.F90
-   PUBLIC   dyn_zdf_init  !  routine called by opa.F90
-
-   INTEGER  ::   nzdf = 0   ! type vertical diffusion algorithm used, defined from ln_zdf... namlist logicals
+   PUBLIC   dyn_zdf   !  routine called by step.F90
+
+   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
@@ -57,22 +58,18 @@
       IF( nn_timing == 1 )   CALL timing_start('dyn_zdf')
       !
-      !                                          ! set time step
+      !                             !* 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
+      IF( l_trddyn )   THEN         !* temporary save of ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
          ztrdu(:,:,:) = ua(:,:,:)
          ztrdv(:,:,:) = va(:,:,:)
       ENDIF
-
-      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      !
-      CASE ( 0 )   ;   CALL dyn_zdf_exp( kt, r2dt )      ! explicit scheme
-      CASE ( 1 )   ;   CALL dyn_zdf_imp( kt, r2dt )      ! implicit scheme
-      !
-      END SELECT
-
+      !                             !* 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(:,:,:)
@@ -90,36 +87,279 @@
 
 
-   SUBROUTINE dyn_zdf_init
+   SUBROUTINE dyn_zdf_imp( kt, p2dt )
       !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE dyn_zdf_init  ***
+      !!                  ***  ROUTINE dyn_zdf_imp  ***
+      !!                   
+      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
+      !!              together with the Leap-Frog time stepping using an 
+      !!              implicit scheme.
       !!
-      !! ** Purpose :   initialization of the vertical diffusion scheme
+      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
+      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
+      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
+      !!               - 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] ]
+      !!      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)
       !!
-      !! ** Method  :   implicit (euler backward) scheme (default)
-      !!                explicit (time-splitting) scheme if ln_zdfexp=T
+      !! ** 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
       !!----------------------------------------------------------------------
-      USE zdftke
-      USE zdfgls
-      !!----------------------------------------------------------------------
-      !
-      ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module)
-      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
-      ELSE                   ;   nzdf = 1           ! use implicit scheme
-      ENDIF
-      !
-      ! Force implicit schemes
-      IF( ln_zdftke .OR. ln_zdfgls   )   nzdf = 1   ! TKE or GLS physics
-      IF( ln_dynldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
-      IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
-      !
-      IF(lwp) THEN                                  ! Print the choice
-         WRITE(numout,*)
-         WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
-      ENDIF
-      !
-   END SUBROUTINE dyn_zdf_init
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_imp')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+         !
+         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
+         ELSE                   ;    r_vvl = 1._wp
+         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
+         DO jk = 1, jpkm1
+            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
+               &          + p2dt * 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
+      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
+         DO jk = 1, jpkm1        ! remove barotropic velocities
+            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
+            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
+         END DO
+         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
+            END DO
+         END DO
+         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
+            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
+               END DO
+            END DO
+         END IF
+      ENDIF
+      !
+      !              !==  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
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zwi(ji,jj,1) = 0._wp
+            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+         END DO
+      END DO
+
+      ! Matrix inversion starting from the first level
+      !-----------------------------------------------------------------------
+      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+      !
+      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+      !        (        ...               )( ...  ) ( ...  )
+      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+      !
+      !   m is decomposed in the product of an upper and a lower triangular matrix
+      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+      !   The solution (the after velocity) is in ua
+      !-----------------------------------------------------------------------
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
+         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) )   &
+               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 1, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      !              !==  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
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            zwi(ji,jj,1) = 0._wp
+            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+         END DO
+      END DO
+
+      ! Matrix inversion
+      !-----------------------------------------------------------------------
+      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+      !
+      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+      !        (        ...               )( ...  ) ( ...  )
+      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+      !
+      !   m is decomposed in the product of an upper and lower triangular matrix
+      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+      !   The solution (after velocity) is in 2d array va
+      !-----------------------------------------------------------------------
+      !
+      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
+         DO jj = 2, jpjm1   
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
+         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) )   &
+               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
+         END DO
+      END DO
+      DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 1, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp')
+      !
+   END SUBROUTINE dyn_zdf_imp
 
    !!==============================================================================
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_exp.F90	(revision 8143)
+++ 	(revision )
@@ -1,149 +1,0 @@
-MODULE dynzdf_exp
-   !!==============================================================================
-   !!                     ***  MODULE  dynzdf_exp  ***
-   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
-   !!==============================================================================
-   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
-   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
-   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
-   !!            3.7  !  2015-11  (J. Chanut) output velocities instead of trends
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dyn_zdf_exp   : update the momentum trend with the vertical diffusion using a split-explicit scheme
-   !!                   and perform the Leap-Frog time integration.
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE dom_oce        ! ocean space and time domain
-   USE phycst         ! physical constants
-   USE zdf_oce        ! ocean vertical physics
-   USE dynadv   , ONLY: ln_dynadv_vec ! Momentum advection form
-   USE sbc_oce        ! surface boundary condition: ocean
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dyn_zdf_exp   ! called by step.F90
-   
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dyn_zdf_exp( kt, p2dt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_zdf_exp  ***
-      !!                   
-      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
-      !!              and perform the Leap-Frog time stepping.
-      !!
-      !! ** Method  : - Split-explicit forward time stepping.
-      !!      The vertical mixing of momentum is given by:
-      !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
-      !!      Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
-      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
-      !!      Add this trend to the general trend ua :
-      !!         ua = ua + dz( avmu dz(u) )
-      !!              - Leap-Frog time stepping (Asselin filter will be applied in dyn_nxt) 
-      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
-      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
-      !!
-      !! ** Action : - (ua,va) after velocity
-      !!---------------------------------------------------------------------
-      INTEGER , INTENT(in) ::   kt     ! ocean time-step index
-      REAL(wp), INTENT(in) ::   p2dt   ! time-step 
-      !
-      INTEGER  ::   ji, jj, jk, jl     ! dummy loop indices
-      REAL(wp) ::   zlavmr, zua, zva   ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwx, zwy, zwz, zww
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('dyn_zdf_exp')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zwx, zwy, zwz, zww ) 
-      !
-      IF( kt == nit000 .AND. lwp ) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion - explicit operator'
-         WRITE(numout,*) '~~~~~~~~~~~ '
-      ENDIF
-      !
-      !                 !==  vertical mixing trend  ==!
-      !
-      zlavmr = 1. / REAL( nn_zdfexp )
-      !
-      DO jj = 2, jpjm1           ! Surface boundary condition
-         DO ji = 2, jpim1
-            zwy(ji,jj,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * r1_rau0
-            zww(ji,jj,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_rau0
-         END DO  
-      END DO  
-      DO jk = 1, jpk             ! Initialization of x, z and contingently trends array
-         DO jj = 2, jpjm1 
-            DO ji = 2, jpim1
-               zwx(ji,jj,jk) = ub(ji,jj,jk)
-               zwz(ji,jj,jk) = vb(ji,jj,jk)
-            END DO  
-         END DO  
-      END DO  
-      !
-      DO jl = 1, nn_zdfexp       ! Time splitting loop
-         !
-         DO jk = 2, jpk                ! First vertical derivative
-            DO jj = 2, jpjm1 
-               DO ji = 2, jpim1
-                  zwy(ji,jj,jk) = avmu(ji,jj,jk) * ( zwx(ji,jj,jk-1) - zwx(ji,jj,jk) ) / e3uw_n(ji,jj,jk) 
-                  zww(ji,jj,jk) = avmv(ji,jj,jk) * ( zwz(ji,jj,jk-1) - zwz(ji,jj,jk) ) / e3vw_n(ji,jj,jk)
-               END DO  
-            END DO  
-         END DO  
-         DO jk = 1, jpkm1              ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
-            DO jj = 2, jpjm1 
-               DO ji = 2, jpim1
-                  zua = zlavmr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
-                  zva = zlavmr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-                  va(ji,jj,jk) = va(ji,jj,jk) + zva
-                  !
-                  zwx(ji,jj,jk) = zwx(ji,jj,jk) + p2dt * zua * umask(ji,jj,jk)
-                  zwz(ji,jj,jk) = zwz(ji,jj,jk) + p2dt * zva * vmask(ji,jj,jk)
-               END DO  
-            END DO  
-         END DO  
-      END DO                     ! End of time splitting
-      !
-      !
-      !                 !==  Leap-Frog time integration  ==!
-      !
-      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
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = (          e3u_b(:,:,jk) * ub(:,:,jk)    &
-               &           + p2dt * 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
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwx, zwy, zwz, zww ) 
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_exp')
-      !
-   END SUBROUTINE dyn_zdf_exp
-
-   !!==============================================================================
-END MODULE dynzdf_exp
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 8143)
+++ 	(revision )
@@ -1,324 +1,0 @@
-MODULE dynzdf_imp
-   !!======================================================================
-   !!                    ***  MODULE  dynzdf_imp  ***
-   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend, implicit scheme
-   !!======================================================================
-   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
-   !!            8.0  !  1997-05  (G. Madec)  vertical component of isopycnal
-   !!   NEMO     0.5  !  2002-08  (G. Madec)  F90: Free form and module
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
-   !!            3.4  !  2012-01  (H. Liu) Semi-implicit bottom friction
-   !!            4.0  !  2017-05  (G. Madec)  drag coef. defined at t-point (zdfdrg.F90)
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   dyn_zdf_imp   : compute the vertical diffusion using a implicit time-stepping
-   !!                   together with the Leap-Frog time integration.
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers
-   USE phycst         ! physical constants
-   USE dom_oce        ! ocean space and time domain
-   USE domvvl         ! variable volume
-   USE sbc_oce        ! surface boundary condition: ocean
-   USE dynadv   , ONLY: ln_dynadv_vec ! Momentum advection form
-   USE zdf_oce        ! ocean vertical physics
-   USE zdfdrg         ! vertical physics: top/bottom drag coef.
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   dyn_zdf_imp   ! called by step.F90
-
-   REAL(wp) ::  r_vvl     ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE dyn_zdf_imp( kt, p2dt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE dyn_zdf_imp  ***
-      !!                   
-      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
-      !!              together with the Leap-Frog time stepping using an 
-      !!              implicit scheme.
-      !!
-      !! ** Method  :  - Leap-Frog time stepping on all trends but the vertical mixing
-      !!         ua =         ub + 2*dt *       ua             vector form or linear free surf.
-      !!         ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a   otherwise
-      !!               - 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] ]
-      !!      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 
-      !!---------------------------------------------------------------------
-      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
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-         !
-         If( ln_linssh ) THEN   ;    r_vvl = 0._wp    ! non-linear free surface indicator
-         ELSE                   ;    r_vvl = 1._wp
-         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
-         DO jk = 1, jpkm1
-            ua(:,:,jk) = (         e3u_b(:,:,jk) * ub(:,:,jk)  &
-               &          + p2dt * 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
-      IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
-         DO jk = 1, jpkm1        ! remove barotropic velocities
-            ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)
-            va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)
-         END DO
-         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
-            END DO
-         END DO
-         IF( ln_isfcav ) THEN    ! Ocean cavities (ISF)
-            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
-               END DO
-            END DO
-         END IF
-      ENDIF
-      !
-      !              !==  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
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwi(ji,jj,1) = 0._wp
-            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-         END DO
-      END DO
-
-      ! Matrix inversion starting from the first level
-      !-----------------------------------------------------------------------
-      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
-      !
-      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
-      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
-      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
-      !        (        ...               )( ...  ) ( ...  )
-      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
-      !
-      !   m is decomposed in the product of an upper and a lower triangular matrix
-      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
-      !   The solution (the after velocity) is in ua
-      !-----------------------------------------------------------------------
-      !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      !
-      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
-         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) )   &
-               &                                      / ( ze3ua * rau0 ) * umask(ji,jj,1) 
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      !
-      DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  ==!
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      !              !==  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
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zwi(ji,jj,1) = 0._wp
-            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
-         END DO
-      END DO
-
-      ! Matrix inversion
-      !-----------------------------------------------------------------------
-      !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
-      !
-      !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
-      !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
-      !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
-      !        (        ...               )( ...  ) ( ...  )
-      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
-      !
-      !   m is decomposed in the product of an upper and lower triangular matrix
-      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
-      !   The solution (after velocity) is in 2d array va
-      !-----------------------------------------------------------------------
-      !
-      DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  ==
-         DO jj = 2, jpjm1   
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      !
-      DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  ==!
-         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) )   &
-               &                                      / ( ze3va * rau0 ) * vmask(ji,jj,1) 
-         END DO
-      END DO
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
-            END DO
-         END DO
-      END DO
-      !
-      DO jj = 2, jpjm1        !==  third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  ==!
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
-         END DO
-      END DO
-      DO jk = jpk-2, 1, -1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-               va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('dyn_zdf_imp')
-      !
-   END SUBROUTINE dyn_zdf_imp
-
-   !!==============================================================================
-END MODULE dynzdf_imp
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 8160)
@@ -4,11 +4,11 @@
    !! Ocean active tracers:  vertical component of the tracer mixing trend
    !!==============================================================================
-   !! History :  1.0  ! 2005-11  (G. Madec)  Original code
-   !!            3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
+   !! History :  1.0  !  2005-11  (G. Madec)  Original code
+   !!            3.0  !  2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!            4.0  !  2017-06  (G. Madec)  remove explict time-stepping option
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   tra_zdf       : Update the tracer trend with the vertical diffusion
-   !!   tra_zdf_init  : initialisation of the computation
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers variables
@@ -20,6 +20,4 @@
    USE ldftra         ! lateral diffusion: eddy diffusivity
    USE ldfslp         ! lateral diffusion: iso-neutral slope 
-   USE trazdf_exp     ! vertical diffusion: explicit (tra_zdf_exp routine)
-   USE trazdf_imp     ! vertical diffusion: implicit (tra_zdf_imp routine)
    USE trd_oce        ! trends: ocean variables
    USE trdtra         ! trends: tracer trend manager
@@ -29,5 +27,4 @@
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory allocation
    USE timing         ! Timing
 
@@ -35,8 +32,6 @@
    PRIVATE
 
-   PUBLIC   tra_zdf        ! routine called by step.F90
-   PUBLIC   tra_zdf_init   ! routine called by nemogcm.F90
-
-   INTEGER ::   nzdf = 0   ! type vertical diffusion algorithm used (defined from ln_zdf...  namlist logicals)
+   PUBLIC   tra_zdf       ! called by step.F90
+   PUBLIC   tra_zdf_imp   ! called by trczdf.F90
 
    !! * Substitutions
@@ -58,5 +53,5 @@
       !
       INTEGER  ::   jk                   ! Dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
       !!---------------------------------------------------------------------
       !
@@ -69,14 +64,13 @@
       ENDIF
       !
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+      IF( l_trdtra )   THEN                  !* Save ta and sa trends
+         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
       ENDIF
       !
-      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )    ;    CALL tra_zdf_exp( kt, nit000, 'TRA', r2dt, nn_zdfexp, tsb, tsa, jpts )  !   explicit scheme 
-      CASE ( 1 )    ;    CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt,            tsb, tsa, jpts )  !   implicit scheme 
-      END SELECT
+      !                                      !* compute lateral mixing trend and add it to the general trend
+      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) 
+
 !!gm WHY here !   and I don't like that !
       ! DRAKKAR SSS control {
@@ -97,5 +91,5 @@
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
+         DEALLOCATE( ztrdt , ztrds )
       ENDIF
       !                                          ! print mean trends (used for debugging)
@@ -107,37 +101,161 @@
    END SUBROUTINE tra_zdf
 
-
-   SUBROUTINE tra_zdf_init
+ 
+   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 
       !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE tra_zdf_init  ***
-      !!
-      !! ** Purpose :   Choose the vertical mixing scheme
-      !!
-      !! ** Method  :   Set nzdf from ln_zdfexp
-      !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T)
-      !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F)
-      !!      NB: rotation of lateral mixing operator or TKE & GLS schemes,
-      !!          an implicit scheme is required.
-      !!----------------------------------------------------------------------
-      !
-      ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module)
-      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
-      ELSE                   ;   nzdf = 1           ! use implicit scheme
-      ENDIF
-      !
-      ! Force implicit schemes when absolutely needed
-      IF( .NOT.ln_zdfcst             )   nzdf = 1   ! turbulent closure schemes
-      IF( ln_traldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
-      IF( ln_traldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
-      !
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme'
-         WRITE(numout,*) '~~~~~~~~~~~'
-         IF( nzdf ==  0 )   WRITE(numout,*) '      ===>>   Explicit time-splitting scheme'
-         IF( nzdf ==  1 )   WRITE(numout,*) '      ===>>   Implicit (euler backward) scheme'
-      ENDIF
-      !
-   END SUBROUTINE tra_zdf_init
+      !!                  ***  ROUTINE tra_zdf_imp  ***
+      !!
+      !! ** Purpose :   Compute the after tracer through a implicit computation
+      !!     of the vertical tracer diffusion (including the vertical component 
+      !!     of lateral mixing (only for 2nd order operator, for fourth order 
+      !!     it is already computed and add to the general trend in traldf) 
+      !!
+      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
+      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
+      !!      It is computed using a backward time scheme (t=after field)
+      !!      which provide directly the after tracer field.
+      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers
+      !!      Surface and bottom boundary conditions: no diffusive flux on
+      !!      both tracers (bottom, applied through the masked field avt).
+      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
+      !!
+      !! ** Action  : - pta  becomes the after tracer
+      !!---------------------------------------------------------------------
+      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
+      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
+      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field
+      !
+      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::  zrhs             ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zwi, zwt, zwd, zws
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp')
+      !
+      IF( kt == kit000 )  THEN
+         IF(lwp)WRITE(numout,*)
+         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
+         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
+      ENDIF
+      !                                               ! ============= !
+      DO jn = 1, kjpt                                 !  tracer loop  !
+         !                                            ! ============= !
+         !  Matrix construction
+         ! --------------------
+         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
+         !
+         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   &
+            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
+            !
+            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
+            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
+            ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk)
+            ENDIF
+            zwt(:,:,1) = 0._wp
+            !
+            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
+               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)  
+                        END DO
+                     END DO
+                  END DO
+               ELSE                          ! standard or triad iso-neutral operator
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
+                        END DO
+                     END DO
+                  END DO
+               ENDIF
+            ENDIF
+            !
+            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+!!gm BUG  I think, use e3w_a instead of e3w_n, not sure of that
+                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
+                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
+                     zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                 END DO
+               END DO
+            END DO
+            !
+            !! Matrix inversion from the first level
+            !!----------------------------------------------------------------------
+            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
+            !
+            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
+            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
+            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
+            !        (        ...               )( ...  ) ( ...  )
+            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
+            !
+            !   m is decomposed in the product of an upper and lower triangular matrix.
+            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
+            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal
+            !   and "superior" (above diagonal) components of the tridiagonal system.
+            !   The solution will be in the 4d array pta.
+            !   The 3d array zwt is used as a work space array.
+            !   En route to the solution pta is used a to evaluate the rhs and then 
+            !   used as a work space array: its value is modified.
+            !
+            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
+               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction)
+                  zwt(ji,jj,1) = zwd(ji,jj,1)
+               END DO
+            END DO
+            DO jk = 2, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF 
+         !         
+         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
+            DO ji = fs_2, fs_jpim1
+               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side
+                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
+               END DO
+            END DO
+         END DO
+         !
+         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
+            DO ji = fs_2, fs_jpim1
+               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+            END DO
+         END DO
+         DO jk = jpk-2, 1, -1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
+                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !                                            ! ================= !
+      END DO                                          !  end tracer loop  !
+      !                                               ! ================= !
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp')
+      !
+   END SUBROUTINE tra_zdf_imp
 
    !!==============================================================================
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 8143)
+++ 	(revision )
@@ -1,152 +1,0 @@
-MODULE trazdf_exp
-   !!==============================================================================
-   !!                    ***  MODULE  trazdf_exp  ***
-   !! Ocean  tracers:  vertical component of the tracer mixing trend using
-   !!                  a split-explicit time-stepping 
-   !!==============================================================================
-   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
-   !!            7.0  !  1991-11  (G. Madec)
-   !!                 !  1992-06  (M. Imbard)  correction on tracer trend loops
-   !!                 !  1996-01  (G. Madec)  statement function for e3
-   !!                 !  1997-05  (G. Madec)  vertical component of isopycnal
-   !!                 !  1997-07  (G. Madec)  geopotential diffusion in s-coord
-   !!                 !  2000-08  (G. Madec)  double diffusive mixing
-   !!   NEMO     1.0  !  2002-08  (G. Madec)  F90: Free form and module
-   !!             -   !  2004-08  (C. Talandier) New trends organisation
-   !!             -   !  2005-11  (G. Madec)  New organisation
-   !!            3.0  !  2008-04  (G. Madec)  leap-frog time stepping done in trazdf
-   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
-   !!----------------------------------------------------------------------
-
-   !!----------------------------------------------------------------------
-   !!   tra_zdf_exp   : compute the tracer the vertical diffusion trend using a
-   !!                   split-explicit time stepping and provide the after tracer
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and active tracers 
-   USE dom_oce        ! ocean space and time domain 
-   USE domvvl         ! variable volume levels
-   USE zdf_oce        ! ocean vertical physics
-   USE zdfddm         ! ocean vertical physics: double diffusion
-   USE trc_oce        ! share passive tracers/Ocean variables
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   tra_zdf_exp   ! routine called by step.F90
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE tra_zdf_exp( kt, kit000, cdtype, p2dt, ksts,   &
-      &                                        ptb , pta , kjpt )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_zdf_exp  ***
-      !!                   
-      !! ** Purpose :   Compute the after tracer fields due to the vertical
-      !!      tracer mixing alone, and then due to the whole tracer trend.
-      !!
-      !! ** Method  : - The after tracer fields due to the vertical diffusion
-      !!      of tracers alone is given by:
-      !!                ztb = ptb + p2dt difft
-      !!      where difft = dz( avt dz(ptb) ) = 1/e3t dk+1( avt/e3w dk(ptb) )
-      !!           (if ln_zdfddm=T use avs on salinity and passive tracers instead of avt)
-      !!      difft is evaluated with an Euler split-explit scheme using a
-      !!      no flux boundary condition at both surface and bottomi boundaries.
-      !!      (N.B. bottom condition is applied through the masked field avt).
-      !!              - the after tracer fields due to the whole trend is 
-      !!      obtained in leap-frog environment applied on thickness weighted tracer by :
-      !!          pta = [ ptb*e3tb + e3tn*( ztb - ptb + p2dt pta ) ] / e3tn
-      !!
-      !! ** Action : - after tracer fields pta
-      !!---------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
-      INTEGER                              , INTENT(in   ) ::   ksts     ! number of sub-time step
-      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field 
-      !
-      INTEGER  ::  ji, jj, jk, jn, jl   ! dummy loop indices
-      REAL(wp) ::  z1_ksts, ze3tr       ! local scalars
-      REAL(wp) ::  ztra, ze3tb    !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztb, zwf
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_exp')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   ztb, zwf ) 
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_zdf_exp : explicit vertical mixing on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-
-      ! Initializations
-      ! ---------------
-      z1_ksts = 1._wp / REAL( ksts, wp )
-      zwf(:,:, 1 ) = 0._wp    ! no flux at the surface and at bottom level
-      zwf(:,:,jpk) = 0._wp
-      !
-      !
-      DO jn = 1, kjpt         !==  loop over tracers  ==!
-         !
-         ztb(:,:,:) = ptb(:,:,:,jn)    ! initial before value for tracer
-         ! 
-         DO jl = 1, ksts         !==  Split-explicit loop  ==!
-            !              
-            DO jk = 2, jpk             ! 1st vertical derivative (w-flux)
-               DO jj = 2, jpjm1 
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN  ! temperature : use of avt
-                        zwf(ji,jj,jk) =   avt(ji,jj,jk) * ( ztb(ji,jj,jk-1) - ztb(ji,jj,jk) ) / e3w_b(ji,jj,jk)
-                     ELSE                                           ! salinity or pass. tracer : use of avs
-                        zwf(ji,jj,jk) = avs(ji,jj,jk) * ( ztb(ji,jj,jk-1) - ztb(ji,jj,jk) ) / e3w_b(ji,jj,jk)
-                     END IF
-                  END DO
-               END DO
-            END DO
-            !
-            DO jk = 1, jpkm1           ! 2nd vertical derivative   ==> tracer at kt+l*2*rdt/nn_zdfexp
-               DO jj = 2, jpjm1 
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ztb(ji,jj,jk) = ztb(ji,jj,jk) + p2dt * ( zwf(ji,jj,jk) - zwf(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
-                  END DO
-               END DO
-            END DO
-            !
-         END DO                  ! end sub-time stepping
-
-         DO jk = 1, jpkm1        !==  After tracer due to all trends
-            DO jj = 2, jpjm1 
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ze3tb = e3t_b(ji,jj,jk) / e3t_n(ji,jj,jk)
-                  ztra  = ( ztb(ji,jj,jk) - ptb(ji,jj,jk,jn) ) + p2dt * pta(ji,jj,jk,jn)  ! total trend * 2dt 
-                  pta(ji,jj,jk,jn) = ( ze3tb * ptb(ji,jj,jk,jn) + ztra ) * tmask(ji,jj,jk)    ! after tracer
-               END DO
-            END DO
-         END DO
-         !
-      END DO                     ! end of tracer loop
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   ztb, zwf ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_exp')
-      !
-   END SUBROUTINE tra_zdf_exp
-
-   !!==============================================================================
-END MODULE trazdf_exp
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 8143)
+++ 	(revision )
@@ -1,216 +1,0 @@
-MODULE trazdf_imp
-   !!======================================================================
-   !!                 ***  MODULE  trazdf_imp  ***
-   !! Ocean  tracers:  vertical component of the tracer mixing trend
-   !!======================================================================
-   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
-   !!            7.0  !  1991-11  (G. Madec)
-   !!                 !  1992-06  (M. Imbard) correction on tracer trend loops
-   !!                 !  1996-01  (G. Madec) statement function for e3
-   !!                 !  1997-05  (G. Madec) vertical component of isopycnal
-   !!                 !  1997-07  (G. Madec) geopotential diffusion in s-coord
-   !!                 !  2000-08  (G. Madec) double diffusive mixing
-   !!   NEMO     1.0  !  2002-08  (G. Madec) F90: Free form and module
-   !!            2.0  !  2006-11  (G. Madec) New step reorganisation
-   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends
-   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
-   !!             -   !  2011-02  (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition
-   !!            3.7  !  2015-11  (G. Madec, A. Coward)  non linear free surface by default 
-   !!----------------------------------------------------------------------
-  
-   !!----------------------------------------------------------------------
-   !!   tra_zdf_imp   : Update the tracer trend with vertical mixing, nad compute the after tracer field
-   !!----------------------------------------------------------------------
-   USE oce            ! ocean dynamics and tracers variables
-   USE dom_oce        ! ocean space and time domain variables 
-   USE zdf_oce        ! ocean vertical physics variables
-   USE trc_oce        ! share passive tracers/ocean variables
-   USE domvvl         ! variable volume
-   USE ldftra         ! lateral mixing type
-   USE ldfslp         ! lateral physics: slope of diffusion
-   USE zdfddm         ! ocean vertical physics: double diffusion
-   USE traldf_triad   ! active tracers: Method of Stabilizing Correction
-   !
-   USE in_out_manager ! I/O manager
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE lib_mpp        ! MPP library
-   USE wrk_nemo       ! Memory Allocation
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   tra_zdf_imp   !  routine called by step.F90
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
- 
-   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_zdf_imp  ***
-      !!
-      !! ** Purpose :   Compute the after tracer through a implicit computation
-      !!     of the vertical tracer diffusion (including the vertical component 
-      !!     of lateral mixing (only for 2nd order operator, for fourth order 
-      !!     it is already computed and add to the general trend in traldf) 
-      !!
-      !! ** Method  :  The vertical diffusion of a tracer ,t , is given by:
-      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
-      !!      It is computed using a backward time scheme (t=after field)
-      !!      which provide directly the after tracer field.
-      !!      If ln_zdfddm=T, use avs for salinity or for passive tracers
-      !!      Surface and bottom boundary conditions: no diffusive flux on
-      !!      both tracers (bottom, applied through the masked field avt).
-      !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
-      !!
-      !! ** Action  : - pta  becomes the after tracer
-      !!---------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
-      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field
-      !
-      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::  zrhs             ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwi, zwt, zwd, zws
-      !!---------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp')
-      !
-      CALL wrk_alloc( jpi,jpj,jpk,   zwi, zwt, zwd, zws ) 
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp)WRITE(numout,*)
-         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
-         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
-      ENDIF
-      !                                               ! ============= !
-      DO jn = 1, kjpt                                 !  tracer loop  !
-         !                                            ! ============= !
-         !  Matrix construction
-         ! --------------------
-         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
-         !
-         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR.   &
-            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
-            !
-            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
-            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt(:,:,2:jpk)
-            ELSE                                            ;   zwt(:,:,2:jpk) = avs(:,:,2:jpk)
-            ENDIF
-            zwt(:,:,1) = 0._wp
-            !
-            IF( l_ldfslp ) THEN            ! isoneutral diffusion: add the contribution 
-               IF( ln_traldf_msc  ) THEN     ! MSC iso-neutral operator 
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)  
-                        END DO
-                     END DO
-                  END DO
-               ELSE                          ! standard or triad iso-neutral operator
-                  DO jk = 2, jpkm1
-                     DO jj = 2, jpjm1
-                        DO ji = fs_2, fs_jpim1   ! vector opt.
-                           zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
-                        END DO
-                     END DO
-                  END DO
-               ENDIF
-            ENDIF
-            !
-            ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-!!gm BUG  I think, use e3w_a instead of e3w_n
-                     zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
-                     zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
-                     zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
-                 END DO
-               END DO
-            END DO
-            !
-            !! Matrix inversion from the first level
-            !!----------------------------------------------------------------------
-            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
-            !
-            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
-            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
-            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
-            !        (        ...               )( ...  ) ( ...  )
-            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
-            !
-            !   m is decomposed in the product of an upper and lower triangular matrix.
-            !   The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
-            !   Suffices i,s and d indicate "inferior" (below diagonal), diagonal
-            !   and "superior" (above diagonal) components of the tridiagonal system.
-            !   The solution will be in the 4d array pta.
-            !   The 3d array zwt is used as a work space array.
-            !   En route to the solution pta is used a to evaluate the rhs and then 
-            !   used as a work space array: its value is modified.
-            !
-            DO jj = 2, jpjm1        !* 1st recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
-               DO ji = fs_2, fs_jpim1            ! done one for all passive tracers (so included in the IF instruction)
-                  zwt(ji,jj,1) = zwd(ji,jj,1)
-               END DO
-            END DO
-            DO jk = 2, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1
-                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
-                  END DO
-               END DO
-            END DO
-            !
-         ENDIF 
-         !         
-         DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
-            DO ji = fs_2, fs_jpim1
-               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)
-            END DO
-         END DO
-         DO jk = 2, jpkm1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side
-                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
-               END DO
-            END DO
-         END DO
-         !
-         DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
-            DO ji = fs_2, fs_jpim1
-               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
-            END DO
-         END DO
-         DO jk = jpk-2, 1, -1
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
-                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         !                                            ! ================= !
-      END DO                                          !  end tracer loop  !
-      !                                               ! ================= !
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk,   zwi, zwt, zwd, zws ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp')
-      !
-   END SUBROUTINE tra_zdf_imp
-
-   !!==============================================================================
-END MODULE trazdf_imp
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 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 8160)
@@ -37,7 +37,4 @@
    LOGICAL , PUBLIC ::   ln_zdfswm   !: surface  wave-induced mixing flag
    LOGICAL , PUBLIC ::   ln_zdfiwm   !: internal wave-induced mixing flag
-   !                             ! time-stepping
-   LOGICAL , PUBLIC ::   ln_zdfexp   !: explicit vertical diffusion scheme flag
-   INTEGER , PUBLIC ::      nn_zdfexp   !: number of sub-time step (explicit time stepping)
    !                             ! coefficients 
    REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfdrg.F90	(revision 8160)
@@ -190,6 +190,4 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'zdf_drg_init: Choose ONE type of drag coef in namdrg' )
       !
-      IF( ln_drgimp .AND. ln_zdfexp )   CALL ctl_stop( 'zdf_drg_init: ln_drgimp=T requires ln_zdfexp' )
-
       !
       !                     !==  BOTTOM drag setting  ==!   (applied at seafloor)
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90	(revision 8160)
@@ -327,5 +327,5 @@
       zdiag(:,:,2) = zdiag(:,:,2) +  zd_lw(:,:,2) ! Remove zd_lw from zdiag
       zd_lw(:,:,2) = 0._wp
-      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
+      zkar (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) ))
       zflxs(:,:)   = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) &
           &                    * (  ( zhsro(:,:)+gdept_n(:,:,1) ) / zhsro(:,:)  )**(1.5_wp*ra_sf)
@@ -556,5 +556,5 @@
          !
          ! One level below
-         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
+         zkar    (:,:)   = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) )))
          zdep    (:,:)   = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:)
          psi     (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1)
@@ -577,5 +577,5 @@
          !
          ! Set psi vertical flux at the surface:
-         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
+         zkar (:,:)   = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope
          zdep (:,:)   = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)
          zflxs(:,:)   = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp)
@@ -841,7 +841,4 @@
       !!
       !!----------------------------------------------------------------------
-      USE dynzdf_exp
-      USE trazdf_exp
-      !
       INTEGER ::   jk    ! dummy loop indices
       INTEGER ::   ios   ! Local integer output status for namelist read
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8160)
@@ -73,5 +73,4 @@
          &             ln_zdfswm,                                    &     ! surface  wave-induced mixing
          &             ln_zdfiwm,                                    &     ! internal  -      -      -
-         &             ln_zdfexp, nn_zdfexp,                         &     ! time-stepping
          &             rn_avm0, rn_avt0, nn_avb, nn_havtb                  ! coefficients
       !!----------------------------------------------------------------------
@@ -110,7 +109,4 @@
          WRITE(numout,*) '         surface  wave (Qiao et al 2010)         ln_zdfswm = ', ln_zdfswm                                          ! surface wave induced mixing
          WRITE(numout,*) '         internal wave (de Lavergne et al 2017)  ln_zdfiwm = ', ln_zdfiwm
-         WRITE(numout,*) '      time-steping scheme'
-         WRITE(numout,*) '         time splitting (T) / implicit (F)       ln_zdfexp = ', ln_zdfexp
-         WRITE(numout,*) '         number of sub-time step (ln_zdfexp=T)   nn_zdfexp = ', nn_zdfexp
          WRITE(numout,*) '      coefficients : '
          WRITE(numout,*) '         vertical eddy viscosity                 rn_avm0   = ', rn_avm0
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 8160)
@@ -471,5 +471,4 @@
                             CALL tra_adv_init      ! horizontal & vertical advection
                             CALL tra_ldf_init      ! lateral mixing
-                            CALL tra_zdf_init      ! vertical mixing and after tracer fields
 
       !                                      ! Dynamics
@@ -479,5 +478,4 @@
                             CALL dyn_ldf_init      ! lateral mixing
                             CALL dyn_hpg_init      ! horizontal gradient of Hydrostatic pressure
-                            CALL dyn_zdf_init      ! vertical diffusion
                             CALL dyn_spg_init      ! surface pressure gradient
 
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 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8160)
@@ -344,3 +344,4 @@
    END SUBROUTINE stp
    
+   !!======================================================================
 END MODULE step
Index: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
===================================================================
--- branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 8143)
+++ branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90	(revision 8160)
@@ -17,7 +17,8 @@
    USE oce_trc       ! ocean dynamics and active tracers
    USE trd_oce       ! trends: ocean variables
-   USE trazdf_exp    ! vertical diffusion: explicit (tra_zdf_exp     routine)
-   USE trazdf_imp    ! vertical diffusion: implicit (tra_zdf_imp     routine)
+   USE trazdf        ! tracer: vertical diffusion
+!!gm do we really need this ?
    USE trcldf        ! passive tracers: lateral diffusion
+!!gm
    USE trdtra        ! trends manager: tracers 
    USE prtctl_trc    ! Print control
