Index: trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90
===================================================================
--- trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 2528)
@@ -1,540 +1,100 @@
 MODULE trcbbl
-   !!======================================================================
+  !!======================================================================
    !!                       ***  MODULE  trcbbl  ***
    !! Ocean passive tracers physics :  advective and/or diffusive bottom boundary 
    !!                                  layer scheme
    !!======================================================================
-   !! History :  8.0  !  96-06  (L. Mortier)  Original code
-   !!            8.0  !  97-11  (G. Madec)  Optimization
-   !!            8.5  !  02-08  (G. Madec)  free form + modules
-   !!            9.0  !  04-03  (C. Ethe)   Adaptation for passive tracers
-   !!                 !  07-02  (C. Deltel)  Diagnose ML trends for passive tracers
+   !!==============================================================================
+   !! History :  OPA  !  1996-06  (L. Mortier)  Original code
+   !!            8.0  !  1997-11  (G. Madec)    Optimization
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  free form + modules
+   !!             -   !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl
+   !!            3.3  !  2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization 
+   !!             -   !  2010-04  (G. Madec)  Campin & Goosse advective bbl 
+   !!             -   !  2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
    !!----------------------------------------------------------------------
-#if  defined key_top && ( defined key_trcbbl_dif   ||   defined key_trcbbl_adv ) && ! defined key_c1d
+#if  defined key_top &&  defined key_trabbl 
    !!----------------------------------------------------------------------
-   !!   'key_trcbbl_dif'   or            diffusive bottom boundary layer
-   !!   'key_trcbbl_adv'                 advective bottom boundary layer
+   !!   'key_trabbl                      diffusive or/and adevective bottom boundary layer
    !!----------------------------------------------------------------------
-   !!   trc_bbl_dif  : update the passive tracer trends due to the bottom
-   !!                  boundary layer (diffusive only)
-   !!   trc_bbl_adv  : update the passive tracer trends due to the bottom
-   !!                  boundary layer (advective and/or diffusive)
+   !!    trc_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
    !!----------------------------------------------------------------------
    USE oce_trc             ! ocean dynamics and active tracers variables
    USE trc                 ! ocean passive tracers variables
-   USE trctrp_lec          ! passive tracers transport
+   USE trcnam_trp      ! passive tracers transport namelist variables
+   USE trabbl              ! 
    USE prtctl_trc          ! Print control for debbuging
-   USE eosbn2
-   USE lbclnk
-   USE trdmld_trc
-   USE trdmld_trc_oce     
+   USE trdmod_oce
+   USE trdtra
 
-   IMPLICIT NONE
-   PRIVATE
+   PUBLIC   trc_bbl       !  routine called by step.F90
 
-   PUBLIC trc_bbl_dif    ! routine called by step.F90
-   PUBLIC trc_bbl_adv    ! routine called by step.F90
-
-# if defined key_trcbbl_dif
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_dif = .TRUE.   !: diffusive bottom boundary layer flag
-# else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_dif = .FALSE.  !: diffusive bottom boundary layer flag
-# endif
-
-# if defined key_trcbbl_adv
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_adv = .TRUE.   !: advective bottom boundary layer flag
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   u_trc_bbl  !: veloc. involved in the advective BBL
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   v_trc_bbl  !: veloc. involved in the advective BBL
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   w_trc_bbl  !: vertic. increment of veloc. due to adv. BBL
-   !                                                        !  only affect tracer vertical advection
-# else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_adv = .FALSE.  !: advective bottom boundary layer flag
-# endif
-
-   INTEGER, DIMENSION(jpi,jpj) ::   mbkt, mbku, mbkv
 
    !! * Substitutions
 #  include "top_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   TOP 1.0 , LOCEAN-IPSL (2005) 
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/TOP_SRC/TRP/trcbbl.F90,v 1.12 2006/09/12 11:10:13 opalod Exp $ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
 CONTAINS
 
-   SUBROUTINE trc_bbl_dif( kt )
+
+   SUBROUTINE trc_bbl( kt )
       !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE trc_bbl_dif  ***
+      !!                  ***  ROUTINE bbl  ***
+      !!                   
+      !! ** Purpose :   Compute the before tracer (t & s) trend associated 
+      !!     with the bottom boundary layer and add it to the general trend
+      !!     of tracer equations.
       !!
