Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90	(revision 5758)
@@ -6,5 +6,7 @@
    !! History :  2.0  !  2005-11  (G. Madec)  Original code
    !!            3.3  !  2010-09  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
-   !!            4.0  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation
+   !!            3.6  !  2011-06  (G. Madec)  Addition of Mixed Layer Eddy parameterisation
+   !!            4.0  !  2014-05  (G. Madec)  Add 2nd/4th order cases for CEN and FCT schemes 
+   !!             -   !  2014-12  (G. Madec) suppression of cross land advection option
    !!----------------------------------------------------------------------
 
@@ -22,8 +24,8 @@
    USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine)
    USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine)
-   USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine)
    USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine)
    USE cla             ! cross land advection      (cla_traadv     routine)
-   USE ldftra_oce      ! lateral diffusion coefficient on tracers
+   USE ldftra          ! lateral diffusion coefficient on tracers
+   USE ldfslp          ! Lateral diffusion: slopes of neutral surfaces
    !
    USE in_out_manager  ! I/O manager
@@ -74,5 +76,4 @@
       !! ** Method  : - Update (ua,va) with the advection term following nadv
       !!----------------------------------------------------------------------
-      !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
       !
@@ -84,4 +85,5 @@
       !
       CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
+      !
       !                                          ! set time step
       IF( neuler == 0 .AND. kt == nit000 ) THEN     ! at nit000
@@ -95,7 +97,7 @@
       !                                               !==  effective transport  ==!
       DO jk = 1, jpkm1
-         zun(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
-         zvn(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)
-         zwn(:,:,jk) = e1t(:,:) * e2t(:,:)      * wn(:,:,jk)
+         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)                  ! eulerian transport only
+         zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)
+         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
       END DO
       !
@@ -109,14 +111,14 @@
       zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
       !
-      IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &
-         &              CALL tra_adv_eiv( kt, nit000, zun, zvn, zwn, 'TRA' )    ! add the eiv transport (if necessary)
-      !
-      IF( ln_mle    )   CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' )    ! add the mle transport (if necessary)
-      !
-      CALL iom_put( "uocetr_eff", zun )                                         ! output effective transport      
+      IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad )   &
+         &              CALL ldf_eiv_trp( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the eiv transport (if necessary)
+      !
+      IF( ln_mle    )   CALL tra_adv_mle( kt, nit000, zun, zvn, zwn, 'TRA' )   ! add the mle transport (if necessary)
+      !
+      CALL iom_put( "uocetr_eff", zun )                                        ! output effective transport      
       CALL iom_put( "vocetr_eff", zvn )
       CALL iom_put( "wocetr_eff", zwn )
       !
-      IF( ln_diaptr )   CALL dia_ptr( zvn )                                     ! diagnose the effective MSF 
+      IF( ln_diaptr )   CALL dia_ptr( zvn )                                    ! diagnose the effective MSF 
       !
    
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_eiv.F90	(revision 5737)
+++ 	(revision )
@@ -1,220 +1,0 @@
-MODULE traadv_eiv
-   !!======================================================================
-   !!                    ***  MODULE  traadv_eiv  ***
-   !! Ocean tracers:  advection trend - eddy induced velocity
-   !!======================================================================
-   !! History :  1.0  !  2005-11 (G. Madec)  Original code, from traldf and zdf _iso
-   !!            3.3  !  2010-05 (C. Ethe, G. Madec)  merge TRC-TRA 
-   !!----------------------------------------------------------------------
-#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
-   !!----------------------------------------------------------------------
-   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
-   USE iom
-   USE trc_oce         ! share passive tracers/Ocean variables
-# if defined key_diaeiv
-   USE phycst          ! physical constants
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-# endif  
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   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"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )
-      !!----------------------------------------------------------------------
-      !!                  ***  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
-      !!----------------------------------------------------------------------
-      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)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pun      ! in : 3 ocean velocity components 
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pvn      ! out: 3 ocean velocity components
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pwn      ! increased by the eiv
-      !!
-      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zuwk, zuwk1, zuwi, zuwi1   ! local scalars
-      REAL(wp) ::   zvwk, zvwk1, zvwj, zvwj1   !   -      -
-# if defined key_diaeiv 
-      REAL(wp) ::   zztmp                      ! local scalar
-# endif  
-      REAL(wp), POINTER, DIMENSION(:,:) :: zu_eiv, zv_eiv, zw_eiv, z2d
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_adv_eiv')
-      !
-# if defined key_diaeiv 
-      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
-# else
-      CALL wrk_alloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
-# endif
-
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_adv_eiv : eddy induced advection on ', cdtype,' :'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   add to velocity fields the eiv component'
-# if defined key_diaeiv 
-         IF( cdtype == 'TRA') THEN
-            u_eiv(:,:,:) = 0.e0
-            v_eiv(:,:,:) = 0.e0
-            w_eiv(:,:,:) = 0.e0
-         END IF
-# endif
-      ENDIF
-
-      zu_eiv(:,:) = 0.e0   ;   zv_eiv(:,:) = 0.e0   ;    zw_eiv(:,:) = 0.e0  
-      
-                                                    ! =================
-      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(ji,jj) = 0.5 * umask(ji,jj,jk) * ( zuwk - zuwk1 ) 
-               zv_eiv(ji,jj) = 0.5 * vmask(ji,jj,jk) * ( zvwk - zvwk1 ) 
-   
-               pun(ji,jj,jk) = pun(ji,jj,jk) + e2u(ji,jj) * zu_eiv(ji,jj)
-               pvn(ji,jj,jk) = pvn(ji,jj,jk) + e1v(ji,jj) * zv_eiv(ji,jj)
-            END DO
-         END DO
-# if defined key_diaeiv 
-         IF( cdtype == 'TRA') THEN
-            u_eiv(:,:,jk) = zu_eiv(:,:) / fse3u(:,:,jk)
-            v_eiv(:,:,jk) = zv_eiv(:,:) / fse3v(:,:,jk)
-         END IF
-# endif
-         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(ji,jj) = - 0.5 * tmask(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj ) 
-# else
-                  zuwi  = ( wslpi(ji,jj,jk) + wslpi(ji-1,jj,jk) ) * e2u(ji-1,jj) * umask(ji-1,jj,jk)
-                  zuwi1 = ( wslpi(ji,jj,jk) + wslpi(ji+1,jj,jk) ) * e2u(ji  ,jj) * umask(ji  ,jj,jk)
-                  zvwj  = ( wslpj(ji,jj,jk) + wslpj(ji,jj-1,jk) ) * e1v(ji,jj-1) * vmask(ji,jj-1,jk)
-                  zvwj1 = ( wslpj(ji,jj,jk) + wslpj(ji,jj+1,jk) ) * e1v(ji  ,jj) * vmask(ji  ,jj,jk)
-
-                  zw_eiv(ji,jj) = - 0.5 * tmask(ji,jj,jk) * fsaeiw(ji,jj,jk) * ( zuwi1 - zuwi + zvwj1 - zvwj )
-# endif
-                  pwn(ji,jj,jk) = pwn(ji,jj,jk) + zw_eiv(ji,jj)
-               END DO
-            END DO
-# if defined key_diaeiv 
-            IF( cdtype == 'TRA')  w_eiv(:,:,jk) = zw_eiv(:,:) / ( e1e2t(:,:) )
-# endif
-         ENDIF
-         !                                          ! =================
-      END DO                                        !    End of slab  
-      !                                             ! =================
-
-# if defined key_diaeiv 
-      IF( cdtype == 'TRA') THEN
-         CALL iom_put( "uoce_eiv", u_eiv )    ! i-eiv current
-         CALL iom_put( "voce_eiv", v_eiv )    ! j-eiv current
-         CALL iom_put( "woce_eiv", w_eiv )    ! vert. eiv current
-         IF( iom_use('ueiv_heattr') ) THEN
-            zztmp = 0.5 * rau0 * rcp 
-            z2d(:,:) = 0.e0 
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + u_eiv(ji,jj,jk) &
-                       &         * (tsn(ji,jj,jk,jp_tem)+tsn(ji+1,jj,jk,jp_tem)) * e2u(ji,jj) * fse3u(ji,jj,jk) 
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( z2d, 'U', -1. )
-            CALL iom_put( "ueiv_heattr", zztmp * z2d )                  ! heat transport in i-direction
-         ENDIF
-            
-         IF( iom_use('veiv_heattr') ) THEN
-            zztmp = 0.5 * rau0 * rcp 
-            z2d(:,:) = 0.e0 
-            DO jk = 1, jpkm1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1   ! vector opt.
-                     z2d(ji,jj) = z2d(ji,jj) + v_eiv(ji,jj,jk) &
-                     &           * (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * e1v(ji,jj) * fse3v(ji,jj,jk) 
-                  END DO
-               END DO
-            END DO
-            CALL lbc_lnk( z2d, 'V', -1. )
-            CALL iom_put( "veiv_heattr", zztmp * z2d )                  !  heat transport in i-direction
-         ENDIF
-    END IF
-# endif  
-      ! 
-# if defined key_diaeiv 
-      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv, z2d )
-# else
-      CALL wrk_dealloc( jpi, jpj, zu_eiv, zv_eiv, zw_eiv )
-# endif
-      !
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_adv_eiv')
-      !
-    END SUBROUTINE tra_adv_eiv
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :             No rotation of the lateral mixing tensor
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE tra_adv_eiv( kt, kit000, pun, pvn, pwn, cdtype )              ! Empty routine
-      INTEGER  ::   kt    
-      INTEGER  ::   kit000    
-      CHARACTER(len=3) ::   cdtype
-      REAL, DIMENSION(:,:,:) ::   pun, pvn, pwn
-      WRITE(*,*) 'tra_adv_eiv: You should not have seen this print! error?',   &
-         &    kt, cdtype, pun(1,1,1), pvn(1,1,1), pwn(1,1,1)
-   END SUBROUTINE tra_adv_eiv
-#endif
-
-   !!==============================================================================
-END MODULE traadv_eiv
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 5758)
@@ -4,6 +4,9 @@
    !! Ocean Active tracers : lateral diffusive trends 
    !!=====================================================================
-   !! History :  9.0  ! 2005-11 (G. Madec)  Original code
-   !!       NEMO 3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
+   !! History :  9.0  ! 2005-11  (G. Madec)  Original code
+   !!  NEMO      3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
+   !!            3.7  ! 2013-12  (G. Madec) remove the optional computation from T & S anomaly profiles and traldf_bilapg
+   !!             -   ! 2013-12  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
+   !!             -   ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of lateral diffusive operators
    !!----------------------------------------------------------------------
 
@@ -11,25 +14,23 @@
    !!   tra_ldf      : update the tracer trend with the lateral diffusion
    !!   tra_ldf_init : initialization, namelist read, and parameters control
-   !!       ldf_ano  : compute lateral diffusion for constant T-S profiles
-   !!----------------------------------------------------------------------
-   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_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine)
-   USE traldf_lap      ! lateral mixing               (tra_ldf_lap routine)
-   USE trd_oce         ! trends: ocean variables
-   USE trdtra          ! trends manager: tracers 
+   !!----------------------------------------------------------------------
+   USE oce           ! ocean dynamics and tracers
+   USE dom_oce       ! ocean space and time domain
+   USE phycst        ! physical constants
+   USE ldftra        ! lateral diffusion: eddy diffusivity & EIV coeff.
+   USE ldfslp        ! lateral diffusion: iso-neutral slope
+   USE traldf_lap    ! lateral diffusion: laplacian iso-level            operator  (tra_ldf_lap   routine)
+   USE traldf_iso    ! lateral diffusion: laplacian iso-neutral standard operator  (tra_ldf_iso   routine)
+   USE traldf_triad  ! lateral diffusion: laplacian iso-neutral triad    operator  (tra_ldf_triad routine)
+   USE traldf_blp    ! lateral diffusion (iso-level lap/blp)                       (tra_ldf_lap   routine)
+   USE trd_oce       ! trends: ocean variables
+   USE trdtra        ! ocean active tracers 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)
-   USE wrk_nemo        ! Memory allocation
-   USE timing          ! Timing
+   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)
+   USE wrk_nemo       ! Memory allocation
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -37,16 +38,13 @@
 
    PUBLIC   tra_ldf        ! called by step.F90 
-   PUBLIC   tra_ldf_init   ! called by opa.F90 
+   PUBLIC   tra_ldf_init   ! called by nemogcm.F90 
    !
    INTEGER ::   nldf = 0   ! type of lateral diffusion used defined from ln_traldf_... namlist logicals)
-
-   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   t0_ldf, s0_ldf   !: lateral diffusion trends of T & S for a cst profile
-   !                                                               !  (key_traldf_ano only)
-
+   
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -65,67 +63,39 @@
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start('tra_ldf')
-      !
-      rldf = 1     ! For active tracers the 
-
+      IF( nn_timing == 1 )   CALL timing_start('tra_ldf')
+      !
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
-         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
+         CALL wrk_alloc( jpi,jpj,jpk,   ztrdt, ztrds ) 
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
          ztrds(:,:,:) = tsa(:,:,:,jp_sal)
       ENDIF
-
-      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
-      CASE ( 0 )   ;   CALL tra_ldf_lap     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
-                               &                                   tsb, tsa, jpts        )  ! iso-level laplacian
-      CASE ( 1 )                                                                              ! rotated laplacian
-         IF( ln_traldf_grif ) THEN                                                          
-                       CALL tra_ldf_iso_grif( kt, nit000,'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )      ! Griffies operator
-         ELSE                                                                                
-                       CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
-                               &                                  tsb, tsa, jpts, ahtb0 )      ! Madec operator
-         ENDIF
-      CASE ( 2 )   ;   CALL tra_ldf_bilap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
-                               &                                   tsb, tsa, jpts        )  ! iso-level bilaplacian
-      CASE ( 3 )   ;   CALL tra_ldf_bilapg  ( kt, nit000, 'TRA',             tsb, tsa, jpts        )  ! s-coord. geopot. bilap.
+      !
+      SELECT CASE ( nldf )                     !* compute lateral mixing trend and add it to the general trend
+      !
+      CASE ( n_lap   )                                   ! laplacian: iso-level operator
+         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb,      tsa, jpts,  1   )
          !
-      CASE ( -1 )                                ! esopa: test all possibility with control print
-         CALL tra_ldf_lap   ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
-         &                                       tsb, tsa, jpts        ) 
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf0 - Ta: ', mask1=tmask,               &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         IF( ln_traldf_grif ) THEN
-            CALL tra_ldf_iso_grif( kt, nit000, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )
-         ELSE
-            CALL tra_ldf_iso     ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
-            &                                               tsb, tsa, jpts, ahtb0 )  
-         ENDIF
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf1 - Ta: ', mask1=tmask,               &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_bilap ( kt, nit000, 'TRA', gtsu, gtsv, gtui, gtvi,        &
-         &                                       tsb, tsa, jpts        ) 
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf2 - Ta: ', mask1=tmask,               &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
-         CALL tra_ldf_bilapg( kt, nit000, 'TRA',             tsb, tsa, jpts        ) 
-         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf3 - Ta: ', mask1=tmask,               &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      CASE ( n_lap_i )                                   ! laplacian: standard iso-neutral operator (Madec)
+         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
+         !
+      CASE ( n_lap_it )                                  ! laplacian: triad iso-neutral operator (griffies)
+         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts,  1   )
+         !
+      CASE ( n_blp , n_blp_i , n_blp_it )                ! bilaplacian: iso-level & iso-neutral operators
+         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb      , tsa, jpts, nldf )
       END SELECT
 
