Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2371)
@@ -422,5 +422,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =  1000.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 2371)
@@ -422,5 +422,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =  1000.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2371)
@@ -430,4 +430,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =  1000.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2371)
@@ -467,5 +467,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =  2000.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 2371)
@@ -467,5 +467,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =  2000.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2371)
@@ -467,5 +467,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =  2000.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/CONFIG/POMME/EXP00/namelist
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 2371)
@@ -467,5 +467,6 @@
    ln_traldf_hor    =  .false.  !     horizontal (geopotential)         (require "key_ldfslp" when ln_sco=T)
    ln_traldf_iso    =  .true.   !     iso-neutral                       (require "key_ldfslp")
-   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")
+   ln_traldf_grif   =  .false.  !     griffies skew flux formulation    (require "key_ldfslp")   ! UNDER TEST, DO NOT USE
+   ln_traldf_gdia   =  .false.  !     griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE
                            !  Coefficient
    rn_aht_0         =   300.    !     horizontal eddy diffusivity for tracers [m2/s]
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 2371)
@@ -87,46 +87,83 @@
       ! 1. Compute lateral diffusive coefficient 
       ! ----------------------------------------
-
-      DO jk = 1, jpk
+      IF (ln_traldf_grif) THEN
+         DO jk = 1, jpk
 #  if defined key_vectopt_loop  