-      !! ** Purpose :   Compute the before tracer trend associated 
-      !!      with the bottom boundary layer and add it to the general trend 
-      !!      of tracer equations. The bottom boundary layer is supposed to be
-      !!      a purely diffusive bottom boundary layer.
-      !!
-      !! ** Method  :   When the product grad( rho) * grad(h) < 0 (where grad 
-      !!      is an along bottom slope gradient) an additional lateral diffu-
-      !!      sive trend along the bottom slope is added to the general tracer
-      !!      trend, otherwise nothing is done.
-      !!      Second order operator (laplacian type) with variable coefficient
-      !!      computed as follow for temperature (idem on s): 
-      !!         difft = 1/(e1t*e2t*e3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
-      !!                                 + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
-      !!      where ztb is a 2D array: the bottom ocean temperature and ahtb
-      !!      is a time and space varying diffusive coefficient defined by:
-      !!         ahbt = zahbp    if grad(rho).grad(h) < 0
-      !!              = 0.       otherwise.
-      !!      Note that grad(.) is the along bottom slope gradient. grad(rho)
-      !!      is evaluated using the local density (i.e. referenced at the
-      !!      local depth). Typical value of ahbt is 2000 m2/s (equivalent to
-      !!      a downslope velocity of 20 cm/s if the condition for slope
-      !!      convection is satified)
-      !!      Add this before trend to the general trend tra of the 
-      !!      botton ocean tracer point:
-      !!         tra = tra + difft
-      !!
-      !! ** Action  : - update tra at the bottom level with the bottom
-      !!                boundary layer trend
-      !!
-      !! References :
-      !!     Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
-      !!----------------------------------------------------------------------
-      USE oce, ONLY :   ztrtrd => ua                      ! use ua as 3D workspace   
-      !!
-      INTEGER, INTENT( in ) ::   kt                         ! ocean time-step
-      INTEGER ::   ji, jj, jn                               ! dummy loop indices
-      INTEGER ::   ik, iku, ikv
-      INTEGER ::   ii0, ii1, ij0, ij1                       ! temporary integers
-      INTEGER ::   iku1, iku2, ikv1, ikv2                   ! temporary intergers
-      REAL(wp) ::  ze3u, ze3v                              ! temporary scalars
-      REAL(wp) ::  zbtr, ztra
-#if ! defined key_off_tra
-      REAL(wp) ::   zgdrho, zalbet, zsign, zt, zs, zh
-      REAL(wp), DIMENSION(jpi,jpj) ::   zki, zkj
-#endif
-      REAL(wp), DIMENSION(jpi,jpj) ::   zkx, zky  ! temporary workspace arrays
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztnb, zsnb, zdep
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztrb, zahu, zahv
-
+      !!----------------------------------------------------------------------  
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
       CHARACTER (len=22) :: charout
-      REAL(wp) ::   fsalbt, pft, pfs, pfh                   ! statement function
-      !!----------------------------------------------------------------------
-      ! ratio alpha/beta
-      ! ================
-      !  fsalbt: ratio of thermal over saline expension coefficients
-      !       pft :  potential temperature in degrees celcius
-      !       pfs :  salinity anomaly (s-35) in psu
-      !       pfh :  depth in meters
-
-      fsalbt( pft, pfs, pfh ) =                                              &
-         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
-                                   - 0.203814e-03 ) * pft                    &
-                                   + 0.170907e-01 ) * pft                    &
-                                   + 0.665157e-01                            &
-         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
-         +  ( ( - 0.302285e-13 * pfh                                         &
-                - 0.251520e-11 * pfs                                         &
-                + 0.512857e-12 * pft * pft          ) * pfh                  &
-                                     - 0.164759e-06   * pfs                  &
-             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
-                                     + 0.380374e-04 ) * pfh   
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   ztrtrd
       !!----------------------------------------------------------------------
 
