Index: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90
===================================================================
--- branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 4750)
+++ branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90	(revision 5965)
@@ -12,4 +12,5 @@
    !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
    !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
+   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
    !!----------------------------------------------------------------------
 #if   defined key_trabbl   ||   defined key_esopa
@@ -28,7 +29,8 @@
    USE phycst         ! physical constant
    USE eosbn2         ! equation of state
-   USE trdmod_oce     ! trends: ocean variables
+   USE trd_oce     ! trends: ocean variables
    USE trdtra         ! trends: active tracers
-   USE iom            ! IOM server
+   !
+   USE iom            ! IOM library               
    USE in_out_manager ! I/O manager
    USE lbclnk         ! ocean lateral boundary conditions
@@ -36,5 +38,5 @@
    USE wrk_nemo       ! Memory Allocation
    USE timing         ! Timing
-
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -57,5 +59,5 @@
    REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
 
-   LOGICAL , PUBLIC ::   l_bbl                  !: flag to compute bbl diffu. flux coef and transport
+   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
@@ -84,5 +86,5 @@
          &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
          &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
-         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , STAT= tra_bbl_alloc                )
+         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
          !
       IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
@@ -104,5 +106,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step
-      !!
+      !
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
       !!----------------------------------------------------------------------
@@ -110,5 +112,5 @@
       IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
       !
-      IF( l_trdtra )   THEN                        !* Save ta and sa trends
+      IF( l_trdtra )   THEN                         !* Save ta and sa trends
          CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
@@ -116,12 +118,12 @@
       ENDIF
 
-      IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
-
-      IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl
+      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
+
+      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
          !
          CALL tra_bbl_dif( tsb, tsa, jpts )
          IF( ln_ctl )  &
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
+            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
          CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
@@ -131,10 +133,10 @@
       END IF
 
-      IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl
+      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
          !
          CALL tra_bbl_adv( tsb, tsa, jpts )
          IF(ln_ctl)   &
          CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
-         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
+            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
          CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
@@ -147,6 +149,6 @@
          ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
          ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds )
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
          CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
       ENDIF
@@ -164,6 +166,5 @@
       !!                advection terms.
       !!
-      !! ** Method  :
-      !!        * diffusive bbl (nn_bbl_ldf=1) :
+      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
       !!        When the product grad( rho) * grad(h) < 0 (where grad is an
       !!      along bottom slope gradient) an additional lateral 2nd order
@@ -179,5 +180,4 @@
       !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
       !!----------------------------------------------------------------------
-      !
       INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
       REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
@@ -196,24 +196,14 @@
       DO jn = 1, kjpt                                     ! tracer loop
          !                                                ! ===========
-#  if defined key_vectopt_loop
-         DO jj = 1, 1   ! vector opt. (forced unrolling)
-            DO ji = 1, jpij
-#else
          DO jj = 1, jpj
             DO ji = 1, jpi
-#endif
-               ik = mbkt(ji,jj)                        ! bottom T-level index
-               zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S
+               ik = mbkt(ji,jj)                              ! bottom T-level index
+               zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S
             END DO
          END DO
-         !                                                ! Compute the trend
-#  if defined key_vectopt_loop
-         DO jj = 1, 1   ! vector opt. (forced unrolling)
-            DO ji = jpi+1, jpij-jpi-1
-#  else
-         DO jj = 2, jpjm1
+         !               
+         DO jj = 2, jpjm1                                    ! Compute the trend
             DO ji = 2, jpim1
-#  endif
-               ik = mbkt(ji,jj)                            ! bottom T-level index
+               ik = mbkt(ji,jj)                              ! bottom T-level index
                zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
@@ -264,11 +254,6 @@
       DO jn = 1, kjpt                                            ! tracer loop
          !                                                       ! ===========
-# if defined key_vectopt_loop
-         DO jj = 1, 1
-            DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling)
-# else
          DO jj = 1, jpjm1
             DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
-# endif
                IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
                   ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
@@ -333,6 +318,5 @@
       !!                advection terms.
       !!
-      !! ** Method  :
-      !!        * diffusive bbl (nn_bbl_ldf=1) :
+      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
       !!        When the product grad( rho) * grad(h) < 0 (where grad is an
       !!      along bottom slope gradient) an additional lateral 2nd order
@@ -342,5 +326,5 @@
       !!      a downslope velocity of 20 cm/s if the condition for slope
       !!      convection is satified)
-      !!        * advective bbl (nn_bbl_adv=1 or 2) :
+      !!              * advective bbl (nn_bbl_adv=1 or 2) :
       !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
       !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
@@ -353,56 +337,20 @@
       !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
       !!----------------------------------------------------------------------
-      !
       INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER         , INTENT(in   ) ::   kit000          ! first time step index
+      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
       CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
       !!
       INTEGER  ::   ji, jj                    ! dummy loop indices
       INTEGER  ::   ik                        ! local integers