-!CDIR NOVERRCHK 
-         DO ji = 1, jpij   ! vector opt.
-            ! Take the max of N^2 and zero then take the vertical sum
-            ! of the square root of the resulting N^2 ( required to compute
-            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
-            zn2 = MAX( rn2b(ji,1,jk), 0.e0 )
-            zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
-            ! Compute elements required for the inverse time scale of baroclinic
-            ! eddies using the isopycnal slopes calculated in ldfslp.F :
-            ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-            ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
+            !CDIR NOVERRCHK 
+            DO ji = 1, jpij   ! vector opt.
+               ! Take the max of N^2 and zero then take the vertical sum
+               ! of the square root of the resulting N^2 ( required to compute
+               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
+               zn2 = MAX( rn2b(ji,1,jk), 0.e0 )
+               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
+               ! Compute elements required for the inverse time scale of baroclinic
+               ! eddies using the isopycnal slopes calculated in ldfslp.F :
+               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
+               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
+               zah(ji,1) = zah(ji,1) + zn2 * wslp2(ji,1,jk) * ze3w
+               zhw(ji,1) = zhw(ji,1) + ze3w
+            END DO
+#  else
+            DO jj = 2, jpjm1
+               !CDIR NOVERRCHK 
+               DO ji = 2, jpim1
+                  ! Take the max of N^2 and zero then take the vertical sum 
+                  ! of the square root of the resulting N^2 ( required to compute 
+                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
+                  zn2 = MAX( rn2b(ji,jj,jk), 0.e0 )
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
+                  ! Compute elements required for the inverse time scale of baroclinic
+                  ! eddies using the isopycnal slopes calculated in ldfslp.F : 
+                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
+                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                  zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
+                  zhw(ji,jj) = zhw(ji,jj) + ze3w
+               END DO
+            END DO
+#  endif
+         END DO
+      ELSE
+         DO jk = 1, jpk
+#  if defined key_vectopt_loop  
+            !CDIR NOVERRCHK 
+            DO ji = 1, jpij   ! vector opt.
+               ! Take the max of N^2 and zero then take the vertical sum
+               ! of the square root of the resulting N^2 ( required to compute
+               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
+               zn2 = MAX( rn2b(ji,1,jk), 0.e0 )
+               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
+               ! Compute elements required for the inverse time scale of baroclinic
+               ! eddies using the isopycnal slopes calculated in ldfslp.F :
+               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
+               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
                zah(ji,1) = zah(ji,1) + zn2   &
-                              * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    &
-                                + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   &
-                              * ze3w
-            zhw(ji,1) = zhw(ji,1) + ze3w
-         END DO
+                    * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    &
+                    + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   &
+                    * ze3w
+               zhw(ji,1) = zhw(ji,1) + ze3w
+            END DO
 #  else
-         DO jj = 2, jpjm1
-!CDIR NOVERRCHK 
-            DO ji = 2, jpim1
-               ! Take the max of N^2 and zero then take the vertical sum 
-               ! of the square root of the resulting N^2 ( required to compute 
-               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
-               zn2 = MAX( rn2b(ji,jj,jk), 0.e0 )
-               zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
-               ! Compute elements required for the inverse time scale of baroclinic
-               ! eddies using the isopycnal slopes calculated in ldfslp.F : 
-               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
-               ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
-               zah(ji,jj) = zah(ji,jj) + zn2   &
-                              * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    &
-                                + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  &
-                              * ze3w
-               zhw(ji,jj) = zhw(ji,jj) + ze3w
-            END DO 
-         END DO 
+            DO jj = 2, jpjm1
+               !CDIR NOVERRCHK 
+               DO ji = 2, jpim1
+                  ! Take the max of N^2 and zero then take the vertical sum 
+                  ! of the square root of the resulting N^2 ( required to compute 
+                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 
+                  zn2 = MAX( rn2b(ji,jj,jk), 0.e0 )
+                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
+                  ! Compute elements required for the inverse time scale of baroclinic
+                  ! eddies using the isopycnal slopes calculated in ldfslp.F : 
+                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
+                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
+                  zah(ji,jj) = zah(ji,jj) + zn2   &
+                       * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    &
+                       + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  &
+                       * ze3w
+                  zhw(ji,jj) = zhw(ji,jj) + ze3w
+               END DO
+            END DO
 #  endif
-      END DO 
+         END DO
+      END IF
 
       DO jj = 2, jpjm1
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90	(revision 2371)
@@ -9,5 +9,5 @@
    !!   NEMO     0.5  ! 2002-10  (G. Madec)  Free form, F90
    !!            1.0  ! 2005-10  (A. Beckmann)  correction for s-coordinates
-   !!            3.3  ! 2006-10  (C. Harris, G. Nurser)  add ldf_slp_grif (Griffies operator)
+   !!            3.3  ! 2010-10  (G. Nurser, C. Harris, G. Madec)  add Griffies operator
    !!----------------------------------------------------------------------
 #if   defined key_ldfslp   ||   defined key_esopa
@@ -15,15 +15,16 @@
    !!   'key_ldfslp'                      Rotation of lateral mixing tensor
    !!----------------------------------------------------------------------
-   !!   ldf_slp      : compute the slopes of neutral surface
-   !!   ldf_slp_mxl  : compute the slopes of iso-neutral surface
+   !!   ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)
+   !!   ldf_slp      : calculates the slopes of neutral surface   (Madec operator)
+   !!   ldf_slp_mxl  : calculates the slopes at the base of the mixed layer (Madec operator)
    !!   ldf_slp_init : initialization of the slopes computation
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE ldftra_oce
-   USE ldfdyn_oce
+   USE ldftra_oce      ! lateral diffusion: traceur
+   USE ldfdyn_oce      ! lateral diffusion: dynamics
    USE phycst          ! physical constants
    USE zdfmxl          ! mixed layer depth
-   USE eosbn2
+   USE eosbn2          ! equation of states
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
@@ -33,20 +34,23 @@
    PRIVATE
 
-   PUBLIC   ldf_slp         ! routine called by step.F90
-   PUBLIC   ldf_slp_init    ! routine called by opa.F90
-   PUBLIC   ldf_slp_grif    ! routine called by step.F90
-
-   LOGICAL , PUBLIC, PARAMETER              ::   lk_ldfslp = .TRUE.   !: slopes flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   uslp, wslpi          !: i_slope at U- and W-points
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   vslp, wslpj          !: j-slope at V- and W-points
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   wslp2                !: wslp**2 from Griffies quarter cells
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   alpha, beta          !: alpha,beta at T points
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   tfw,sfw,ftu,fsu,ftv,fsv,ftud,fsud,ftvd,fsvd
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psix_eiv
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   psiy_eiv
-   
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   omlmask           ! mask of the surface mixed layer at T-pt
-   REAL(wp), DIMENSION(jpi,jpj)     ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
-   REAL(wp), DIMENSION(jpi,jpj)     ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
+   PUBLIC   ldf_slp        ! routine called by step.F90
+   PUBLIC   ldf_slp_grif   ! routine called by step.F90
+   PUBLIC   ldf_slp_init   ! routine called by opa.F90
+
+   LOGICAL , PUBLIC, PARAMETER ::   lk_ldfslp = .TRUE.     !: slopes flag
+   !                                                                             !! Madec operator
+   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   uslp, wslpi          !: i_slope at U- and W-points
+   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   vslp, wslpj          !: j-slope at V- and W-points
+   !                                                                             !! Griffies operator
+   REAL(wp), PUBLIC, DIMENSION(:,:,:)    , ALLOCATABLE ::   wslp2                !: wslp**2 from Griffies quarter cells
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE ::   triadi_g, triadj_g   !: skew flux  slopes relative to geopotentials 
+   REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE ::   triadi  , triadj     !: isoneutral slopes relative to model-coordinate
+
+   !                                                              !! Madec operator
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   omlmask           ! mask of the surface mixed layer at T-pt
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   uslpml, wslpiml   ! i_slope at U- and W-points just below the mixed layer
+   REAL(wp), DIMENSION(:,:)  , ALLOCATABLE ::   vslpml, wslpjml   ! j_slope at V- and W-points just below the mixed layer
+
+   REAL(wp) ::   repsln = 1.e-25_wp       ! tiny value used as minium of di(rho), dj(rho) and dk(rho)
 
    !! * Substitutions
@@ -67,5 +71,5 @@
       !! 
       !! ** Purpose :   Compute the slopes of neutral surface (slope of isopycnal
-      !!              surfaces referenced locally) ('key_traldfiso').
+      !!              surfaces referenced locally) (ln_traldf_iso=T).
       !! 
       !! ** Method  :   The slope in the i-direction is computed at U- and 
@@ -347,447 +351,222 @@
       !
    END SUBROUTINE ldf_slp
-
+   
 
    SUBROUTINE ldf_slp_grif ( kt )
-     !!----------------------------------------------------------------------
-     !!                 ***  ROUTINE ldf_slp_grif  ***
-     !!
-     !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope
-     !!      of iso-pycnal surfaces referenced locally) ('key_traldfiso')
-     !!      at W-points using the Griffies quarter-cells.  Also calculates
-     !!      alpha and beta at T-points for use in the Griffies isopycnal
-     !!      scheme.
-     !!
-     !! ** Method  :   The slope in the i-direction is computed at U- and
-     !!      W-points (uslp, wslpi) and the slope in the j-direction is
-     !!      computed at V- and W-points (vslp, wslpj).
-     !!
-     !! ** Action : - alpha, beta
-     !!               wslp2 squared slope of neutral surfaces at w-points.
-     !!----------------------------------------------------------------------
-     USE oce, zdit  => ua   ! use ua as workspace
-     USE oce, zdis  => va   ! use va as workspace
-     USE oce, zdjt  => ta   ! use ta as workspace
-     USE oce, zdjs  => sa   ! use sa as workspace
-     !!
-     INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
-     !!
-     INTEGER  ::   ji, jj, jk, ip, jp, kp  ! dummy loop indices
-     INTEGER  ::   iku, ikv                ! local integer
-     REAL(wp) ::   zt, zs, zh, zt2, zsp5, zp1t1           ! local scalars
-     REAL(wp) ::     zdenr, zrhotmp, zdndt, zdddt         !   -      -
-     REAL(wp) ::     zdnds, zddds, znum, zden             !   -      -
-     REAL(wp) ::     zslope, za_sxe, zslopec, zdsloper    !   -      -
-     REAL(wp) ::     zfact, zepsln, zatempw,zatempu,zatempv    !   -      -
-     REAL(wp) ::     ze1ur, ze2vr, ze3wr, zdxt, zdxs, zdyt, zdys, zdzt, zdzs, zvolf 
-     REAL(wp) ::     zr_slpmax, zdxrho, zdyrho, zabs_dzrho
-     REAL(wp), DIMENSION(jpi,jpj,jpk,0:1,0:1) ::   zsx,zsy
-     REAL(wp), DIMENSION(jpi,jpj    ,0:1,0:1) ::   zsx_ml_base, zsy_ml_base
-     REAL(wp), DIMENSION(jpi,jpj,jpk)         ::   zdkt, zdks
-     REAL(wp), DIMENSION(jpi,jpj) ::   zr_ml_basew
-     !!----------------------------------------------------------------------
-
-     ! 0. Local constant initialization
-     ! --------------------------------
-     zr_slpmax = 1.0_wp/slpmax
-
-     ! zslopec=0.004
-     ! zdsloper=1000.0
-     zepsln=1e-25
-
-     SELECT CASE ( nn_eos )
-
-     CASE ( 0 )               ! Jackett and McDougall (1994) formulation
-
-        !                                                ! ===============
-        DO jk = 1, jpk                                   ! Horizontal slab
-           !                                             ! ===============
-           DO jj = 1, jpjm1
-              DO ji = 1, fs_jpim1
-                 zt = tb(ji,jj,jk)          ! potential temperature
-                 zs = sb(ji,jj,jk) - 35.0   ! salinity anomaly (s-35)
-                 zh = fsdept(ji,jj,jk)      ! depth in meters
-
-                 beta(ji,jj,jk) = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &
-                      &                            - 0.301985e-05 ) * zt      &
-                      &   + 0.785567e-03                                      &
-                      &   + (     0.515032e-08 * zs                           &
-                      &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
-                      &   +(  (   0.121551e-17 * zh                           &
-                      &         - 0.602281e-15 * zs                           &
-                      &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
-                      &                             + 0.408195e-10   * zs     &
-                      &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
-                      &                             - 0.121555e-07 ) * zh
-
-                 alpha(ji,jj,jk) = - beta(ji,jj,jk) *                       &
-                      &     (((( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &
-                      &                               - 0.203814e-03 ) * zt   &
-                      &                               + 0.170907e-01 ) * zt   &
-                      &   + 0.665157e-01                                      &
-                      &   +     ( - 0.678662e-05 * zs                         &
-                      &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
-                      &   +   ( ( - 0.302285e-13 * zh                         &
-                      &           - 0.251520e-11 * zs                         &
-                      &           + 0.512857e-12 * zt * zt           ) * zh   &
-                      &           - 0.164759e-06 * zs                         &
-                      &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
-                      &                               + 0.380374e-04 ) * zh)
-
-              END DO
-           END DO
-        END DO
-
-     CASE ( 1 )
-
-        alpha(:,:,:)=-rn_alpha
-        beta(:,:,:)=0.0
-
-     CASE ( 2 )
-
-        alpha(:,:,:)=-rn_alpha
-        beta (:,:,:)=rn_beta
-
-     CASE ( 3 )
-
-        DO jk =1, jpk
-           DO jj = 1, jpjm1
-              DO ji = 1, fs_jpim1
-                 zt = tb(ji,jj,jk)
-                 zs = sb(ji,jj,jk)
-                 zh = fsdept(ji,jj,jk)
-                 zt2 = zt**2
-                 zsp5 = sqrt(ABS(zs))
-                 zp1t1=zh*zt
-                 znum=9.99843699e+02+zt*(7.35212840e+00+zt*(-5.45928211e-02+3.98476704e-04*zt)) &
-                      +zs*(2.96938239e+00-7.23268813e-03*zt+2.12382341e-03*zs)  &
-                      +zh*(1.04004591e-02+1.03970529e-07*zt2+5.18761880e-06*zs+ &
-                      zh*(-3.24041825e-08-1.23869360e-11*zt2))
-                 zden=1.00000000e+00+zt*(7.28606739e-03+zt*(-4.60835542e-05+zt*(3.68390573e-07+zt*1.80809186e-10))) &
-                      +zs*(2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(4.76534122e-06+1.63410736e-09*zt2)) &
-                      + zh*(5.30848875e-06+zh*zt*(-3.03175128e-16*zt2-1.27934137e-17*zh))
-                 zdenr=1.0/zden
-                 zrhotmp=znum*zdenr
-                 zdndt=7.35212840e+00+zt*(-1.091856422e-01+1.195430112e-03*zt)-7.23268813e-03*zs &
-                      +zp1t1*(2.07941058e-07-2.4773872e-11*zh)
-                 zdddt=7.28606739e-03+zt*(-9.21671084e-05+zt*(1.105171719e-06+7.23236744e-10*zt))  &
-                      +zs*(-9.27062484e-06+zt*(-5.35030929e-10*zt+3.26821472e-09*zsp5))  &
-                      +zh*zh*(-9.09525384e-16*zt2-1.27934137e-17*zh)
-                 zdnds=2.96938239e+00-7.23268813e-03*zt+2*2.12382341e-03*zs+5.18761880e-06*zh
-                 zddds=2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(7.14801183e-06+2.45116104e-09*zt2)
-                 alpha(ji,jj,jk)=(zdndt-zrhotmp*zdddt)*zdenr
-                 beta(ji,jj,jk)=zdenr*(zdnds-zrhotmp*zddds)
-
-              END DO
-           END DO
-        END DO
-
-     CASE DEFAULT
-
-        IF(lwp) WRITE(numout,cform_err)
-        IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
-        nstop = nstop + 1
-
-     END SELECT
-
-     CALL lbc_lnk( alpha, 'T', 1. )
-     CALL lbc_lnk( beta, 'T', 1. )
-
-
-     ! ---------------------------------------
-     ! 1. Horizontal tracer gradients at T-level jk
-     ! ---------------------------------------
-     DO jk = 1, jpkm1
-        DO jj = 1, jpjm1
-           DO ji = 1, fs_jpim1   ! vector opt.
-              ! i-gradient of T and S at jj
-              zdit (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk)
-              zdis (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk)
-              ! j-gradient of T and S at jj
-              zdjt (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk)
-              zdjs (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk)
-           END DO
-        END DO
-     END DO
-
-     IF( ln_zps ) THEN      ! partial steps correction at the last level
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE ldf_slp_grif  ***
+      !!
+      !! ** Purpose :   Compute the squared slopes of neutral surfaces (slope
+      !!      of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T)
+      !!      at W-points using the Griffies quarter-cells.  
+      !!
+      !! ** Method  :   calculates alpha and beta at T-points 
+      !!
+      !! ** Action : - triadi_g, triadj_g   T-pts i- and j-slope triads relative to geopot. (used for eiv)
+      !!             - triadi , triadj    T-pts i- and j-slope triads relative to model-coordinate
+      !!             - wslp2              squared slope of neutral surfaces at w-points.
+      !!----------------------------------------------------------------------
+      USE oce,   zdit  => ua   ! use ua as workspace
+      USE oce,   zdis  => va   ! use va as workspace
+      USE oce,   zdjt  => ta   ! use ta as workspace
+      USE oce,   zdjs  => sa   ! use sa as workspace
+      !!
+      INTEGER, INTENT( in ) ::   kt         ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk, jl, ip, jp, kp  ! dummy loop indices
+      INTEGER  ::   iku, ikv                ! temporary integer
+      REAL(wp) ::   zfacti, zfactj, zatempw,zatempu,zatempv   ! local scalars
+      REAL(wp) ::   zbu, zbv, zbti, zbtj
+      REAL(wp) ::   zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim
+      REAL(wp) ::   zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim
+      REAL(wp) ::   zdzrho_raw
+      REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) ::   zdzrho, zdyrho, zdxrho     ! Horizontal and vertical density gradients
+      REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) ::   zti_mlb, ztj_mlb
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     ::   zdkt, zdks
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     ::   zalpha, zbeta       ! alpha, beta at T points, at depth fsgdept
+      REAL(wp), DIMENSION(jpi,jpj)         ::   z1_mlbw
+      !!----------------------------------------------------------------------
+
+      !--------------------------------!
+      !  Some preliminary calculation  !
+      !--------------------------------!
+      !
+      CALL eos_alpbet( tsb, zalpha, zbeta )     !==  before thermal and haline expension coeff. at T-points  ==!
+      !
+      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) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk)    ! i-gradient of T and S at jj
+               zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk)
+               zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk)    ! j-gradient of T and S at jj
+               zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      IF( ln_zps ) THEN                               ! partial steps: correction at the last level
 # if defined key_vectopt_loop
-     DO jj = 1, 1
-        DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+         DO jj = 1, 1
+            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
 # else
