Index: /trunk/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 408)
+++ /trunk/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 409)
@@ -16,9 +16,10 @@
    !!----------------------------------------------------------------------
    !! * Modules used
-   USE oce             ! ocean dynamics and active tracers
-   USE dom_oce         ! ocean space and time domain
-   USE trdmod_oce      ! ocean variables trends
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
+   USE oce                  ! ocean dynamics and active tracers
+   USE dom_oce              ! ocean space and time domain
+   USE trdmod_oce           ! ocean variables trends
+   USE eosbn2 , ONLY : neos ! type of equation of state
+   USE in_out_manager       ! I/O manager
+   USE prtctl               ! Print control
 
    IMPLICIT NONE
@@ -30,9 +31,14 @@
 
    !! * Shared module variables
+   REAL(wp), PUBLIC ::            &  !!: * bbl namelist *
+      atrbbl = 1.e+3                  !: lateral coeff. for BBL scheme (m2/s)
+#if defined key_trabbl_dif
    LOGICAL, PUBLIC, PARAMETER ::   &  !:
       lk_trabbl_dif = .TRUE.          !: diffusive bottom boundary layer flag
-   REAL(wp), PUBLIC ::            &  !!: * bbl namelist *
-      atrbbl = 1.e+3                  !: lateral coeff. for bottom boundary 
-      !                               !  layer scheme (m2/s)
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   &  !:
+      lk_trabbl_dif = .FALSE.         !: diffusive bottom boundary layer flag
+#endif
+
 # if defined key_trabbl_adv
    LOGICAL, PUBLIC, PARAMETER ::    &  !:
@@ -228,4 +234,8 @@
       ! multiplied by the slope of the ocean bottom
 
+      SELECT CASE ( neos )
+
+      CASE ( 0   )               ! 0 :Jackett and McDougall (1994) formulation
+
 #  if defined key_vectopt_loop   &&   ! defined key_autotasking
       jj = 1
@@ -276,4 +286,54 @@
       END DO
 
+      CASE ( 1 )               ! Linear formulation function of temperature only
+                               ! 
+#  if defined key_vectopt_loop   &&   ! defined key_autotasking
+      jj = 1
+      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+#  else
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+#  endif
+            ! local 'density/temperature' gradient along i-bathymetric slope
+            zgdrho =  ztnb(ji+1,jj) - ztnb(ji,jj) 
+            ! sign of local i-gradient of density multiplied by the i-slope
+            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
+            zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
+#  if ! defined key_vectopt_loop   ||   defined key_autotasking
+         END DO
+#  endif
+      END DO
+
+#  if defined key_vectopt_loop   &&   ! defined key_autotasking
+      jj = 1
+      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
+#  else
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+#  endif
+            ! local density gradient along j-bathymetric slope
+            zgdrho =  ztnb(ji,jj+1) - ztnb(ji,jj) 
+            ! sign of local j-gradient of density multiplied by the j-slope
+            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
+            zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
+#  if ! defined key_vectopt_loop   ||   defined key_autotasking
+         END DO
+#  endif
+      END DO
+
+      CASE ( 2 )               ! Linear formulation function of temperature and salinity
+
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rbeta * S - ralpha * T )'
+         IF(lwp) WRITE(numout,*) '          bbl not implented: easy to do it '
+         nstop = nstop + 1
+
+      CASE DEFAULT
+
+         IF(lwp) WRITE(numout,cform_err)
+         IF(lwp) WRITE(numout,*) '          bad flag value for neos = ', neos
+         nstop = nstop + 1
+
+      END SELECT
 
       ! 2. Additional second order diffusive trends
@@ -420,6 +480,9 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'tra_bbl_init : * Diffusive Bottom Boundary Layer'
+         WRITE(numout,*) 'tra_bbl_init : '
          WRITE(numout,*) '~~~~~~~~~~~~'
+         IF (lk_trabbl_dif ) THEN
+            WRITE(numout,*) '               * Diffusive Bottom Boundary Layer'
+         ENDIF 
          IF( lk_trabbl_adv ) THEN
             WRITE(numout,*) '               * Advective Bottom Boundary Layer'
Index: /trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90
===================================================================
--- /trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90	(revision 408)
+++ /trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90	(revision 409)
@@ -73,6 +73,5 @@
          zgdrho, zbtr, zta, zsa            !    "         " 
       REAL(wp), DIMENSION(jpi,jpj) ::   &
