Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbc.F90	(revision 10985)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE tra_bbc( kt, Kmm, Krhs )
+   SUBROUTINE tra_bbc( kt, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbc  ***
@@ -73,6 +73,7 @@
       !!              Emile-Geay and Madec, 2009, Ocean Science.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
+      INTEGER,                                   INTENT(in   ) :: kt         ! ocean time-step index
+      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts        ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj    ! dummy loop indices
@@ -84,22 +85,22 @@
       IF( l_trdtra )   THEN         ! Save the input temperature trend
          ALLOCATE( ztrdt(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
       ENDIF
       !                             !  Add the geothermal trend on temperature
       DO jj = 2, jpjm1
          DO ji = 2, jpim1
-            ts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = ts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
+            pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) = pts(ji,jj,mbkt(ji,jj),jp_tem,Krhs) + qgh_trd0(ji,jj) / e3t(ji,jj,mbkt(ji,jj),Kmm)
          END DO
       END DO
       !
-      CALL lbc_lnk( 'trabbc', ts(:,:,:,jp_tem,Krhs) , 'T', 1. )
+      CALL lbc_lnk( 'trabbc', pts(:,:,:,jp_tem,Krhs) , 'T', 1. )
       !
       IF( l_trdtra ) THEN        ! Send the trend for diagnostics
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbc, ztrdt )
          DEALLOCATE( ztrdt )
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
       IF( ln_timing )   CALL timing_stop('tra_bbc')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90	(revision 10985)
@@ -89,5 +89,5 @@
 
 
-   SUBROUTINE tra_bbl( kt, Kbb, Kmm, Krhs )
+   SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE bbl  ***
@@ -101,6 +101,7 @@
       !!              is added to the general tracer trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
-      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! time level indices
+      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
@@ -111,6 +112,6 @@
       IF( l_trdtra )   THEN                         !* Save the T-S input trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
-         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
       ENDIF
 
@@ -119,8 +120,8 @@
       IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
          !
-         CALL tra_bbl_dif( ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts, Kmm )
+         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
          IF( ln_ctl )  &
-         CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
-            &          tab3d_2=ts(:,:,:,jp_sal,Krhs), 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. )
@@ -132,8 +133,8 @@
       IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
          !
-         CALL tra_bbl_adv( ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts, Kmm )
+         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
          IF(ln_ctl)   &
-         CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
-            &          tab3d_2=ts(:,:,:,jp_sal,Krhs), 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. )
@@ -144,6 +145,6 @@
 
       IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
-         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs) - 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 )
@@ -156,5 +157,5 @@
 
 
-   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt, Kmm )
+   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_bbl_dif  ***
@@ -172,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.
@@ -178,6 +179,6 @@
       !!----------------------------------------------------------------------
       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
       !
@@ -193,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
@@ -200,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(ji,jj,ik,Kmm)
+               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
@@ -214,5 +215,5 @@
 
 
-   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt, Kmm )
+   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE trc_bbl  ***
@@ -230,6 +231,6 @@
       !!----------------------------------------------------------------------
       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
       !
@@ -254,16 +255,16 @@
                   !                                               ! up  -slope T-point (shelf bottom point)
                   zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
-                  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
+                  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(iid,jj,jk,Kmm)
-                     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
+                     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(iid,jj,ikud,Kmm)
-                  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
+                  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
                !
@@ -276,16 +277,16 @@
                   ! up  -slope T-point (shelf bottom point)
                   zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
-                  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
+                  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(ji,ijd,jk,Kmm)
-                     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
+                     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(ji,ijd,ikvd,Kmm)
-                  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
+                  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
@@ -348,9 +349,9 @@
          DO ji = 1, jpi
             ik = mbkt(ji,jj)                             ! bottom T-level index