-           DO jj = 1, jpjm1
-              DO ji = 1, jpim1
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
 # endif
-                 ! last ocean level
-                 iku  = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
-                 ikv  = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
-                 ! i-gradient of T and S
-                 zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem)
-                 zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal)
-                 ! j-gradient of T and S
-                 zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem)
-                 zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal)
-              END DO
-           END DO
-        ENDIF
-
-        ! ---------------------------------------
-        ! 1. Vertical tracer gradient at w-level jk
-        ! ---------------------------------------
-        DO jk = 2, jpk
-           zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk)
-           zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk)
-        END DO
-
-        zdkt(:,:,1) = 0.0
-        zdks(:,:,1) = 0.0
-
-        ! ---------------------------------------
-        ! Depth of the w-point below ML base
-        ! ---------------------------------------
-        DO jj = 1, jpj
-           DO ji = 1, jpi
-              jk = nmln(ji,jj)
-              zr_ml_basew(ji,jj)=1.0/fsdepw(ji,jj,jk+1)
-           END DO
-        END DO
-
-
-        wslp2(:,:,:)=0.0
-        tfw(:,:,:) = 0.0
-        sfw(:,:,:) = 0.0
-        ftu(:,:,:) = 0.0
-        fsu(:,:,:) = 0.0
-        ftv(:,:,:) = 0.0
-        fsv(:,:,:) = 0.0
-        ftud(:,:,:) = 0.0
-        fsud(:,:,:) = 0.0
-        ftvd(:,:,:) = 0.0
-        fsvd(:,:,:) = 0.0
-        psix_eiv(:,:,:) = 0.0
-        psiy_eiv(:,:,:) = 0.0
-
-        ! ----------------------------------------------------------------------
-        ! x-z plane
-        ! ----------------------------------------------------------------------
-
-        ! calculate limited triad x-slopes zsx in interior (1=<jk=<jpk-1)
-        DO ip=0,1
-           DO kp=0,1
-
-              DO jk = 1, jpkm1
-                 DO jj = 1, jpjm1
-                    DO ji = 1, fs_jpim1   ! vector opt.
-
-                       ze1ur=1.0/e1u(ji,jj)
-                       zdxt=zdit(ji,jj,jk)*ze1ur
-                       zdxs=zdis(ji,jj,jk)*ze1ur
-
-                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp)
-                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr
-                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr
-                       ! Calculate the horizontal and vertical density gradient
-                       zdxrho = alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs
-                       zabs_dzrho = ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)+zepsln
-                       ! Limit by slpmax, and mask by psi-point
-                       zsx(ji+ip,jj,jk,1-ip,kp) = umask(ji,jj,jk+kp) &
-                            & *zdxrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdxrho))
-                    END DO
-                 END DO
-              END DO
-
-           END DO
-        END DO
-
-        ! calculate slope of triad immediately below mixed-layer base
-        DO ip=0,1
-           DO kp=0,1
-              DO jj = 1, jpjm1
-                 DO ji = 1, fs_jpim1
-                    jk = nmln(ji+ip,jj)
-                    zsx_ml_base(ji+ip,jj,1-ip,kp)=zsx(ji+ip,jj,jk+1-kp,1-ip,kp)
-                 END DO
-              END DO
-           END DO
-        END DO
-
-        ! Below ML use limited zsx as is
-        ! Inside ML replace by linearly reducing zsx_ml_base towards surface
-        DO ip=0,1
-           DO kp=0,1
-
-              DO jk = 1, jpkm1
-
-                 DO jj = 1, jpjm1
-
-                    DO ji = 1, fs_jpim1   ! vector opt.
-                       ! k index of uppermost point(s) of triad is jk+kp-1
-                       ! must be .ge. nmln(ji,jj) for zfact=1.
-                       !                   otherwise  zfact=0.
-                       zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj))
-                       zsx(ji+ip,jj,jk,1-ip,kp) = zfact*zsx(ji+ip,jj,jk,1-ip,kp) + &
-                            & (1.0_wp-zfact)*(fsdepw(ji+ip,jj,jk+kp)*zr_ml_basew(ji+ip,jj))*zsx_ml_base(ji+ip,jj,1-ip,kp) 
-                    END DO
-
-                 END DO
-
-              END DO
-           END DO
-        END DO
-
-        ! Use zsx to calculate fluxes and save averaged slopes psix_eiv at psi-points
-        DO ip=0,1
-           DO kp=0,1
-
-              DO jk = 1, jpkm1
-
-                 DO jj = 1, jpjm1
-
-                    DO ji = 1, fs_jpim1   ! vector opt.
-
-                       ze1ur=1.0/e1u(ji,jj)
-                       zdxt=zdit(ji,jj,jk)*ze1ur
-                       zdxs=zdis(ji,jj,jk)*ze1ur
-
-                       ze3wr=1.0/fse3w(ji+ip,jj,jk+kp)
-                       zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr
-                       zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr
-                       zslope=zsx(ji+ip,jj,jk,1-ip,kp)
-
-                       zvolf = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk)
-
-                       ftu(ji,jj,jk)= ftu(ji,jj,jk)+zslope*zdzt*zvolf*ze1ur
-                       fsu(ji,jj,jk)= fsu(ji,jj,jk)+zslope*zdzs*zvolf*ze1ur
-                       ftud(ji,jj,jk)=ftud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxt*zvolf*ze1ur
-                       fsud(ji,jj,jk)=fsud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxs*zvolf*ze1ur
-                       tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxt
-                       sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxs
-                       wslp2(ji+ip,jj,jk+kp)=wslp2(ji+ip,jj,jk+kp)+ &
-                            & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zslope)**2
-                       psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zslope
-
-                    END DO
-                 END DO
-
-              END DO
-           END DO
-        END DO
-
-        ! ----------------------------------------------------------------------
-        ! y-z plane
-        ! ----------------------------------------------------------------------
-
-        ! calculate limited triad y-slopes zsy in interior (1=<jk=<jpk-1)
-        DO jp=0,1
-           DO kp=0,1
-
-              DO jk = 1, jpkm1
-                 DO jj = 1, jpjm1
-                    DO ji = 1, fs_jpim1   ! vector opt.
-
-                       ze2vr=1.0/e2v(ji,jj)
-                       zdyt=zdjt(ji,jj,jk)*ze2vr
-                       zdys=zdjs(ji,jj,jk)*ze2vr
-
-                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp)
-                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr
-                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr
-                       ! Calculate the horizontal and vertical density gradient
-                       zdyrho = alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys
-                       zabs_dzrho = ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)+zepsln
-                       ! Limit by slpmax, and mask by psi-point
-                       zsy(ji,jj+jp,jk,1-jp,kp) = vmask(ji,jj,jk+kp) &
-                            & *zdyrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdyrho))
-                    END DO
-                 END DO
-              END DO
-
-           END DO
-        END DO
-
-        ! calculate slope of triad immediately below mixed-layer base
-        DO jp=0,1
-           DO kp=0,1
-              DO jj = 1, jpjm1
-                 DO ji = 1, fs_jpim1
-                    jk = nmln(ji,jj+jp)
-                    zsy_ml_base(ji,jj+jp,1-jp,kp)=zsy(ji,jj+jp,jk+1-kp,1-jp,kp)
-                 END DO
-              END DO
-           END DO
-        END DO
-
-        ! Below ML use limited zsy as is
-        ! Inside ML replace by linearly reducing zsy_ml_base towards surface
-        DO jp=0,1
-           DO kp=0,1
-              DO jk = 1, jpkm1
-                 DO jj = 1, jpjm1
-                    DO ji = 1, fs_jpim1   ! vector opt.
-                       ! k index of uppermost point(s) of triad is jk+kp-1
-                       ! must be .ge. nmln(ji,jj) for zfact=1.
-                       !                   otherwise  zfact=0.
-                       zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp))
-                       zsy(ji,jj+jp,jk,1-jp,kp) = zfact*zsy(ji,jj+jp,jk,1-jp,kp) + &
-                           & (1.0_wp-zfact)*(fsdepw(ji,jj+jp,jk+kp)*zr_ml_basew(ji,jj+jp))*zsy_ml_base(ji,jj+jp,1-jp,kp) 
-                    END DO
-
-                 END DO
-
-              END DO
-           END DO
-        END DO
-
-        ! Use zsy to calculate fluxes and save averaged slopes psiy_eiv at psi-points
-        DO jp=0,1
-           DO kp=0,1
-              DO jk = 1, jpkm1
-                 DO jj = 1, jpjm1
-                    DO ji = 1, fs_jpim1   ! vector opt.
-
-                       ze2vr=1.0/e2v(ji,jj)
-                       zdyt=zdjt(ji,jj,jk)*ze2vr
-                       zdys=zdjs(ji,jj,jk)*ze2vr
-
-                       ze3wr=1.0/fse3w(ji,jj+jp,jk+kp)
-                       zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr
-                       zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr
-                       zslope=zsy(ji,jj+jp,jk,1-jp,kp)
-
-                       zvolf = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk)
-
-                       ftv(ji,jj,jk)= ftv(ji,jj,jk)+zslope*zdzt*zvolf*ze2vr
-                       fsv(ji,jj,jk)= fsv(ji,jj,jk)+zslope*zdzs*zvolf*ze2vr
-                       ftvd(ji,jj,jk)=ftvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdyt*zvolf*ze2vr
-                       fsvd(ji,jj,jk)=fsvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdys*zvolf*ze2vr
-                       tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdyt
-                       sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdys
-                       wslp2(ji,jj+jp,jk+kp)=wslp2(ji,jj+jp,jk+kp)+ &
-                            & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zslope)**2
-                       psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zslope
-
-                    END DO
-                 END DO
-              END DO
-           END DO
-        END DO
-
-        tfw(:,:,1)=0.e0
-        sfw(:,:,1)=0.e0
-        wslp2(:,:,1)=0.e0
-
-        CALL lbc_lnk( wslp2, 'W', 1. )
-        CALL lbc_lnk( tfw, 'W', 1. )
-        CALL lbc_lnk( sfw, 'W', 1. )
-        CALL lbc_lnk( ftu, 'U', -1. )
-        CALL lbc_lnk( fsu, 'U', -1. )
-        CALL lbc_lnk( ftv, 'V', -1. )
-        CALL lbc_lnk( fsv, 'V', -1. )
-        CALL lbc_lnk( ftud, 'U', -1. )
-        CALL lbc_lnk( fsud, 'U', -1. )
-        CALL lbc_lnk( ftvd, 'V', -1. )
-        CALL lbc_lnk( fsvd, 'V', -1. )
-        CALL lbc_lnk( psix_eiv, 'U', -1. )
-        CALL lbc_lnk( psiy_eiv, 'V', -1. )
+               iku = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1, 2 )    ! last ocean level
+               ikv = MAX( MIN( mbathy(ji,jj), mbathy(ji,jj+1  ) ) - 1, 2 )
+               zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem)                          ! i-gradient of T and S
+               zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal)
+               zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem)                          ! j-gradient of T and S
+               zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal)
+            END DO
+         END DO
+      ENDIF
+      !
+      zdkt(:,:,1) = 0._wp                    !==  before vertical T & S gradient at w-level  ==!
+      zdks(:,:,1) = 0._wp
+      DO jk = 2, jpk
+         zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk)
+         zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk)
+      END DO
+      !
+      !
+      DO jl = 0, 1                           !==  density i-, j-, and k-gradients  ==!
+         ip = jl   ;   jp = jl         ! guaranteed nonzero gradients ( absolute value larger than repsln)
+         DO jk = 1, jpkm1                          ! done each pair of triad
+            DO jj = 1, jpjm1                       ! NB: not masked due to the minimum value set
+               DO ji = 1, fs_jpim1   ! vector opt. 
+                  zdxrho_raw = ( zalpha(ji+ip,jj   ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj   ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj)
+                  zdyrho_raw = ( zalpha(ji   ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji   ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj)
+                  zdxrho(ji+ip,jj   ,jk,1-ip) = SIGN( MAX(   repsln, ABS( zdxrho_raw ) ), zdxrho_raw )    ! keep the sign
+                  zdyrho(ji   ,jj+jp,jk,1-jp) = SIGN( MAX(   repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
+               END DO
+            END DO
+         END DO
+      END DO
+     DO kp = 0, 1                           !==  density i-, j-, and k-gradients  ==!
+         DO jk = 1, jpkm1                          ! done each pair of triad
+            DO jj = 1, jpj                       ! NB: not masked due to the minimum value set
+               DO ji = 1, jpi   ! vector opt. 
+                   zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) )   &
+                     &       / fse3w(ji,jj,jk+kp)
+                  zdzrho(ji   ,jj   ,jk,  kp) =     - MIN( - repsln,      zdzrho_raw )                    ! force zdzrho >= repsln
+               END DO
+            END DO
+         END DO
+      END DO
+      !
+      DO jj = 1, jpj                         !==  Reciprocal depth of the w-point below ML base  ==!
+         DO ji = 1, jpi                            ! i.e. 1 / (hmld+e3t(nmln))  where hmld=depw(nmln)
+            jk = MIN( nmln(ji,jj), mbathy(ji,jj) - 1 ) + 1     ! MIN in case ML depth is the ocean depth
+            z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)
+         END DO
+      END DO
+      !
+      !                                      !==  intialisations to zero  ==!
+      !
+      wslp2  (:,:,:)     = 0._wp                                           ! wslp2 will be cumulated 3D field set to zero
+      triadi_g(:,:,1,:,:) = 0._wp   ;   triadi_g(:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero
+      triadj_g(:,:,1,:,:) = 0._wp   ;   triadj_g(:,:,jpk,:,:) = 0._wp
+!!gm _iso set to zero missing
+      triadi (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp        ! set surface and bottom slope to zero
+      triadj (:,:,1,:,:) = 0._wp   ;   triadj (:,:,jpk,:,:) = 0._wp
+      
+      !-------------------------------------!
+      !  Triads just below the Mixed Layer  !
+      !-------------------------------------!
+      !
+      DO jl = 0, 1               ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
+         DO kp = 0, 1            ! with only the slope-max limit   and   MASKED 
+            DO jj = 1, jpjm1
+               DO ji = 1, fs_jpim1
+                  ip = jl   ;   jp = jl
+                  jk = MIN( nmln(ji+ip,jj), mbathy(ji+ip,jj) - 1 ) + 1     ! ML level+1 (MIN in case ML depth is the ocean depth)
+                  zti_g_raw = (  zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp)      &
+                     &      + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj)  ) * umask(ji,jj,jk)
+                  jk = MIN( nmln(ji,jj+jp), mbathy(ji,jj+jp) - 1 ) + 1
+                  ztj_g_raw = (  zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp)      &
+                     &      + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj)  ) * vmask(ji,jj,jk)
+                  zti_mlb(ji+ip,jj   ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
+                  ztj_mlb(ji   ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
+               END DO
+            END DO
+         END DO
+      END DO
+
+      !-------------------------------------!
+      !  Triads with surface limits         !
+      !-------------------------------------!
+      !
+      DO kp = 0, 1                        ! k-index of triads
+         DO jl = 0, 1
+            ip = jl   ;   jp = jl         ! i- and j-indices of triads (i-k and j-k planes)
+            DO jk = 1, jpkm1
+               DO jj = 1, jpjm1
+                  DO ji = 1, fs_jpim1   ! vector opt.
+                     !
+                     ! Calculate slope relative to geopotentials used for GM skew fluxes
+                     ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth)
+                     ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point
+                     ! masked by umask taken at the level of dz(rho)
+                     !
+                     ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti)
+                     !
+                     zti_raw   = zdxrho(ji+ip,jj   ,jk,1-ip) / zdzrho(ji+ip,jj   ,jk,kp)                   ! unmasked
+                     ztj_raw   = zdyrho(ji   ,jj+jp,jk,1-jp) / zdzrho(ji   ,jj+jp,jk,kp)
+                     zti_coord = ( fsdept(ji+1,jj  ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
+                     ztj_coord = ( fsdept(ji  ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
+                  ! unmasked
+                     zti_g_raw = zti_raw + zti_coord      ! ref to geopot surfaces
+                     ztj_g_raw = ztj_raw + ztj_coord
+                     zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
+                     ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
+                     !
+                     ! Below  ML use limited zti_g as is
+                     ! Inside ML replace by linearly reducing sx_mlb towards surface
+                     !
+                     zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp )  ! k index of uppermost point(s) of triad is jk+kp-1
+                     zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp )  ! must be .ge. nmln(ji,jj) for zfact=1
+                     !                                                          !                   otherwise  zfact=0
+                     zti_g_lim =           zfacti   * zti_g_lim                       &
+                        &      + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp)   &
+                        &                           * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj)
+                     ztj_g_lim =           zfactj   * ztj_g_lim                       &
+                        &      + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp)   &
+                        &                           * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp)                   ! masked
+                     !
+                     triadi_g(ji+ip,jj   ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp)
+                     triadj_g(ji   ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp)
+                     !
+                     ! Get coefficients of isoneutral diffusion tensor
+                     ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients)
+                     ! 2. We require that isoneutral diffusion  gives no vertical buoyancy flux
+                     !     i.e. 33 term = (real slope* 31, 13 terms)
+                     ! To do this, retain limited sx**2  in vertical flux, but divide by real slope for 13/31 terms
+                     ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2
+                     !
+                     zti_lim  = zti_g_lim - zti_coord    ! remove the coordinate slope  ==> relative to coordinate surfaces
+                     ztj_lim  = ztj_g_lim - ztj_coord                                 
+                     zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp)      ! square of limited slopes            ! masked <<==
+                     ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp)
+                     !
+                     zbu = e1u(ji    ,jj) * e2u(ji   ,jj) * fse3u(ji   ,jj,jk   )
+                     zbv = e1v(ji    ,jj) * e2v(ji   ,jj) * fse3v(ji   ,jj,jk   )
+                     zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp)
+                     zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp)
+                     !
+                     triadi(ji+ip,jj   ,jk,1-ip,kp) = zti_lim2 / zti_raw                                          ! masked
+                     triadj(ji   ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw
+                     !
+!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers....  ==> to be checked
+                     wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2             ! masked
+                     wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2
+                  END DO
+               END DO
+            END DO
+         END DO
+      END DO
+      !
+      wslp2(:,:,1) = 0._wp                ! force the surface wslp to zero
+      
+      CALL lbc_lnk( wslp2, 'W', 1. )      ! lateral boundary confition on wslp2 only   ==>>> gm : necessary ? to be checked
       !
    END SUBROUTINE ldf_slp_grif
@@ -1008,52 +787,67 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk   ! dummy loop indices
+      INTEGER ::   ierr         ! local integer
       !!----------------------------------------------------------------------
       
       IF(lwp) THEN    
          WRITE(numout,*)
-         WRITE(numout,*) 'ldf_slp : direction of lateral mixing'
-         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing'
+         WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
-
-      ! Direction of lateral diffusion (tracers and/or momentum)
-      ! ------------------------------
-      ! set the slope to zero (even in s-coordinates)
-
-      uslp (:,:,:) = 0.e0
-      vslp (:,:,:) = 0.e0
-      wslpi(:,:,:) = 0.e0
-      wslpj(:,:,:) = 0.e0
-
-      uslpml (:,:) = 0.e0
-      vslpml (:,:) = 0.e0
-      wslpiml(:,:) = 0.e0
-      wslpjml(:,:) = 0.e0
-
-      IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
-         IF(lwp) THEN
-            WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
+      
+      IF( ln_traldf_grif ) THEN        ! Griffies operator : triad of slopes
+         ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr )
+         ALLOCATE( triadi  (jpi,jpj,jpk,0:1,0:1) , triadj  (jpi,jpj,jpk,0:1,0:1)                      , STAT=ierr )
+         IF( ierr > 0 ) THEN
+            CALL ctl_stop( 'ldf_slp_init : unable to allocate Griffies operator slope ' )   ;   RETURN
          ENDIF