-#if defined key_traldf_ano
-      tsa(:,:,:,jp_tem) = tsa(:,:,:,jp_tem) - t0_ldf(:,:,:)      ! anomaly: substract the reference diffusivity
-      tsa(:,:,:,jp_sal) = tsa(:,:,:,jp_sal) - s0_ldf(:,:,:)
-#endif
-
-      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
+      IF( l_trdtra )   THEN                    !* save the horizontal diffusive trends for further diagnostics
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_ldf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_ldf, ztrds )
-         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
-      ENDIF
-      !                                          ! print mean trends (used for debugging)
+         CALL wrk_dealloc( jpi,jpj,jpk,   ztrdt, ztrds ) 
+      ENDIF
+      !                                        !* print mean trends (used for debugging)
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' ldf  - Ta: ', mask1=tmask,               &
          &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf')
+      IF( nn_timing == 1 )   CALL timing_stop('tra_ldf')
       !
    END SUBROUTINE tra_ldf
@@ -139,17 +109,9 @@
       !!
       !! ** Method  :   set nldf from the namtra_ldf 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
-      !!----------------------------------------------------------------------
-      INTEGER ::   ioptio, ierr         ! temporary integers 
-      !!----------------------------------------------------------------------
-
-      !  Define the lateral mixing oparator for tracers
-      ! ===============================================
-    
-      IF(lwp) THEN                    ! Namelist print
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio, ierr   ! temporary integers 
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN                     ! Namelist print
          WRITE(numout,*)
          WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator'
@@ -159,204 +121,79 @@
          WRITE(numout,*)
       ENDIF
-
-      !                               ! control the input
+      !                                ! control the input
       ioptio = 0
-      IF( ln_traldf_lap   )   ioptio = ioptio + 1
-      IF( ln_traldf_bilap )   ioptio = ioptio + 1
-      IF( ioptio >  1 )   CALL ctl_stop( '          use ONE or NONE of the 2 lap/bilap operator type on tracer' )
-      IF( ioptio == 0 )   nldf = -2   ! No lateral diffusion
+      IF( ln_traldf_lap )   ioptio = ioptio + 1
+      IF( ln_traldf_blp )   ioptio = ioptio + 1
+      IF( ioptio >  1   )   CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' )
+      IF( ioptio == 0   )   nldf = n_no_ldf     ! No lateral diffusion
       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 )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
-
-      ! defined the type of lateral diffusion from ln_traldf_... logicals
-      ! CAUTION : nldf = 1 is used in trazdf_imp, change it carefully
+      IF( ln_traldf_lev )   ioptio = ioptio + 1
+      IF( ln_traldf_hor )   ioptio = ioptio + 1
+      IF( ln_traldf_iso )   ioptio = ioptio + 1
+      IF( ioptio >  1 )   CALL ctl_stop( 'tra_ldf_init: use only ONE direction (level/hor/iso)' )
+      !
+      !                                ! 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             ! zps-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             ! s-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             ! zps-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             ! s-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( nldf == 3 )   CALL ctl_warn( 'geopotential bilaplacian tracer diffusion in s-coords not thoroughly tested' )
+      IF( ln_traldf_lap ) THEN         ! laplacian operator
+         IF ( ln_zco ) THEN               ! z-coordinate
+            IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_hor   )   nldf = n_lap     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard  (   rotation)
+            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad     (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN               ! z-coordinate with partial step
+            IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed 
+            IF ( ln_traldf_hor   )   nldf = n_lap     ! horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation)
+            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation)
+         ENDIF
+         IF ( ln_sco ) THEN               ! s-coordinate
+            IF ( ln_traldf_lev   )   nldf = n_lap     ! iso-level  (no rotation)
+            IF ( ln_traldf_hor   )   nldf = n_lap_it  ! horizontal (   rotation)       !!gm   a checker....
+            IF ( ln_traldf_iso   )   nldf = n_lap_i   ! iso-neutral: standard (rotation)
+            IF ( ln_traldf_triad )   nldf = n_lap_it  ! iso-neutral: triad    (rotation)
+         ENDIF
+      ENDIF
+      !
+      IF( ln_traldf_blp ) THEN         ! bilaplacian operator
+         IF ( ln_zco ) THEN               ! z-coordinate
+            IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_hor   )   nldf = n_blp     ! iso-level = horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
+            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+         ENDIF
+         IF ( ln_zps ) THEN               ! z-coordinate with partial step
+            IF ( ln_traldf_lev   )   ierr = 1         ! iso-level not allowed 
+            IF ( ln_traldf_hor   )   nldf = n_blp     ! horizontal (no rotation)
+            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
+            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+         ENDIF
+         IF ( ln_sco ) THEN               ! s-coordinate
+            IF ( ln_traldf_lev   )   nldf = n_blp     ! iso-level  (no rotation)
+            IF ( ln_traldf_hor   )   nldf = n_blp_it  ! horizontal (   rotation)       !!gm   a checker....
+            IF ( ln_traldf_iso   )   nldf = n_blp_i   ! iso-neutral: standard (rotation)
+            IF ( ln_traldf_triad )   nldf = n_blp_it  ! iso-neutral: triad    (rotation)
+         ENDIF
+      ENDIF
+      !
       IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
-      IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
-      IF( lk_traldf_eiv .AND. .NOT.ln_traldf_iso )   &
-           CALL ctl_stop( '          eddy induced velocity on tracers',   &
-           &              ' the eddy induced velocity on tracers requires isopycnal laplacian diffusion' )
-      IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
-         IF( .NOT.lk_ldfslp )   CALL ctl_stop( '          the rotation of the diffusive tensor require key_ldfslp' )
-         l_traldf_rot = .TRUE.                 ! needed for trazdf_imp
-      ENDIF
-
-      IF( lk_esopa ) THEN
-         IF(lwp) WRITE(numout,*) '          esopa control: use all lateral physics options'
-         nldf = -1
-      ENDIF
-
+      IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) )                                    &
+           &            CALL ctl_stop( '          eddy induced velocity on tracers requires isopycnal',    &
+           &                                                                    ' laplacian diffusion' )
+      IF(  nldf == n_lap_i .OR. nldf == n_lap_it .OR. &
+         & nldf == n_blp_i .OR. nldf == n_blp_it  )   l_ldfslp = .TRUE.    ! slope of neutral surfaces required 
+      !
       IF(lwp) THEN
          WRITE(numout,*)
-         IF( nldf == -2 )   WRITE(numout,*) '          NO lateral diffusion'
-         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
+         IF( nldf == n_no_ldf )   WRITE(numout,*) '          NO lateral diffusion'
+         IF( nldf == n_lap    )   WRITE(numout,*) '          laplacian iso-level operator'
+         IF( nldf == n_lap_i  )   WRITE(numout,*) '          Rotated laplacian operator (standard)'
+         IF( nldf == n_lap_it )   WRITE(numout,*) '          Rotated laplacian operator (triad)'
+         IF( nldf == n_blp    )   WRITE(numout,*) '          bilaplacian iso-level operator'
+         IF( nldf == n_blp_i  )   WRITE(numout,*) '          Rotated bilaplacian operator (standard)'
+         IF( nldf == n_blp_it )   WRITE(numout,*) '          Rotated bilaplacian operator (triad)'
+      ENDIF
       !
    END SUBROUTINE tra_ldf_init