-            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb)    ! bottom before T and S
+            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(ji,jj,ik,Kmm)              ! bottom T-level reference depth
-            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)          ! bottom velocity
+            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
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tradmp.F90	(revision 10985)
@@ -72,5 +72,5 @@
 
 
-   SUBROUTINE tra_dmp( kt, Kbb, Kmm, Krhs )
+   SUBROUTINE tra_dmp( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tra_dmp  ***
@@ -90,6 +90,7 @@
       !! ** Action  : - tsa: tracer trends updated with the damping trend
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
+      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index
+      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
       !
       INTEGER ::   ji, jj, jk, jn   ! dummy loop indices
@@ -102,5 +103,5 @@
       IF( l_trdtra )   THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) ) 
-         ztrdts(:,:,:,:) = ts(:,:,:,:,Krhs) 
+         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) 
       ENDIF
       !                           !==  input T-S data at kt  ==!
@@ -114,5 +115,6 @@
                DO jj = 2, jpjm1
                   DO ji = fs_2, fs_jpim1   ! vector opt.
-                     ts(ji,jj,jk,jn,Krhs) = ts(ji,jj,jk,jn,Krhs) + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - ts(ji,jj,jk,jn,Kbb) )
+                     pts(ji,jj,jk,jn,Krhs) = pts(ji,jj,jk,jn,Krhs)           &
+                        &                  + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jn) - pts(ji,jj,jk,jn,Kbb) )
                   END DO
                END DO
@@ -125,8 +127,8 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,Kbb) )
-                     ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,Kbb) )
+                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
                   ENDIF
                END DO
@@ -139,8 +141,8 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
-                     ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - ts(ji,jj,jk,jp_tem,Kbb) )
-                     ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs)   &
-                        &                 + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - ts(ji,jj,jk,jp_sal,Kbb) )
+                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_tem) - pts(ji,jj,jk,jp_tem,Kbb) )
+                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)   &
+                        &                      + resto(ji,jj,jk) * ( zts_dta(ji,jj,jk,jp_sal) - pts(ji,jj,jk,jp_sal,Kbb) )
                   ENDIF
                END DO
@@ -151,5 +153,5 @@
       !
       IF( l_trdtra )   THEN       ! trend diagnostic
-         ztrdts(:,:,:,:) = ts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
+         ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_dmp, ztrdts(:,:,:,jp_tem) )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_dmp, ztrdts(:,:,:,jp_sal) )
@@ -157,6 +159,6 @@
       ENDIF
       !                           ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' dmp  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_dmp')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranpc.F90	(revision 10985)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE tra_npc( kt, Kmm, Krhs )
+   SUBROUTINE tra_npc( kt, Kmm, Krhs, pts, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tranpc  ***
@@ -58,6 +58,7 @@
       !! References :     Madec, et al., 1991, JPO, 21, 9, 1349-1371.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
+      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index
+      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -67,10 +68,10 @@
       REAL(wp) ::   zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z
       REAL(wp) ::   zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw, z1_r2dt
-      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp      ! acceptance criteria for neutrality (N2==0)
-      REAL(wp), DIMENSION(        jpk     ) ::   zvn2         ! vertical profile of N2 at 1 given point...
-      REAL(wp), DIMENSION(        jpk,jpts) ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
-      REAL(wp), DIMENSION(jpi,jpj,jpk     ) ::   zn2          ! N^2 
-      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) ::   zab          ! alpha and beta
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds   ! 3D workspace
+      REAL(wp), PARAMETER ::   zn2_zero = 1.e-14_wp             ! acceptance criteria for neutrality (N2==0)
+      REAL(wp), DIMENSION(        jpk     )   ::   zvn2         ! vertical profile of N2 at 1 given point...
+      REAL(wp), DIMENSION(        jpk,jpts)   ::   zvts, zvab   ! vertical profile of T & S , and  alpha & betaat 1 given point
+      REAL(wp), DIMENSION(jpi,jpj,jpk     )   ::   zn2          ! N^2 
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts)   ::   zab          ! alpha and beta
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds ! 3D workspace
       !
       LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is
@@ -85,6 +86,6 @@
          IF( l_trdtra )   THEN                    !* Save initial after fields
             ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-            ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) 
-            ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
+            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 
+            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
          ENDIF
          !
@@ -96,6 +97,6 @@
          ENDIF
          !