-
-         ! geopotential diffusion in s-coordinates on tracers and/or momentum
-         ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
-         ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
-
-         ! set the slope of diffusion to the slope of s-surfaces
-         !      ( c a u t i o n : minus sign as fsdep has positive value )
-         DO jk = 1, jpk
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1   ! vector opt.
-                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
-                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
-                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
-                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+         !
+         IF( ln_dynldf_iso )   CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
+         !
+         IF( ( ln_traldf_hor .AND. ln_dynldf_hor ) .AND. ln_sco )   &
+            &     CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator ',   &
+            &                                              'in s-coordinate not supported' )
+         !
+      ELSE                             ! Madec operator : slopes at u-, v-, and w-points
+         ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) ,                &
+            &   omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj)   , vslpml(jpi,jpj)    , wslpiml(jpi,jpj)   , wslpjml(jpi,jpj) , STAT=ierr )
+         IF( ierr > 0 ) THEN
+            CALL ctl_stop( 'ldf_slp_init : unable to allocate Madec operator slope ' )   ;   RETURN
+         ENDIF
+
+         ! Direction of lateral diffusion (tracers and/or momentum)
+         ! ------------------------------
+         uslp (:,:,:) = 0._wp   ;   uslpml (:,:) = 0._wp      ! set the slope to zero (even in s-coordinates)
+         vslp (:,:,:) = 0._wp   ;   vslpml (:,:) = 0._wp
+         wslpi(:,:,:) = 0._wp   ;   wslpiml(:,:) = 0._wp
+         wslpj(:,:,:) = 0._wp   ;   wslpjml(:,:) = 0._wp
+
+!!gm I no longer understand this.....
+         IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN
+            IF(lwp) THEN
+               WRITE(numout,*) '          Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
+            ENDIF
+
+            ! geopotential diffusion in s-coordinates on tracers and/or momentum
+            ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
+            ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
+
+            ! set the slope of diffusion to the slope of s-surfaces
+            !      ( c a u t i o n : minus sign as fsdep has positive value )
+            DO jk = 1, jpk
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1   ! vector opt.
+                     uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
+                     vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
+                     wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
+                     wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
+                  END DO
                END DO
             END DO