-
-#if defined key_traldf_ano
-   !!----------------------------------------------------------------------
-   !!   'key_traldf_ano'               T & S lateral diffusion on anomalies
-   !!----------------------------------------------------------------------
-
-   SUBROUTINE ldf_ano
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE ldf_ano  ***
-      !!
-      !! ** Purpose :   initializations of 
-      !!----------------------------------------------------------------------
-      !
-      USE zdf_oce         ! vertical mixing
-      USE trazdf          ! vertical mixing: double diffusion
-      USE zdfddm          ! vertical mixing: double diffusion
-      !
-      INTEGER  ::   jk              ! Dummy loop indice
-      INTEGER  ::   ierr            ! local integer
-      LOGICAL  ::   llsave          ! local logical
-      REAL(wp) ::   zt0, zs0, z12   ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_ref, zs_ref, ztb, zsb, zavt     
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ldf_ano')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, zt_ref, zs_ref, ztb, zsb, zavt ) 
-      !
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'tra:ldf_ano : lateral diffusion acting on anomalies'
-         WRITE(numout,*) '~~~~~~~~~~~'
-      ENDIF
-
-      !                              ! allocate trabbl arrays
-      ALLOCATE( t0_ldf(jpi,jpj,jpk) , s0_ldf(jpi,jpj,jpk) , STAT=ierr )
-      IF( lk_mpp    )   CALL mpp_sum( ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'ldf_ano: unable to allocate arrays' )
-
-      ! 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
-
-      ! Initialisation of gtui/gtvi in case of no cavity
-      IF ( .NOT. ln_isfcav ) THEN
-         gtui(:,:,:) = 0.0_wp
-         gtvi(:,:,:) = 0.0_wp
-      END IF
-      !                                        ! 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   (:,:,:) = tsb (:,:,:,jp_tem)
-      zsb   (:,:,:) = tsb (:,:,:,jp_sal)
-      ua    (:,:,:) = tsa (:,:,:,jp_tem)
-      va    (:,:,:) = tsa (:,:,:,jp_sal)
-      zavt  (:,:,:) = avt(:,:,:)
-      IF( lk_zdfddm ) THEN CALL ctl_stop( ' key_traldf_ano with key_zdfddm not implemented' )
-      ! set tb, sb to reference values and avr to zero
-      tsb (:,:,:,jp_tem) = zt_ref(:,:,:)
-      tsb (:,:,:,jp_sal) = zs_ref(:,:,:)
-      tsa (:,:,:,jp_tem) = 0.e0
-      tsa (:,:,:,jp_sal) = 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(:,:,:) = tsa(:,:,:,jp_tem)
-         s0_ldf(:,:,:) = tsa(:,:,:,jp_sal)
-      ELSE
-         DO jk = 1, jpkm1
-            t0_ldf(:,:,jk) = ( tsa(:,:,jk,jp_tem) - tsb(:,:,jk,jp_tem) ) / ( z12 *rdttra(jk) )
-            s0_ldf(:,:,jk) = ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / ( z12 *rdttra(jk) )
-         END DO
-      ENDIF
-      tsb(:,:,:,jp_tem) = ztb (:,:,:)
-      tsb(:,:,:,jp_sal) = zsb (:,:,:)
-      tsa(:,:,:,jp_tem) = ua  (:,:,:)
-      tsa(:,:,:,jp_sal) = va  (:,:,:)
-      avt(:,:,:)        = zavt(:,:,:)
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, zt_ref, zs_ref, ztb, zsb, zavt ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('ldf_ano')
-      !
-   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
 
    !!======================================================================
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilap.F90	(revision 5737)
+++ 	(revision )
@@ -1,189 +1,0 @@
-MODULE traldf_bilap
-   !!==============================================================================
-   !!                   ***  MODULE  traldf_bilap  ***
-   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
-   !!==============================================================================
-   !! History :  OPA  !  1991-11  (G. Madec)  Original code
-   !!                 !  1993-03  (M. Guyon)  symetrical conditions
-   !!                 !  1995-11  (G. Madec)  suppress volumetric scale factors
-   !!                 !  1996-01  (G. Madec)  statement function for e3
-   !!                 !  1996-01  (M. Imbard)  mpp exchange
-   !!                 !  1997-07  (G. Madec)  optimization, and ahtt
-   !!            8.5  !  2002-08  (G. Madec)  F90: Free form and module
-   !!   NEMO     1.0  !  2004-08  (C. Talandier) New trends organization
-   !!             -   !  2005-11  (G. Madec)  zps or sco as default option
-   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA 
-   !!==============================================================================
-
-   !!----------------------------------------------------------------------
-   !!   tra_ldf_bilap : update the tracer trend with the horizontal diffusion
-   !!                   using a iso-level biharmonic operator
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE ldftra_oce      ! ocean tracer   lateral physics
-   USE in_out_manager  ! I/O manager
-   USE ldfslp          ! iso-neutral slopes 
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE diaptr          ! poleward transport diagnostics
-   USE trc_oce         ! share passive tracers/Ocean variables
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo       ! Memory Allocation
-   USE timing         ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   tra_ldf_bilap   ! routine called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
-#  include "ldfeiv_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_ldf_bilap( kt, kit000, cdtype, pgu, pgv,            &
-      &                                          pgui, pgvi,          &
-      &                                  ptb, pta, kjpt )  
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_ldf_bilap  ***
-      !!
-      !! ** Purpose :   Compute the before horizontal tracer diffusive 
-      !!      trend and add it to the general trend of tracer equation.
-      !!
-      !! ** Method  :   4th order diffusive operator along model level surfaces 
-      !!      evaluated using before fields (forward time scheme). The hor.
-      !!      diffusive trends  is given by:
-      !!      Laplacian of tb:
-      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
-      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
-      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
-      !!        zlt   = ahtt * zlt
-      !!        call to lbc_lnk
-      !!      Bilaplacian (laplacian of zlt):
-      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
-      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
-      !!
-      !!      Add this trend to the general trend
-      !!         (pta) = (pta) + ( difft )
-      !!
-      !! ** Action : - Update pta arrays with the before iso-level
-      !!               biharmonic mixing trend.
-      !!----------------------------------------------------------------------
-      USE oce     , ONLY:   ztu  => ua       , ztv  => va                           ! (ua,va) used as workspace
-      !!
-      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), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at pstep levels
-      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        ! tracer trend 
-      !!
-      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::  zbtr, ztra       ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zeeu, zeev, zlt
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start( 'tra_ldf_bilap')
-      !
-      CALL wrk_alloc( jpi, jpj, zeeu, zeev, zlt ) 
-      !
-
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap : iso-level biharmonic operator on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
-      ENDIF
-      !                                                          ! ===========
-      DO jn = 1, kjpt                                            ! tracer loop
-         !                                                       ! ===========
-         !                                               
-         DO jk = 1, jpkm1                                        ! Horizontal slab
-            !                                             
-            !                          !==  Initialization of metric arrays (for z- or s-coordinates)  ==!
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zeeu(ji,jj) = e2_e1u(ji,jj) * fse3u_n(ji,jj,jk) * umask(ji,jj,jk)
-                  zeev(ji,jj) = e1_e2v(ji,jj) * fse3v_n(ji,jj,jk) * vmask(ji,jj,jk)
-               END DO
-            END DO
-            !                          !==  Laplacian  ==!
-            !
-            DO jj = 1, jpjm1                 ! First derivative (gradient)
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
-                  ztv(ji,jj,jk) = zeev(ji,jj) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
-               END DO
-            END DO
-            !
-            IF( ln_zps ) THEN                ! set gradient at partial step level (last ocean level)
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     IF( mbku(ji,jj) == jk )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgu(ji,jj,jn)
-                     IF( mbkv(ji,jj) == jk )  ztv(ji,jj,jk) = zeev(ji,jj) * pgv(ji,jj,jn)
-                  END DO
-               END DO
-            ENDIF
-            ! (ISH)
-            IF( ln_zps .AND. ln_isfcav ) THEN ! set gradient at partial step level (first ocean level in a cavity)
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1
-                     IF( miku(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeeu(ji,jj) * pgui(ji,jj,jn)
-                     IF( mikv(ji,jj) == MAX(jk,2) )  ztu(ji,jj,jk) = zeev(ji,jj) * pgvi(ji,jj,jn)
-                  END DO
-               END DO
-            ENDIF
-            !
-            DO jj = 2, jpjm1                 ! Second derivative (divergence) time the eddy diffusivity coefficient
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  zlt(ji,jj) = fsahtt(ji,jj,jk) * zbtr * (   ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
-                     &                                     + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)   )
-               END DO
-            END DO
-            CALL lbc_lnk( zlt, 'T', 1. )     ! Lateral boundary conditions (unchanged sgn)
-
-            !                          !==  Bilaplacian  ==!
-            !
-            DO jj = 1, jpjm1                 ! third derivative (gradient)
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  ztu(ji,jj,jk) = zeeu(ji,jj) * ( zlt(ji+1,jj  ) - zlt(ji,jj) )
-                  ztv(ji,jj,jk) = zeev(ji,jj) * ( zlt(ji  ,jj+1) - zlt(ji,jj) )
-               END DO
-            END DO
-            DO jj = 2, jpjm1                 ! fourth derivative (divergence) and add to the general tracer trend
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ! horizontal diffusive trends
-                  zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
-                  ! add it to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
-               END DO
-            END DO
-            !                                             
-         END DO                                           ! Horizontal slab
-         !                                                
-         ! "zonal" mean lateral diffusive heat and salt transport
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
-           IF( jn == jp_tem )  htr_ldf(:) = ptr_sj( ztv(:,:,:) )
-           IF( jn == jp_sal )  str_ldf(:) = ptr_sj( ztv(:,:,:) )
-         ENDIF
-         !                                                ! ===========
-      END DO                                              ! tracer loop
-      !                                                   ! ===========
-      IF( nn_timing == 1 )  CALL timing_stop( 'tra_ldf_bilap')
-      !
-      CALL wrk_dealloc( jpi, jpj, zeeu, zeev, zlt ) 
-      !
-   END SUBROUTINE tra_ldf_bilap
-
-   !!==============================================================================
-END MODULE traldf_bilap
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_bilapg.F90	(revision 5737)
+++ 	(revision )
@@ -1,363 +1,0 @@
-MODULE traldf_bilapg
-   !!==============================================================================
-   !!                       ***  MODULE  traldf_bilapg  ***
-   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
-   !!==============================================================================
-   !! History : 8.0   ! 1997-07  (G. Madec)  Original code
-   !!   NEMO    1.0   ! 2002-08  (G. Madec)  F90: Free form and module
-   !!           3.3   ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
-   !!==============================================================================
-#if defined key_ldfslp   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_ldfslp'                  rotation of the lateral mixing tensor
-   !!----------------------------------------------------------------------
-   !!   tra_ldf_bilapg : update the tracer trend with the horizontal diffusion
-   !!                    using an horizontal biharmonic operator in s-coordinate 
-   !!   ldfght         :  ???
-   !!----------------------------------------------------------------------
-   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 in_out_manager  ! I/O manager
-   USE ldfslp          ! iso-neutral slopes available
-   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
-   USE diaptr          ! poleward transport diagnostics 
-   USE trc_oce         ! share passive tracers/Ocean variables
-   USE lib_mpp         ! MPP library
-   USE wrk_nemo        ! Memory Allocation
-   USE timing          ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   tra_ldf_bilapg   ! routine called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
-#  include "ldfeiv_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE tra_ldf_bilapg  ***
-      !!                    
-      !! ** Purpose :   Compute the before horizontal tracer diffusive 
-      !!      trend and add it to the general trend of tracer equation.
-      !!
-      !! ** Method  :   The lateral diffusive trends is provided by a 4th order
-      !!      operator rotated along geopotential surfaces. It is computed
-      !!      using before fields (forward in time) and geopotential slopes
-      !!      computed in routine inildf.
-      !!         -1- compute the geopotential harmonic operator applied to
-      !!        ptb and multiply it by the eddy diffusivity coefficient
-      !!       (done by a call to ldfght routine, result in wk1 arrays).
-      !!      Applied the domain lateral boundary conditions by call to lbc_lnk
-      !!         -2- compute the geopotential harmonic operator applied to
-      !!      wk1 by a second call to ldfght routine (result in wk2)
-      !!      arrays).
-      !!         -3- Add this trend to the general trend 
-      !!            pta = pta + wk2
-      !!
-      !! ** Action : - Update pta arrays with the before geopotential 
-      !!               biharmonic mixing trend.
-      !!----------------------------------------------------------------------
-      !
-      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), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend 
-      !
-      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zwk1, zwk2 
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_bilapg')
-      !
-      CALL wrk_alloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 ) 
-      !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_bilapg : horizontal biharmonic operator in s-coordinate on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
-      ENDIF
-
-      ! 1. Laplacian of ptb * aht
-      ! ----------------------------- 
-      CALL ldfght( kt, cdtype, ptb, zwk1, kjpt, 1 )      ! rotated harmonic operator applied to ptb and multiply by aht 
-      !                                                 ! output in wk1 
-      !
-      DO jn = 1, kjpt
-         CALL lbc_lnk( zwk1(:,:,:,jn) , 'T', 1. )        ! Lateral boundary conditions on wk1   (unchanged sign)
-      END DO
-
-      ! 2. Bilaplacian of ptb
-      ! -------------------------
-      CALL ldfght( kt, cdtype, zwk1, zwk2, kjpt, 2 )      ! rotated harmonic operator applied to wk1 ; output in wk2
-
-
-      ! 3. Update the tracer trends                    (j-slab :   2, jpj-1)
-      ! ---------------------------
-      DO jn = 1, kjpt
-         DO jj = 2, jpjm1
-            DO jk = 1, jpkm1
-               DO ji = 2, jpim1
-                  ! add it to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zwk2(ji,jj,jk,jn)
-               END DO
-            END DO
-         END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi, jpj, jpk, kjpt, zwk1, zwk2 ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_bilapg')
-      !
-   END SUBROUTINE tra_ldf_bilapg
-
-
-   SUBROUTINE ldfght ( kt, cdtype, pt, plt, kjpt, kaht )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE ldfght  ***
-      !!          
-      !! ** Purpose :   Apply a geopotential harmonic operator to (pt) and 
-      !!      multiply it by the eddy diffusivity coefficient (if kaht=1).
-      !!      Routine only used in s-coordinates (l_sco=T) with bilaplacian
-      !!      operator (ln_traldf_bilap=T) acting along geopotential surfaces
-      !!      (ln_traldf_hor).
-      !!
-      !! ** Method  :   The harmonic operator rotated along geopotential 
-      !!      surfaces is applied to (pt) using the slopes of geopotential
-      !!      surfaces computed in inildf routine. The result is provided in
-      !!      (plt,pls) arrays. It is computed in 2 steps:
-      !!
-      !!      First step: horizontal part of the operator. It is computed on
-      !!      ==========  pt as follows (idem on ps)
-      !!      horizontal fluxes :
-      !!         zftu = e2u*e3u/e1u di[ pt ] - e2u*uslp dk[ mi(mk(pt)) ]
-      !!         zftv = e1v*e3v/e2v dj[ pt ] - e1v*vslp dk[ mj(mk(pt)) ]
-      !!      take the horizontal divergence of the fluxes (no divided by
-      !!      the volume element :
-      !!         plt  = di-1[ zftu ] +  dj-1[ zftv ]
-      !!
-      !!      Second step: vertical part of the operator. It is computed on
-      !!      ===========  pt as follows (idem on ps)
-      !!      vertical fluxes :
-      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pt ]
-      !!              -     e2t     *       wslpi        di[ mi(mk(pt)) ]
-      !!              -     e1t     *       wslpj        dj[ mj(mk(pt)) ]
-      !!      take the vertical divergence of the fluxes add it to the hori-
-      !!      zontal component, divide the result by the volume element and
-      !!      if kaht=1, multiply by the eddy diffusivity coefficient:
-      !!         plt  = aht / (e1t*e2t*e3t) { plt + dk[ zftw ] }
-      !!      else:
-      !!         plt  =  1  / (e1t*e2t*e3t) { plt + dk[ zftw ] }
-      !!
-      !!----------------------------------------------------------------------
-      USE oce     , ONLY:   zftv => ua       ! ua used as workspace
-      !
-      INTEGER         , INTENT(in )                              ::  kt      ! ocean time-step index
-      CHARACTER(len=3), INTENT(in )                              ::  cdtype  ! =TRA or TRC (tracer indicator) 
-      INTEGER         , INTENT(in )                              ::  kjpt    !: dimension of 
-      REAL(wp)        , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) ::  pt      ! tracer fields ( before for 1st call
-      !                                                         ! and laplacian of these fields for 2nd call. 
-      REAL(wp)        , INTENT(out), DIMENSION(jpi,jpj,jpk,kjpt) ::  plt     !: partial harmonic operator applied to  pt  components except
-      !                                                             !: second order vertical derivative term  
-      INTEGER         , INTENT(in )                              ::  kaht    !: =1 multiply the laplacian by the eddy diffusivity coeff.
-      !                                                             !: =2 no multiplication 
-      !!
-      INTEGER ::   ji, jj, jk,jn          ! dummy loop indices
-      !                                   ! temporary scalars
-      REAL(wp) ::  zabe1, zabe2, zmku, zmkv 
-      REAL(wp) ::  zbtr, ztah, ztav
-      REAL(wp) ::  zcof0, zcof1, zcof2, zcof3, zcof4
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zftu, zdkt, zdk1t
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zftw, zdit, zdjt, zdj1t
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('ldfght')
-      !
-      CALL wrk_alloc( jpi, jpj, zftu, zdkt, zdk1t ) 
-      CALL wrk_alloc( jpi, jpk, zftw, zdit, zdjt, zdj1t ) 
-      !
-      DO jn = 1, kjpt
-         !                               ! ********** !   ! ===============
-         DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
-            !                            ! ********** !   ! ===============
-            
-            ! I.1 Vertical gradient of pt and ps at level jk and jk+1
-            ! -------------------------------------------------------
-            !     surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
-            
-            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
-            IF( jk == 1 ) THEN
-               zdkt(:,:) = zdk1t(:,:)
-            ELSE
-               zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
-            ENDIF
-
-
-            ! I.2 Horizontal fluxes
-            ! ---------------------
-            
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zabe1 = e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
-                  zabe2 = e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
-                  
-                  zmku = 1./MAX( tmask(ji+1,jj,jk  )+tmask(ji,jj,jk+1)   &
-                     &          +tmask(ji+1,jj,jk+1)+tmask(ji,jj,jk  ),1. )
-                  zmkv = 1./MAX( tmask(ji,jj+1,jk  )+tmask(ji,jj,jk+1)   &
-                     &          +tmask(ji,jj+1,jk+1)+tmask(ji,jj,jk  ),1. )
-                  
-                  zcof1 = -e2u(ji,jj) * uslp(ji,jj,jk) * zmku
-                  zcof2 = -e1v(ji,jj) * vslp(ji,jj,jk) * zmkv
-                  
-                  zftu(ji,jj)= umask(ji,jj,jk) *   &
-                     &         (  zabe1 *( pt   (ji+1,jj,jk,jn) - pt(ji,jj,jk,jn) )   &
-                     &          + zcof1 *( zdkt (ji+1,jj) + zdk1t(ji,jj)     &
-                     &                    +zdk1t(ji+1,jj) + zdkt (ji,jj) )  )
-                  
-                  zftv(ji,jj,jk)= vmask(ji,jj,jk) *   &
-                     &         (  zabe2 *( pt(ji,jj+1,jk,jn) - pt(ji,jj,jk,jn) )   &
-                     &          + zcof2 *( zdkt (ji,jj+1) + zdk1t(ji,jj)     &
-                     &                    +zdk1t(ji,jj+1) + zdkt (ji,jj) )  )                  
-               END DO
-            END DO
-
-
-            ! I.3 Second derivative (divergence) (not divided by the volume)
-            ! ---------------------
-            
-            DO jj = 2 , jpjm1
-               DO ji = 2 , jpim1
-                  ztah = zftu(ji,jj) - zftu(ji-1,jj) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)
-                  plt(ji,jj,jk,jn) = ztah
-               END DO
-            END DO
-            !                                             ! ===============
-         END DO                                           !   End of slab
-         !                                                ! ===============
-         ! "Poleward" diffusive heat or salt transport
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( kaht == 2 ) ) THEN
-            ! note sign is reversed to give down-gradient diffusive transports (#1043)
-            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) )
-            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) )
-         ENDIF
-
-         !                             ! ************ !   ! ===============
-         DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
-            !                          ! ************ !   ! ===============
-            
-            ! II.1 horizontal tracer gradient
-            ! -------------------------------
-            
-            DO jk = 1, jpk
-               DO ji = 1, jpim1
-                  zdit (ji,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj  ,jk,jn) ) * umask(ji,jj  ,jk)
-                  zdjt (ji,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj  ,jk,jn) ) * vmask(ji,jj  ,jk)
-                  zdj1t(ji,jk) = ( pt(ji  ,jj  ,jk,jn) - pt(ji,jj-1,jk,jn) ) * vmask(ji,jj-1,jk)
-               END DO
-            END DO
-
-
-            ! II.2 Vertical fluxes
-            ! --------------------
-            
-            ! Surface and bottom vertical fluxes set to zero
-            zftw(:, 1 ) = 0.e0
-            zftw(:,jpk) = 0.e0
-            
-            ! interior (2=<jk=<jpk-1)
-            DO jk = 2, jpkm1
-               DO ji = 2, jpim1
-                  zcof0 = e1e2t(ji,jj) / fse3w_n(ji,jj,jk)   &
-                     &     * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)        &
-                     &        + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
-                  
-                  zmku = 1./MAX(  umask(ji  ,jj,jk-1)+umask(ji-1,jj,jk)   &
-                     &           +umask(ji-1,jj,jk-1)+umask(ji  ,jj,jk), 1. )
-                  
-                  zmkv = 1./MAX(  vmask(ji,jj  ,jk-1)+vmask(ji,jj-1,jk)   &
-                     &           +vmask(ji,jj-1,jk-1)+vmask(ji,jj  ,jk), 1. )
-                  
-                  zcof3 = - e2t(ji,jj) * wslpi (ji,jj,jk) * zmku
-                  zcof4 = - e1t(ji,jj) * wslpj (ji,jj,jk) * zmkv
-                  
-                  zftw(ji,jk) = tmask(ji,jj,jk) *   &
-                     &    (  zcof0 * ( pt   (ji,jj,jk-1,jn) - pt   (ji  ,jj,jk,jn) )   &
-                     &     + zcof3 * ( zdit (ji   ,jk-1   ) + zdit (ji-1,jk      )     &
-                     &                +zdit (ji-1 ,jk-1   ) + zdit (ji  ,jk      ) )   &
-                     &     + zcof4 * ( zdjt (ji   ,jk-1   ) + zdj1t(ji  ,jk)     &
-                     &                +zdj1t(ji   ,jk-1   ) + zdjt (ji  ,jk      ) )  )                 
-               END DO
-            END DO
-
-
-            ! II.3 Divergence of vertical fluxes added to the horizontal divergence
-            ! ---------------------------------------------------------------------
-            
-            IF( kaht == 1 ) THEN
-               ! multiply the laplacian by the eddy diffusivity coefficient
-               DO jk = 1, jpkm1
-                  DO ji = 2, jpim1
-                     ! eddy coef. divided by the volume element
-                     zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                     ! vertical divergence
-                     ztav = fsahtt(ji,jj,jk) * ( zftw(ji,jk) - zftw(ji,jk+1) )
-                     ! harmonic operator applied to (pt,ps) and multiply by aht
-                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
-                  END DO
-               END DO
-            ELSEIF( kaht == 2 ) THEN
-               ! second call, no multiplication
-               DO jk = 1, jpkm1
-                  DO ji = 2, jpim1
-                     ! inverse of the volume element
-                     zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                     ! vertical divergence
-                     ztav = zftw(ji,jk) - zftw(ji,jk+1)
-                     ! harmonic operator applied to (pt,ps) 
-                     plt(ji,jj,jk,jn) = ( plt(ji,jj,jk,jn) + ztav ) * zbtr
-                  END DO
-               END DO
-            ELSE
-               IF(lwp) WRITE(numout,*) ' ldfght: kaht= 1 or 2, here =', kaht
-               IF(lwp) WRITE(numout,*) '         We stop'
-               STOP 'ldfght'
-            ENDIF
-            !                                             ! ===============
-         END DO                                           !   End of slab
-         !                                                ! ===============
-      END DO
-      !
-      CALL wrk_dealloc( jpi, jpj, zftu, zdkt, zdk1t ) 
-      CALL wrk_dealloc( jpi, jpk, zftw, zdit, zdjt, zdj1t ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('ldfght')
-      !
-   END SUBROUTINE ldfght
-
-#else  
-   !!----------------------------------------------------------------------
-   !!   Dummy module :             NO rotation of the lateral mixing tensor
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE tra_ldf_bilapg( kt, kit000, cdtype, ptb, pta, kjpt )      ! Empty routine
-      INTEGER :: kt, kit000
-      CHARACTER(len=3) ::   cdtype
-      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
-      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', &
-        &         kt, kit000, cdtype, ptb(1,1,1,1), pta(1,1,1,1), kjpt
-   END SUBROUTINE tra_ldf_bilapg
-#endif
-
-   !!==============================================================================
-END MODULE traldf_bilapg
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_blp.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_blp.F90	(revision 5758)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_blp.F90	(revision 5758)
@@ -0,0 +1,147 @@
+MODULE traldf_blp
+   !!======================================================================
+   !!                   ***  MODULE  traldf_blp  ***
+   !! Ocean  tracers:  bilaplacian lateral mixing trend on tracers
+   !!======================================================================
+   !! History :  3.7  ! 2015-09  (G. Madec)  original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   traldf_blp   : update the tracer trend with the bilaplacian lateral mixing trend
+   !!----------------------------------------------------------------------
+   USE oce           ! ocean dynamics and active tracers
+   USE dom_oce       ! ocean space and time domain
+   USE phycst        ! physical constants
+   USE trc_oce       ! share passive tracers/Ocean variables
+   USE zdf_oce       ! ocean vertical physics
+   USE ldftra        ! lateral physics: eddy diffusivity
+   USE ldfslp        ! lateral physics: iso-neutral slopes
+   USE traldf_lap    ! lateral diffusion (Standard operator)         (tra_ldf_lap   routine)
+   USE traldf_iso    ! lateral diffusion (Standard operator)         (tra_ldf_iso   routine)
+   USE traldf_triad  ! lateral diffusion (Standard operator)         (tra_ldf_triad routine)
+   USE diaptr        ! poleward transport diagnostics
+   USE zpshde        ! partial step: hor. derivative     (zps_hde routine)
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
+   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_ldf_blp         ! routine called by traldf.F90
+
+   !                                    ! Flag to control the type of lateral diffusive operator
+   INTEGER, PARAMETER, PUBLIC ::   n_no_ldf = 00                      ! without operator (i.e. no lateral diffusive trend)
+   !               !!      laplacian    !    bilaplacian   !
+   INTEGER, PARAMETER, PUBLIC ::   n_lap    = 10   ,   n_blp    = 20  ! iso-level operator
+   INTEGER, PARAMETER, PUBLIC ::   n_lap_i  = 11   ,   n_blp_i  = 21  ! standard iso-neutral or geopotential operator 
+   INTEGER, PARAMETER, PUBLIC ::   n_lap_it = 12   ,   n_blp_it = 22  ! triad    iso-neutral or geopotential operator
+
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
+   !! $Id: traldf_iso_grif.F90 5147 2015-03-13 10:01:32Z cetlod $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+      &                                                    pgui, pgvi,   &
+      &                                                    ptb , pta , kjpt, kldf )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tra_ldf_blp  ***
+      !!                    
+      !! ** Purpose :   Compute the before lateral tracer diffusive 
+      !!      trend and add it to the general trend of tracer equation.
+      !!
+      !! ** Method  :   The lateral diffusive trends is provided by a bilaplacian
+      !!      operator applied to before field (forward in time).
+      !!      It is computed by two successive calls to laplacian routine
+      !!
+      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
+      !!----------------------------------------------------------------------
+      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   ) ::   kldf       ! type of operator used
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
+      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        ! tracer trend
+      !
+      INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zlap         ! laplacian at t-point
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zglu, zglv   ! bottom GRADh of the laplacian (u- and v-points)
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zgui, zgvi   ! top    GRADh of the laplacian (u- and v-points)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_blp')
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,kjpt,   zlap ) 
+      CALL wrk_alloc( jpi,jpj,    kjpt,   zglu, zglv, zgui, zgvi ) 
+      !
+      IF( kt == kit000 .AND. lwp )  THEN
+         WRITE(numout,*)
+         SELECT CASE ( kldf )
+         CASE ( n_blp    )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-level   bilaplacian operator on ', cdtype
+         CASE ( n_blp_i  )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'
+         CASE ( n_blp_it )   ;   WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'
+         END SELECT
+         WRITE(numout,*) '~~~~~~~~~~~'
+      ENDIF
+
+      zlap(:,:,:,:) = 0._wp
+      !
+      SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==!
+      !
+      CASE ( n_blp    )                ! iso-level bilaplacian
+         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1 )
+         !
+      CASE ( n_blp_i  )                ! rotated   bilaplacian : standard operator (Madec)
+         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+         !
+      CASE ( n_blp_it )                ! rotated  bilaplacian : triad operator (griffies)
+         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 )
+         !
+      END SELECT
+      !
+      DO jn = 1, kjpt
+         CALL lbc_lnk( zlap(:,:,:,jn) , 'T', 1. )     ! Lateral boundary conditions (unchanged sign)
+      END DO
+      !                                               ! Partial top/bottom cell: GRADh( zlap )  
+      IF( ln_isfcav .AND. ln_zps ) THEN   ;   CALL zps_hde_isf( kt, kjpt, zlap, zglu, zglv, zgui, zgvi )  ! both top & bottom
+      ELSEIF(             ln_zps ) THEN   ;   CALL zps_hde    ( kt, kjpt, zlap, zglu, zglv )              ! only bottom 
+      ENDIF
+      !
+      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==!
+      !
+      CASE ( n_blp    )                ! iso-level bilaplacian
+         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2 )
+         !
+      CASE ( n_blp_i  )                ! rotated   bilaplacian : standard operator (Madec)
+         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+         !
+      CASE ( n_blp_it )                ! rotated  bilaplacian : triad operator (griffies)
+         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 )
+         !
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,kjpt,   zlap ) 
+      CALL wrk_dealloc( jpi,jpj    ,kjpt,   zglu, zglv, zgui, zgvi ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_blp')
+      !
+   END SUBROUTINE tra_ldf_blp
+
+   !!==============================================================================
+END MODULE traldf_blp
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso.F90	(revision 5758)
@@ -4,18 +4,16 @@
    !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
    !!======================================================================
-   !! History :  OPA  !  1994-08  (G. Madec, M. Imbard)
-   !!            8.0  !  1997-05  (G. Madec)  split into traldf and trazdf
-   !!            NEMO !  2002-08  (G. Madec)  Free form, F90
-   !!            1.0  !  2005-11  (G. Madec)  merge traldf and trazdf :-)
-   !!            3.3  !  2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
+   !! History :  OPA  ! 1994-08  (G. Madec, M. Imbard)
+   !!            8.0  ! 1997-05  (G. Madec)  split into traldf and trazdf
+   !!            NEMO ! 2002-08  (G. Madec)  Free form, F90
+   !!            1.0  ! 2005-11  (G. Madec)  merge traldf and trazdf :-)
+   !!            3.3  ! 2010-09  (C. Ethe, G. Madec) Merge TRA-TRC
+   !!            3.7  ! 2014-01  (G. Madec, S. Masson)  restructuration/simplification of aht/aeiv specification
+   !!             -   ! 2014-02  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
    !!----------------------------------------------------------------------