-
-      IF( kt == nittrc000 )   CALL trc_bbl_init
-
-
-      ! 0. 2D fields of bottom temperature and salinity, and bottom slope
-      ! -----------------------------------------------------------------
-      ! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
-
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-#  endif
-            ik = mbkt(ji,jj)                              ! index of the bottom ocean T-level
-            ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1)   ! masked now T and S at ocean bottom 
-            zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1)
-            zdep(ji,jj) = fsdept(ji,jj,ik)                ! depth of the ocean bottom T-level
-#  if ! defined key_vectopt_loop
-         END DO
-#  endif
-      END DO
-
-      IF( ln_zps ) THEN      ! partial steps correction
-
-#   if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#   else
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-#   endif
-               iku1 = MAX( mbathy(ji+1,jj  )-1, 1 )
-               iku2 = MAX( mbathy(ji  ,jj  )-1, 1 )
-               ikv1 = MAX( mbathy(ji  ,jj+1)-1, 1 )
-               ikv2 = MAX( mbathy(ji  ,jj  )-1, 1 )
-               ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) ) 
-               ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) ) 
-               zahu(ji,jj) = atrcbbl * e2u(ji,jj) * ze3u / e1u(ji,jj) * umask(ji,jj,1)
-               zahv(ji,jj) = atrcbbl * e1v(ji,jj) * ze3v / e2v(ji,jj) * vmask(ji,jj,1)
-#   if ! defined key_vectopt_loop
-            END DO
-#   endif
-         END DO
-      ELSE                  ! z-coordinate - full steps or s-coordinate
-#   if defined key_vectopt_loop
-         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) = atrcbbl * e2u(ji,jj) * fse3u(ji,jj,iku) / e1u(ji,jj) * umask(ji,jj,1)
-               zahv(ji,jj) = atrcbbl * e1v(ji,jj) * fse3v(ji,jj,ikv) / e2v(ji,jj) * vmask(ji,jj,1)
-#   if ! defined key_vectopt_loop
-            END DO
-#   endif
-         END DO
+      IF( .NOT. lk_offline ) THEN
+         CALL bbl( kt, 'TRC' )         ! Online coupling with dynamics  : Computation of bbl coef and bbl transport
+         l_bbl = .FALSE.               ! Offline coupling with dynamics : Read bbl coef and bbl transport from input files
       ENDIF
 