-         END DO
-         ! Lateral boundary conditions on the slopes
-         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
-         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
-      ENDIF
-      !
+            ! Lateral boundary conditions on the slopes
+            CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
+            CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
+         ENDIF
+      ENDIF      !
    END SUBROUTINE ldf_slp_init
 
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90	(revision 2371)
@@ -67,5 +67,5 @@
       NAMELIST/namtra_ldf/ ln_traldf_lap  , ln_traldf_bilap,                  &
          &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   &
-         &                 ln_traldf_grif ,                                   &
+         &                 ln_traldf_grif , ln_traldf_gdia,                   &
          &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0,       &
          &                 rn_slpmax
@@ -82,15 +82,22 @@
          WRITE(numout,*) 'ldf_tra_init : lateral tracer physics'
          WRITE(numout,*) '~~~~~~~~~~~~ '
-         WRITE(numout,*) '   Namelist namtra_ldf : lateral mixing coefficients'
+         WRITE(numout,*) '   Namelist namtra_ldf : lateral mixing parameters (type, direction, coefficients)'
          WRITE(numout,*) '      laplacian operator            ln_traldf_lap   = ', ln_traldf_lap
          WRITE(numout,*) '      bilaplacian operator          ln_traldf_bilap = ', ln_traldf_bilap
-         WRITE(numout,*) '      griffies    operator          ln_traldf_grif  = ', ln_traldf_grif 
+         WRITE(numout,*) '      iso-level                     ln_traldf_level = ', ln_traldf_level
+         WRITE(numout,*) '      horizontal (geopotential)     ln_traldf_hor   = ', ln_traldf_hor
+         WRITE(numout,*) '      iso-neutral                   ln_traldf_iso   = ', ln_traldf_iso
+         WRITE(numout,*) '      iso-neutral (Griffies)        ln_traldf_grif  = ', ln_traldf_grif
+         WRITE(numout,*) '      Griffies strmfn diagnostics   ln_traldf_gdia  = ', ln_traldf_gdia
          WRITE(numout,*) '      lateral eddy diffusivity      rn_aht_0        = ', rn_aht_0
          WRITE(numout,*) '      background hor. diffusivity   rn_ahtb_0       = ', rn_ahtb_0
          WRITE(numout,*) '      eddy induced velocity coef.   rn_aeiv_0       = ', rn_aeiv_0
+         WRITE(numout,*) '      maximum isoppycnal slope      rn_slpmax       = ', rn_slpmax
+         WRITE(numout,*) '   + griffies operator internal controls not set via the namelist (experimental): '
+         WRITE(numout,*) '      calculate triads twice        l_triad_iso     = ', l_triad_iso
+         WRITE(numout,*) '      no Shapiro filter             l_no_smooth     = ', l_no_smooth
          WRITE(numout,*)
       ENDIF
 
-      slpmax = rn_slpmax
       !                                ! convert DOCTOR namelist names into OLD names
       aht0  = rn_aht_0
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90	(revision 2371)
@@ -21,4 +21,5 @@
    LOGICAL , PUBLIC ::   ln_traldf_iso   = .TRUE.    !: iso-neutral direction
    LOGICAL , PUBLIC ::   ln_traldf_grif  = .FALSE.   !: griffies skew flux
+   LOGICAL , PUBLIC ::   ln_traldf_gdia  = .FALSE.   !: griffies skew flux streamfunction diagnostics
    REAL(wp), PUBLIC ::   rn_aht_0        = 2000._wp  !: lateral eddy diffusivity (m2/s)
    REAL(wp), PUBLIC ::   rn_ahtb_0       =    0._wp  !: lateral background eddy diffusivity (m2/s)
@@ -27,5 +28,6 @@
 
    REAL(wp), PUBLIC ::   aht0, ahtb0, aeiv0         !!: OLD namelist names
-   REAL(wp), PUBLIC ::   slpmax                     !: slope limit 
+   LOGICAL , PUBLIC ::   l_triad_iso     = .FALSE.   !: calculate triads twice
+   LOGICAL , PUBLIC ::   l_no_smooth     = .FALSE.   !: no Shapiro smoothing
 
 #if defined key_traldf_c3d
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 2371)
@@ -17,4 +17,5 @@
    !!            3.0  ! 2006-08  (G. Madec)  add tfreez function
    !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!             -   ! 2010-10  (G. Nurser, G. Madec)  add eos_alpbet used in ldfslp
    !!----------------------------------------------------------------------
 
@@ -26,4 +27,5 @@
    !!   eos_insitu_2d  : Compute the in situ density for 2d fields
    !!   eos_bn2        : Compute the Brunt-Vaisala frequency
+   !!   eos_alpbet     : calculates the in situ thermal and haline expansion coeff.
    !!   tfreez         : Compute the surface freezing temperature
    !!   eos_init       : set eos parameters (namelist)
@@ -38,5 +40,5 @@
    PRIVATE
 
-   !! * Interface 
+   !                   !! * Interface 
    INTERFACE eos
       MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
@@ -49,12 +51,13 @@
    PUBLIC   eos_init   ! called by istate module
    PUBLIC   bn2        ! called by step module
+   PUBLIC   eos_alpbet ! called by ldfslp module
    PUBLIC   tfreez     ! called by sbcice_... modules
 
-   !                                        !!* Namelist (nameos) *
-   INTEGER , PUBLIC ::   nn_eos   = 0        !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
-   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4   !: thermal expension coeff. (linear equation of state)
-   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4   !: saline  expension coeff. (linear equation of state)
-
-   REAL(wp), PUBLIC ::   ralpbet           !: alpha / beta ratio
+   !                                          !!* Namelist (nameos) *
+   INTEGER , PUBLIC ::   nn_eos   = 0         !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
+   REAL(wp), PUBLIC ::   rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state)
+   REAL(wp), PUBLIC ::   rn_beta  = 7.7e-4_wp !: saline  expension coeff. (linear equation of state)
+
+   REAL(wp), PUBLIC ::   ralpbet              !: alpha / beta ratio
    
    !! * Substitutions
@@ -64,7 +67,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -135,11 +137,11 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4)*zt   &
-                     &      -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594
+                  zr1= ( ( ( ( 6.536332e-9_wp  *zt - 1.120083e-6_wp )*zt + 1.001685e-4_wp )*zt   &
+                     &        -9.095290e-3_wp )*zt + 6.793952e-2_wp )*zt +  999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt        &
-                     &                   -4.0899e-3 ) *zt+0.824493
-                  zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3
-                  zr4= 4.8314e-4
+                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt        &
+                     &                      -4.0899e-3_wp ) *zt+0.824493_wp
+                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zr4= 4.8314e-4_wp
                   !
                   ! potential volumic mass (reference to the surface)
@@ -147,20 +149,20 @@
                   !
                   ! add the compression terms
-                  ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
-                  zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
+                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
                   zb = zbw + ze * zs
                   !
-                  zd = -2.042967e-2
-                  zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
-                  zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788
+                  zd = -2.042967e-2_wp
+                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
+                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
                   za = ( zd*zsr + zc ) *zs + zaw
                   !
-                  zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545
-                  za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077
-                  zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6
+                  zb1=   (-0.1909078_wp*zt+7.390729_wp )        *zt-55.87545_wp
+                  za1= ( ( 2.326469e-3_wp*zt+1.553190_wp)       *zt-65.00517_wp ) *zt+1044.077_wp
+                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
                   zk0= ( zb1*zsr + za1 )*zs + zkw
                   !
                   ! masked in situ density anomaly
-                  prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
+                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
                      &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
                END DO
@@ -170,5 +172,5 @@
       CASE( 1 )                !==  Linear formulation function of temperature only  ==!
          DO jk = 1, jpkm1
-            prd(:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
+            prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
          END DO
          !
@@ -258,11 +260,11 @@
                   !
                   ! compute volumic mass pure water at atm pressure
-                  zr1= ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   &
-                     &                       -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594
+                  zr1= ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
+                     &                          -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
                   ! seawater volumic mass atm pressure
-                  zr2= ( ( ( 5.3875e-9*zt-8.2467e-7 ) *zt+7.6438e-5 ) *zt   &
-                     &                                   -4.0899e-3 ) *zt+0.824493
-                  zr3= ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3
-                  zr4= 4.8314e-4
+                  zr2= ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp ) *zt+7.6438e-5_wp ) *zt   &
+                     &                                         -4.0899e-3_wp ) *zt+0.824493_wp
+                  zr3= ( -1.6546e-6_wp*zt+1.0227e-4_wp )    *zt-5.72466e-3_wp
+                  zr4= 4.8314e-4_wp
                   !
                   ! potential volumic mass (reference to the surface)