-#if   defined key_ldfslp   ||   defined key_esopa
+
    !!----------------------------------------------------------------------
-   !!   'key_ldfslp'               slope of the lateral diffusive direction
-   !!----------------------------------------------------------------------
-   !!   tra_ldf_iso  : update the tracer trend with the horizontal 
-   !!                  component of a iso-neutral laplacian operator
-   !!                  and with the vertical part of
-   !!                  the isopycnal or geopotential s-coord. operator 
+   !!   tra_ldf_iso  : update the tracer trend with the horizontal component of a iso-neutral laplacian operator
+   !!                  and with the vertical part of the isopycnal or geopotential s-coord. operator 
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and active tracers
@@ -23,7 +21,8 @@
    USE trc_oce         ! share passive tracers/Ocean variables
    USE zdf_oce         ! ocean vertical physics
-   USE ldftra_oce      ! ocean active tracers: lateral physics
+   USE ldftra          ! lateral diffusion: tracer eddy coefficients
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
+   !
    USE in_out_manager  ! I/O manager
    USE iom             ! I/O library
@@ -40,8 +39,7 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -49,7 +47,7 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pgu, pgv,              &
-      &                                pgui, pgvi,                    &
-      &                                ptb, pta, kjpt, pahtb0 )
+  SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+      &                                                   pgui, pgvi,   &
+      &                                       ptb , ptbb, pta , kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_iso  ***
@@ -66,14 +64,15 @@
       !!
       !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
-      !!      ========    with partial cell update if ln_zps=T.
+      !!      ========    with partial cell update if ln_zps=T
+      !!                  with top     cell update if ln_isfcav
       !!
       !!      2nd part :  horizontal fluxes of the lateral mixing operator
       !!      ========    
-      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
-      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
-      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
-      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ]
+      !!         zftu =  pahu e2u*e3u/e1u di[ tb ]
+      !!               - pahu e2u*uslp    dk[ mi(mk(tb)) ]
+      !!         zftv =  pahv e1v*e3v/e2v dj[ tb ]
+      !!               - pahv e2u*vslp    dk[ mj(mk(tb)) ]
       !!      take the horizontal divergence of the fluxes:
-      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
+      !!         difft = 1/(e1e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
       !!      Add this trend to the general trend (ta,sa):
       !!         ta = ta + difft
@@ -82,8 +81,8 @@
       !!      ========  (excluding the vertical flux proportional to dk[t] )
       !!      vertical fluxes associated with the rotated lateral mixing:
-      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
-      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
+      !!         zftw = - {  mi(mk(pahu)) * e2t*wslpi di[ mi(mk(tb)) ]
+      !!                   + mj(mk(pahv)) * e1t*wslpj dj[ mj(mk(tb)) ]  }
       !!      take the horizontal divergence of the fluxes:
-      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
+      !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
       !!      Add this trend to the general trend (ta,sa):
       !!         pta = pta + difft
@@ -91,35 +90,116 @@
       !! ** Action :   Update pta arrays with the before rotated diffusion
       !!----------------------------------------------------------------------
-      USE oce     , ONLY:   zftu => ua       , zftv  => va         ! (ua,va) used as workspace
-      !
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000          ! first 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), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv    ! tracer gradient at pstep levels
-      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgui, pgvi   ! tracer gradient at pstep levels
-      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        ! tracer trend 
-      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef
+      INTEGER                              , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
-      INTEGER  ::  ikt
-      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars
-      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
-      REAL(wp) ::  zcoef0, zbtr, ztra            !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:  ) ::  z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zdkt, zdk1t, zdit, zdjt, ztfw 
+      INTEGER  ::  ierr             ! local integer
+      REAL(wp) ::  zmsku, zahu_w, zabe1, zcof1, zcoef3   ! local scalars
+      REAL(wp) ::  zmskv, zahv_w, zabe2, zcof2, zcoef4   !   -      -
+      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt   !   -      -
+#if defined key_diaar5
+      REAL(wp) ::   zztmp   ! local scalar
+#endif
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdkt, zdk1t, z2d
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdit, zdjt, zftu, zftv, ztfw 
       !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso')
       !
-      CALL wrk_alloc( jpi, jpj,      z2d ) 
-      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
-      !
-
+      CALL wrk_alloc( jpi,jpj,       zdkt, zdk1t, z2d ) 
+      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt , zftu, zftv, ztfw  ) 
+      !
       IF( kt == kit000 )  THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+         !
+         akz     (:,:,:) = 0._wp      
+         ah_wslp2(:,:,:) = 0._wp
+      ENDIF
+      !
+      !                                               ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler)
+      ELSE                                        ;   z2dt = 2.* rdttra(1)      !             (Leapfrog)
+      ENDIF
+      z1_2dt = 1._wp / z2dt
+      !
+      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
+      ELSE                    ;   zsign = -1._wp
+      ENDIF
+         
+         
+      !!----------------------------------------------------------------------
+      !!   0 - calculate  ah_wslp2 and akz
+      !!----------------------------------------------------------------------
+      !
+      IF( kpass == 1 ) THEN                  !==  first pass only  ==!
+         !
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  !
+                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
+                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
+                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
+                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
+                     !
+                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
+                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
+                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
+                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
+                     !
+                  ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
+                     &               + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
+            DO jk = 2, jpkm1
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     akz(ji,jj,jk) = 0.25_wp * (                                                                     &
+                        &              ( pahu(ji  ,jj,jk) + pahu(ji  ,jj,jk-1) ) / ( e1u(ji  ,jj) * e1u(ji  ,jj) )   &
+                        &            + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) )   &
+                        &            + ( pahv(ji,jj  ,jk) + pahv(ji,jj  ,jk-1) ) / ( e2v(ji,jj  ) * e2v(ji,jj  ) )   &
+                        &            + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) )   )
+                  END DO
+               END DO
+            END DO
+            !
+            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
+               DO jk = 2, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, fs_jpim1
+                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )  )
+                     END DO
+                  END DO
+               END DO
+            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
+               DO jk = 2, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, fs_jpim1
+                        ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
+                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
+                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
+                     END DO
+                  END DO
+               END DO
+           ENDIF
+           !
+         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
+            akz(:,:,:) = ah_wslp2(:,:,:)      
+         ENDIF
       ENDIF
       !
@@ -131,7 +211,7 @@
          !!   I - masked horizontal derivative 
          !!----------------------------------------------------------------------
-         !!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient....
-         zdit (1,:,:) = 0.e0     ;     zdit (jpi,:,:) = 0.e0
-         zdjt (1,:,:) = 0.e0     ;     zdjt (jpi,:,:) = 0.e0
+!!gm : bug.... why (x,:,:)?   (1,jpj,:) and (jpi,1,:) should be sufficient....
+         zdit (1,:,:) = 0._wp     ;     zdit (jpi,:,:) = 0._wp
+         zdjt (1,:,:) = 0._wp     ;     zdjt (jpi,:,:) = 0._wp
          !!end
 
@@ -145,61 +225,50 @@
             END DO
          END DO