-         CALL eos_rab( ts(:,:,:,:,Krhs), zab, Kmm )         ! after alpha and beta (given on T-points)
-         CALL bn2    ( ts(:,:,:,:,Krhs), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
+         CALL eos_rab( pts(:,:,:,:,Kaa), zab, Kmm )         ! after alpha and beta (given on T-points)
+         CALL bn2    ( pts(:,:,:,:,Kaa), zab, zn2, Kmm )    ! after Brunt-Vaisala  (given on W-points)
          !
          inpcc = 0
@@ -106,6 +107,6 @@
                IF( tmask(ji,jj,2) == 1 ) THEN      ! At least 2 ocean points
                   !                                     ! consider one ocean column 
-                  zvts(:,jp_tem) = ts(ji,jj,:,jp_tem,Krhs)      ! temperature
-                  zvts(:,jp_sal) = ts(ji,jj,:,jp_sal,Krhs)      ! salinity
+                  zvts(:,jp_tem) = pts(ji,jj,:,jp_tem,Kaa)      ! temperature
+                  zvts(:,jp_sal) = pts(ji,jj,:,jp_sal,Kaa)      ! salinity
                   !
                   zvab(:,jp_tem)  = zab(ji,jj,:,jp_tem)     ! Alpha 
@@ -287,7 +288,7 @@
                   END DO ! DO WHILE ( .NOT. l_column_treated )
 
-                  !! Updating tsa:
-                  ts(ji,jj,:,jp_tem,Krhs) = zvts(:,jp_tem)
-                  ts(ji,jj,:,jp_sal,Krhs) = zvts(:,jp_sal)
+                  !! Updating pts:
+                  pts(ji,jj,:,jp_tem,Kaa) = zvts(:,jp_tem)
+                  pts(ji,jj,:,jp_sal,Kaa) = zvts(:,jp_sal)
 
                   !! LB:  Potentially some other global variable beside theta and S can be treated here
@@ -303,6 +304,6 @@
          IF( l_trdtra ) THEN         ! send the Non penetrative mixing trends for diagnostic
             z1_r2dt = 1._wp / (2._wp * rdt)
-            ztrdt(:,:,:) = ( ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) ) * z1_r2dt
-            ztrds(:,:,:) = ( ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) ) * z1_r2dt
+            ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_r2dt
+            ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_r2dt
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt )
             CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds )
@@ -310,5 +311,5 @@
          ENDIF
          !
-         CALL lbc_lnk_multi( 'tranpc', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. )
+         CALL lbc_lnk_multi( 'tranpc', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
          !
          IF( lwp .AND. l_LB_debug ) THEN
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traqsr.F90	(revision 10985)
@@ -75,5 +75,5 @@
 CONTAINS
 
-   SUBROUTINE tra_qsr( kt, Kmm, Krhs )
+   SUBROUTINE tra_qsr( kt, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_qsr  ***
@@ -101,6 +101,7 @@
       !!              Morel, A. et Berthon, JF, 1989, Limnol Oceanogr 34(8), 1545-1562
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ocean time-step
-      INTEGER, INTENT(in) ::   Kmm, Krhs     ! time level indices
+      INTEGER,                                   INTENT(in   ) :: kt            ! ocean time-step
+      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs     ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts           ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk               ! dummy loop indices
@@ -127,5 +128,5 @@
       IF( l_trdtra ) THEN      ! trends diagnostic: save the input temperature trend
          ALLOCATE( ztrdt(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
       ENDIF
       !
@@ -262,6 +263,6 @@
          DO jj = 2, jpjm1        !-----------------------------!
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)   &
-                  &                 + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
+               pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)   &
+                  &                      + z1_2 * ( qsr_hc_b(ji,jj,jk) + qsr_hc(ji,jj,jk) ) / e3t(ji,jj,jk,Kmm)
             END DO
          END DO
@@ -296,10 +297,10 @@
       !
       IF( l_trdtra ) THEN     ! qsr tracers trends saved for diagnostics
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_qsr, ztrdt )
          DEALLOCATE( ztrdt ) 
       ENDIF
       !                       ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' qsr  - Ta: ', mask1=tmask, clinfo3='tra-ta' )
       !
       IF( ln_timing )   CALL timing_stop('tra_qsr')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90	(revision 10985)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE tra_sbc ( kt, Kmm, Krhs )