@@ -273,20 +275,20 @@
                   !
                   ! add the compression terms
-                  ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
-                  zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
+                  ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+                  zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
                   zb = zbw + ze * zs
                   !
-                  zd = -2.042967e-2
-                  zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
-                  zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt - 4.721788
+                  zd = -2.042967e-2_wp
+                  zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
+                  zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt - 4.721788_wp
                   za = ( zd*zsr + zc ) *zs + zaw
                   !
-                  zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545
-                  za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077
-                  zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt + 2098.925 ) *zt+190925.6
+                  zb1=   (  -0.1909078_wp  *zt+7.390729_wp    ) *zt-55.87545_wp
+                  za1= ( (   2.326469e-3_wp*zt+1.553190_wp    ) *zt-65.00517_wp ) *zt + 1044.077_wp
+                  zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp ) *zt-30.41638_wp ) *zt + 2098.925_wp ) *zt+190925.6_wp
                   zk0= ( zb1*zsr + za1 )*zs + zkw
                   !
                   ! masked in situ density anomaly
-                  prd(ji,jj,jk) = (  zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
+                  prd(ji,jj,jk) = (  zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  )    &
                      &             - rau0  ) * zrau0r * tmask(ji,jj,jk)
                END DO
@@ -296,6 +298,6 @@
       CASE( 1 )                !==  Linear formulation = F( temperature )  ==!
          DO jk = 1, jpkm1
-            prd  (:,:,jk) = ( 0.0285 - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk)
+            prd  (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_sal) )        * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0_wp   +            prd (:,:,jk)       ) * rau0 * tmask(:,:,jk)
          END DO
          !
@@ -303,5 +305,5 @@
          DO jk = 1, jpkm1
             prd  (:,:,jk) = ( rn_beta  * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) )        * tmask(:,:,jk)
-            prhop(:,:,jk) = ( 1.e0     + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk)
+            prhop(:,:,jk) = ( 1.e0_wp  + prd (:,:,jk) )                                       * rau0 * tmask(:,:,jk)
          END DO
          !
@@ -382,11 +384,11 @@
                !
                ! compute volumic mass pure water at atm pressure
-               zr1 = ( ( ( ( 6.536332e-9*zt-1.120083e-6 )*zt+1.001685e-4 )*zt   &
-                  &                        -9.095290e-3 )*zt+6.793952e-2 )*zt+999.842594
+               zr1 = ( ( ( ( 6.536332e-9_wp*zt-1.120083e-6_wp )*zt+1.001685e-4_wp )*zt   &
+                  &                        -9.095290e-3_wp )*zt+6.793952e-2_wp )*zt+999.842594_wp
                ! seawater volumic mass atm pressure
-               zr2 = ( ( ( 5.3875e-9*zt-8.2467e-7 )*zt+7.6438e-5 ) *zt   &
-                  &                                   -4.0899e-3 ) *zt+0.824493
-               zr3 = ( -1.6546e-6*zt+1.0227e-4 ) *zt-5.72466e-3
-               zr4 = 4.8314e-4
+               zr2 = ( ( ( 5.3875e-9_wp*zt-8.2467e-7_wp )*zt+7.6438e-5_wp ) *zt   &
+                  &                                   -4.0899e-3_wp ) *zt+0.824493_wp
+               zr3 = ( -1.6546e-6_wp*zt+1.0227e-4_wp ) *zt-5.72466e-3_wp
+               zr4 = 4.8314e-4_wp
                !
                ! potential volumic mass (reference to the surface)
@@ -394,21 +396,21 @@
                !
                ! add the compression terms
-               ze = ( -3.508914e-8*zt-1.248266e-8 ) *zt-2.595994e-6
-               zbw= (  1.296821e-6*zt-5.782165e-9 ) *zt+1.045941e-4
+               ze = ( -3.508914e-8_wp*zt-1.248266e-8_wp ) *zt-2.595994e-6_wp
+               zbw= (  1.296821e-6_wp*zt-5.782165e-9_wp ) *zt+1.045941e-4_wp
                zb = zbw + ze * zs
                !
-               zd = -2.042967e-2
-               zc =   (-7.267926e-5*zt+2.598241e-3 ) *zt+0.1571896
-               zaw= ( ( 5.939910e-6*zt+2.512549e-3 ) *zt-0.1028859 ) *zt -4.721788
+               zd =    -2.042967e-2_wp
+               zc =   (-7.267926e-5_wp*zt+2.598241e-3_wp ) *zt+0.1571896_wp
+               zaw= ( ( 5.939910e-6_wp*zt+2.512549e-3_wp ) *zt-0.1028859_wp ) *zt -4.721788_wp
                za = ( zd*zsr + zc ) *zs + zaw
                !
-               zb1=   (-0.1909078*zt+7.390729 ) *zt-55.87545
-               za1= ( ( 2.326469e-3*zt+1.553190)*zt-65.00517 ) *zt+1044.077
-               zkw= ( ( (-1.361629e-4*zt-1.852732e-2 ) *zt-30.41638 ) *zt   &
-                  &                                       +2098.925 ) *zt+190925.6
+               zb1=     (-0.1909078_wp  *zt+7.390729_wp      ) *zt-55.87545_wp
+               za1=   ( ( 2.326469e-3_wp*zt+1.553190_wp      ) *zt-65.00517_wp ) *zt+1044.077_wp
+               zkw= ( ( (-1.361629e-4_wp*zt-1.852732e-2_wp   ) *zt-30.41638_wp ) *zt   &
+                  &                             +2098.925_wp ) *zt+190925.6_wp
                zk0= ( zb1*zsr + za1 )*zs + zkw
                !
                ! masked in situ density anomaly
-               prd(ji,jj) = ( zrhop / (  1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
+               prd(ji,jj) = ( zrhop / (  1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) )  ) - rau0 ) / rau0 * zmask
             END DO
          END DO
@@ -417,5 +419,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               prd(ji,jj) = ( 0.0285 - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
+               prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -490,28 +492,28 @@
                   zh = fsdepw(ji,jj,jk)                                     ! depth in meters  at w-point
                   !
-                  zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt   &   ! ratio alpha/beta
-                     &                               - 0.203814e-03 ) * zt   &
-                     &                               + 0.170907e-01 ) * zt   &
-                     &   + 0.665157e-01                                      &
-                     &   +     ( - 0.678662e-05 * zs                         &
-                     &           - 0.846960e-04 * zt + 0.378110e-02 ) * zs   &
-                     &   +   ( ( - 0.302285e-13 * zh                         &
-                     &           - 0.251520e-11 * zs                         &
-                     &           + 0.512857e-12 * zt * zt           ) * zh   &
-                     &           - 0.164759e-06 * zs                         &
-                     &        +(   0.791325e-08 * zt - 0.933746e-06 ) * zt   &
-                     &                               + 0.380374e-04 ) * zh
+                  zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &   ! ratio alpha/beta
+                     &                                  - 0.203814e-03_wp ) * zt   &
+                     &                                  + 0.170907e-01_wp ) * zt   &
+                     &   +         0.665157e-01_wp                                 &
+                     &   +     ( - 0.678662e-05_wp * zs                            &
+                     &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
+                     &   +   ( ( - 0.302285e-13_wp * zh                            &
+                     &           - 0.251520e-11_wp * zs                            &
+                     &           + 0.512857e-12_wp * zt * zt              ) * zh   &
+                     &           - 0.164759e-06_wp * zs                            &
+                     &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
+                     &                                  + 0.380374e-04_wp ) * zh
                      !
-                  zbeta  = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt      &   ! beta
-                     &                            - 0.301985e-05 ) * zt      &
-                     &   + 0.785567e-03                                      &
-                     &   + (     0.515032e-08 * zs                           &
-                     &         + 0.788212e-08 * zt - 0.356603e-06 ) * zs     &
-                     &   +(  (   0.121551e-17 * zh                           &
-                     &         - 0.602281e-15 * zs                           &
-                     &         - 0.175379e-14 * zt + 0.176621e-12 ) * zh     &
-                     &                             + 0.408195e-10   * zs     &
-                     &     + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt     &
-                     &                             - 0.121555e-07 ) * zh
+                  zbeta  = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt      &   ! beta
+                     &                               - 0.301985e-05_wp ) * zt      &
+                     &   +       0.785567e-03_wp                                   &
+                     &   + (     0.515032e-08_wp * zs                              &
+                     &         + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs     &
+                     &   + ( (   0.121551e-17_wp * zh                              &
+                     &         - 0.602281e-15_wp * zs                              &
+                     &         - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh     &
+                     &                                + 0.408195e-10_wp   * zs     &
+                     &     + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt     &
+                     &                                - 0.121555e-07_wp ) * zh
                      !
                   pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk)           &   ! N^2 
@@ -521,5 +523,5 @@
                   !                                                         !!bug **** caution a traiter zds=dk[S]= 0 !!!!
                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )                    ! Rrau = (alpha / beta) (dk[t] / dk[s])
-                  IF ( ABS( zds) <= 1.e-20 ) zds = 1.e-20
+                  IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
                   rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
 #endif
@@ -544,5 +546,5 @@
                DO ji = 1, jpi
                   zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )  
-                  IF ( ABS( zds ) <= 1.e-20 ) zds = 1.e-20
+                  IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp
                   rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
                END DO