-
-         ! partial cell correction
-         IF( ln_zps ) THEN      ! partial steps correction at the last ocean level 
-            DO jj = 1, jpjm1
+         IF( ln_zps ) THEN      ! botton and surface ocean correction of the horizontal gradient
+            DO jj = 1, jpjm1              ! bottom correction (partial bottom cell)
                DO ji = 1, fs_jpim1   ! vector opt.
-! IF useless if zpshde defines pgu everywhere
+!!gm  the following anonymous remark is to considered: ! IF useless if zpshde defines pgu everywhere
                   zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)          
                   zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
                END DO
             END DO
+            IF( ln_isfcav ) THEN      ! first wet level beneath a cavity
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     IF( miku(ji,jj) > 1 )   zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)          
+                     IF( mikv(ji,jj) > 1 )   zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
+                  END DO
+               END DO
+            ENDIF
          ENDIF
-         IF( ln_zps .AND. ln_isfcav ) THEN      ! partial steps correction at the first wet level beneath a cavity
-            DO jj = 1, jpjm1
+
+         !!----------------------------------------------------------------------
+         !!   II - horizontal trend  (full)
+         !!----------------------------------------------------------------------
+         !
+         DO jk = 1, jpkm1                                 ! Horizontal slab
+            !
+            !                             !== Vertical tracer gradient
+            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
+            !
+            IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
+            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
+            ENDIF
+!!gm I don't understand why we should need this.... since wmask is used instead of tmask
+!         IF ( ln_isfcav ) THEN
+!            DO jj = 1, jpj
+!               DO ji = 1, jpi   ! vector opt.
+!                  ikt = mikt(ji,jj) ! surface level
+!                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)
+!                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)
+!               END DO
+!            END DO
+!         END IF
+!!gm 
+
+            DO jj = 1 , jpjm1            !==  Horizontal fluxes
                DO ji = 1, fs_jpim1   ! vector opt.
-                  IF (miku(ji,jj) > 1) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn)          
-                  IF (mikv(ji,jj) > 1) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn)     
-               END DO
-            END DO
-         END IF
-
-         !!----------------------------------------------------------------------
-         !!   II - horizontal trend  (full)
-         !!----------------------------------------------------------------------
-!!!!!!!!!!CDIR PARALLEL DO PRIVATE( zdk1t ) 
-            ! 1. Vertical tracer gradient at level jk and jk+1
-            ! ------------------------------------------------
-         ! 
-         ! interior value 
-         DO jk = 2, jpkm1               
-            DO jj = 1, jpj
-               DO ji = 1, jpi   ! vector opt.
-                  zdk1t(ji,jj,jk) = ( ptb(ji,jj,jk,jn  ) - ptb(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1)
-                  !
-                  zdkt(ji,jj,jk)  = ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn  ) ) * wmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
-         zdk1t(:,:,1) = ( ptb(:,:,1,jn  ) - ptb(:,:,2,jn) ) * wmask(:,:,2)
-         zdkt (:,:,1) = zdk1t(:,:,1)
-         IF ( ln_isfcav ) THEN
-            DO jj = 1, jpj
-               DO ji = 1, jpi   ! vector opt.
-                  ikt = mikt(ji,jj) ! surface level
-                  zdk1t(ji,jj,ikt) = ( ptb(ji,jj,ikt,jn  ) - ptb(ji,jj,ikt+1,jn) ) * wmask(ji,jj,ikt+1)
-                  zdkt (ji,jj,ikt) = zdk1t(ji,jj,ikt)
-               END DO
-            END DO
-         END IF
-
-         ! 2. Horizontal fluxes
-         ! --------------------   
-         DO jk = 1, jpkm1
-            DO jj = 1 , jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = ( fsahtu(ji,jj,jk) + pahtb0 ) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
-                  zabe2 = ( fsahtv(ji,jj,jk) + pahtb0 ) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
+                  zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
+                  zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
                   !
                   zmsku = 1. / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
@@ -209,73 +278,31 @@
                      &             + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1. )
                   !
-                  zcof1 = - fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
-                  zcof2 = - fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
+                  zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
+                  zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
                   !
                   zftu(ji,jj,jk ) = (  zabe1 * zdit(ji,jj,jk)   &
-                     &              + zcof1 * (  zdkt (ji+1,jj,jk) + zdk1t(ji,jj,jk)      &
-                     &                         + zdk1t(ji+1,jj,jk) + zdkt (ji,jj,jk)  )  ) * umask(ji,jj,jk)
+                     &               + zcof1 * (  zdkt (ji+1,jj) + zdk1t(ji,jj)      &
+                     &                          + zdk1t(ji+1,jj) + zdkt (ji,jj)  )  ) * umask(ji,jj,jk)
                   zftv(ji,jj,jk) = (  zabe2 * zdjt(ji,jj,jk)   &
-                     &              + zcof2 * (  zdkt (ji,jj+1,jk) + zdk1t(ji,jj,jk)      &
-                     &                         + zdk1t(ji,jj+1,jk) + zdkt (ji,jj,jk)  )  ) * vmask(ji,jj,jk)                  
-               END DO
-            END DO
-
-            ! II.4 Second derivative (divergence) and add to the general trend
-            ! ----------------------------------------------------------------
-            DO jj = 2 , jpjm1
+                     &               + zcof2 * (  zdkt (ji,jj+1) + zdk1t(ji,jj)      &
+                     &                          + zdk1t(ji,jj+1) + zdkt (ji,jj)  )  ) * vmask(ji,jj,jk)                  
+               END DO
+            END DO
+            !
+            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ztra = zbtr * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
-               END DO
-            END DO
-            !                                          ! ===============
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
+                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
+                     &                                        / (  e1e2t(ji,jj) * fse3t(ji,jj,jk)  )
+               END DO
+            END DO
          END DO                                        !   End of slab  
-         !                                             ! ===============
-         !
-         ! "Poleward" diffusive heat or salt transports (T-S case only)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-            ! note sign is reversed to give down-gradient diffusive transports (#1043)
-            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) )
-            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) )
-         ENDIF
- 
-         IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
-           !
-           IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-               z2d(:,:) = 0._wp 
-               DO jk = 1, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
-                     END DO
-                  END DO
-               END DO
-               z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
-               CALL lbc_lnk( z2d, 'U', -1. )
-               CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
-               !
-               z2d(:,:) = 0._wp 
-               DO jk = 1, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
-                     END DO
-                  END DO
-               END DO
-               z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
-               CALL lbc_lnk( z2d, 'V', -1. )
-               CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
-            END IF
-            !
-         ENDIF
-
-         !!----------------------------------------------------------------------
-         !!   III - vertical trend of T & S (extra diagonal terms only)
-         !!----------------------------------------------------------------------
-         
-         ! Local constant initialization
-         ! -----------------------------
-         ztfw(1,:,:) = 0.e0     ;     ztfw(jpi,:,:) = 0.e0
+
+
+         !!----------------------------------------------------------------------
+         !!   III - vertical trend (full)
+         !!----------------------------------------------------------------------
+         
+         ztfw(1,:,:) = 0._wp     ;     ztfw(jpi,:,:) = 0._wp
          
          ! Vertical fluxes
@@ -283,5 +310,5 @@
          
          ! Surface and bottom vertical fluxes set to zero
-         ztfw(:,:, 1 ) = 0.e0      ;      ztfw(:,:,jpk) = 0.e0
+         ztfw(:,:, 1 ) = 0._wp      ;      ztfw(:,:,jpk) = 0._wp
          
          ! interior (2=<jk=<jpk-1)