-         zki, zkj, zkw, zkx, zky, zkz,  &  ! temporary workspace arrays
-         ztnb, zsnb, zdep, ztbb, zsbb,  &  !    "                  "
+         ztnb, zsnb, zdep, ztbb, zsbb,  &  ! temporary workspace arrays
          zahu, zahv                        !    "                  "
       REAL(wp), DIMENSION(jpi,jpj) ::   &  ! temporary workspace arrays
@@ -160,25 +159,4 @@
        CALL lbc_lnk( zsnb, 'T', 1. )    ;   CALL lbc_lnk( zsbb, 'T', 1. )
 
-      ! Conditional diffusion along the slope in the bottom boundary layer
-
-#if defined key_trabbl_dif
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-# else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-# endif
-            iku = mbku(ji,jj)
-            ikv = mbkv(ji,jj)
-            zahu(ji,jj) = atrbbl*e2u(ji,jj)*fse3u(ji,jj,iku)/e1u(ji,jj) * umask(ji,jj,1)
-            zahv(ji,jj) = atrbbl*e1v(ji,jj)*fse3v(ji,jj,ikv)/e2v(ji,jj) * vmask(ji,jj,1)
-# if ! defined key_vectopt_loop   ||   defined key_autotasking
-         END DO
-# endif
-      END DO
-#endif
-
-
       ! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
       ! --------------------------------------------