@@ -560,4 +562,90 @@
 
 
+   SUBROUTINE eos_alpbet( pts, palph, pbeta )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE ldf_slp_grif  ***
+      !!
+      !! ** Purpose :   Calculates the thermal and haline expansion coefficients at T-points. 
+      !!
+      !! ** Method  :   calculates alpha and beta at T-points 
+      !!       * nn_eos = 0  : UNESCO sea water properties
+      !!         The brunt-vaisala frequency is computed using the polynomial
+      !!      polynomial expression of McDougall (1987):
+      !!            N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
+      !!      If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
+      !!      computed and used in zdfddm module :
+      !!              Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
+      !!       * nn_eos = 1  : linear equation of state (temperature only)
+      !!            N^2 = grav * rn_alpha * dk[ t ]/e3w
+      !!       * nn_eos = 2  : linear equation of state (temperature & salinity)
+      !!            N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
+      !!       * nn_eos = 3  : Jackett JAOT 2003 ???
+      !!
+      !! ** Action  : - palph, pbeta : thermal and haline expansion coeff. at T-point
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in   ) ::   pts            ! pot. temperature & salinity
+      REAL(wp), DIMENSION(jpi,jpj,jpk)     , INTENT(  out) ::   palph, pbeta   ! thermal & haline expansion coeff.
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zt, zs, zh   ! local scalars 
+      !!----------------------------------------------------------------------
+      !
+      SELECT CASE ( nn_eos )
+      !
+      CASE ( 0 )               ! Jackett and McDougall (1994) formulation
+         DO jk = 1, jpk
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zt = pts(ji,jj,jk,jp_tem)           ! potential temperature
+                  zs = pts(ji,jj,jk,jp_sal) - 35._wp  ! salinity anomaly (s-35)
+                  zh = fsdept(ji,jj,jk)              ! depth in meters 
+                  !
+                  pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt   &
+                     &                                        - 0.301985e-05_wp ) * zt   &
+                     &           + 0.785567e-03_wp                                       &
+                     &           + (     0.515032e-08_wp * zs                            &
+                     &                 + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs   &
+                     &           + ( (   0.121551e-17_wp * zh                            &
+                     &                 - 0.602281e-15_wp * zs                            &
+                     &                 - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh   &
+                     &                                        + 0.408195e-10_wp   * zs   &
+                     &             + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt   &
+                     &                                        - 0.121555e-07_wp ) * zh
+                     !
+                  palph(ji,jj,jk) = - pbeta(ji,jj,jk) *                             &
+                      &     ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt   &
+                      &                                  - 0.203814e-03_wp ) * zt   &
+                      &                                  + 0.170907e-01_wp ) * zt   &
+                      &   + 0.665157e-01_wp                                         &
+                      &   +     ( - 0.678662e-05_wp * zs                            &
+                      &           - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs   &
+                      &   +   ( ( - 0.302285e-13_wp * zh                            &
+                      &           - 0.251520e-11_wp * zs                            &
+                      &           + 0.512857e-12_wp * zt * zt              ) * zh   &
+                      &           - 0.164759e-06_wp * zs                            &
+                      &        +(   0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt   &
+                      &                                  + 0.380374e-04_wp ) * zh)
+               END DO
+            END DO
+         END DO
+         !
+      CASE ( 1 )
+         palph(:,:,:) = - rn_alpha
+         pbeta(:,:,:) =   0._wp
+         !
+      CASE ( 2 )
+         palph(:,:,:) = - rn_alpha
+         pbeta(:,:,:) =   rn_beta
+         !
+      CASE DEFAULT
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          bad flag value for nn_eos = ', nn_eos
+         nstop = nstop + 1
+         !
+      END SELECT
+      !
+   END SUBROUTINE eos_alpbet
+
+
    FUNCTION tfreez( psal ) RESULT( ptf )
       !!----------------------------------------------------------------------
@@ -576,6 +664,6 @@
       !!----------------------------------------------------------------------
       !
-      ptf(:,:) = ( - 0.0575 + 1.710523e-3 * SQRT( psal(:,:) )   &
-         &                  - 2.154996e-4 *       psal(:,:)   ) * psal(:,:)
+      ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) )   &
+         &                     - 2.154996e-4_wp *       psal(:,:)   ) * psal(:,:)
       !
    END FUNCTION tfreez
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90	(revision 2371)
@@ -73,5 +73,6 @@
       CASE ( 1 )    
          IF ( ln_traldf_grif ) THEN
-            CALL tra_ldf_iso_grif    ( kt )           ! Griffies quarter-cell formulation
+
+            CALL tra_ldf_iso_grif    ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0  )           ! Griffies quarter-cell formulation
          ELSE
             CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  ! rotated laplacian
@@ -85,5 +86,5 @@
          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
          IF ( ln_traldf_grif ) THEN
-            CALL tra_ldf_iso_grif    ( kt )
+            CALL tra_ldf_iso_grif    ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0  )
          ELSE
             CALL tra_ldf_iso   ( kt, 'TRA', gtsu, gtsv, tsb, tsa, jpts, ahtb0 )  
@@ -133,7 +134,4 @@
       INTEGER ::   ioptio, ierr         ! temporary integers 
 !      
-!     NAMELIST/namtra_ldf/ ln_traldf_lap  , ln_traldf_bilap,                  &
-!        &                 ln_traldf_level, ln_traldf_hor  , ln_traldf_iso,   &
-!        &                 rn_aht_0       , rn_ahtb_0      , rn_aeiv_0
       !!----------------------------------------------------------------------
 
@@ -141,18 +139,11 @@
       ! ===============================================
     
-!     REWIND( numnam )                ! Namelist namtra_ldf already read in ldftra module
-!     READ  ( numnam, namtra_ldf )    
-
       IF(lwp) THEN                    ! Namelist print
          WRITE(numout,*)
          WRITE(numout,*) 'tra_ldf_init : lateral tracer diffusive operator'
          WRITE(numout,*) '~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namtra_ldf : set lateral mixing parameters (type, direction, coefficients)'
-         WRITE(numout,*) '      laplacian operator          ln_traldf_lap   = ', ln_traldf_lap
-         WRITE(numout,*) '      bilaplacian operator        ln_traldf_bilap = ', ln_traldf_bilap
-         WRITE(numout,*) '      iso-level                   ln_traldf_level = ', ln_traldf_level
-         WRITE(numout,*) '      horizontal (geopotential)   ln_traldf_hor   = ', ln_traldf_hor
-         WRITE(numout,*) '      iso-neutral                 ln_traldf_iso   = ', ln_traldf_iso
-         WRITE(numout,*) '      iso-neutral (Griffies)      ln_traldf_grif  = ', ln_traldf_grif
+         WRITE(numout,*) '   Namelist namtra_ldf already read in ldftra module'
+         WRITE(numout,*) '   see ldf_tra_init report for lateral mixing parameters'
+         WRITE(numout,*)
       ENDIF
 
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90	(revision 2371)
@@ -1,8 +1,9 @@
 MODULE traldf_iso_grif
-   !!==============================================================================
+   !!======================================================================
    !!                   ***  MODULE  traldf_iso_grif  ***
-   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
-   !!==============================================================================
-   !! History :   9.0  !  06-10  (C. Harris)
+   !! 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
@@ -14,17 +15,18 @@
    !!                       and with the vertical part of
    !!                       the isopycnal or geopotential s-coord. operator 
-   !!                       vector optimization, use k-j-i loops.
-   !!----------------------------------------------------------------------
-
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and active tracers
    USE dom_oce         ! ocean space and time domain
    USE ldftra_oce      ! ocean active tracers: lateral physics
-   USE trdmod          ! ocean active tracers trends 
-   USE trdmod_oce      ! ocean variables trends
    USE zdf_oce         ! ocean vertical physics
    USE in_out_manager  ! I/O manager
+   USE iom             !
    USE ldfslp          ! iso-neutral slopes
    USE diaptr          ! poleward transport diagnostics
-   USE prtctl          ! Print control
+   USE trc_oce         ! share passive tracers/Ocean variables
+#if defined key_diaar5
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+#endif
 
    IMPLICIT NONE
@@ -32,11 +34,14 @@
 
    PUBLIC tra_ldf_iso_grif   ! routine called by traldf.F90
+
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   psix_eiv
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   psiy_eiv
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   ah_wslp2
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
 #  include "ldftra_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
 #  include "ldfeiv_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -47,146 +52,304 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_iso_grif( kt )