@@ -289,13 +316,17 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zcoef0 = - fsahtw(ji,jj,jk) * wmask(ji,jj,jk)
-                  !
-                  zmsku = 1./MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
-                     &            + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.  )
-                  zmskv = 1./MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
-                     &            + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.  )
-                  !
-                  zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)
-                  zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
+                  !
+                  zmsku = tmask(ji,jj,jk) / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)          &
+                     &                           + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk) , 1._wp  )
+                  zmskv = tmask(ji,jj,jk) / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)          &
+                     &                           + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk) , 1._wp  )
+                     !
+                  zahu_w = (   pahu(ji  ,jj,jk-1) + pahu(ji-1,jj,jk)    &
+                     &       + pahu(ji-1,jj,jk-1) + pahu(ji  ,jj,jk)  ) * zmsku
+                  zahv_w = (   pahv(ji,jj  ,jk-1) + pahv(ji,jj-1,jk)    &
+                     &       + pahv(ji,jj-1,jk-1) + pahv(ji,jj  ,jk)  ) * zmskv
+                     !
+                  zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk)   !wslpi & j are already w-masked
+                  zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk)
                   !
                   ztfw(ji,jj,jk) = zcoef3 * (   zdit(ji  ,jj  ,jk-1) + zdit(ji-1,jj  ,jk)      &
@@ -306,41 +337,105 @@
             END DO
          END DO
-         
-         
-         ! I.5 Divergence of vertical fluxes added to the general tracer trend
-         ! -------------------------------------------------------------------
-         DO jk = 1, jpkm1
+         !
+         !                                !==  add the vertical 33 flux  ==!
+         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
+            DO jk = 2, jpkm1       
+               DO jj = 1, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   &
+                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
+                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                  END DO
+               END DO
+            END DO
+            !
+         ELSE                                   ! bilaplacian 
+            SELECT CASE( kpass )
+            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
+               DO jk = 2, jpkm1 
+                  DO jj = 1, jpjm1
+                     DO ji = fs_2, fs_jpim1
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
+                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
+                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) * tmask(ji,jj,jk) / fse3w(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO 
+            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+               DO jk = 2, jpkm1 
+                  DO jj = 1, jpjm1
+                     DO ji = fs_2, fs_jpim1
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
+                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                     END DO
+                  END DO
+               END DO
+            END SELECT
+         ENDIF
+         !         
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1.0 / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ztra = (  ztfw(ji,jj,jk) - ztfw(ji,jj,jk+1)  ) * zbtr
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
+                     &                                        / (  e1e2t(ji,jj) * fse3t_n(ji,jj,jk)  )
                END DO
             END DO
          END DO
          !
-      END DO
-      !
-      CALL wrk_dealloc( jpi, jpj, z2d ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw, zdkt, zdk1t ) 
+         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
+             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
+            !
+            !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
+            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
+               ! note sign is reversed to give down-gradient diffusive transports (#1043)
+               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( -zftv(:,:,:) )
+               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( -zftv(:,:,:) )
+            ENDIF
+            !
+            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
+              !
+              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
+                  z2d(:,:) = zftu(ji,jj,1) 
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
+                        END DO
+                     END DO
+                  END DO
+!!gm CAUTION I think there is an error of sign when using BLP operator....
+!!gm         a multiplication by zsign is required (to be checked twice !)
+                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
+                  CALL lbc_lnk( z2d, 'U', -1. )
+                  CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
+                  !
+                  z2d(:,:) = zftv(ji,jj,1) 
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
+                        END DO
+                     END DO
+                  END DO
+                  z2d(:,:) = - rau0_rcp * z2d(:,:)     ! note sign is reversed to give down-gradient diffusive transports (#1043)
+                  CALL lbc_lnk( z2d, 'V', -1. )
+                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
+               END IF
+               !
+            ENDIF
+            !
+         ENDIF                                                    !== end pass selection  ==!
+         !
+         !                                                        ! ===============
+      END DO                                                      ! end tracer loop
+      !                                                           ! ===============
+      !
+      CALL wrk_dealloc( jpi, jpj,      zdkt, zdk1t, z2d ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt , zftu, zftv, ztfw  ) 
       !
       IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso')
       !
    END SUBROUTINE tra_ldf_iso
-
-#else
-   !!----------------------------------------------------------------------
-   !!   default option :   Dummy code   NO rotation of the diffusive tensor
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE tra_ldf_iso( kt, kit000,cdtype, pgu, pgv, pgui, pgvi, ptb, pta, kjpt, pahtb0 )      ! Empty routine
-      INTEGER:: kt, kit000
-      CHARACTER(len=3) ::   cdtype
-      REAL, DIMENSION(:,:,:) ::   pgu, pgv, pgui, pgvi    ! tracer gradient at pstep levels
-      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
-      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt, kit000, cdtype,   &
-         &                       pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0
-   END SUBROUTINE tra_ldf_iso
-#endif
 
    !!==============================================================================
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 5737)
+++ 	(revision )
@@ -1,450 +1,0 @@
-MODULE traldf_iso_grif
-   !!======================================================================
-   !!                   ***  MODULE  traldf_iso_grif  ***
-   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
-   !!======================================================================
-   !! History : 3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)
-   !!                !          Griffies operator version adapted from traldf_iso.F90
-   !!----------------------------------------------------------------------
-#if   defined key_ldfslp   ||   defined key_esopa
-   !!----------------------------------------------------------------------
-   !!   'key_ldfslp'               slope of the lateral diffusive direction
-   !!----------------------------------------------------------------------
-   !!   tra_ldf_iso_grif  : update the tracer trend with the horizontal component
-   !!                       of the Griffies iso-neutral laplacian operator
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE phycst          ! physical constants
-   USE trc_oce         ! share passive tracers/Ocean variables
-   USE zdf_oce         ! ocean vertical physics
-   USE ldftra_oce      ! ocean active tracers: lateral physics
-   USE ldfslp          ! iso-neutral slopes
-   USE diaptr          ! poleward transport diagnostics
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O library
-   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_ldf_iso_grif   ! routine called by traldf.F90
-
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   psix_eiv, psiy_eiv   !: eiv stream function (diag only)
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   ah_wslp2             !: aeiv*w-slope^2
-   REAL(wp),         DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d               !: vertical tracer gradient at 2 levels
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-#  include "ldfeiv_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-  SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv,              &
-       &                                   ptb, pta, kjpt, pahtb0 )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_ldf_iso_grif  ***
-      !!
-      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
-      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
-      !!      add it to the general trend of tracer equation.
-      !!
-      !! ** Method  :   The horizontal component of the lateral diffusive trends
-      !!      is provided by a 2nd order operator rotated along neural or geopo-
-      !!      tential surfaces to which an eddy induced advection can be added
-      !!      It is computed using before fields (forward in time) and isopyc-
-      !!      nal or geopotential slopes computed in routine ldfslp.
-      !!
-      !!      1st part :  masked horizontal derivative of T  ( di[ t ] )
-      !!      ========    with partial cell update if ln_zps=T.
-      !!
-      !!      2nd part :  horizontal fluxes of the lateral mixing operator
-      !!      ========
-      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
-      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
-      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
-      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ]
-      !!      take the horizontal divergence of the fluxes:
-      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
-      !!      Add this trend to the general trend (ta,sa):
-      !!         ta = ta + difft
-      !!
-      !!      3rd part: vertical trends of the lateral mixing operator
-      !!      ========  (excluding the vertical flux proportional to dk[t] )
-      !!      vertical fluxes associated with the rotated lateral mixing:
-      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
-      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
-      !!      take the horizontal divergence of the fluxes:
-      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
-      !!      Add this trend to the general trend (ta,sa):
-      !!         pta = pta + difft
-      !!
-      !! ** Action :   Update pta arrays with the before rotated diffusion
-      !!----------------------------------------------------------------------
-      USE oce     , ONLY:   zftu => ua       , zftv => va            ! (ua,va) used as 3D workspace
-      !
-      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), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
-      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        ! tracer trend
-      REAL(wp)                             , INTENT(in   ) ::   pahtb0     ! background diffusion coef
-      !
-      INTEGER  ::  ji, jj, jk,jn   ! dummy loop indices
-      INTEGER  ::  ip,jp,kp        ! dummy loop indices
-      INTEGER  ::  ierr            ! temporary integer
-      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars
-      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
-      REAL(wp) ::  zcoef0, zbtr                  !   -      -
-      !
-      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv
-      REAL(wp) ::   ze1ur, zdxt, ze2vr, ze3wr, zdyt, zdzt
-      REAL(wp) ::   zah, zah_slp, zaei_slp
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, ztfw 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
-      !!----------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_iso_grif')
-      !
-      CALL wrk_alloc( jpi, jpj,      z2d ) 
-      CALL wrk_alloc( jpi, jpj, jpk, zdit, zdjt, ztfw  ) 
-      !
-
-      IF( kt == kit000 .AND. .NOT.ALLOCATED(ah_wslp2) )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
-         ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt3d(jpi,jpj,0:1), STAT=ierr )
-         IF( lk_mpp   )   CALL mpp_sum ( ierr )
-         IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate arrays')
-         IF( ln_traldf_gdia ) THEN
-            IF (.NOT. ALLOCATED(psix_eiv))THEN
-                ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
-                IF( lk_mpp   )   CALL mpp_sum ( ierr )
-                IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate diagnostics')
-            ENDIF
-         ENDIF
-      ENDIF
-
-      !!----------------------------------------------------------------------
-      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv
-      !!----------------------------------------------------------------------
-
-      !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads
-
-      ah_wslp2(:,:,:) = 0._wp
-      IF( ln_traldf_gdia ) THEN
-         psix_eiv(:,:,:) = 0._wp
-         psiy_eiv(:,:,:) = 0._wp
-      ENDIF
-
-      DO ip = 0, 1
-         DO kp = 0, 1
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1
-                     ze1ur = 1._wp / e1u(ji,jj)
-                     ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
-                     zbu   = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                     zah   = fsahtu(ji,jj,jk)                                  ! fsaht(ji+ip,jj,jk)
-                     zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
-                     ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
-                     ! (do this by *adding* gradient of depth)
-                     zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp)
-                     zslope2 = zslope2 *zslope2
-                     ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp)    &
-                        &                     + zah * ( zbu * ze3wr / ( e1t(ji+ip,jj) * e2t(ji+ip,jj) ) ) * zslope2
-                     IF( ln_traldf_gdia ) THEN
-                        zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew           ! fsaeit(ji+ip,jj,jk)*zslope_skew
-                        psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp
-                     ENDIF
-                  END DO
-               END DO
-            END DO
-         END DO
-      END DO
-      !
-      DO jp = 0, 1
-         DO kp = 0, 1
-            DO jk = 1, jpkm1
-               DO jj = 1, jpjm1
-                  DO ji=1,fs_jpim1
-                     ze2vr = 1._wp / e2v(ji,jj)
-                     ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)
-                     zbv   = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                     zah   = fsahtv(ji,jj,jk)                                  ! fsaht(ji,jj+jp,jk)
-                     zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
-                     ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
-                     !    (do this by *adding* gradient of depth)
-                     zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp)
-                     zslope2 = zslope2 * zslope2
-                     ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp)   &
-                        &                     + zah * ( zbv * ze3wr / ( e1t(ji,jj+jp) * e2t(ji,jj+jp) ) ) * zslope2
-                     IF( ln_traldf_gdia ) THEN
-                        zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew           ! fsaeit(ji,jj+jp,jk)*zslope_skew
-                        psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp
-                     ENDIF
-                  END DO
-               END DO
-            END DO
-         END DO
-      END DO
-      !
-      IF( iom_use("uoce_eiv") .OR. iom_use("voce_eiv") .OR. iom_use("woce_eiv") )  THEN
-         !
-         IF( ln_traldf_gdia .AND. cdtype == 'TRA' ) THEN
-            CALL wrk_alloc( jpi , jpj , jpk  , zw3d )
-            DO jk=1,jpkm1
-               zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
-            END DO
-            zw3d(:,:,jpk) = 0._wp
-            CALL iom_put( "uoce_eiv", zw3d )    ! i-eiv current
-
-            DO jk=1,jpk-1
-               zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
-            END DO
-            zw3d(:,:,jpk) = 0._wp
-            CALL iom_put( "voce_eiv", zw3d )    ! j-eiv current
-
-            DO jk=1,jpk-1
-               DO jj = 2, jpjm1
-                  DO ji = fs_2, fs_jpim1  ! vector opt.
-                     zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2t(ji,jj) + &
-                          &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1t(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
-                  END DO
-               END DO
-            END DO
-            zw3d(:,:,jpk) = 0._wp
-            CALL iom_put( "woce_eiv", zw3d )    ! vert. eiv current
-            CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
-         ENDIF
-         !
-      ENDIF
-      !                                                          ! ===========
-      DO jn = 1, kjpt                                            ! tracer loop
-         !                                                       ! ===========
-         ! Zero fluxes for each tracer
-         ztfw(:,:,:) = 0._wp
-         zftu(:,:,:) = 0._wp
-         zftv(:,:,:) = 0._wp
-         !
-         DO jk = 1, jpkm1                          !==  before lateral T & S gradients at T-level jk  ==!
-            DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-         IF( ln_zps.and.l_grad_zps ) THEN              ! partial steps: correction at the last level
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1
-                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
-                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
-               END DO
-            END DO
-         ENDIF
-
-         !!----------------------------------------------------------------------
-         !!   II - horizontal trend  (full)
-         !!----------------------------------------------------------------------
-         !
-         DO jk = 1, jpkm1
-            !
-            !                    !==  Vertical tracer gradient at level jk and jk+1
-            zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
-            !
-            !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
-            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
-            ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
-            ENDIF
-
-
-            IF (ln_botmix_grif) THEN
-               DO ip = 0, 1              !==  Horizontal & vertical fluxes
-                  DO kp = 0, 1
-                     DO jj = 1, jpjm1
-                        DO ji = 1, fs_jpim1
-                           ze1ur = 1._wp / e1u(ji,jj)
-                           zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
-                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
-                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
-                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
-
-                           zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                           ! ln_botmix_grif is .T. don't mask zah for bottom half cells
-                           zah = fsahtu(ji,jj,jk)   !*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)           ===>>  ????
-                           zah_slp  = zah * zslope_iso
-                           zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew    !fsaeit(ji+ip,jj,jk)*zslope_skew
-                           zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
-                        END DO
-                     END DO
-                  END DO
-               END DO
-
-               DO jp = 0, 1
-                  DO kp = 0, 1
-                     DO jj = 1, jpjm1
-                        DO ji = 1, fs_jpim1
-                           ze2vr = 1._wp / e2v(ji,jj)
-                           zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
-                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
-                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
-                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                           ! ln_botmix_grif is .T. don't mask zah for bottom half cells
-                           zah = fsahtv(ji,jj,jk)        !*vmask(ji,jj,jk+kp)  ! fsaht(ji,jj+jp,jk)
-                           zah_slp = zah * zslope_iso
-                           zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew
-                           zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
-                        END DO
-                     END DO
-                  END DO
-               END DO
-            ELSE
-               DO ip = 0, 1              !==  Horizontal & vertical fluxes
-                  DO kp = 0, 1
-                     DO jj = 1, jpjm1
-                        DO ji = 1, fs_jpim1
-                           ze1ur = 1._wp / e1u(ji,jj)
-                           zdxt  = zdit(ji,jj,jk) * ze1ur
-                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
-                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
-                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
-                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
-
-                           zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
-                           ! ln_botmix_grif is .F. mask zah for bottom half cells
-                           zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! fsaht(ji+ip,jj,jk)   ===>>  ????
-                           zah_slp  = zah * zslope_iso
-                           zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew        ! fsaeit(ji+ip,jj,jk)*zslope_skew
-                           zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
-                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
-                        END DO
-                     END DO
-                  END DO
-               END DO
-
-               DO jp = 0, 1
-                  DO kp = 0, 1
-                     DO jj = 1, jpjm1
-                        DO ji = 1, fs_jpim1
-                           ze2vr = 1._wp / e2v(ji,jj)
-                           zdyt  = zdjt(ji,jj,jk) * ze2vr
-                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
-                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
-                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
-                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
-                           zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
-                           ! ln_botmix_grif is .F. mask zah for bottom half cells
-                           zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! fsaht(ji,jj+jp,jk)
-                           zah_slp = zah * zslope_iso
-                           zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew
-                           zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
-                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
-                        END DO
-                     END DO
-                  END DO
-               END DO
-            END IF
-            !                          !==  divergence and add to the general trend  ==!
-            DO jj = 2 , jpjm1
-               DO ji = fs_2, fs_jpim1  ! vector opt.
-                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (   zftu(ji-1,jj,jk) - zftu(ji,jj,jk)   &
-                     &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )
-               END DO
-            END DO
-            !
-         END DO
-         !
-         DO jk = 1, jpkm1              !== Divergence of vertical fluxes added to the general tracer trend
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-                     &                                / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
-               END DO
-            END DO
-         END DO
-         !
-         !                             ! "Poleward" diffusive heat or salt transports (T-S case only)
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names
-            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) )
-         ENDIF
-
-         IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
-           !
-           IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
-               z2d(:,:) = 0._wp 
-               DO jk = 1, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
-                     END DO
-                  END DO
-               END DO
-               z2d(:,:) = rau0_rcp * z2d(:,:) 
-               CALL lbc_lnk( z2d, 'U', -1. )
-               CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
-               !
-               z2d(:,:) = 0._wp 
-               DO jk = 1, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
-                     END DO
-                  END DO
-               END DO
-               z2d(:,:) = rau0_rcp * z2d(:,:)     
-               CALL lbc_lnk( z2d, 'V', -1. )
-               CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
-            END IF
-            !
-         ENDIF
-         !
-      END DO
-      !
-      CALL wrk_dealloc( jpi, jpj,      z2d ) 
-      CALL wrk_dealloc( jpi, jpj, jpk, zdit, zdjt, ztfw  ) 
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_iso_grif')
-      !
-  END SUBROUTINE tra_ldf_iso_grif
-
-#else
-   !!----------------------------------------------------------------------
-   !!   default option :   Dummy code   NO rotation of the diffusive tensor
-   !!----------------------------------------------------------------------
-   REAL, PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   psix_eiv, psiy_eiv   !: eiv stream function (diag only)
-CONTAINS
-   SUBROUTINE tra_ldf_iso_grif( kt, kit000, cdtype, pgu, pgv,              &
-       &                                   ptb, pta, kjpt, pahtb0 )
-      CHARACTER(len=3) ::   cdtype
-      INTEGER          ::   kit000     ! first time step index
-      REAL, DIMENSION(:,:,:) ::   pgu, pgv   ! tracer gradient at pstep levels
-      REAL, DIMENSION(:,:,:,:) ::   ptb, pta
-      WRITE(*,*) 'tra_ldf_iso_grif: You should not have seen this print! error?', kt, cdtype,    &
-         &                  pgu(1,1,1), pgv(1,1,1), ptb(1,1,1,1), pta(1,1,1,1), kjpt, pahtb0
-   END SUBROUTINE tra_ldf_iso_grif
-#endif
-
-   !!==============================================================================
-END MODULE traldf_iso_grif
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90	(revision 5758)
@@ -2,40 +2,44 @@
    !!==============================================================================
    !!                       ***  MODULE  traldf_lap  ***
-   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
+   !! Ocean tracers:  lateral diffusivity trend  (laplacian and bilaplacian)
    !!==============================================================================
-   !! History :  OPA  !  87-06  (P. Andrich, D. L Hostis)  Original code
-   !!                 !  91-11  (G. Madec)
-   !!                 !  95-11  (G. Madec)  suppress volumetric scale factors
-   !!                 !  96-01  (G. Madec)  statement function for e3
-   !!            NEMO !  02-06  (G. Madec)  F90: Free form and module
-   !!            1.0  !  04-08  (C. Talandier) New trends organization
-   !!                 !  05-11  (G. Madec)  add zps case
-   !!            3.0  !  10-06  (C. Ethe, G. Madec) Merge TRA-TRC
+   !! History :  OPA  ! 1987-06  (P. Andrich, D. L Hostis)  Original code
+   !!                 ! 1991-11  (G. Madec)
+   !!                 ! 1995-11  (G. Madec)  suppress volumetric scale factors
+   !!                 ! 1996-01  (G. Madec)  statement function for e3
+   !!            NEMO ! 2002-06  (G. Madec)  F90: Free form and module
+   !!            1.0  ! 2004-08  (C. Talandier) New trends organization
+   !!                 ! 2005-11  (G. Madec)  add zps case
+   !!            3.0  ! 2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
+   !!            3.7  ! 2014-01  (G. Madec, S. Masson) re-entrant laplacian 
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   tra_ldf_lap  : update the tracer trend with the horizontal diffusion
-   !!                 using a iso-level harmonic (laplacien) operator.
+   !!   tra_ldf_lap : update the tracer trend with the lateral diffusion : iso-level laplacian operator
+   !!   tra_ldf_blp : update the tracer trend with the lateral diffusion : iso-level bilaplacian operator
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
-   USE ldftra_oce      ! ocean active tracers: lateral physics
-   USE in_out_manager  ! I/O manager
+   USE ldftra          ! lateral physics: eddy diffusivity
    USE diaptr          ! poleward transport diagnostics
    USE trc_oce         ! share passive tracers/Ocean variables
-   USE lib_mpp         ! MPP library
+   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
+   !
+   USE in_out_manager  ! I/O manager
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! distribued memory computing library
    USE timing          ! Timing
+   USE wrk_nemo        ! Memory allocation
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tra_ldf_lap   ! routine called by step.F90
+   PUBLIC   tra_ldf_lap   ! routine called by traldf.F90
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -43,7 +47,7 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pgu , pgv ,    &
-      &                                        pgui, pgvi,    &
-      &                                ptb, pta, kjpt ) 
+   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+      &                                                    pgui, pgvi,   &
+      &                                        ptb , pta , kjpt, kpass ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_lap  ***
@@ -55,6 +59,6 @@
       !!      fields (forward time scheme). The horizontal diffusive trends of 
       !!      the tracer is given by:
-      !!          difft = 1/(e1t*e2t*e3t) {  di-1[ aht e2u*e3u/e1u di(tb) ]
-      !!                                   + dj-1[ aht e1v*e3v/e2v dj(tb) ] }
+      !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
+      !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
       !!      Add this trend to the general tracer trend pta :
       !!          pta = pta + difft
@@ -63,106 +67,100 @@
       !!                harmonic mixing trend.
       !!----------------------------------------------------------------------
-      USE oce, ONLY:   ztu => ua , ztv => va  ! (ua,va) used as workspace
-      !
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000          ! first 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   ) ::   kpass      ! =1/2 first or second passage
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
       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        ! tracer trend 
       !
-      INTEGER  ::   ji, jj, jk, jn       ! dummy loop indices
-      INTEGER  ::   iku, ikv, ierr       ! local integers
-      REAL(wp) ::   zabe1, zabe2, zbtr   ! local scalars
+      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
+      REAL(wp) ::   zsign            ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztu, ztv, zaheeu, zaheev
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 ) CALL timing_start('tra_ldf_lap')
+      IF( nn_timing == 1 )   CALL timing_start('tra_ldf_lap')
       !