@@ -204,5 +182,4 @@
       !   ... sign of local i-gradient of density multiplied by the i-slope
           zsign = sign( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
-          zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
 
           zsigna= sign(0.5, zunb(ji,jj)*(  zdep(ji+1,jj) - zdep(ji,jj) ))
@@ -225,5 +202,4 @@
       !   ... sign of local j-gradient of density multiplied by the j-slope
           zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
-          zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
 
           zsigna= sign(0.5, zvnb(ji,jj)*(zdep(ji,jj+1) - zdep(ji,jj) ) )
@@ -234,9 +210,28 @@
 
       CASE ( 1 )               ! Linear formulation function of temperature only
-
-         IF(lwp) WRITE(numout,cform_err)
-         IF(lwp) WRITE(numout,*) '          use of linear eos rho(T,S) = rau0 * ( rbeta * S - ralpha * T )'
-         IF(lwp) WRITE(numout,*) '          bbl not implented: easy to do it '
-         nstop = nstop + 1
+                               !
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            ! local 'density/temperature' gradient along i-bathymetric slope
+            zgdrho =  ztnb(ji+1,jj) - ztnb(ji,jj)
+            ! sign of local i-gradient of density multiplied by the i-slope
+            zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
+
+            zsigna= sign(0.5, zunb(ji,jj)*(  zdep(ji+1,jj) - zdep(ji,jj) ))
+            zalphax(ji,jj)=(0.5+zsigna)*(0.5-zsign)*umask(ji,jj,1)
+         END DO
+      END DO
+
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+            ! local density gradient along j-bathymetric slope
+            zgdrho =  ztnb(ji,jj+1) - ztnb(ji,jj)
+            ! sign of local j-gradient of density multiplied by the j-slope
+            zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
+
+            zsigna= sign(0.5, zvnb(ji,jj)*(zdep(ji,jj+1) - zdep(ji,jj) ) )
+            zalphay(ji,jj)=(0.5+zsigna)*(0.5-zsign)*vmask(ji,jj,1)
+         END DO
+      END DO
 
       CASE ( 2 )               ! Linear formulation function of temperature and salinity
@@ -285,93 +280,4 @@
       ! lateral boundary conditions on u_bbl and v_bbl   (changed sign)
        CALL lbc_lnk( u_bbl, 'U', -1. )   ;   CALL lbc_lnk( v_bbl, 'V', -1. )
-
-
-
-#if defined key_trabbl_dif
-      ! 4. Additional second order diffusive trends
-      ! -------------------------------------------
-
-      ! ... first derivative (gradient)
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vertor opt.
-            zkx(ji,jj) = zki(ji,jj)*( ztbb(ji+1,jj) - ztbb(ji,jj) )
-            zkz(ji,jj) = zki(ji,jj)*( zsbb(ji+1,jj) - zsbb(ji,jj) )
-
-            zky(ji,jj) = zkj(ji,jj)*( ztbb(ji,jj+1) - ztbb(ji,jj) )
-            zkw(ji,jj) = zkj(ji,jj)*( zsbb(ji,jj+1) - zsbb(ji,jj) )
-         END DO
-      END DO
-
-      IF( cp_cfg == "orca" ) THEN   
-         SELECT CASE ( jp_cfg )
-            !                                        ! =======================
-            CASE ( 2 )                               !  ORCA_R2 configuration
-               !                                     ! =======================
-               ! Gibraltar enhancement of BBL
-               zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zkx( mi0(139):mi1(140) , mj0(102):mj1(102) )
-               zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zky( mi0(139):mi1(140) , mj0(102):mj1(102) )
-   
-               ! Red Sea enhancement of BBL
-               zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zkx( mi0(161):mi1(162) , mj0(88):mj1(88) )
-               zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zky( mi0(161):mi1(162) , mj0(88):mj1(88) )
-   
-               !                                     ! =======================
-            CASE ( 4 )                               !  ORCA_R4 configuration
-               !                                     ! =======================
-               ! Gibraltar enhancement of BBL
-               zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zkx( mi0(70):mi1(71) , mj0(52):mj1(52) )
-               zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zky( mi0(70):mi1(71) , mj0(52):mj1(52) )
- 
-         END SELECT
- 
-      ENDIF
-
-      ! ... second derivative (divergence) and add to the general tracer trend
-
-# if defined key_vectopt_loop   &&   ! defined key_autotasking
-      jj = 1
-      DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-# endif
-            ik = mbkt(ji,jj)
-            zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
-            zta = (  zkx(ji,jj) - zkx(ji-1,jj  )   & 
-               &   + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
-            zsa = (  zkz(ji,jj) - zkz(ji-1,jj  )   &
-               &   + zkw(ji,jj) - zkw(ji  ,jj-1)  ) * zbtr
-            ta(ji,jj,ik) = ta(ji,jj,ik) + zta
-            sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
-#if ! defined key_vectopt_loop   ||   defined key_autotasking
-         END DO
-#endif
-      END DO
-
-      ! save the trends for diagnostic
-      ! BBL lateral diffusion tracers trends
-      IF( l_trdtra )   THEN
-#  if defined key_vectopt_loop   &&   ! defined key_autotasking
-         jj = 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-#  else
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-#  endif
-            ik = mbkt(ji,jj)
-            tldfbbl(ji,jj) = ta(ji,jj,ik) - ztdta(ji,jj,ik)
-            sldfbbl(ji,jj) = sa(ji,jj,ik) - ztdsa(ji,jj,ik)
-#  if ! defined key_vectopt_loop   ||   defined key_autotasking
-            END DO
-#  endif
-         END DO
-
-         ! save the new ta & sa trends
-         ztdta(:,:,:) = ta(:,:,:)
-         ztdsa(:,:,:) = sa(:,:,:)
-
-      ENDIF
-
-#endif
 
       ! 5. Along sigma advective trend
@@ -398,10 +304,10 @@
             zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztbb(ji  ,jj  )   &
                &          +( zfui - ABS( zfui ) ) * ztbb(ji+1,jj  ) ) * 0.5
-            zwy(ji,jj) = ( ( zfui + ABS( zfvj ) ) * ztbb(ji  ,jj  )   &
-               &          +( zfui - ABS( zfvj ) ) * ztbb(ji  ,jj+1) ) * 0.5
+            zwy(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * ztbb(ji  ,jj  )   &
+               &          +( zfvj - ABS( zfvj ) ) * ztbb(ji  ,jj+1) ) * 0.5
             zww(ji,jj) = ( ( zfui + ABS( zfui ) ) * zsbb(ji  ,jj  )   &
                &          +( zfui - ABS( zfui ) ) * zsbb(ji+1,jj  ) ) * 0.5
-            zwz(ji,jj) = ( ( zfui + ABS( zfvj ) ) * zsbb(ji  ,jj  )   &
-               &          +( zfui - ABS( zfvj ) ) * zsbb(ji  ,jj+1) ) * 0.5
+            zwz(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * zsbb(ji  ,jj  )   &
+               &          +( zfvj - ABS( zfvj ) ) * zsbb(ji  ,jj+1) ) * 0.5
 #if ! defined key_vectopt_loop   ||   defined key_autotasking
          END DO
@@ -462,7 +368,7 @@
       DO jk= 1, jpkm1
          DO jj=1, jpjm1
-            DO ji = 1, fs_jpim1   ! vertor opt.
-               zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk)
-               zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk)
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk) * fse3u(ji,jj,jk)
+               zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk) * fse3v(ji,jj,jk)
             END DO
          END DO
@@ -471,5 +377,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zbt = e1t(ji,jj) * e2t(ji,jj)
+               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
                zhdivn(ji,jj,jk) = (  zwu(ji,jj) - zwu(ji-1,jj  )   &
                                    + zwv(ji,jj) - zwv(ji  ,jj-1)  ) / zbt