+  SUBROUTINE tra_ldf_iso_grif( kt, 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,   zftu => ua   ! use ua as workspace
+    USE oce,   zftv => va   ! use va 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       ! 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  ::  iku, ikv        ! temporary integer
+    INTEGER  ::  ierr            ! temporary integer
+    REAL(wp) ::  zmsku, zabe1, zcof1, zcoef3   ! local scalars
+    REAL(wp) ::  zmskv, zabe2, zcof2, zcoef4   !   -      -
+    REAL(wp) ::  zcoef0, zbtr                  !   -      -
+    REAL(wp), DIMENSION(jpi,jpj,0:1) ::   zdkt               ! 2D+1 workspace
+    REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, ztfw   ! 3D workspace
+      !
+    REAL(wp) ::   zslope_skew,zslope_iso,zslope2, zbu, zbv
+    REAL(wp) ::   ze1ur,zdxt,ze2vr,ze3wr,zdyt,zdzt
+    REAL(wp) ::   ah,zah_slp,zaei_slp
+#if defined key_diaar5
+    REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                ! 2D workspace
+    REAL(wp)                         ::   zztmp              ! local scalar
+#endif
       !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_ldf_iso_grif  ***
-      !!
-      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive 
-      !!      trend for a laplacian tensor (except 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 neutral or geopo-
-      !!      tential surfaces to which an eddy induced advection can be added
-      !!      It is computed using before fields (forward in time)
-      !!
-      !!      1st part :  masked horizontal derivative of T & S ( di[ t ] )
-      !!      ========    with partial cell update if ln_zps=T.
-      !!
-      !!      2nd part :  horizontal fluxes of the lateral mixing operator
-      !!      ========    
-      !!      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:
-
-      !!      take the horizontal divergence of the fluxes:
-      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
-      !!      Add this trend to the general trend (ta,sa):
-      !!         ta = ta + difft
-      !!
-      !! ** Action :
-      !!         Update (ta,sa) arrays with the before rotated diffusion trend
-      !!            (except the dk[ dk[.] ] term)
-      !!
+
+    IF( kt == nit000 )  THEN
+       IF(lwp) WRITE(numout,*)
+       IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype
+       IF(lwp) WRITE(numout,*) '                   WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL'
+       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+       ALLOCATE( ah_wslp2(jpi,jpj,jpk) , STAT=ierr )
+       IF( ierr > 0 ) THEN
+          CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator ah_wslp2 ' )   ;   RETURN
+       ENDIF
+       IF( ln_traldf_gdia ) THEN
+          ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
+          IF( ierr > 0 ) THEN
+             CALL ctl_stop( 'tra_ldf_iso_grif : unable to allocate Griffies operator diagnostics ' )   ;   RETURN
+          ENDIF
+       ENDIF
+    ENDIF
+
+      !                                               
       !!----------------------------------------------------------------------
-      USE oce           , zftv => ua   ! use ua as workspace
-      USE oce           , zfsv => va   ! use va as workspace
-      !!
-      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
-      !!
-      INTEGER  ::   ji, jj, jk, ip, jp, kp    ! dummy loop indices
-      INTEGER  ::   iku, ikv                  ! temporary integer
-      REAL(wp) ::   zatempu, zdx, zta         ! temporary scalars
-      REAL(wp) ::   zatempv, zdy, zsa         !    "         "
-      REAL(wp) ::   zslopec, zdsloper, zepsln !    "         "
-      REAL(wp) ::   zsxe, za_sxe, zfact       !    "         "
-      REAL(wp) ::   zbtr                      !    "         "
-      REAL(wp), DIMENSION(2) ::   zdelta   ! 1D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zftu   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zfsu   !    "           "
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdit, zdjt, zdkt     ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdis, zdjs, zdks     !  "      "
-
-
-
+      !!   0 - calculate  ah_wslp2, psix_eiv, psiy_eiv 
       !!----------------------------------------------------------------------
 
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~'
-      ENDIF
-
-      IF ( .NOT. lk_traldf_eiv ) THEN
-        fsaeiu(:,:,:)=0.0
-        fsaeiv(:,:,:)=0.0
-        fsaeiw(:,:,:)=0.0
-      ENDIF
-
-      DO jk=1,jpkm1
-      DO jj=1,jpjm1
-      DO ji=1,fs_jpim1
-               ftu(ji,jj,jk)=ftud(ji,jj,jk)+ftu(ji,jj,jk)*(fsahtu(ji,jj,jk)-fsaeiu(ji,jj,jk))
-               fsu(ji,jj,jk)=fsud(ji,jj,jk)+fsu(ji,jj,jk)*(fsahtu(ji,jj,jk)-fsaeiu(ji,jj,jk))
-               ftv(ji,jj,jk)=ftvd(ji,jj,jk)+ftv(ji,jj,jk)*(fsahtv(ji,jj,jk)-fsaeiv(ji,jj,jk))
-               fsv(ji,jj,jk)=fsvd(ji,jj,jk)+fsv(ji,jj,jk)*(fsahtv(ji,jj,jk)-fsaeiv(ji,jj,jk))
+!!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads
+
+    ah_wslp2(:,:,:) = 0.0
+    IF( ln_traldf_gdia ) THEN
+       psix_eiv(:,:,:) = 0.0
+       psiy_eiv(:,:,:) = 0.0
+    ENDIF
+
+    DO ip=0,1
+       DO kp=0,1
+          DO jk=1,jpkm1
+             DO jj=1,jpjm1
+                DO ji=1,fs_jpim1
+                   ze3wr=1.0_wp/fse3w(ji+ip,jj,jk+kp)
+                   zbu = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk)
+                   ah = fsahtu(ji,jj,jk)                                  !  fsaht(ji+ip,jj,jk)
+                   zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp)
+                   zslope2=(zslope_skew &
+                        & - umask(ji,jj,jk+kp)*(fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk))*ze1ur &
+                        & )**2
+                   ah_wslp2(ji+ip,jj,jk+kp)=ah_wslp2(ji+ip,jj,jk+kp) + &
+                        & ah*((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
+                   ze3wr=1.0_wp/fse3w(ji,jj+jp,jk+kp)
+                   zbv = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk)
+                   ah = fsahtu(ji,jj,jk)!fsaht(ji,jj+jp,jk)
+                   zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp)
+                   zslope2=(zslope_skew - vmask(ji,jj,jk+kp)*(fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk))*ze2vr &
+                        & )**2
+                   ah_wslp2(ji,jj+jp,jk+kp)=ah_wslp2(ji,jj+jp,jk+kp) + &
+                        & ah*((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
-      END DO
-      END DO
-
-      DO jk = 1, jpkm1
-
-         ! II.4 Second derivative (divergence) and add to the general trend
-         ! ----------------------------------------------------------------
-         DO jj = 2 , jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zbtr= 1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
-               zta = zbtr * ( ftu(ji,jj,jk) - ftu(ji-1,jj,jk) + ftv(ji,jj,jk) - ftv(ji,jj-1,jk)  )
-               zsa = zbtr * ( fsu(ji,jj,jk) - fsu(ji-1,jj,jk) + fsv(ji,jj,jk) - fsv(ji,jj-1,jk)  )
-               ta (ji,jj,jk) = ta (ji,jj,jk) + zta
-               sa (ji,jj,jk) = sa (ji,jj,jk) + zsa
+
+      !
+      !                                                          ! ===========
+      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
-         !                                          ! ===============
-      END DO                                        !   End of slab  
-      !                                             ! ===============
-
-      IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN   ! Poleward diffusive heat and salt transports
-         pht_ldf(:) = ptr_vj( zftv(:,:,:) )
-         pst_ldf(:) = ptr_vj( zfsv(:,:,:) )
-      ENDIF
-
-      !!----------------------------------------------------------------------
-      !!   III - vertical trend of T & S (extra diagonal terms only)
-      !!----------------------------------------------------------------------
-
-      ! I.5 Divergence of vertical fluxes added to the general tracer trend
-      ! -------------------------------------------------------------------
-
-      DO jk=1,jpk
-      DO jj=2,jpjm1
-      DO ji=fs_2,fs_jpim1
-               tfw(ji,jj,jk)=tfw(ji,jj,jk)*(fsahtw(ji,jj,jk)+fsaeiw(ji,jj,jk))
-               sfw(ji,jj,jk)=sfw(ji,jj,jk)*(fsahtw(ji,jj,jk)+fsaeiw(ji,jj,jk))
-      END DO
-      END DO
-      END DO
-
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zbtr =  1. / ( e1t(ji,jj)*e2t(ji,jj)*fse3t(ji,jj,jk) )
-               zta  = (  tfw(ji,jj,jk) - tfw(ji,jj,jk+1)  ) * zbtr
-               zsa  = (  sfw(ji,jj,jk) - sfw(ji,jj,jk+1)  ) * zbtr
-               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
-               sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
+         IF( ln_zps ) THEN                               ! partial steps: correction at the last level
+# if defined key_vectopt_loop
+            DO jj = 1, 1
+               DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+# else
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1
+# endif
+                  iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1                ! last level
+                  ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
+                  zdit(ji,jj,iku) = pgu(ji,jj,jn)          
+                  zdjt(ji,jj,ikv) = pgv(ji,jj,jn)      
+               END DO
             END DO
-         END DO
-      END DO
-      !
-      DO jk=1,jpk
-      DO jj=2,jpjm1
-      DO ji=fs_2,fs_jpim1
-       psix_eiv(ji,jj,jk) = psix_eiv(ji,jj,jk)*fsaeiu(ji,jj,jk)
-       psiy_eiv(ji,jj,jk) = psiy_eiv(ji,jj,jk)*fsaeiv(ji,jj,jk)
-      END DO
-      END DO
-      END DO
-
-   END SUBROUTINE tra_ldf_iso_grif
+         ENDIF
+
+         !!----------------------------------------------------------------------
+         !!   II - horizontal trend  (full)
+         !!----------------------------------------------------------------------
+         !
+         DO jk = 1, jpkm1
+            !
+            !                    !==  Vertical tracer gradient at level jk and jk+1
+            zdkt(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
+            !
+            !                          ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
+            IF( jk == 1 ) THEN   ;   zdkt(:,:,0) = zdkt(:,:,1)
+            ELSE                 ;   zdkt(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
+            ENDIF
+
+            IF( .NOT. l_triad_iso ) THEN
+               triadi = triadi_g
+               triadj = triadj_g
+            ENDIF
+
+            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  = zdkt(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)
+                        ah = fsahtu(ji,jj,jk)!*umask(ji,jj,jk+kp)         !fsaht(ji+ip,jj,jk)             ===>>  ????
+                        zah_slp = ah*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) - (ah*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 = zdkt(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)
+                        ah    = fsahtv(ji,jj,jk)!*vmask(ji,jj,jk+kp)         !fsaht(ji,jj+jp,jk)
+                        zah_slp = ah * 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) - (ah*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
+
+          !                        !==  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 .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
+         IF( jn == jp_tem)   pht_ldf(:) = ptr_vj( zftv(:,:,:) )        ! 3.3  names
+         IF( jn == jp_sal)   pst_ldf(:) = ptr_vj( zftv(:,:,:) )
+       ENDIF
+
+#if defined key_diaar5
+       IF( cdtype == 'TRA' .AND. jn == jp_tem  ) 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) + zztmp * zftu(ji,jj,jk)       &
+                        &       * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) * e1u(ji,jj) * fse3u(ji,jj,jk) 
+                END DO
+             END DO
+          END DO
+          CALL lbc_lnk( z2d, 'U', -1. )
+          CALL iom_put( "udiff_heattr", z2d )                  ! heat transport in i-direction
+          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) + zztmp * zftv(ji,jj,jk)       &
+                        &       * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) * e2v(ji,jj) * fse3v(ji,jj,jk) 
+                END DO
+             END DO
+          END DO
+          CALL lbc_lnk( z2d, 'V', -1. )
+          CALL iom_put( "vdiff_heattr", z2d )                  !  heat transport in i-direction
+       END IF
+#endif
+       !
+    END DO
+    !
+  END SUBROUTINE tra_ldf_iso_grif
 
 #else
@@ -196,5 +359,5 @@
 CONTAINS
    SUBROUTINE tra_ldf_iso_grif( kt )               ! Empty routine
-      WRITE(*,*) 'tra_ldf_iso_grif: You should not have seen this print! error?', kt
+      WRITE(*,*) 'tra_ldf_iso: You should not have seen this print! error?', kt
    END SUBROUTINE tra_ldf_iso_grif
 #endif
Index: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
===================================================================
--- branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2370)
+++ branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90	(revision 2371)
@@ -91,4 +91,5 @@
       USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
       USE oce    , ONLY :   zws   => va   ! va  -          -
+      USE traldf_iso_grif    , ONLY :   ah_wslp2
       !! 
       INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
@@ -145,5 +146,6 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
+                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
+                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
                         zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
                      END DO
@@ -201,5 +203,6 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
+                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
+                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
                         zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
                      END DO
