Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/trabbl.F90	(revision 11947)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/trabbl.F90	(revision 11949)
@@ -89,5 +89,5 @@
 
 
-   SUBROUTINE tra_bbl( kt )
+   SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -101,5 +101,7 @@
       !!              is added to the general tracer trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
+      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step
+      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
@@ -110,16 +112,16 @@
       IF( l_trdtra )   THEN                         !* Save the T-S input trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
-      ENDIF
-
-      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
+      ENDIF
+
+      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* 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 )
+         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
          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' )
+         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
+            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
          CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
@@ -131,8 +133,8 @@
       IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
          !
-         CALL tra_bbl_adv( tsb, tsa, jpts )
+         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
          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' )
+         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
+            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
          ! lateral boundary conditions ; just need for outputs
          CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
@@ -143,8 +145,8 @@
 
       IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
          DEALLOCATE( ztrdt, ztrds )
       ENDIF
@@ -155,5 +157,5 @@
 
 
-   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
+   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbl_dif  ***
@@ -171,5 +173,5 @@
       !!      convection is satified)
       !!
-      !! ** Action  :   pta   increased by the bbl diffusive trend
+      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
       !!
       !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
@@ -177,6 +179,7 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
+      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
       !
       INTEGER  ::   ji, jj, jn   ! dummy loop indices
@@ -191,5 +194,5 @@
             DO ji = 1, jpi
                ik = mbkt(ji,jj)                             ! bottom T-level index
-               zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S
+               zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
             END DO
          END DO
@@ -198,10 +201,10 @@
             DO ji = 2, jpim1
                ik = mbkt(ji,jj)                            ! bottom T-level index
-               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  &
-                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
-                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
-                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
-                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
-                  &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik)
+               pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
+                  &                + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
+                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
+                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
+                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
+                  &                * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
             END DO
          END DO
@@ -212,5 +215,5 @@
 
 
-   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
+   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_bbl  ***
@@ -228,6 +231,7 @@
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
+      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
       !
       INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
@@ -250,17 +254,17 @@
                   !
                   !                                               ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus)
-                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
-                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
+                  zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
+                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
+                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
                   !
                   DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk)
-                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
-                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
+                     zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
+                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
+                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
                   END DO
                   !
-                  zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud)
-                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
-                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
+                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
+                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
+                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
                ENDIF
                !
@@ -272,17 +276,17 @@
                   !
                   ! up  -slope T-point (shelf bottom point)
-                  zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs)
-                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
-                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
+                  zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
+                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
+                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
                   !
                   DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
-                     zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk)
-                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
-                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
+                     zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
+                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
+                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
                   END DO
                   !                                               ! down-slope T-point (deep bottom point)
-                  zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd)
-                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
-                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
+                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
+                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
+                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
                ENDIF
             END DO
@@ -295,5 +299,5 @@
 
 
-   SUBROUTINE bbl( kt, kit000, cdtype )
+   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -324,4 +328,5 @@
       INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
       CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
       !
       INTEGER  ::   ji, jj                    ! dummy loop indices
@@ -344,14 +349,14 @@
          DO ji = 1, jpi
             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)
+            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
+            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
             !
-            zdep(ji,jj) = gdept_n(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))
+            zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
+            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
+            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
          END DO
       END DO
       !
-      CALL eos_rab( zts, zdep, zab )
+      CALL eos_rab( zts, zdep, zab, Kmm )
       !
       !                                   !-------------------!