-#if defined key_off_tra
-      !!=====================================================================
-      !!                I. OFFLINE VERSION OF DIFFUSIVE BBL
-      !!=====================================================================
-      
-      ! 1. Criteria of additional bottom diffusivity : grad(rho).grad(h) < 0
-      ! --------------------------------------------------------------------
-      
-      !    Only used in the online version of diffusive BBL (see below)
-      
-      ! 2. Additional second order diffusive trends
-      ! -------------------------------------------
-      !                                                          ! ===========
-      DO jn = 1, jptra                                           ! tracer loop
-         !                                                       ! ===========       
-
-         IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)
-     
-         ! first derivative (gradient)         
-#  if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij   ! vector opt. (forced unrolling)
-#  else
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-#  endif
-               ik = mbkt(ji,jj) 
-               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
-#  if ! defined key_vectopt_loop
-            END DO
-#  endif
-         END DO
-
-#  if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-#  endif
-               zkx(ji,jj) = bblx(ji,jj) * zahu(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
-               zky(ji,jj) = bbly(ji,jj) * zahv(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
-#  if ! defined key_vectopt_loop
-            END DO
-#  endif
-         END DO
-
-#else
-      !!=====================================================================
-      !!               II. ONLINE VERSION OF DIFFUSIVE BBL
-      !!=====================================================================
-
-      ! 1. Criteria of additional bottom diffusivity : grad(rho).grad(h) < 0
-      ! --------------------------------------------------------------------
-      ! Sign of the local density gradient along the i- and j-slopes
-      ! multiplied by the slope of the ocean bottom
-      SELECT CASE ( nn_eos )
-
-      CASE ( 0 )                 ! Jackett and McDougall (1994) formulation
-         
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-#  endif
-            ! temperature, salinity anomalie and depth
-            zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
-            zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
-            zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
-            ! masked ratio alpha/beta
-            zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
-            ! local density gradient along i-bathymetric slope
-            zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) )   &
-                   -          ( zsnb(ji+1,jj) - zsnb(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
-         END DO
-#  endif
-      END DO
-
-#  if defined key_vectopt_loop
-      jj = 1
-      DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-#  endif
-            ! temperature, salinity anomalie and depth
-            zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
-            zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
-            zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
-            ! masked ratio alpha/beta
-            zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
-            ! local density gradient along j-bathymetric slope
-            zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) )   &
-                   -          ( zsnb(ji,jj+1) - zsnb(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
-         END DO
-#  endif
-      END DO
-      
-      CASE ( 1 )                 ! Linear formulation function of temperature only
-
-#  if defined key_vectopt_loop
-      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 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
-         END DO
-#  endif
-      END DO
-
-#  if defined key_vectopt_loop
-      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
-         END DO
-#  endif
-      END DO
-
-      CASE ( 2 )                 ! Linear formulation function of temperature and salinity
-
-      DO jj = 1, jpjm1
-        DO ji = 1, fs_jpim1   ! vector opt.
-            ! local density gradient along i-bathymetric slope
-            zgdrho = - ( rn_beta*( zsnb(ji+1,jj) - zsnb(ji,jj) )   &
-                     -  rn_alpha*( 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)
-        END DO
-      END DO
-
-      DO jj = 1, jpjm1
-        DO ji = 1, fs_jpim1   ! vector opt.
-            ! local density gradient along j-bathymetric slope
-            zgdrho = - ( rn_beta*( zsnb(ji,jj+1) - zsnb(ji,jj) )   &
-                     -  rn_alpha*( 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)
-         END DO
-      END DO
-
-      CASE DEFAULT
-
-         WRITE(ctmp1,*) '          bad flag value for nn_eos = ', nn_eos
-         CALL ctl_stop( ctmp1 )
-
-      END SELECT
-      
-      ! 2. Additional second order diffusive trends
-      ! -------------------------------------------
-      !                                                          ! ===========
-      DO jn = 1, jptra                                           ! tracer loop
-         !                                                       ! ===========
-         IF( l_trdtrc )   ztrtrd(:,:,:) = tra(:,:,:,jn)
-
-         ! first derivative (gradient)
-#  if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij   ! vector opt. (forced unrolling)
-#  else
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-#  endif
-               ik = mbkt(ji,jj)
-               ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1)
-#  if ! defined key_vectopt_loop
-            END DO
-#  endif
-         END DO
-#  if defined key_vectopt_loop
-         jj = 1
-         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling)
-#  else
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-#  endif
-               zkx(ji,jj) = zki(ji,jj) * ( ztrb(ji+1,jj) - ztrb(ji,jj) )
-               zky(ji,jj) = zkj(ji,jj) * ( ztrb(ji,jj+1) - ztrb(ji,jj) )
-#  if ! defined key_vectopt_loop
-            END DO
-#  endif
-         END DO
-#endif
-
-      !!=====================================================================
-      !!     III. COMMON CODE FOR OFFLINE/ONLINE VERSIONS OF DIFFUSIVE BBL
-      !!=====================================================================
-
-         IF( cp_cfg == "orca" ) THEN
-            
-            SELECT CASE ( jp_cfg )
-               !                                        ! =======================
-            CASE ( 2 )                                  !  ORCA_R2 configuration
-               !                                        ! =======================
-               ! Gibraltar enhancement of BBL
-               ij0 = 102   ;   ij1 = 102
-               ii0 = 139   ;   ii1 = 140  
-               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
-               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
-               
-               ! Red Sea enhancement of BBL
-               ij0 =  88   ;   ij1 =  88
-               ii0 = 161   ;   ii1 = 162
-               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
-               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 10.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
-               
-               !                                        ! =======================
-            CASE ( 4 )                                  !  ORCA_R4 configuration
-               !                                        ! =======================
-               ! Gibraltar enhancement of BBL
-               ij0 =  52   ;   ij1 =  52
-               ii0 =  70   ;   ii1 =  71  
-               zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zkx( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
-               zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0 * zky( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) )
-               
-            END SELECT
-            
-         ENDIF
-         
-         ! second derivative (divergence) and add to the general tracer trend
-#  if defined key_vectopt_loop
-         jj = 1
-         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
-#  else
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-#  endif
-               ik = MAX( mbathy(ji,jj)-1, 1 )
-               zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
-               ztra = (  zkx(ji,jj) - zkx(ji-1,jj  )    &
-                  &    + zky(ji,jj) - zky(ji  ,jj-1)  ) * zbtr
-               tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
-#  if ! defined key_vectopt_loop
-            END DO
-#  endif
-         END DO
-
-         ! save the trends for diagnostic
-         IF( l_trdtrc ) THEN
-            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
-            IF (luttrd(jn)) CALL trd_mod_trc( ztrtrd, jn, jptrc_trd_bbl, kt )
-         END IF
-         !                                                       ! ===========
-      END DO                                                     ! tracer loop
-      !                                                          ! ===========
-
-      IF( ln_ctl ) THEN    ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('bbl - dif')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+      IF( l_trdtrc )  THEN
+         ALLOCATE( ztrtrd(jpi,jpj,jpk,jptra) )   ! temporary save of trends
+         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
       ENDIF
 
-   END SUBROUTINE trc_bbl_dif
+      !* Diffusive bbl :
+      IF( nn_bbl_ldf == 1 ) THEN
+         !
+         CALL tra_bbl_dif( trb, tra, jptra )  
+         IF( ln_ctl )   THEN
+            WRITE(charout, FMT="(' bbl_dif')")  ;  CALL prt_ctl_trc_info(charout)
+            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         ENDIF
+         !
+      END IF
 
-# if defined key_trcbbl_adv
-   !!----------------------------------------------------------------------
-   !!   'key_trcbbl_adv'                    advective bottom boundary layer
-   !!----------------------------------------------------------------------
-#  include "trcbbl_adv.h90"
-# else
-   !!----------------------------------------------------------------------
-   !!   Default option :                 NO advective bottom boundary layer
-   !!----------------------------------------------------------------------
-   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_bbl_adv
-# endif
+      !* Advective bbl : bbl upstream advective trends added to the tracer trends
+      IF( nn_bbl_adv /= 0 ) THEN
+         !
+         CALL tra_bbl_adv( trb, tra, jptra )  
+         IF( ln_ctl )   THEN
+            WRITE(charout, FMT="(' bbl_adv')")  ;  CALL prt_ctl_trc_info(charout)
+            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         ENDIF
+         !
+      END IF
 
-   SUBROUTINE trc_bbl_init
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE trc_bbl_init  ***
-      !!
-      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj
-      REAL(wp),  DIMENSION(jpi,jpj) ::   zmbk  
-      !!----------------------------------------------------------------------
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )   ! vertical index of the bottom ocean T-level
-         END DO
-      END DO
-      
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-            mbku(ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )
-            mbkv(ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
-         END DO
-      END DO
-
-      zmbk(:,:) = FLOAT( mbku (:,:) )   
-      CALL lbc_lnk(zmbk,'U',1.)
-      mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
-	
-      zmbk(:,:) = FLOAT( mbkv (:,:) )   
-      CALL lbc_lnk(zmbk,'V',1.)
-      mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
-
-# if defined key_trcbbl_adv
-      w_trc_bbl(:,:,:) = 0.e0    ! initialisation of w_trc_bbl to zero
-# endif
-
-   END SUBROUTINE trc_bbl_init
+      IF( l_trdtrc )   THEN                      ! save the horizontal diffusive trends for further diagnostics
+        DO jn = 1, jptra
+           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+           CALL trd_tra( kt, 'TRC', jn, jptra_trd_ldf, ztrtrd(:,:,:,jn) )
+        END DO
+        DEALLOCATE( ztrtrd )
+      ENDIF
+      !
+   END SUBROUTINE trc_bbl
 
 #else
@@ -542,15 +102,8 @@
    !!   Dummy module :                      No bottom boundary layer scheme
    !!----------------------------------------------------------------------
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_dif = .FALSE.   !: diff bbl flag
-   LOGICAL, PUBLIC, PARAMETER ::   lk_trcbbl_adv = .FALSE.   !: adv  bbl flag
 CONTAINS
-   SUBROUTINE trc_bbl_dif (kt )              ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'trc_bbl_dif: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_bbl_dif
-   SUBROUTINE trc_bbl_adv (kt )              ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'trc_bbl_adv: You should not have seen this print! error?', kt
-   END SUBROUTINE trc_bbl_adv
+   SUBROUTINE trc_bbl( kt )              ! Empty routine
+      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_bbl
 #endif
 