-      IF( kt == kit000 )  THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+      IF( kt == nit000 .AND. lwp )  THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'tra_ldf_lap : iso-level laplacian diffusion on ', cdtype, ', pass=', kpass
+         WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
-
-      !                                                          ! =========== !
-      DO jn = 1, kjpt                                            ! tracer loop !
-         !                                                       ! =========== !    
-         DO jk = 1, jpkm1                                            ! slab loop
-            !                                           
-            ! 1. First derivative (gradient)
-            ! -------------------
+      !
+      CALL wrk_alloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
+      !
+      !                                !==  Initialization of metric arrays used for all tracers  ==!
+      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
+      ELSE                    ;   zsign = -1._wp
+      ENDIF
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)   !!gm   * umask(ji,jj,jk) pah masked!
+               zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)   !!gm   * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      !                             ! =========== !
+      DO jn = 1, kjpt               ! tracer loop !
+         !                          ! =========== !    
+         !                               
+         DO jk = 1, jpkm1              !== First derivative (gradient)  ==!
             DO jj = 1, jpjm1
-               DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e2_e1u(ji,jj) * fse3u_n(ji,jj,jk)
-                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * fse3v_n(ji,jj,jk)
-                  ztu(ji,jj,jk) = zabe1 * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
-                  ztv(ji,jj,jk) = zabe2 * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
+               DO ji = 1, fs_jpim1
+                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
+                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
                END DO
             END DO
-            IF( ln_zps ) THEN      ! set gradient at partial step level for the last ocean cell
+         END DO  
+         IF( ln_zps ) THEN                ! set gradient at bottom/top ocean level
+            DO jj = 1, jpjm1                    ! bottom
+               DO ji = 1, fs_jpim1
+                  ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn)
+                  ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn)
+               END DO
+            END DO  
+            IF( ln_isfcav ) THEN                ! top in ocean cavities only
                DO jj = 1, jpjm1
                   DO ji = 1, fs_jpim1   ! vector opt.
-                     ! last level
-                     iku = mbku(ji,jj)
-                     ikv = mbkv(ji,jj)
-                     IF( iku == jk ) THEN
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e2_e1u(ji,jj) * fse3u_n(ji,jj,iku)
-                        ztu(ji,jj,jk) = zabe1 * pgu(ji,jj,jn)
-                     ENDIF
-                     IF( ikv == jk ) THEN
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e1_e2v(ji,jj) * fse3v_n(ji,jj,ikv)
-                        ztv(ji,jj,jk) = zabe2 * pgv(ji,jj,jn)
-                     ENDIF
+                     IF( miku(ji,jj) > 1 )   ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 
+                     IF( mikv(ji,jj) > 1 )   ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) 
                   END DO
                END DO
             ENDIF
-            ! (ISH)
-            IF( ln_zps .AND. ln_isfcav ) THEN      ! set gradient at partial step level for the first ocean cell
-                                                   ! into a cavity
-               DO jj = 1, jpjm1
-                  DO ji = 1, fs_jpim1   ! vector opt.
-                     ! ice shelf level level MAX(2,jk) => only where ice shelf
-                     iku = miku(ji,jj) 
-                     ikv = mikv(ji,jj) 
-                     IF( iku == MAX(2,jk) ) THEN 
-                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e2_e1u(ji,jj) * fse3u_n(ji,jj,iku) 
-                        ztu(ji,jj,jk) = zabe1 * pgui(ji,jj,jn) 
-                     ENDIF 
-                     IF( ikv == MAX(2,jk) ) THEN 
-                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e1_e2v(ji,jj) * fse3v_n(ji,jj,ikv) 
-                        ztv(ji,jj,jk) = zabe2 * pgvi(ji,jj,jn) 
-                     END IF 
-                  END DO
-               END DO
-            ENDIF
-         
-         
-            ! 2. Second derivative (divergence) added to the general tracer trends
-            ! ---------------------------------------------------------------------
+         ENDIF
+         !
+         DO jk = 1, jpkm1              !== Second derivative (divergence) added to the general tracer trends  ==!
             DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
-                  ! horizontal diffusive trends added to the general tracer trends
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
-                     &                                          + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
+               DO ji = fs_2, fs_jpim1
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
+                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
+                     &                                / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) )
                END DO
             END DO
-            !
-         END DO                                             !  End of slab  
+         END DO  
          !
-         ! "Poleward" diffusive heat or salt transports
-         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
-            IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( ztv(:,:,:) )
-            IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( ztv(:,:,:) )
+         !                             !== "Poleward" diffusive heat or salt transports  ==!
+         IF( ( kpass == 1 .AND. .NOT.ln_traldf_blp ) .OR.  &     !==  first pass only (  laplacian)  ==!
+             ( kpass == 2 .AND.      ln_traldf_blp ) ) THEN      !==  2nd   pass only (bilaplacian)  ==!
+            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
+               IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( -ztv(:,:,:) )
+               IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( -ztv(:,:,:) )
+            ENDIF
          ENDIF
-         !                                                  ! ==================
-      END DO                                                ! end of tracer loop
-      !                                                     ! ==================
-      IF( nn_timing == 1 ) CALL timing_stop('tra_ldf_lap')
+         !                          ! ==================
+      END DO                        ! end of tracer loop
+      !                             ! ==================
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk,   ztu, ztv, zaheeu, zaheev ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_lap')
       !
    END SUBROUTINE tra_ldf_lap
-
+   
    !!==============================================================================
 END MODULE traldf_lap
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 5758)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90	(revision 5758)
@@ -0,0 +1,470 @@
+MODULE traldf_triad
+   !!======================================================================
+   !!                   ***  MODULE  traldf_triad  ***
+   !! Ocean  tracers:  horizontal component of the lateral tracer mixing trend
+   !!======================================================================
+   !! History :  3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  Griffies operator (original code)
+   !!            3.7  ! 2013-12  (F. Lemarie, G. Madec)  triad operator (Griffies) + Method of Stabilizing Correction
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and active tracers
+   USE dom_oce         ! ocean space and time domain
+   USE phycst          ! physical constants
+   USE trc_oce         ! share passive tracers/Ocean variables
+   USE zdf_oce         ! ocean vertical physics
+   USE ldftra          ! lateral physics: eddy diffusivity
+   USE ldfslp          ! lateral physics: iso-neutral slopes
+   USE traldf_iso      ! lateral diffusion (Madec operator)         (tra_ldf_iso routine)
+   USE diaptr          ! poleward transport diagnostics
+   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
+   !
+   USE in_out_manager  ! I/O manager
+   USE iom             ! I/O library
+   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_ldf_triad   ! routine called by traldf.F90
+
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE ::   zdkt3d   !: vertical tracer gradient at 2 levels
+
+   !! * Substitutions
+#  include "domzgr_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_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
+      &                                                     pgui, pgvi,   &
+      &                                         ptb , ptbb, pta , kjpt, kpass )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_ldf_triad  ***
+      !!
+      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
+      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
+      !!      add it to the general trend of tracer equation.
+      !!
+      !! ** Method  :   The horizontal component of the lateral diffusive trends
+      !!      is provided by a 2nd order operator rotated along neural or geopo-
+      !!      tential surfaces to which an eddy induced advection can be added
+      !!      It is computed using before fields (forward in time) and isopyc-
+      !!      nal or geopotential slopes computed in routine ldfslp.
+      !!
+      !!      see documentation for the desciption
+      !!
+      !! ** Action :   pta   updated with the before rotated diffusion
+      !!               ah_wslp2 ....
+      !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp)
+      !!----------------------------------------------------------------------
+      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   ) ::   kpass      ! =1/2 first or second passage
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(in   ) ::   pahu, pahv ! eddy diffusivity at u- and v-points  [m2/s]
+      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      !
+      INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
+      INTEGER  ::  ip,jp,kp         ! dummy loop indices
+      INTEGER  ::  ierr            ! local integer
+      REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3          ! local scalars
+      REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4          !   -      -
+      REAL(wp) ::  zcoef0, ze3w_2, zsign, z2dt, z1_2dt  !   -      -
+      !
+      REAL(wp) ::   zslope_skew, zslope_iso, zslope2, zbu, zbv
+      REAL(wp) ::   ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt
+      REAL(wp) ::   zah, zah_slp, zaei_slp
+#if defined key_diaar5
+      REAL(wp) ::   zztmp              ! local scalar
+#endif
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: z2d                                            ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw   ! 3D     -
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_ldf_triad')
+      !
+      CALL wrk_alloc( jpi,jpj,       z2d ) 
+      CALL wrk_alloc( jpi,jpj,jpk,   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw  ) 
+      !
+      IF( .NOT.ALLOCATED(zdkt3d) )  THEN
+         ALLOCATE( zdkt3d(jpi,jpj,0:1) , STAT=ierr )
+         IF( lk_mpp   )   CALL mpp_sum ( ierr )
+         IF( ierr > 0 )   CALL ctl_stop('STOP', 'tra_ldf_triad: unable to allocate arrays')
+      ENDIF
+     !
+      IF( kpass == 1 .AND. kt == kit000 )  THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_ldf_triad : rotated laplacian diffusion operator on ', cdtype
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+      !
+      !                                               ! set time step size (Euler/Leapfrog)
+      IF( neuler == 0 .AND. kt == kit000 ) THEN   ;   z2dt =     rdttra(1)      ! at nit000   (Euler)
+      ELSE                                        ;   z2dt = 2.* rdttra(1)      !             (Leapfrog)
+      ENDIF
+      z1_2dt = 1._wp / z2dt
+      !
+      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign (eddy diffusivity >0)
+      ELSE                    ;   zsign = -1._wp
+      ENDIF
+                  
+      !!----------------------------------------------------------------------
+      !!   0 - calculate  ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw
+      !!----------------------------------------------------------------------
+      !
+      IF( kpass == 1 ) THEN         !==  first pass only  and whatever the tracer is  ==!
+         !
+         akz     (:,:,:) = 0._wp      
+         ah_wslp2(:,:,:) = 0._wp
+         IF( ln_ldfeiv_dia ) THEN
+            zpsi_uw(:,:,:) = 0._wp
+            zpsi_vw(:,:,:) = 0._wp
+         ENDIF
+         !
+         DO ip = 0, 1                            ! i-k triads
+            DO kp = 0, 1
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, fs_jpim1
+                        ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                        zbu   = e1e2u(ji,jj) * fse3u(ji,jj,jk)
+                        zah   = 0.25_wp * pahu(ji,jj,jk)
+                        zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
+                        ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth)
+                        zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                        zslope2 = zslope2 *zslope2
+                        ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2
+                        akz     (ji+ip,jj,jk+kp) = akz     (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj)       &
+                           &                                                      * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)
+                        !
+                       IF( ln_ldfeiv_dia )   zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp)   &
+                           &                                       + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew
+                     END DO
+                  END DO
+               END DO
+            END DO
+         END DO
+         !
+         DO jp = 0, 1                            ! j-k triads 
+            DO kp = 0, 1
+               DO jk = 1, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, fs_jpim1
+                        ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)
+                        zbv   = e1e2v(ji,jj) * fse3v(ji,jj,jk)
+                        zah   = 0.25_wp * pahv(ji,jj,jk)
+                        zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
+                        ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces
+                        !    (do this by *adding* gradient of depth)
+                        zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                        zslope2 = zslope2 * zslope2
+                        ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2
+                        akz     (ji,jj+jp,jk+kp) = akz     (ji,jj+jp,jk+kp) + zah * r1_e2v(ji,jj)     &
+                           &                                                      * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)
+                        !
+                        IF( ln_ldfeiv_dia )   zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp)   &
+                           &                                       + 0.25 * aeiv(ji,jj,jk) * e1v(ji,jj) * zslope_skew
+                     END DO
+                  END DO
+               END DO
+            END DO
+         END DO
+         !
+         IF( ln_traldf_msc ) THEN                ! stabilizing vertical diffusivity coefficient
+            !
+            IF( ln_traldf_blp ) THEN                ! bilaplacian operator
+               DO jk = 2, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, fs_jpim1
+                        akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk)   &
+                           &          * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )  )
+                     END DO
+                  END DO
+               END DO
+            ELSEIF( ln_traldf_lap ) THEN              ! laplacian operator
+               DO jk = 2, jpkm1
+                  DO jj = 1, jpjm1
+                     DO ji = 1, fs_jpim1
+                        ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
+                        zcoef0 = z2dt * (  akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2  )
+                        akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt
+                     END DO
+                  END DO
+               END DO
+           ENDIF
+           !
+         ELSE                                    ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit
+            akz(:,:,:) = ah_wslp2(:,:,:)      
+         ENDIF
+         !
+         IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' )   CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )
+         !
+      ENDIF                                  !==  end 1st pass only  ==!
+      !
+      !                                                           ! ===========
+      DO jn = 1, kjpt                                             ! tracer loop
+         !                                                        ! ===========
+         ! Zero fluxes for each tracer
+!!gm  this should probably be done outside the jn loop
+         ztfw(:,:,:) = 0._wp
+         zftu(:,:,:) = 0._wp
+         zftv(:,:,:) = 0._wp
+         !
+         DO jk = 1, jpkm1        !==  before lateral T & S gradients at T-level jk  ==!
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         IF( ln_zps .AND. l_grad_zps ) THEN    ! partial steps: correction at top/bottom ocean level
+            DO jj = 1, jpjm1                       ! bottom level
+               DO ji = 1, fs_jpim1   ! vector opt.
+                  zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn)
+                  zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn)
+               END DO
+            END DO
+            IF( ln_isfcav ) THEN                   ! top level (ocean cavities only)
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     IF( miku(ji,jj)  > 1 )   zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 
+                     IF( mikv(ji,jj)  > 1 )   zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) 
+                  END DO
+               END DO
+            ENDIF
+         ENDIF
+
+         !!----------------------------------------------------------------------
+         !!   II - horizontal trend  (full)
+         !!----------------------------------------------------------------------
+         !
+         DO jk = 1, jpkm1
+            !
+            !                    !==  Vertical tracer gradient at level jk and jk+1
+            zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
+            !
+            !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
+            IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
+            ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
+            ENDIF
+            !
+            zaei_slp = 0._wp
+            !
+            IF( ln_botmix_triad ) THEN
+               DO ip = 0, 1              !==  Horizontal & vertical fluxes
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze1ur = r1_e1u(ji,jj)
+                           zdxt  = zdit(ji,jj,jk) * ze1ur
+                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
+                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
+                           zslope_iso  = triadi  (ji+ip,jj,jk,1-ip,kp)
+
+                           zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)
+                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????   ahu is masked....
+                           zah = pahu(ji,jj,jk)
+                           zah_slp  = zah * zslope_iso
+                           IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew
+                           zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
+                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt                 * zbu * ze3wr
+                        END DO
+                     END DO
+                  END DO
+               END DO
+
+               DO jp = 0, 1
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze2vr = r1_e2v(ji,jj)
+                           zdyt  = zdjt(ji,jj,jk) * ze2vr
+                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
+                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
+                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
+                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)
+                           ! ln_botmix_triad is .T. don't mask zah for bottom half cells    !!gm ?????  ahv is masked...
+                           zah = pahv(ji,jj,jk)
+                           zah_slp = zah * zslope_iso
+                           IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew
+                           zftv(ji,jj   ,jk   ) = zftv(ji,jj   ,jk   ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
+                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt                * zbv * ze3wr
+                        END DO
+                     END DO
+                  END DO
+               END DO
+
+            ELSE
+
+               DO ip = 0, 1               !==  Horizontal & vertical fluxes
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze1ur = r1_e1u(ji,jj)
+                           zdxt  = zdit(ji,jj,jk) * ze1ur
+                           ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)
+                           zdzt  = zdkt3d(ji+ip,jj,kp) * ze3wr
+                           zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
+                           zslope_iso  = triadi(ji+ip,jj,jk,1-ip,kp)
+
+                           zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)
+                           ! ln_botmix_triad is .F. mask zah for bottom half cells
+                           zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp)         ! pahu(ji+ip,jj,jk)   ===>>  ????
+                           zah_slp  = zah * zslope_iso
+                           IF( ln_ldfeiv )   zaei_slp = aeiu(ji,jj,jk) * zslope_skew        ! fsaeit(ji+ip,jj,jk)*zslope_skew
+                           zftu(ji   ,jj,jk   ) = zftu(ji   ,jj,jk   ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur
+                           ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr
+                        END DO
+                     END DO
+                  END DO
+               END DO
+
+               DO jp = 0, 1
+                  DO kp = 0, 1
+                     DO jj = 1, jpjm1
+                        DO ji = 1, fs_jpim1
+                           ze2vr = r1_e2v(ji,jj)
+                           zdyt  = zdjt(ji,jj,jk) * ze2vr
+                           ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)
+                           zdzt  = zdkt3d(ji,jj+jp,kp) * ze3wr
+                           zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
+                           zslope_iso  = triadj(ji,jj+jp,jk,1-jp,kp)
+                           zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)
+                           ! ln_botmix_triad is .F. mask zah for bottom half cells
+                           zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp)         ! pahv(ji,jj+jp,jk)  ????
+                           zah_slp = zah * zslope_iso
+                           IF( ln_ldfeiv )   zaei_slp = aeiv(ji,jj,jk) * zslope_skew        ! fsaeit(ji,jj+jp,jk)*zslope_skew
+                           zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr
+                           ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr
+                        END DO
+                     END DO
+                  END DO
+               END DO
+            ENDIF
+            !                             !==  horizontal divergence and add to the general trend  ==!
+            DO jj = 2 , jpjm1
+               DO ji = fs_2, fs_jpim1  ! vector opt.
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
+                     &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
+                     &                                        / (  e1e2t(ji,jj) * fse3t(ji,jj,jk)  )
+               END DO
+            END DO
+            !
+         END DO
+         !
+         !                                !==  add the vertical 33 flux  ==!
+         IF( ln_traldf_lap ) THEN               ! laplacian case: eddy coef = ah_wslp2 - akz
+            DO jk = 2, jpkm1       
+               DO jj = 1, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)   &
+                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
+                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                  END DO
+               END DO
+            END DO
+         ELSE                                   ! bilaplacian 
+            SELECT CASE( kpass )
+            CASE(  1  )                            ! 1st pass : eddy coef = ah_wslp2
+               DO jk = 2, jpkm1 
+                  DO jj = 1, jpjm1
+                     DO ji = fs_2, fs_jpim1
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)             &
+                           &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                     END DO
+                  END DO
+               END DO 
+            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+               DO jk = 2, jpkm1 
+                  DO jj = 1, jpjm1
+                     DO ji = fs_2, fs_jpim1
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk)                      &
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
+                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                     END DO
+                  END DO
+               END DO
+            END SELECT 
+         ENDIF
+         !
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1  ! vector opt.
+                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
+                     &                                        / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         !
+         IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR.  &     !==  first pass only (  laplacian)  ==!
+             ( kpass == 2 .AND. ln_traldf_blp ) ) THEN      !==  2nd   pass      (bilaplacian)  ==!
+            !
+            !                          ! "Poleward" diffusive heat or salt transports (T-S case only)
+            IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
+               IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )        ! 3.3  names
+               IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) )
+            ENDIF
+            !
+            IF( iom_use("udiff_heattr") .OR. iom_use("vdiff_heattr") ) THEN
+              !
+              IF( cdtype == 'TRA' .AND. jn == jp_tem  ) THEN
+                  z2d(:,:) = zftu(ji,jj,1) 
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 
+                        END DO
+                     END DO
+                  END DO
+                  z2d(:,:) = rau0_rcp * z2d(:,:) 
+                  CALL lbc_lnk( z2d, 'U', -1. )
+                  CALL iom_put( "udiff_heattr", z2d )                  ! heat i-transport
+                  !
+                  z2d(:,:) = zftv(ji,jj,1) 
+                  DO jk = 2, jpkm1
+                     DO jj = 2, jpjm1
+                        DO ji = fs_2, fs_jpim1   ! vector opt.
+                           z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 
+                        END DO
+                     END DO
+                  END DO
+                  z2d(:,:) = rau0_rcp * z2d(:,:)     
+                  CALL lbc_lnk( z2d, 'V', -1. )
+                  CALL iom_put( "vdiff_heattr", z2d )                  !  heat j-transport
+               ENDIF
+               !
+            ENDIF
+            !
+         ENDIF                                                    !== end pass selection  ==!
+         !
+         !                                                        ! ===============
+      END DO                                                      ! end tracer loop
+      !                                                           ! ===============
+      !
+      CALL wrk_dealloc( jpi,jpj,       z2d ) 
+      CALL wrk_dealloc( jpi,jpj,jpk,   zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw  ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_ldf_triad')
+      !
+   END SUBROUTINE tra_ldf_triad
+
+   !!==============================================================================
+END MODULE traldf_triad
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 5758)
@@ -37,5 +37,6 @@
    USE traqsr          ! penetrative solar radiation (needed for nksr)
    USE phycst          ! physical constant