+   SUBROUTINE tra_sbc ( kt, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_sbc  ***
@@ -72,6 +72,7 @@
       !!              - send trends to trdtra module for further diagnostics(l_trdtra=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kmm, Krhs  ! time level indices
+      INTEGER,                                   INTENT(in   ) :: kt         ! ocean time-step index
+      INTEGER,                                   INTENT(in   ) :: Kmm, Krhs  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts        ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices  
@@ -91,6 +92,6 @@
       IF( l_trdtra ) THEN                    !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs)
-         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
       ENDIF
       !
@@ -132,10 +133,10 @@
          DO jj = 2, jpj                         !==>> add concentration/dilution effect due to constant volume cell
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
-               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
+               sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm)
+               sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rau0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm)
             END DO
          END DO                                 !==>> output c./d. term
-         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * ts(:,:,1,jp_tem,Kmm) )
-         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * ts(:,:,1,jp_sal,Kmm) )
+         IF( iom_use('emp_x_sst') )   CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) )
+         IF( iom_use('emp_x_sss') )   CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) )
       ENDIF
       !
@@ -143,5 +144,5 @@
          DO jj = 2, jpj
             DO ji = fs_2, fs_jpim1   ! vector opt.  
-               ts(ji,jj,1,jn,Krhs) = ts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm)
+               pts(ji,jj,1,jn,Krhs) = pts(ji,jj,1,jn,Krhs) + zfact * ( sbc_tsc_b(ji,jj,jn) + sbc_tsc(ji,jj,jn) ) / e3t(ji,jj,1,Kmm)
             END DO
          END DO
@@ -174,14 +175,14 @@
                DO jk = ikt, ikb - 1
                ! compute trend
-                  ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)                                                &
-                     &           + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
-                     &           * r1_hisf_tbl(ji,jj)
+                  pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)                                      &
+                     &                      + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )  &
+                     &                      * r1_hisf_tbl(ji,jj)
                END DO
    
                ! level partially include in ice shelf boundary layer 
                ! compute trend
-               ts(ji,jj,ikb,jp_tem,Krhs) = ts(ji,jj,ikb,jp_tem,Krhs)                                                 &
-                  &              + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )             &
-                  &              * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
+               pts(ji,jj,ikb,jp_tem,Krhs) = pts(ji,jj,ikb,jp_tem,Krhs)                                       &
+                  &                       + zfact * ( risf_tsc_b(ji,jj,jp_tem) + risf_tsc(ji,jj,jp_tem) )    &
+                  &                       * r1_hisf_tbl(ji,jj) * ralpha(ji,jj)
 
             END DO
@@ -200,8 +201,8 @@
                   zdep = zfact / h_rnf(ji,jj)
                   DO jk = 1, nk_rnf(ji,jj)
-                                        ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs)                                 &
-                                           &                 +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
-                     IF( ln_rnf_sal )   ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs)                                 &
-                                           &                 +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
+                                        pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs)                                  &
+                                           &                      +  ( rnf_tsc_b(ji,jj,jp_tem) + rnf_tsc(ji,jj,jp_tem) ) * zdep
+                     IF( ln_rnf_sal )   pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs)                                  &
+                                           &                      +  ( rnf_tsc_b(ji,jj,jp_sal) + rnf_tsc(ji,jj,jp_sal) ) * zdep 
                   END DO
                ENDIF
@@ -210,6 +211,6 @@
       ENDIF
 
-      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*ts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst
-      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*ts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss
+      IF( iom_use('rnf_x_sst') )   CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) )   ! runoff term on sst
+      IF( iom_use('rnf_x_sss') )   CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) )   ! runoff term on sss
 
 #if defined key_asminc
@@ -225,6 +226,6 @@
                DO ji = fs_2, fs_jpim1
                   ztim = ssh_iau(ji,jj) / e3t(ji,jj,1,Kmm)