-      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
-      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
-      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
-      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
-      !!
-      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
-      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep
-      !!----------------------- zv_bbl-----------------------------------------------
-      ! 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
-      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
-      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
-         ( ( ( -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
-      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
-         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
-                                 - 0.301985e-05 ) * pft                      &
-                                 + 0.785567e-03                              &
-         + (     0.515032e-08 * pfs                                          &
-               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
-               +(  (   0.121551e-17 * pfh                                    &
-                     - 0.602281e-15 * pfs                                    &
-                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
-                                          + 0.408195e-10   * pfs             &
-                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
-                                          - 0.121555e-07 ) * pfh
-      !!----------------------------------------------------------------------
-
+      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
+      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
+      REAL(wp) ::   za, zb, zgdrho            ! local scalars
+      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
+      !!----------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start( 'bbl')
       !
-      CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )
-      !
-
       IF( kt == kit000 )  THEN
          IF(lwp)  WRITE(numout,*)
@@ -410,50 +358,43 @@
          IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
       ENDIF
-
-      !                                        !* bottom temperature, salinity, velocity and depth
-#if defined key_vectopt_loop
-      DO jj = 1, 1   ! vector opt. (forced unrolling)
-         DO ji = 1, jpij
-#else
+      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
       DO jj = 1, jpj
          DO ji = 1, jpi
-#endif
-            ik = mbkt(ji,jj)                        ! bottom T-level index
-            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
-            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
-            zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth
+            ik = mbkt(ji,jj)                             ! bottom T-level index
+            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
+            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
             !
-            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
-            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
+            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
+            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
+            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
          END DO
       END DO
-
+      !
+      CALL eos_rab( zts, zdep, zab )
+      !
       !                                   !-------------------!
       IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
          !                                !-------------------!
          DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
-            DO ji = 1, jpim1
-               !                                                ! i-direction
-               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
-               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
-               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
-               !                                                         ! masked bbl i-gradient of density
-               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
-                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
+            DO ji = 1, fs_jpim1   ! vector opt.
+               !                                                   ! i-direction
+               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
+               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
+               !                                                         ! 2*masked bottom density gradient
+               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
+                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
                !
-               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
-               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
+               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
+               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
                !
-               !                                                ! j-direction
-               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
-               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
-               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
-               !                                                         ! masked bbl j-gradient of density
-               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
-                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
+               !                                                   ! j-direction
+               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
+               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
+               !                                                         ! 2*masked bottom density gradient
+               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
+                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
                !
-               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
+               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
                ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
-               !
             END DO
          END DO
@@ -469,29 +410,27 @@
             DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
                DO ji = 1, fs_jpim1   ! vector opt.
-                  !                                               ! i-direction
-                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
-                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
-                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
-                  !                                                           ! masked bbl i-gradient of density
-                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
-                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
-                  !
-                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
-                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
-                  !
-                  !                                                           ! bbl velocity
+                  !                                                  ! i-direction
+                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
+                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
+                  !                                                          ! 2*masked bottom density gradient 
+                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
+                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
+                  !
+                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
+                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
+                  !
+                  !                                                          ! bbl velocity
                   utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
                   !
-                  !                                               ! j-direction
-                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
-                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
-                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
-                  !                                                           ! masked bbl j-gradient of density
-                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
-                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
-                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
-                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
-                  !
-                  !                                                           ! bbl velocity
+                  !                                                  ! j-direction
+                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
+                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
+                  !                                                          ! 2*masked bottom density gradient
+                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
+                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
+                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
+                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
+                  !
+                  !                                                          ! bbl transport
                   vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
                END DO
@@ -502,36 +441,38 @@
             DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
                DO ji = 1, fs_jpim1   ! vector opt.
-                  !                                         ! i-direction
+                  !                                                  ! i-direction
                   ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
-                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
-                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
-                  !
-                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
-                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
-                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
-                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
-                  zgdrho =    fsbeta( zt, zs, zh )                                    &
-                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    &
-                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
-                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
-                  !
-                  !                                             ! bbl transport (down-slope direction)
+                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
+                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
+                  !
+                  ikud = mbku_d(ji,jj)
+                  ikus = mbku(ji,jj)
+                  !
+                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
+                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
+                  !                                                          !   masked bottom density gradient
+                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
+                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
+                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
+                  !
+                  !                                                          ! bbl transport (down-slope direction)
                   utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
                   !
-                  !                                         ! j-direction
+                  !                                                  ! j-direction
                   !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
-                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
-                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
-                  !
-                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
-                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
-                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
-                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
-                  zgdrho =    fsbeta( zt, zs, zh )                                    &
-                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    &
-                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
-                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
-                  !
-                  !                                             ! bbl transport (down-slope direction)
+                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
+                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
+                  !
+                  ikvd = mbkv_d(ji,jj)
+                  ikvs = mbkv(ji,jj)
+                  !
+                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
+                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
+                  !                                                          !   masked bottom density gradient
+                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
+                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
+                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
+                  !
+                  !                                                          ! bbl transport (down-slope direction)
                   vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
                END DO
@@ -541,6 +482,4 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )
-      !
       IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
       !
@@ -558,6 +497,6 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj               ! dummy loop indices
-      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
+      INTEGER ::   ios                  !   -      -
       REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
       !!
@@ -598,6 +537,4 @@
       IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
 
-      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
-
       !                             !* vertical index of  "deep" bottom u- and v-points
       DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
@@ -607,10 +544,10 @@
          END DO
       END DO
-      ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
+      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
       zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
       zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
 
-                                     !* sign of grad(H) at u- and v-points
-      mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
+                                        !* sign of grad(H) at u- and v-points
+      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
       DO jj = 1, jpjm1
          DO ji = 1, jpim1
@@ -621,5 +558,5 @@
 
       DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
-         DO ji = 1, jpim1           ! minimum of top & bottom e3u_0 (e3v_0)
+         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
             e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
             e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