-   USE ldftra_oce      ! lateral physics on tracers
+   USE ldftra          ! lateral physics on tracers
+   USE ldfslp
    USE bdy_oce         ! BDY open boundary condition variables
    USE bdytra          ! open boundary condition (bdy_tra routine)
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90	(revision 5758)
@@ -19,7 +19,10 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE dynspg_oce
+   !
+   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 ldftra_oce      ! ocean active tracers: lateral physics
+   !
    USE trd_oce         ! trends: ocean variables
    USE trdtra          ! trends manager: tracers 
@@ -45,5 +48,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -88,8 +91,10 @@
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       END SELECT
+!!gm WHY here !   and I don't like that !
       ! DRAKKAR SSS control {
       ! JMM avoid negative salinities near river outlet ! Ugly fix
       ! JMM : restore negative salinities to small salinities:
       WHERE ( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp
+!!gm
 
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
@@ -98,6 +103,8 @@
             ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal) - tsb(:,:,jk,jp_sal) ) / r2dtra(jk) ) - ztrds(:,:,jk)
          END DO
+!!gm this should be moved in trdtra.F90 and done on all trends
          CALL lbc_lnk( ztrdt, 'T', 1. )
          CALL lbc_lnk( ztrds, 'T', 1. )
+!!gm
          CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
          CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
@@ -123,10 +130,9 @@
       !!      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 or KPP scheme,
-      !!      the implicit scheme is required.
+      !!      NB: rotation of lateral mixing operator or TKE & GLS schemes,
+      !!          an implicit scheme is required.
       !!----------------------------------------------------------------------
       USE zdftke
       USE zdfgls
-      USE zdfkpp
       !!----------------------------------------------------------------------
 
@@ -137,9 +143,9 @@
 
       ! Force implicit schemes
-      IF( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp )   nzdf = 1      ! TKE, GLS 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
+      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',   &
-            &                         ' TKE or KPP scheme, the implicit scheme is required, set ln_zdfexp = .false.' )
+            &                         ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' )
 
       ! Test: esopa
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 5758)
@@ -19,22 +19,21 @@
   
    !!----------------------------------------------------------------------
-   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical  
-   !!                 part of the mixing tensor.
-   !!----------------------------------------------------------------------
-   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_oce      ! ocean active tracers: lateral physics
-   USE ldftra          ! lateral mixing type
-   USE ldfslp          ! lateral physics: slope of diffusion
-   USE zdfddm          ! ocean vertical physics: double diffusion
-   USE traldf_iso_grif ! active tracers: Griffies operator
-   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
+   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical part of the mixing tensor.
+   !!----------------------------------------------------------------------
+   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_iso_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
@@ -47,9 +46,8 @@
    !! * Substitutions
 #  include "domzgr_substitute.h90"
-#  include "ldftra_substitute.h90"
 #  include "zdfddm_substitute.h90"
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -120,32 +118,26 @@
             ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
             ENDIF
-            DO jj=1, jpj
-               DO ji=1, jpi
-                  zwt(ji,jj,1) = 0._wp
-               END DO
-            END DO
-!
-#if defined key_ldfslp
-            ! isoneutral diffusion: add the contribution 
-            IF( ln_traldf_grif    ) THEN     ! Griffies isoneutral diff
-               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)       
+            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
-               END DO
-            ELSE IF( l_traldf_rot ) THEN     ! standard isoneutral diff
-               DO jk = 2, jpkm1
-                  DO jj = 2, jpjm1
-                     DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       &
-                           &                          * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
-                           &                             + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
+               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
-               END DO
+               ENDIF
             ENDIF
-#endif
+            !
             ! Diagonal, lower (i), upper (s)  (including the bottom boundary condition since avt is masked)
             DO jk = 1, jpkm1
@@ -202,6 +194,5 @@
                ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1)
                ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1)
-               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn)                     &
-                  &                      + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
+               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
             END DO
          END DO
Index: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90
===================================================================
--- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5737)
+++ branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90	(revision 5758)
@@ -93,16 +93,15 @@
       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom)
       !
-      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
-      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points
-      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
-      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             ! 
+      INTEGER  ::   ji, jj, jn                  ! Dummy loop indices
+      INTEGER  ::   iku, ikv, ikum1, ikvm1      ! partial step level (ocean bottom level) at u- and v-points
+      REAL(wp) ::   ze3wu, ze3wv, zmaxu, zmaxv  ! local scalars
+      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
+      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj             ! 
       !!----------------------------------------------------------------------
       !
-      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde')
-      !
-      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
-      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ;
-      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ;
+      IF( nn_timing == 1 )   CALL timing_start( 'zps_hde')
+      !
+      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp
+      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -149,8 +148,8 @@
          !
       END DO
-
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
+      !                
+      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
+         pgru(:,:) = 0._wp
+         pgrv(:,:) = 0._wp                ! depth of the partial step level
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
@@ -167,12 +166,9 @@
             END DO
          END DO
-
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
+         !
+         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
+         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
+         !
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
                iku = mbku(ji,jj)
@@ -192,11 +188,12 @@
       END IF
       !
-      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde')
+      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde')
       !
    END SUBROUTINE zps_hde
-   !
-   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv,   &
-      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  &
-      &                   pgtui, pgtvi, pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
+
+
+   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi,                                   &
+      &                              prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv ,   &
+      &                                   pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE zps_hde  ***
@@ -245,18 +242,20 @@
       !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points 
       !!----------------------------------------------------------------------
-      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
-      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts 
-      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi  ! hor. grad. of stra at u- & v-pts (ISF)
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi      ! hor. grad of prd at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui, pmrvi      ! hor. sum  of prd at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui, pgzvi      ! hor. grad of z   at u- & v-pts (top)
-      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui, pge3rvi  ! hor. grad of prd weighted by local e3w at u- & v-pts (top)
+      INTEGER                              , INTENT(in   )           ::  kt                ! ocean time-step index
+      INTEGER                              , INTENT(in   )           ::  kjpt              ! number of tracers
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta               ! 4D tracers fields
+      !                                                              !!  u-point ! v-point !
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu    , pgtv    ! bottom GRADh( ptra )  
+      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui   , pgtvi   ! top    GRADh( ptra )
+      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd               ! 3D density anomaly fields
+      !                                                              !!  u-point ! v-point !
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru    , pgrv    ! bottom GRADh( prd  )
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru    , pmrv    ! bottom SUM  ( prd  )
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu    , pgzv    ! bottom GRADh( z    ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru  , pge3rv  ! bottom GRADh( prd  ) weighted by e3w
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui   , pgrvi   ! top    GRADh( prd  ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui   , pmrvi   ! top    SUM  ( prd  ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui   , pgzvi   ! top    GRADh( z    ) 
+      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui , pge3rvi ! top    GRADh( prd  ) weighted by e3w
       !
       INTEGER  ::   ji, jj, jn      ! Dummy loop indices
@@ -269,8 +268,8 @@
       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf')
       !
-      pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
-      pgtui(:,:,:)=0.0_wp ; pgtvi(:,:,:)=0.0_wp ;
-      zti (:,:,:)=0.0_wp ; ztj (:,:,:)=0.0_wp ;
-      zhi (:,:  )=0.0_wp ; zhj (:,:  )=0.0_wp ;
+      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
+      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp
+      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp
+      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp
       !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
@@ -322,11 +321,12 @@
       END DO
 
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
-         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
-         pgzu(:,:)=0.0_wp   ; pgzv(:,:)=0.0_wp ;
-         pmru(:,:)=0.0_wp   ; pmru(:,:)=0.0_wp ;
-         pge3ru(:,:)=0.0_wp ; pge3rv(:,:)=0.0_wp ;
-         DO jj = 1, jpjm1
+      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
+         !
+         pgru  (:,:)=0._wp   ;   pgrv  (:,:) = 0._wp
+         pgzu  (:,:)=0._wp   ;   pgzv  (:,:) = 0._wp 
+         pmru  (:,:)=0._wp   ;   pmru  (:,:) = 0._wp 
+         pge3ru(:,:)=0._wp   ;   pge3rv(:,:) = 0._wp 
+         !
+         DO jj = 1, jpjm1                 ! depth of the partial step level
             DO ji = 1, jpim1
                iku = mbku(ji,jj)
@@ -334,5 +334,5 @@
                ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
                ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
-
+               !
                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
@@ -343,12 +343,9 @@
             END DO
          END DO
-         
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
+         !
+         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
+         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
+
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
                iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
@@ -394,5 +391,7 @@
          !
       END IF
-         ! (ISH)  compute grui and gruvi
+      !
+      !     !==  (ISH)  compute grui and gruvi  ==!
+      !
       DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
          DO jj = 1, jpjm1
@@ -442,12 +441,12 @@
       END DO
 
-      ! horizontal derivative of density anomalies (rd)
-      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
+      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
+         !
          pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
          pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ;
          pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ;
          pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ;
-
-         DO jj = 1, jpjm1
+         !
+         DO jj = 1, jpjm1        ! depth of the partial step level
             DO ji = 1, jpim1
                iku = miku(ji,jj)
@@ -455,5 +454,5 @@
                ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
                ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
-
+               !
                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
                ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
@@ -464,12 +463,9 @@
             END DO
          END DO
-
-         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
-         ! step and store it in  zri, zrj for each  case
-         CALL eos( zti, zhi, zri )  
-         CALL eos( ztj, zhj, zrj )
-
-         ! Gradient of density at the last level 
-         DO jj = 1, jpjm1
+         !
+         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj 
+         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj 
+         !
+         DO jj = 1, jpjm1                 ! Gradient of density at the last level 
             DO ji = 1, jpim1
                iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
@@ -482,6 +478,6 @@
                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1 
                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  &
-                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
-                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
+                    &           * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
+                    &              - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
                ELSE
                  pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
@@ -489,6 +485,6 @@
                  pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
                  pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   &
-                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
-                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
+                    &           * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
+                    &              -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
                ENDIF
                IF( ze3wv >= 0._wp ) THEN
@@ -497,6 +493,6 @@
                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  & 
-                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
-                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
+                     &           * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
+                     &                - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
                                   ! + 2 due to the formulation in density and not in anomalie in hpg sco
                ELSE
@@ -505,6 +501,6 @@
                  pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
                  pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   &
-                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
-                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
+                    &           * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
+                    &              -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
                ENDIF
             END DO
@@ -517,5 +513,5 @@
       END IF  
       !
-      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_isf')
+      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf')
       !
    END SUBROUTINE zps_hde_isf