-                  ts(ji,jj,1,jp_tem,Krhs) = ts(ji,jj,1,jp_tem,Krhs) + ts(ji,jj,1,jp_tem,Kmm) * ztim
-                  ts(ji,jj,1,jp_sal,Krhs) = ts(ji,jj,1,jp_sal,Krhs) + ts(ji,jj,1,jp_sal,Kmm) * ztim
+                  pts(ji,jj,1,jp_tem,Krhs) = pts(ji,jj,1,jp_tem,Krhs) + pts(ji,jj,1,jp_tem,Kmm) * ztim
+                  pts(ji,jj,1,jp_sal,Krhs) = pts(ji,jj,1,jp_sal,Krhs) + pts(ji,jj,1,jp_sal,Kmm) * ztim
                END DO
             END DO
@@ -233,6 +234,6 @@
                DO ji = fs_2, fs_jpim1
                   ztim = ssh_iau(ji,jj) / ( ht_n(ji,jj) + 1. - ssmask(ji, jj) )
-                  ts(ji,jj,:,jp_tem,Krhs) = ts(ji,jj,:,jp_tem,Krhs) + ts(ji,jj,:,jp_tem,Kmm) * ztim
-                  ts(ji,jj,:,jp_sal,Krhs) = ts(ji,jj,:,jp_sal,Krhs) + ts(ji,jj,:,jp_sal,Kmm) * ztim
+                  pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim
+                  pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim
                END DO  
             END DO  
@@ -252,6 +253,6 @@
                DO ji = fs_2, fs_jpim1
                   zdep = 1._wp / e3t(ji,jj,jk,Kmm) 
-                  ts(ji,jj,jk,jp_tem,Krhs) = ts(ji,jj,jk,jp_tem,Krhs) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
-                  ts(ji,jj,jk,jp_sal,Krhs) = ts(ji,jj,jk,jp_sal,Krhs) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
+                  pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) - htsc_iscpl(ji,jj,jk,jp_tem) * zdep
+                  pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) - htsc_iscpl(ji,jj,jk,jp_sal) * zdep  
                END DO  
             END DO  
@@ -260,6 +261,6 @@
 
       IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
-         ztrdt(:,:,:) = ts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
-         ztrds(:,:,:) = ts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_nsr, ztrdt )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_nsr, ztrds )
@@ -267,6 +268,6 @@
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
-         &                       tab3d_2=ts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' sbc  - Ta: ', mask1=tmask,   &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_sbc')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10980)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10985)
@@ -242,10 +242,10 @@
       IF(  lk_asminc .AND. ln_asmiau .AND. &
          & ln_trainc )   CALL tra_asm_inc   ( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment
-                         CALL tra_sbc       ( kstp,      Nnn, Nrhs )  ! surface boundary condition
-      IF( ln_traqsr  )   CALL tra_qsr       ( kstp,      Nnn, Nrhs )  ! penetrative solar radiation qsr
-      IF( ln_trabbc  )   CALL tra_bbc       ( kstp,      Nnn, Nrhs )  ! bottom heat flux
-      IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nbb, Nnn, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
-      IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nbb, Nnn, Nrhs )  ! internal damping trends
-      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp, Nbb, ts, Nrhs )  ! bdy damping trends
+                         CALL tra_sbc       ( kstp,      Nnn, ts, Nrhs )  ! surface boundary condition
+      IF( ln_traqsr  )   CALL tra_qsr       ( kstp,      Nnn, ts, Nrhs )  ! penetrative solar radiation qsr
+      IF( ln_trabbc  )   CALL tra_bbc       ( kstp,      Nnn, ts, Nrhs )  ! bottom heat flux
+      IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nbb, Nnn, ts, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
+      IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nbb, Nnn, ts, Nrhs )  ! internal damping trends
+      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp, Nbb,      ts, Nrhs )  ! bdy damping trends
 #if defined key_agrif
       IF(.NOT. Agrif_Root())  & 
@@ -262,5 +262,5 @@
 !!gm
                          CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer	==> after
-      IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Nnn, Nrhs           )  ! update after fields by non-penetrative convection
+      IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
