Index: trunk/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv.F90	(revision 458)
+++ trunk/NEMO/OPA_SRC/TRA/traadv.F90	(revision 458)
@@ -0,0 +1,199 @@
+MODULE traadv
+   !!==============================================================================
+   !!                       ***  MODULE  traadv  ***
+   !! Ocean active tracers:  advection trend 
+   !!==============================================================================
+   !! History :
+   !!   9.0  !  05-11  (G. Madec)  Original code
+   !!----------------------------------------------------------------------
+   !!   tra_adv      : compute ocean tracer advection trend
+   !!   tra_adv_ctl  : control the different options of advection scheme
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE traadv_cen2     ! 2nd order centered scheme   (tra_adv_cen2 routine)
+   USE traadv_cen2_jki ! 2nd order centered scheme   (tra_adv_cen2 routine)
+   USE traadv_tvd      ! TVD scheme                (tra_adv_tvd    routine)
+   USE traadv_muscl    ! MUSCL scheme              (tra_adv_muscl  routine)
+   USE traadv_muscl2   ! MUSCL2 scheme             (tra_adv_muscl2 routine)
+   USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine)
+   USE trabbl          ! ???
+   USE ldftra_oce      ! ???
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC tra_adv         ! routine called by step module
+ 
+   !! * Share module variables
+   LOGICAL, PUBLIC ::   &
+      ln_traadv_cen2   = .TRUE.  ,   &  ! 2nd order centered scheme flag
+      ln_traadv_tvd    = .FALSE. ,   &  ! TVD scheme flag
+      ln_traadv_muscl  = .FALSE. ,   &  ! MUSCL scheme flag
+      ln_traadv_muscl2 = .FALSE.        ! MUSCL2 scheme flag
+
+   !! * Module variables
+   INTEGER ::   &
+      nadv                              ! choice of the type of advection scheme
+
+   !! * Substitutions
+#   include "domzgr_substitute.h90"
+#   include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tra_adv( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv  ***
+      !!
+      !! ** Purpose :   compute the ocean tracer advection trend.
+      !!
+      !!----------------------------------------------------------------------
+#if ( defined key_trabbl_adv || defined key_traldf_eiv )
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  &  ! temporary arrays
+         &         zun, zvn, zwn
+#else
+      USE oce                , zun => un,  &  ! When no advective bbl, zun == un
+         &                     zvn => vn,  &  !      "        "      , zvn == vn
+         &                     zwn => wn      !      "        "      , zwn == wn
+#endif
+
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 )   CALL tra_adv_ctl                  ! initialisation & control of options
+
+#if defined key_trabbl_adv
+      ! Advective bottom boundary layer                      ! add the bbl velocity
+      ! -------------------------------
+      zun(:,:,:) = un(:,:,:) - u_bbl(:,:,:)
+      zvn(:,:,:) = vn(:,:,:) - v_bbl(:,:,:)
+      zwn(:,:,:) = wn(:,:,:) + w_bbl(:,:,:)
+#endif
+      IF( lk_traldf_eiv ) THEN                               ! add the eiv velocity
+         IF( .NOT. lk_trabbl_adv ) THEN 
+            zun(:,:,:) = un(:,:,:)
+            zvn(:,:,:) = vn(:,:,:)
+            zwn(:,:,:) = wn(:,:,:)
+         ENDIF
+         CALL tra_adv_eiv( kt, zun, zvn, zwn )                  ! compute and add the eiv velocity 
+      ENDIF
+
+      SELECT CASE ( nadv )                       ! compute advection trend and add it to general trend
+      CASE ( -1 )                                       ! esopa: test all possibility with control print
+         CALL tra_adv_cen2    ( kt, zun, zvn, zwn )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf0 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_adv_cen2_jki( kt, zun, zvn, zwn )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf1 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_adv_tvd     ( kt, zun, zvn, zwn )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf2 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_adv_muscl   ( kt, zun, zvn, zwn )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf3 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_adv_muscl2  ( kt, zun, zvn, zwn )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf4 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+      CASE ( 0 )                                       ! 2nd order centered scheme k-j-i loops
+         CALL tra_adv_cen2    ( kt, zun, zvn, zwn )
+      CASE ( 1 )                                       ! 2nd order centered scheme
+         CALL tra_adv_cen2_jki( kt, zun, zvn, zwn )
+      CASE ( 2 )                                       ! TVD scheme
+         CALL tra_adv_tvd    ( kt, zun, zvn, zwn )
+      CASE ( 3 )                                       ! MUSCL  scheme
+         CALL tra_adv_muscl  ( kt, zun, zvn, zwn )
+      CASE ( 4 )                                       ! MUSCL2 scheme
+         CALL tra_adv_muscl2 ( kt, zun, zvn, zwn )
+      END SELECT
+
+      !                                                ! print mean trends (used for debugging)
+!      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' adv  - Ta: ', mask1=tmask,               &
+!         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+   END SUBROUTINE tra_adv
+
+
+   SUBROUTINE tra_adv_ctl
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_ctl  ***
+      !!                
+      !! ** Purpose :   Control the consistency between cpp options for 
+      !!      tracer advection schemes
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio
+      NAMELIST/nam_traadv/ ln_traadv_cen2 , ln_traadv_tvd,   &
+         &                 ln_traadv_muscl, ln_traadv_muscl2
+      !!----------------------------------------------------------------------
+
+      ! Read Namelist nam_traadv : tracer advection scheme
+      ! -------------------------
+      REWIND ( numnam )
+      READ   ( numnam, nam_traadv )
+
+      ! Parameter control and print
+      ! ---------------------------
+      ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_adv_ctl : choice/control of the tracer advection scheme'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '          Namelist nam_tra_adv : chose a advection scheme for tracers'
+         WRITE(numout,*)
+         WRITE(numout,*) '             2nd order advection scheme     ln_traadv_cen2   = ', ln_traadv_cen2
+         WRITE(numout,*) '             TVD advection scheme           ln_traadv_tvd    = ', ln_traadv_tvd
+         WRITE(numout,*) '             MUSCL  advection scheme        ln_traadv_muscl  = ', ln_traadv_muscl
+         WRITE(numout,*) '             MUSCL2 advection scheme        ln_traadv_muscl2 = ', ln_traadv_muscl2
+      ENDIF
+
+      ! Control of Advection scheme options
+      ! -----------------------------------
+      ioptio = 0
+      IF( ln_traadv_cen2   )   ioptio = ioptio + 1
+      IF( ln_traadv_tvd    )   ioptio = ioptio + 1
+      IF( ln_traadv_muscl  )   ioptio = ioptio + 1
+      IF( ln_traadv_muscl2 )   ioptio = ioptio + 1
+
+      IF( .NOT.lk_esopa .AND. ( ioptio > 1 .OR. ioptio == 0 ) ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) ' Choose ONE advection scheme in namelist nam_traadv'
+         nstop = nstop + 1
+      ENDIF
+
+      IF( n_cla == 1 .AND. .NOT. ln_traadv_cen2 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '     cross-land advection only with 2nd order advection scheme'
+         nstop = nstop + 1
+      ENDIF
+
+      ! Set nadv
+      ! --------
+      IF( ln_traadv_cen2   )   nadv = 0
+#if defined key_mpp_omp
+      IF( ln_traadv_cen2   )   nadv = 1
+#endif
+      IF( ln_traadv_tvd    )   nadv = 2
+      IF( ln_traadv_muscl  )   nadv = 3
+      IF( ln_traadv_muscl2 )   nadv = 4
+
+      IF( lk_esopa ) THEN
+         IF(lwp) WRITE(numout,*) '          esopa control: use all lateral physics options'
+         nadv = -1
+      ENDIF
+      IF(lwp) WRITE(numout,*) '             choice of tra_adv_...                nadv   = ', nadv
+
+   END SUBROUTINE tra_adv_ctl
+
+  !!======================================================================
+END MODULE traadv
Index: trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90	(revision 458)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90	(revision 458)
@@ -0,0 +1,361 @@
+MODULE traadv_cen2_jki
+   !!==============================================================================
+   !!                       ***  MODULE  traadv_cen2_jki  ***
+   !! Ocean active tracers:  horizontal & vertical advective trend
+   !!==============================================================================
+   !! History :
+   !!   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv
+   !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
+   !!   9.0  !  04-08  (C. Talandier) New trends organization
+   !!    "   !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!    "   !  06-04  (R. Benshila, G. Madec) Step reorganization
+   !!----------------------------------------------------------------------
+   !!   tra_adv_cen2_jki : update the tracer trend with the horizontal and
+   !!                  vertical advection trends using a seconder order
+   !!                  centered scheme. Auto-tasking case, k-slab for
+   !!                  hor. adv., j-slab for vert. adv. (j-k-i loops)
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE trdmod          ! ocean active tracers trends
+   USE trdmod_oce      ! ocean variables trends
+   USE flxrnf          !
+   USE trabbl          ! advective term in the BBL
+   USE ocfzpt          !
+   USE lib_mpp
+   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+   USE in_out_manager  ! I/O manager
+   USE diaptr          ! poleward transport diagnostics
+   USE dynspg_oce      ! choice/control of key cpp for surface pressure gradient
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC tra_adv_cen2_jki    ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE tra_adv_cen2_jki( kt, pun, pvn, pwn )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE tra_adv_cen2_jki  ***
+      !!              
+      !! ** Purpose :   Compute the now trend due to the advection of tracers
+      !!      and add it to the general trend of passive tracer equations.
+      !!
+      !! ** Method  :   The advection is evaluated by a second order centered
+      !!      scheme using now fields (leap-frog scheme). In specific areas
+      !!      (vicinity of major river mouths, some straits, or where tn is
+      !!      approaching the freezing point) it is mixed with an upstream 
+      !!      scheme for stability reasons.
+      !!         Part 0 : compute the upstream / centered flag
+      !!                  (3D array, zind, defined at T-point (0<zind<1))
+      !!         Part I : horizontal advection
+      !!       * centered flux:
+      !!               zcenu = e2u*e3u  un  mi(tn)
+      !!               zcenv = e1v*e3v  vn  mj(tn)
+      !!       * upstream flux:
+      !!               zupsu = e2u*e3u  un  (tb(i) or tb(i-1) ) [un>0 or <0]
+      !!               zupsv = e1v*e3v  vn  (tb(j) or tb(j-1) ) [vn>0 or <0]
+      !!       * mixed upstream / centered horizontal advection scheme
+      !!               zcofi = max(zind(i+1), zind(i))
+      !!               zcofj = max(zind(j+1), zind(j))
+      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
+      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
+      !!       * horizontal advective trend (divergence of the fluxes)
+      !!               zta = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
+      !!       * Add this trend now to the general trend of tracer (ta,sa):
+      !!              (ta,sa) = (ta,sa) + ( zta , zsa )
+      !!       * trend diagnostic ('key_trdtra' defined): the trend is
+      !!      saved for diagnostics. The trends saved is expressed as
+      !!      Uh.gradh(T), i.e.
+      !!                     save trend = zta + tn divn
+      !!         In addition, the advective trend in the two horizontal direc-
+      !!      tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is
+      !!      equal to (in s-coordinates, and similarly in z-coord.):
+      !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] )
+      !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  }
+      !!         C A U T I O N : the trend saved is the centered trend only.
+      !!      It doesn't take into account the upstream part of the scheme.
+      !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so 
+      !!      they vanish from the expression of the flux and divergence.
+      !!
+      !!         Part II : vertical advection
+      !!      For temperature (idem for salinity) the advective trend is com-
+      !!      puted as follows :
+      !!            zta = 1/e3t dk+1[ zwz ]
+      !!      where the vertical advective flux, zwz, is given by :
+      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
+      !!      with 
+      !!        zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
+      !!        zcenu = centered flux = wn * mk(tn)
+      !!         The surface boundary condition is : 
+      !!      rigid-lid (lk_dynspg_frd = T) : zero advective flux
+      !!      free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1)
+      !!         Add this trend now to the general trend of tracer (ta,sa):
+      !!            (ta,sa) = (ta,sa) + ( zta , zsa )
+      !!         Trend diagnostic ('key_trdtra' defined): the trend is
+      !!      saved for diagnostics. The trends saved is expressed as :
+      !!             save trend =  w.gradz(T) = zta - tn divn.
+      !!
+      !! ** Action :  - update (ta,sa) with the now advective tracer trends
+      !!              - save trends in (ttrdh,ttrd,strdhi,strd) ('key_trdtra')
+      !!
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      USE oce              , zwx => ua,    &  ! use ua as workspace
+         &                   zwy => va        ! use va as workspace
+
+      !! * Arguments
+      INTEGER , INTENT( in ) ::   kt          ! ocean time-step index
+      REAL(wp), INTENT( in ), DIMENSION(jpi,jpj,jpk) ::   &
+         pun, pvn, pwn                        ! now ocean velocity fields
+
+      !! * Local save
+      REAL(wp), DIMENSION(jpi,jpj), SAVE ::   &
+         zbtr2
+
+      !! * Local declarations
+      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
+      REAL(wp) ::                           &
+         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars
+         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
+         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         "
+         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
+         zcofk, zupst, zupss, zcent,        &  !    "         "
+         zcens, zfp_w, zfm_w,               &  !    "         "
+         zcenut, zcenvt, zcenus, zcenvs        !    "         "
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         zwz, zww, zind,                    &  ! temporary workspace arrays
+         ztdta, ztdsa                          !    "         "
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   Auto-tasking case'
+         IF(lwp) WRITE(numout,*)
+         
+         zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
+      ENDIF
+
+      ! Save ta and sa trends
+      IF( l_trdtra )   THEN
+         ztdta(:,:,:) = ta(:,:,:) 
+         ztdsa(:,:,:) = sa(:,:,:) 
+         l_adv = 'ce2'
+      ENDIF
+
+!CDIR PARALLEL DO
+!$OMP PARALLEL DO
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+
+         ! 0. Upstream / centered scheme indicator
+         ! ---------------------------------------
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zind(ji,jj,jk) =  MAX (              &
+                  upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
+                  upsadv(ji,jj)                     &  ! in the vicinity of some straits
+#if defined key_ice_lim
+                  , tmask(ji,jj,jk)                 &  ! half upstream tracer fluxes if tn < ("freezing"+0.1 )
+                  * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) )   &
+#endif
+               )
+            END DO
+         END DO
+
+
+         ! I. Horizontal advective fluxes
+         ! ------------------------------
+         ! Second order centered tracer flux at u and v-points
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               ! upstream indicator
+               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
+               zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
+               ! volume fluxes * 1/2
+#if defined key_zco
+               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
+               zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
+#else
+               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
+               zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
+#endif
+               ! upstream scheme
+               zfp_ui = zfui + ABS( zfui )
+               zfp_vj = zfvj + ABS( zfvj )
+               zfm_ui = zfui - ABS( zfui )
+               zfm_vj = zfvj - ABS( zfvj )
+               zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj  ,jk)
+               zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji  ,jj+1,jk)
+               zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
+               zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
+               ! centered scheme
+               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
+               zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
+               zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
+               zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
+               ! mixed centered / upstream scheme
+               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
+               zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
+               zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
+               zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
+            END DO
+         END DO
+
+         ! 2. Tracer flux divergence at t-point added to the general trend
+         ! ---------------------------------------------------------------
+
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+#if defined key_zco
+               zbtr = zbtr2(ji,jj) 
+#else
+               zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
+#endif
+               ! horizontal advective trends
+               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+                  &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
+               zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
+                  &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
+               ! add it to the general tracer trends
+               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
+               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
+            END DO
+         END DO
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+      ! 3. Save the horizontal advective trends for diagnostic
+      ! ------------------------------------------------------
+      IF( l_trdtra )   THEN
+         ! Recompute the hoizontal advection zta & zsa trends computed 
+         ! at the step 2. above in making the difference between the new 
+         ! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
+         ! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
+         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:) 
+         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
+
+         CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
+
+         ! Save the new ta and sa trends
+         ztdta(:,:,:) = ta(:,:,:) 
+         ztdsa(:,:,:) = sa(:,:,:) 
+
+      ENDIF
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 had  - Ta: ', mask1=tmask, &
+         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+      ! "zonal" mean advective heat and salt transport
+      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         IF( lk_zco ) THEN
+            DO jk = 1, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                    zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
+                    zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         ENDIF
+         pht_adv(:) = ptr_vj( zwy(:,:,:) )
+         pst_adv(:) = ptr_vj( zwz(:,:,:) )
+      ENDIF
+
+      ! II. Vertical advection
+      ! ----------------------
+!CDIR PARALLEL DO
+!$OMP PARALLEL DO
+      !                                                ! ===============
+      DO jj = 2, jpjm1                                 !  Vertical slab
+         !                                             ! ===============
+         ! Bottom value : flux and indicator set to zero
+         zwz (:,jj,jpk) = 0.e0     ;    zww(:,jj,jpk) = 0.e0
+         zind(:,jj,jpk) = 0.e0
+
+         ! Surface value
+         IF( lk_dynspg_rl ) THEN                          ! rigid lid : flux set to zero
+            zwz(:,jj, 1 ) = 0.e0    ;    zww(:,jj, 1 ) = 0.e0
+         ELSE                                             ! free surface-constant volume
+            zwz(:,jj, 1 ) = pwn(:,jj,1) * tn(:,jj,1)
+            zww(:,jj, 1 ) = pwn(:,jj,1) * sn(:,jj,1)
+         ENDIF
+         
+         ! 1. Vertical advective fluxes
+         ! ----------------------------
+         ! Second order centered tracer flux at w-point
+         DO jk = 2, jpk
+            DO ji = 2, jpim1
+               ! upstream indicator
+               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
+               ! velocity * 1/2
+               zhw = 0.5 * pwn(ji,jj,jk)
+               ! upstream scheme
+               zfp_w = zhw + ABS( zhw )
+               zfm_w = zhw - ABS( zhw )
+               zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
+               zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
+               ! centered scheme
+               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
+               zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
+               ! mixed centered / upstream scheme
+               zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
+               zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
+            END DO
+         END DO
+
+         ! 2. Tracer flux divergence at t-point added to the general trend
+         ! -------------------------
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               ze3tr = 1. / fse3t(ji,jj,jk)
+               ! vertical advective trends
+               zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
+               zsa = - ze3tr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) )
+               ! add it to the general tracer trends
+               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
+               sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
+            END DO
+         END DO
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+      ! 3. Save the vertical advective trends for diagnostic
+      ! ----------------------------------------------------
+      IF( l_trdtra )   THEN
+         ! Recompute the vertical advection zta & zsa trends computed 
+         ! at the step 2. above in making the difference between the new 
+         ! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
+         ! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
+         ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:) 
+         ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
+
+         CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
+      ENDIF
+
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' cen2 zad  - Ta: ', mask1=tmask, &
+         &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+   END SUBROUTINE tra_adv_cen2_jki
+
+   !!======================================================================
+END MODULE traadv_cen2_jki
Index: trunk/NEMO/OPA_SRC/TRA/traadv_eiv.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 458)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 458)
@@ -0,0 +1,156 @@
+MODULE traadv_eiv
+   !!==============================================================================
+   !!                    ***  MODULE  traadv_eiv  ***
+   !! Ocean active tracers:  advection trend - eddy induced velocity
+   !!==============================================================================
+   !! History :
+   !!   9.0  !  05-11  (G. Madec)  Original code from traldf & zdf _iso
+   !!----------------------------------------------------------------------
+#if defined key_traldf_eiv   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_traldf_eiv'                  rotation of the lateral mixing tensor
+   !!----------------------------------------------------------------------
+   !!   tra_ldf_iso : update the tracer trend with the horizontal component
+   !!                 of iso neutral laplacian operator or horizontal 
+   !!                 laplacian operator in s-coordinate
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE ldfslp          ! iso-neutral slopes
+   USE in_out_manager  ! I/O manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC tra_adv_eiv  ! routine called by step.F90
+
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "ldftra_substitute.h90"
+#  include "ldfeiv_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_adv_eiv  ***
+      !! 
+      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
+      !!      trend and add it to the general trend of tracer equation.
+      !!
+      !! ** Method  :   The eddy induced advection is computed from the slope
+      !!      of iso-neutral surfaces computed in routine ldf_slp as follows:
+      !!         zu_eiv =  1/(e2u e3u)   dk[ aeiu e2u mi(wslpi) ]
+      !!         zv_eiv =  1/(e1v e3v)   dk[ aeiv e1v mj(wslpj)
+      !!         zw_eiv = -1/(e1t e2t) { di[ aeiu e2u mi(wslpi) ]
+      !!                               + dj[ aeiv e1v mj(wslpj) ] }
+      !!      add the eiv component to the model velocity:
+      !!         p.n = p.n + z._eiv
+      !!
+      !! ** Action  : - add to p.n the eiv component
+      !!
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER , INTENT( in ) ::   kt       ! ocean time-step index
+      REAL(wp), INTENT( inout ), DIMENSION(jpi,jpj,jpk) ::   &
+         pun, pvn, pwn                     ! in : now ocean velocity fields 
+         !                                 ! out: fields increased by the eiv components
+
+      !! * Local declarations
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      REAL(wp) ::                       &
+         zuwk, zuwk1, zuwi, zuwi1,      &  ! temporary scalar
+         zvwk, zvwk1, zvwj, zvwj1,      &  !
+         zu_eiv, zv_eiv, zw_eiv            !
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection :'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
+# if defined key_diaeiv
+         u_eiv(:,:,:) = 0.e0
+         v_eiv(:,:,:) = 0.e0
+         w_eiv(:,:,:) = 0.e0
+# endif
+      ENDIF
+
+      !                                             ! =================
+      DO jk = 1, jpkm1                              !  Horizontal slab
+         !                                          ! =================
+
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zuwk = ( wslpi(ji,jj,jk  ) + wslpi(ji+1,jj,jk  ) ) * fsaeiu(ji,jj,jk  ) * umask(ji,jj,jk  )
+               zuwk1= ( wslpi(ji,jj,jk+1) + wslpi(ji+1,jj,jk+1) ) * fsaeiu(ji,jj,jk+1) * umask(ji,jj,jk+1)
+               zvwk = ( wslpj(ji,jj,jk  ) + wslpj(ji,jj+1,jk  ) ) * fsaeiv(ji,jj,jk  ) * vmask(ji,jj,jk  )
+               zvwk1= ( wslpj(ji,jj,jk+1) + wslpj(ji,jj+1,jk+1) ) * fsaeiv(ji,jj,jk+1) * vmask(ji,jj,jk+1)
+
+               zu_eiv =  0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) / fse3u(ji,jj,jk)
+               zv_eiv =  0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) / fse3v(ji,jj,jk)
+   
+               pun(ji,jj,jk) = pun(ji,jj,jk) + zu_eiv
+               pvn(ji,jj,jk) = pvn(ji,jj,jk) + zv_eiv
+
+# if defined key_diaeiv
+               u_eiv(ji,jj,jk) = zu_eiv
+               v_eiv(ji,jj,jk) = zv_eiv
+# endif
+            END DO
+         END DO
+         IF( jk >=2 ) THEN                             ! jk=1 zw_eiv=0, not computed
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+# if defined key_traldf_c2d || defined key_traldf_c3d
+                  zuwi  = ( wslpi(ji,jj,jk)+wslpi(ji-1,jj,jk) ) * fsaeiu(ji-1,jj,jk) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
+                  zuwi1 = ( wslpi(ji,jj,jk)+wslpi(ji+1,jj,jk) ) * fsaeiu(ji  ,jj,jk) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
+                  zvwj  = ( wslpj(ji,jj,jk)+wslpj(ji,jj-1,jk) ) * fsaeiv(ji,jj-1,jk) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
+                  zvwj1 = ( wslpj(ji,jj,jk)+wslpj(ji,jj+1,jk) ) * fsaeiv(ji,jj  ,jk) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
+  
+                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) / ( e1t(ji,jj)*e2t(ji,jj) )
+# else
+                  zuwki = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
+                  zuwk  = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
+                  zvwki = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
+                  zvwk  = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
+
+                  zw_eiv = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwk - zuwki + zvwk - zvwki )
+                     &                                                / ( e1t(ji,jj)*e2t(ji,jj) )
+# endif
+                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv
+
+# if defined key_diaeiv
+                  w_eiv(ji,jj,jk) = zw_eiv
+# endif
+               END DO
+            END DO
+         ENDIF
+         !                                          ! =================
+      END DO                                        !    End of slab  
+      !                                             ! =================
+
+   END SUBROUTINE tra_adv_eiv
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :             No rotation of the lateral mixing tensor
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE tra_adv_eiv( kt, pun, pvn, pwn )              ! Empty routine
+      INTEGER ::   kt 
+      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
+      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?',   &
+         &       kt, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
+   END SUBROUTINE tra_adv_eiv
+#endif
+
+   !!==============================================================================
+END MODULE traadv_eiv
Index: trunk/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traldf.F90	(revision 458)
+++ trunk/NEMO/OPA_SRC/TRA/traldf.F90	(revision 458)
@@ -0,0 +1,375 @@
+MODULE traldf
+   !!======================================================================
+   !!                       ***  MODULE  traldf  ***
+   !! Ocean Active tracers : lateral diffusive trends 
+   !!=====================================================================
+   !! History :
+   !!   9.0  !  05-11  (G. Madec)  Original code
+   !!----------------------------------------------------------------------
+   !!   tra_ldf     : update the tracer trend with the lateral diffusion
+   !!       ldf_ctl : initialization, namelist read, and parameters control
+   !!       ldf_ano : compute lateral diffusion for constant T-S profiles
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE ldftra_oce      ! ocean tracer   lateral physics
+   USE ldfslp          ! ???
+   USE traldf_bilapg   ! lateral mixing            (tra_ldf_bilapg routine)
+   USE traldf_bilap    ! lateral mixing             (tra_ldf_bilap routine)
+   USE traldf_iso      ! lateral mixing               (tra_ldf_iso routine)
+   USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine)
+   USE trdmod          ! ocean active tracers trends
+   USE trdmod_oce      ! ocean variables trends
+   USE prtctl          ! Print control
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! *  Routine accessibility
+   PUBLIC tra_ldf      ! called by step.F90 
+
+   !! * module variables
+   INTEGER ::                        & 
+      nldf = 0                         ! type of lateral diffusion used
+      !                                ! defined from ln_traldf_... namlist logicals)
+#if defined key_traldf_ano
+   REAL, DIMENSION(jpi,jpj,jpk) ::   &  
+      t0_ldf, s0_ldf                    ! lateral diffusion trends of T & S
+      !                                 ! for a constant vertical profile
+#endif
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!---------------------------------------------------------------------------------
+   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tra_ldf( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_ldf  ***
+      !! 
+      !! ** Purpose :   compute the lateral ocean tracer physics.
+      !!
+      !!----------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
+
+      !! * local declarations
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         ztrdt, ztrds                         ! 3D temporary workspace
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 )   CALL ldf_ctl      ! initialisation & control of options
+
+      IF( l_trdtra )   THEN                      ! temporary save of ta and sa trends
+         ztrdt(:,:,:) = ta(:,:,:) 
+         ztrds(:,:,:) = sa(:,:,:) 
+      ENDIF
+
+      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
+      CASE ( -1 )                                       ! esopa: test all possibility with control print
+         CALL tra_ldf_lap    ( kt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf0 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_ldf_iso    ( kt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf1 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_ldf_bilap  ( kt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf2 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_ldf_bilapg ( kt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf3 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+      CASE ( 0 )                                       ! iso-level laplacian
+         CALL tra_ldf_lap    ( kt )
+      CASE ( 1 )                                       ! rotated laplacian (except dk[ dk[.] ] part)
+         CALL tra_ldf_iso    ( kt )
+      CASE ( 2 )                                       ! iso-level bilaplacian
+         CALL tra_ldf_bilap  ( kt )
+      CASE ( 3 )                                       ! s-coord. horizontal bilaplacian
+         CALL tra_ldf_bilapg ( kt )
+      END SELECT
+
+#if defined key_traldf_ano
+      !                                          ! anomaly: substract the reference diffusivity
+      ta(:,:,:) = ta(:,:,:) - t0_ldf(:,:,:)
+      sa(:,:,:) = sa(:,:,:) - s0_ldf(:,:,:)
+#endif
+
+      !                                          ! save the horizontal diffusive trends for further diagnostics
+      IF( l_trdtra )   THEN
+         ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
+         ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
+         CALL trd_mod( ztrdt, ztrds, jpttdldf, 'TRA', kt )
+      ENDIF
+
+      !                                          ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' ldf  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+   END SUBROUTINE tra_ldf
+
+
+   SUBROUTINE ldf_ctl
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_ctl  ***
+      !! 
+      !! ** Purpose :   Choice of the operator for the lateral tracer diffusion
+      !!
+      !! ** Method  :   set nldf from the nam_traldf logicals
+      !!      nldf == -1   ESOPA test: ALL operators are used
+      !!      nldf ==  0   laplacian operator
+      !!      nldf ==  1   Rotated laplacian operator
+      !!      nldf ==  2   bilaplacian operator
+      !!      nldf ==  3   Rotated bilaplacian
+      !!
+      !!----------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   ioptio, ierr         ! temporary integers 
+!      
+!     NAMELIST/nam_traldf/ ln_traldf_lap  , ln_traldf_bilap,                &
+!        &                 ln_traldf_level, ln_traldf_hor, ln_traldf_iso,   &
+!        &                 aht0, ahtb0, aeiv0
+      !!----------------------------------------------------------------------
+
+      !  Define the lateral mixing oparator for tracers
+      ! ===============================================
+    
+      ! Namelist nam_traldf already read in ldftra module
+!     ! Read Namelist nam_traldf : Lateral physics on tracers
+!     REWIND( numnam )
+!     READ  ( numnam, nam_traldf )
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra:ldf_ctl : lateral tracer diffusive operator'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         WRITE(numout,*) '          Namelist nam_traldf : set lateral mixing parameters (type, direction, coefficients)'
+         WRITE(numout,*) '             laplacian operator          ln_traldf_lap   = ', ln_traldf_lap
+         WRITE(numout,*) '             bilaplacian operator        ln_traldf_bilap = ', ln_traldf_bilap
+         WRITE(numout,*) '             iso-level                   ln_traldf_level = ', ln_traldf_level
+         WRITE(numout,*) '             horizontal (geopotential)   ln_traldf_hor   = ', ln_traldf_hor
+         WRITE(numout,*) '             iso-neutral                 ln_traldf_iso   = ', ln_traldf_iso
+      ENDIF
+
+      ! Parameter control
+
+      ! control the input
+      ioptio = 0
+      IF( ln_traldf_lap   )   ioptio = ioptio + 1
+      IF( ln_traldf_bilap )   ioptio = ioptio + 1
+      IF( ioptio /= 1 )   THEN
+          IF(lwp) WRITE(numout,cform_err)
+          IF(lwp) WRITE(numout,*) '          use ONE of the 2 lap/bilap operator type on tracer'
+          nstop = nstop + 1
+      ENDIF
+      ioptio = 0
+      IF( ln_traldf_level )   ioptio = ioptio + 1
+      IF( ln_traldf_hor   )   ioptio = ioptio + 1
+      IF( ln_traldf_iso   )   ioptio = ioptio + 1
+      IF( ioptio /= 1 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          use only ONE direction (level/hor/iso)'
+         nstop = nstop + 1
+      ENDIF
+
+      ! defined the type of lateral diffusion from ln_traldf_... logicals
+      ierr = 0
+      IF ( ln_traldf_lap ) THEN      ! laplacian operator
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation)
+            IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed
+            IF ( ln_traldf_hor   )   nldf = 0      ! horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! z-coordinate
+            IF ( ln_traldf_level )   nldf = 0      ! iso-level  (no rotation)
+            IF ( ln_traldf_hor   )   nldf = 1      ! horizontal (   rotation)
+            IF ( ln_traldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
+      IF( ln_traldf_bilap ) THEN      ! bilaplacian operator
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation)
+            IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation)
+            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_traldf_level )   ierr = 1      ! iso-level not allowed 
+            IF ( ln_traldf_hor   )   nldf = 2      ! horizontal (no rotation)
+            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! z-coordinate
+            IF ( ln_traldf_level )   nldf = 2      ! iso-level  (no rotation)
+            IF ( ln_traldf_hor   )   nldf = 3      ! horizontal (   rotation)
+            IF ( ln_traldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
+      IF( ierr == 1 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) ' iso-level in z-coordinate - partial step, not allowed'
+         nstop = nstop + 1
+      ENDIF
+      IF( ierr == 2 ) THEN
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) ' isoneutral bilaplacian operator does not exist'
+         nstop = nstop + 1
+      ENDIF
+      IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso ) THEN
+         IF(lwp) WRITE(numout,*) '          eddy induced velocity on tracers'
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) ' the eddy induced velocity on tracers requires isopycnal laplacian diffusion'
+         nstop = nstop + 1
+      ENDIF
+      IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
+         IF( .NOT.lk_ldfslp ) THEN
+            IF(lwp) WRITE(numout,cform_err)
+            IF(lwp) WRITE(numout,*) '          the rotation of the diffusive tensor require key_ldfslp'
+            nstop = nstop + 1
+         ENDIF
+      ENDIF
+
+      IF( lk_esopa ) THEN
+         IF(lwp) WRITE(numout,*) '          esopa control: use all lateral physics options'
+         nldf = -1
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF( nldf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
+         IF( nldf ==  0 )   WRITE(numout,*) '              laplacian operator'
+         IF( nldf ==  1 )   WRITE(numout,*) '              Rotated laplacian operator'
+         IF( nldf ==  2 )   WRITE(numout,*) '              bilaplacian operator'
+         IF( nldf ==  3 )   WRITE(numout,*) '              Rotated bilaplacian'
+      ENDIF
+
+
+      ! Reference T & S diffusivity (if necessary)
+      ! ===========================
+
+      CALL ldf_ano
+
+   END SUBROUTINE ldf_ctl
+
+#if defined key_traldf_ano
+   !!----------------------------------------------------------------------
+   !!   'key_traldf_ano'               T & S lateral diffusion on anomalies
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE ldf_ano
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_ano  ***
+      !!
+      !! ** Purpose :   initializations of 
+      !!
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      USE zdf_oce         ! vertical mixing
+      USE trazdf          ! vertical mixing: double diffusion
+      USE zdfddm          ! vertical mixing: double diffusion
+
+      !! * local declarations
+      INTEGER  ::   jk                         ! Dummy loop indice
+      LOGICAL  ::   llsave                     !
+      REAL(wp) ::   zt0, zs0, z12              ! temporary scalar
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         zt_ref, ztb, zavt,                 &  ! 3D temporary workspace
+         zs_ref, zsb                           ! 3D temporary workspace
+      !!----------------------------------------------------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on anomalies'
+         WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+
+      ! defined the T & S reference profiles
+      ! ------------------------------------
+      zt0 =10.e0                               ! homogeneous ocean
+      zs0 =35.e0
+      zt_ref(:,:,:) = 10.0 * tmask(:,:,:)
+      zs_ref(:,:,:) = 35.0 * tmask(:,:,:)
+      IF(lwp) WRITE(numout,*) '              homogeneous ocean T = ', zt0, ' S = ',zs0
+
+      !                                        ! T & S profile (to be coded +namelist parameter
+
+      ! prepare the ldf computation
+      ! ---------------------------
+      llsave = l_trdtra
+      l_trdtra = .false.      ! desactivate trend computation
+      t0_ldf(:,:,:) = 0.e0
+      s0_ldf(:,:,:) = 0.e0
+      ztb   (:,:,:) = tb (:,:,:)
+      zsb   (:,:,:) = sb (:,:,:)
+      ua    (:,:,:) = ta (:,:,:)
+      va    (:,:,:) = sa (:,:,:)
+      zavt  (:,:,:) = avt(:,:,:)
+      IF( lk_zdfddm ) THEN
+          IF(lwp) WRITE(numout,cform_err)
+          IF(lwp) WRITE(numout,*) ' key_traldf_ano with key_zdfddm not implemented'
+          nstop = nstop + 1
+      ENDIF
+      ! set tb, sb to reference values and avr to zero
+      tb (:,:,:) = zt_ref(:,:,:)
+      sb (:,:,:) = zs_ref(:,:,:)
+      ta (:,:,:) = 0.e0
+      sa (:,:,:) = 0.e0
+      avt(:,:,:) = 0.e0
+
+      ! Compute the ldf trends
+      ! ----------------------
+      CALL tra_ldf( nit000+1 )      ! horizontal components (+1: no more init)
+      CALL tra_zdf( nit000   )      ! vertical component (if necessary nit000 to performed the init)
+
+      ! finalise the computation and recover all arrays
+      ! -----------------------------------------------
+      l_trdtra = llsave
+      z12 = 2.e0
+      IF( neuler == 1)   z12 = 1.e0
+      IF( ln_zdfexp ) THEN      ! ta,sa are the trends
+         t0_ldf(:,:,:) = ta(:,:,:)
+         s0_ldf(:,:,:) = sa(:,:,:)
+      ELSE
+         DO jk = 1, jpkm1
+            t0_ldf(:,:,jk) = ( ta(:,:,jk) - tb(:,:,jk) ) / ( z12 *rdttra(jk) )
+            s0_ldf(:,:,jk) = ( sa(:,:,jk) - tb(:,:,jk) ) / ( z12 *rdttra(jk) )
+         END DO
+      ENDIF
+      tb    (:,:,:) = ztb (:,:,:)
+      sb    (:,:,:) = zsb (:,:,:)
+      ta    (:,:,:) = ua  (:,:,:)
+      sa    (:,:,:) = va  (:,:,:)
+      avt   (:,:,:) = zavt(:,:,:)
+
+   END SUBROUTINE ldf_ano
+
+#else
+   !!----------------------------------------------------------------------
+   !!   default option :   Dummy code   NO T & S background profiles
+   !!----------------------------------------------------------------------
+   SUBROUTINE ldf_ano
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on the full fields'
+         WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+   END SUBROUTINE ldf_ano
+#endif
+
+   !!======================================================================
+END MODULE traldf
Index: trunk/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 458)
+++ trunk/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 458)
@@ -0,0 +1,191 @@
+MODULE trazdf
+   !!==============================================================================
+   !!                 ***  MODULE  trazdf  ***
+   !! Ocean active tracers:  vertical component of the tracer mixing trend
+   !!==============================================================================
+   !! History :
+   !!   9.0  !  05-11  (G. Madec)  Original code
+   !!----------------------------------------------------------------------
+   !!   tra_zdf      : Update the tracer trend with the vertical diffusion
+   !!       zdf_ctl  : ???
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE zdf_oce         ! ocean vertical physics variables
+
+   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine)
+   USE trazdf_imp      ! vertical diffusion: implicit (tra_zdf_imp     routine)
+   USE trazdf_imp_jki  ! vertical diffusion  implicit (tra_zdf_imp_jki routine)
+
+   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE trdmod          ! ocean active tracers trends 
+   USE trdmod_oce      ! ocean variables trends
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC tra_zdf   !  routine called by step.F90
+
+   !! * module variables
+   INTEGER ::                        &
+      nzdf = 0                         ! type vertical diffusion algorithm used
+      !                                ! defined from ln_zdf...  namlist logicals)
+
+   !! * Module variables
+   REAL(wp), DIMENSION(jpk) ::  &
+      r2dt                          ! vertical profile time-step, = 2 rdttra
+      !                             ! except at nit000 (=rdttra) if neuler=0
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+   
+   SUBROUTINE tra_zdf( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_zdf  ***
+      !!
+      !! ** Purpose :   compute the vertical ocean tracer physics.
+      !! ** Method  :
+      !! ** Action  :
+      !!
+      !!---------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+
+      !! * local declarations
+      INTEGER  ::   jk                   ! Dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   &
+         ztrdt, ztrds                         ! 3D temporary workspace
+      !!---------------------------------------------------------------------
+
+      IF( kt == nit000 )   CALL zdf_ctl          ! initialisation & control of options
+
+      !                                          ! set time step
+      IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
+         r2dt(:) =  rdttra(:)                          ! = rdtra (restarting with Euler time stepping)
+      ELSEIF( kt <= nit000 + 1) THEN                ! at nit000 or nit000+1
+         r2dt(:) = 2. * rdttra(:)                      ! = 2 rdttra (leapfrog)
+      ENDIF
+
+
+      IF( l_trdtra )   THEN                      ! temporary save of ta and sa trends
+         ztrdt(:,:,:) = ta(:,:,:)
+         ztrds(:,:,:) = sa(:,:,:)
+      ENDIF
+
+      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
+      CASE ( -1 )                                       ! esopa: test all possibility with control print
+         CALL tra_zdf_exp    ( kt, r2dt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' zdf0 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_zdf_imp    ( kt, r2dt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' zdf1 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+         CALL tra_zdf_imp_jki( kt, r2dt )
+         CALL prt_ctl( tab3d_1=ta, clinfo1=' zdf2 - Ta: ', mask1=tmask,               &
+            &          tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+      CASE ( 0 )                                       ! explicit scheme
+         CALL tra_zdf_exp    ( kt, r2dt )
+         IF( l_trdtra )   THEN                         ! save the vertical diffusive trends for further diagnostics
+            ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
+            ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
+            CALL trd_mod( ztrdt, ztrds, jpttdzdf, 'TRA', kt )
+         ENDIF
+
+      CASE ( 1 )                                       ! implicit scheme (k-j-i loop)
+         CALL tra_zdf_imp    ( kt, r2dt )
+         IF( l_trdtra )   THEN                         ! save the vertical diffusive trends for further diagnostics
+            DO jk = 1, jpkm1
+               ztrdt(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - ztrdt(:,:,jk)
+               ztrds(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - ztrds(:,:,jk)
+            END DO
+            CALL trd_mod( ztrdt, ztrds, jpttdzdf, 'TRA', kt )
+         ENDIF
+
+      CASE ( 2 )                                       ! implicit scheme (j-k-i loop)
+         CALL tra_zdf_imp_jki( kt, r2dt )
+         IF( l_trdtra )   THEN                         ! save the vertical diffusive trends for further diagnostics
+            DO jk = 1, jpkm1
+               ztrdt(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) - ztrdt(:,:,jk)
+               ztrds(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) - ztrds(:,:,jk)
+            END DO
+            CALL trd_mod( ztrdt, ztrds, jpttdzdf, 'TRA', kt )
+         ENDIF
+
+      END SELECT
+
+      !                                                ! print mean trends (used for debugging)
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ta, clinfo1=' zdf  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=sa, clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+
+   END SUBROUTINE tra_zdf
+
+
+   SUBROUTINE zdf_ctl
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE zdf_ctl  ***
+      !!
+      !! ** Purpose :   Choose the vertical mixing scheme
+      !!
+      !! ** Method  :   Set nzdf from ln_zdfexp and 'key_mpp_omp'.
+      !!      nzdf = 0   explicit (time-splitting) scheme (ln_zdfexp=T)
+      !!           = 1   implicit (euler backward) scheme (ln_zdfexp=F)
+      !!           = 2   implicit (euler backward) scheme with j-k-i loops
+      !!                 (ln_zdfexp=T and 'key_mpp_omp')
+      !!      NB: rotation of lateral mixing operator or TKE or KPP scheme,
+      !!      the implicit scheme is required.
+      !!
+      !!----------------------------------------------------------------------
+      !! * Module used
+      USE zdftke
+      USE zdfkpp
+      !!----------------------------------------------------------------------
+
+      !  Define the vertical tracer physics scheme
+      ! ==========================================
+
+      ! Choice from ln_zdfexp already read in namelist in zdfini module
+      IF( ln_zdfexp ) THEN               ! use explicit scheme
+         nzdf = 0
+      ELSE                               ! use implicit scheme
+         nzdf = 1
+      ENDIF
+
+      ! Force implicit schemes
+      IF( lk_zdftke .OR. lk_zdfkpp   )   nzdf = 1      ! TKE or KPP 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
+
+      ! NEC autotasking / OpenMP
+#if defined key_mpp_omp
+      IF( nzdf == 1 )   nzdf = 2                       ! j-k-i loop
+#endif
+
+      ! Test: esopa
+      IF( lk_esopa )    nzdf = -1                      ! All schemes used
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra:zdf_ctl : vertical tracer physics scheme'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
+         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+         IF( nzdf ==  2 )   WRITE(numout,*) '              Implicit (euler backward) scheme with j-k-i loops'
+      ENDIF
+
+   END SUBROUTINE zdf_ctl
+
+   !!==============================================================================
+END MODULE trazdf
Index: trunk/NEMO/OPA_SRC/TRA/trazdf_imp_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/trazdf_imp_jki.F90	(revision 458)
+++ trunk/NEMO/OPA_SRC/TRA/trazdf_imp_jki.F90	(revision 458)
@@ -0,0 +1,255 @@
+MODULE trazdf_imp_jki
+   !!==============================================================================
+   !!                    ***  MODULE  trazdf_imp_jki  ***
+   !! Ocean active tracers:  vertical component of the tracer mixing trend
+   !!==============================================================================
+   !! History :
+   !!   7.0  !  91-11  (G. Madec)  Original code
+   !!        !  92-06  (M. Imbard)  correction on tracer trend loops
+   !!        !  96-01  (G. Madec)  statement function for e3
+   !!        !  97-05  (G. Madec)  vertical component of isopycnal
+   !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord
+   !!        !  00-08  (G. Madec)  double diffusive mixing
+   !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
+   !!   9.0  !  04-08  (C. Talandier) New trends organization
+   !!        !  05-11  (G. Madec)  New step reorganisation
+   !!----------------------------------------------------------------------
+   !!   tra_zdf_imp_jki  : Update the tracer trend with the diagonal vertical
+   !!                 part of the mixing tensor.
+   !!                 use j-k-i loops.
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE ldfslp          ! Make iso-neutral slopes available
+   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE zdf_oce         ! ocean vertical physics 
+   USE zdfddm          ! ocean vertical physics: double diffusion
+   USE trdmod          ! ocean active tracers trends 
+   USE trdmod_oce      ! ocean variables trends
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl          ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC tra_zdf_imp_jki    ! routine called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "ldftra_substitute.h90"
+#  include "zdfddm_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tra_zdf_imp_jki( kt, p2dt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_zdf_imp_jki  ***
+      !!
+      !! ** Purpose :
+      !!     Compute the trend due to the vertical tracer diffusion inclu-
+      !!     ding the vertical component of lateral mixing (only for second
+      !!     order operator, for fourth order it is already computed and
+      !!     add to the general trend in traldf.F) and add it to the general
+      !!     trend of the tracer equations.
+      !!
+      !! ** Method :
+      !!      save avt coef. resulting from vertical physics alone in zavt:
+      !!         zavt = avt
+      !!      update and save in zavt the vertical eddy viscosity coefficient:
+      !!         avt = avt + wslpi^2+wslj^2
+      !!
+      !!      Second part: vertical trend associated with the vertical physics
+      !!      ===========  (including the vertical flux proportional to dk[t]
+      !!                  associated with the lateral mixing, through the
+      !!                  update of avt)
+      !!      The vertical diffusion of tracers (t & s) 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=ta.
+      !!      Surface and bottom boundary conditions: no diffusive flux on
+      !!      both tracers (bottom, applied through the masked field avt).
+      !!      Add this trend to the general trend ta,sa :
+      !!         ta = ta + dz( avt dz(t) )
+      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
+      !!
+      !!      Third part: recover avt resulting from the vertical physics
+      !!      ==========  alone, for further diagnostics (for example to
+      !!                  compute the turbocline depth in zdfmxl.F90).
+      !!         avt = zavt
+      !!         (avs = zavs if lk_zdfddm=T )
+      !!
+      !! ** Action :
+      !!         Update (ta,sa) arrays with the before vertical diffusion trend
+      !!
+      !!---------------------------------------------------------------------
+      !! * Modules used
+      USE oce                ,   &
+# if defined key_zdfddm
+         zavs => va,  &  ! use va as workspace
+# endif
+         zavt => ua      ! use ua as workspace
+
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
+      REAL(wp), DIMENSION(jpk), INTENT( in ) ::   &
+         p2dt                                 ! vertical profile of tracer time-step
+
+      !! * Local declarations
+      INTEGER  ::   ji, jj, jk                ! dummy loop indices
+      INTEGER  ::   ikst, ikenm2, ikstp1      ! temporary integers
+      REAL(wp) ::   zavi, zavip1              ! temporary scalar
+      REAL(wp), DIMENSION(jpi,jpk) ::   &
+         zwd, zws, zwi,         &  ! ???
+         zwx, zwy, zwz, zwt        ! ???
+      !!---------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_zdf_imp_jki : implicit vertical mixing (j-k-i loops)'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
+      ENDIF
+
+!CDIR PARALLEL DO PRIVATE( zwd, zws, zwi, zwx, zwy, zwz, zwt )
+!$OMP PARALLEL DO PRIVATE( zwd, zws, zwi, zwx, zwy, zwz, zwt )
+      !                                                ! ===============
+      DO jj = 2, jpjm1                                 !  Vertical slab
+         !                                             ! ===============
+
+         ! II. Vertical trend associated with the vertical physics
+         ! =======================================================
+         !     (including the vertical flux proportional to dk[t] associated
+         !      with the lateral mixing, through the avt update)
+         !     dk[ avt dk[ (t,s) ] ] diffusive trends
+
+
+         ! II.0 Matrix construction
+         ! ------------------------
+
+         ! Diagonal, inferior, superior
+         ! (including the bottom boundary condition via avt masked)
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+#if defined key_ldfslp
+               zavi   = avt(ji,jj,jk  ) + fsahtw(ji,jj,jk  )*( wslpi(ji,jj,jk  )*wslpi(ji,jj,jk  )   &
+                  &                                           +wslpj(ji,jj,jk  )*wslpj(ji,jj,jk  ) )
+               zavip1 = avt(ji,jj,jk+1) + fsahtw(ji,jj,jk+1)*( wslpi(ji,jj,jk+1)*wslpi(ji,jj,jk+1)   &
+                  &                                           +wslpj(ji,jj,jk+1)*wslpj(ji,jj,jk+1) )
+#else
+               zavi   = avt(ji,jj,jk  )
+               zavip1 = avt(ji,jj,jk+1)
+#endif
+               zwi(ji,jk) = - p2dt(jk) * zavi   / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
+               zws(ji,jk) = - p2dt(jk) * zavip1 / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
+               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
+            END DO
+         END DO
+
+         ! Surface boudary conditions
+         DO ji = 2, jpim1
+            zwi(ji,1) = 0.e0
+            zwd(ji,1) = 1. - zws(ji,1)
+         END DO
+
+
+         ! II.1. Vertical diffusion on t
+         ! ---------------------------
+
+         ! Second member construction
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               zwy(ji,jk) = tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk)
+            END DO
+         END DO
+
+         ! Matrix inversion from the first level
+         ikst = 1
+#   include "zdf.matrixsolver.h90"
+
+         ! Save the masked temperature after in ta
+         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done
+         !                  it will not be done in tranxt)
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               ta(ji,jj,jk) = zwx(ji,jk) * tmask(ji,jj,jk)
+            END DO
+         END DO
+
+
+         ! II.2 Vertical diffusion on s
+         ! ---------------------------
+
+#if defined key_zdfddm
+         ! Rebuild the Matrix as avt /= avs
+
+         ! Diagonal, inferior, superior
+         ! (including the bottom boundary condition via avs masked)
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+#if defined key_ldfslp
+               zavi   = fsavs(ji,jj,jk  ) + fsahtw(ji,jj,jk  )*( wslpi(ji,jj,jk  )*wslpi(ji,jj,jk  )   &
+                  &                                             +wslpj(ji,jj,jk  )*wslpj(ji,jj,jk  ) )
+               zavip1 = fsavs(ji,jj,jk+1) + fsahtw(ji,jj,jk+1)*( wslpi(ji,jj,jk+1)*wslpi(ji,jj,jk+1)   &
+                  &                                             +wslpj(ji,jj,jk+1)*wslpj(ji,jj,jk+1) )
+#else
+               zavi   = fsavs(ji,jj,jk  )
+               zavip1 = fsavs(ji,jj,jk+1)
+#endif
+               zwi(ji,jk) = - p2dt(jk) * zavi   / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
+               zws(ji,jk) = - p2dt(jk) * zavip1 / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
+               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
+            END DO
+         END DO
+
+         ! Surface boudary conditions
+         DO ji = 2, jpim1
+            zwi(ji,1) = 0.e0
+            zwd(ji,1) = 1. - zws(ji,1)
+         END DO
+#endif
+         ! Second member construction
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               zwy(ji,jk) = sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk)
+            END DO
+         END DO
+
+         ! Matrix inversion from the first level
+         ikst = 1
+#   include "zdf.matrixsolver.h90"
+
+         ! Save the masked salinity after in sa
+         ! (c a u t i o n:  salinity not its trend, Leap-frog scheme done
+         !                  it will not be done in tranxt)
+         DO jk = 1, jpkm1
+            DO ji = 2, jpim1
+               sa(ji,jj,jk) = zwx(ji,jk)  * tmask(ji,jj,jk)
+            END DO
+         END DO
+
+
+#if defined key_ldfslp
+         ! III. recover the avt (avs) resulting from vertical physics only
+         ! ===============================================================
+
+         DO jk = 2, jpkm1
+            DO ji = 2, jpim1
+               avt(ji,jj,jk) = zavt(ji,jj,jk)
+# if defined key_zdfddm
+               fsavs(ji,jj,jk) = zavs(ji,jj,jk)
+# endif
+            END DO
+         END DO
+#endif
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+   END SUBROUTINE tra_zdf_imp_jki
+
+   !!==============================================================================
+END MODULE trazdf_imp_jki
