Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/trazdf_exp.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_exp.F90	(revision 8185)
+++ 	(revision )
@@ -1,153 +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 "zdfddm_substitute.h90"
-#  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 lk_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) = fsavs(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_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 8185)
+++ 	(revision )
@@ -1,217 +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 "zdfddm_substitute.h90"
-#  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 lk_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. lk_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) = fsavs(:,:,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_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/eosbn2.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 8185)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/eosbn2.F90	(revision 8367)
@@ -928,5 +928,5 @@
                pn2(ji,jj,jk) = grav * (  zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) )     &
                   &                    - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  )  &
-                  &            / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  &            / e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
             END DO
          END DO
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/trabbl.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 8185)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/trabbl.F90	(revision 8367)
@@ -13,8 +13,7 @@
    !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
    !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
+   !!            4.0  ! 2017-04  (G. Madec)  ln_trabbl namelist variable instead of a CPP key
    !!----------------------------------------------------------------------
-#if   defined key_trabbl
-   !!----------------------------------------------------------------------
-   !!   'key_trabbl'   or                             bottom boundary layer
+
    !!----------------------------------------------------------------------
    !!   tra_bbl_alloc : allocate trabbl arrays
@@ -49,7 +48,6 @@
    PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
 
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
-
    !                                !!* Namelist nambbl *
+   LOGICAL , PUBLIC ::   ln_trabbl   !: bottom boundary layer flag
    INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
    INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
@@ -82,8 +80,8 @@
       !!                ***  FUNCTION tra_bbl_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
-         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
-         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
-         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
+      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
+         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
+         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
+         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
          !
       IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
@@ -111,5 +109,5 @@
       IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
       !
-      IF( l_trdtra )   THEN                         !* Save the input trends
+      IF( l_trdtra )   THEN                         !* Save the T-S input trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
@@ -301,7 +299,8 @@
             !
          END DO
-         !                                                       ! ===========
-      END DO                                                     ! end tracer
-      !                                                          ! ===========
+         !                                                  ! ===========
+      END DO                                                ! end tracer
+      !                                                     ! ===========
+      !
       IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
       !
@@ -498,6 +497,6 @@
       INTEGER ::   ios                  !   -      -
       REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
-      !
-      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
+      !!
+      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
       !!----------------------------------------------------------------------
       !
@@ -519,17 +518,22 @@
          WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist nambbl : set bbl parameters'
-         WRITE(numout,*) '      diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
-         WRITE(numout,*) '      advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
-         WRITE(numout,*) '      diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
-         WRITE(numout,*) '      advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
-      ENDIF
-
+         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
+         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
+      ENDIF
+      IF( .NOT.ln_trabbl )   RETURN
+      !
+      IF(lwp) THEN
+         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
+         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
+         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
+         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
+      ENDIF
+      !
       !                              ! allocate trabbl arrays
       IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
-
+      !
       IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
       IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
-
+      !
       !                             !* vertical index of  "deep" bottom u- and v-points
       DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
@@ -544,5 +548,5 @@
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
       CALL wrk_dealloc( jpi, jpj, zmbk )
-
+      !
       !                                 !* sign of grad(H) at u- and v-points
       mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
@@ -565,22 +569,8 @@
       ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
       ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
-
       !
       IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
       !
    END SUBROUTINE tra_bbl_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                      No bottom boundary layer scheme
-   !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
-CONTAINS
-   SUBROUTINE tra_bbl_init               ! Dummy routine
-   END SUBROUTINE tra_bbl_init
-   SUBROUTINE tra_bbl( kt )              ! Dummy routine
-      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
-   END SUBROUTINE tra_bbl
-#endif
 
    !!======================================================================
Index: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/trazdf.F90
===================================================================
--- branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 8185)
+++ branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/TRA/trazdf.F90	(revision 8367)
@@ -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,14 +32,11 @@
    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
-#  include "zdfddm_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -59,5 +53,5 @@
       !
       INTEGER  ::   jk                   ! Dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
       !!---------------------------------------------------------------------
       !
@@ -65,19 +59,18 @@
       !
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
-         r2dt =  rdt                          ! = rdt (restarting with Euler time stepping)
+         r2dt =  rdt                                     ! = rdt (restarting with Euler time stepping)
       ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
-         r2dt = 2. * rdt                      ! = 2 rdt (leapfrog)
-      ENDIF
-      !
-      IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
+         r2dt = 2. * rdt                                 ! = 2 rdt (leapfrog)
+      ENDIF
+      !
+      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 {
@@ -98,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)
@@ -108,42 +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.
-      !!----------------------------------------------------------------------
-      USE zdftke
-      USE zdfgls
-      !!----------------------------------------------------------------------
-      !
-      ! Choice from ln_zdfexp already read in namelist in zdfini module
-      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
-      ELSE                   ;   nzdf = 1           ! use implicit scheme
-      ENDIF
-      !
-      ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfgls   )   nzdf = 1   ! TKE, or GLS physics
-      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( ln_zdfexp .AND. nzdf == 1 )   CALL ctl_stop( 'tra_zdf : If using the rotation of lateral mixing operator',   &
-            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' )
-            !
-      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
 
    !!==============================================================================
