Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10979)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf.F90	(revision 10980)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf( kt, Kbb, Kmm, Krhs )
+   SUBROUTINE tra_ldf( kt, Kbb, Kmm, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf  ***
@@ -53,6 +53,7 @@
       !! ** Purpose :   compute the lateral ocean tracer physics.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kbb, Kmm, Krhs  ! ocean time level indices
+      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step index
+      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! ocean 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
@@ -63,22 +64,22 @@
       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
       !
       SELECT CASE ( nldf_tra )                 !* compute lateral mixing trend and add it to the general trend
       CASE ( np_lap   )                                  ! laplacian: iso-level operator
-         CALL tra_ldf_lap  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs),                   jpts,  1, Kmm )
+         CALL tra_ldf_lap  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),                   jpts,  1 )
       CASE ( np_lap_i )                                  ! laplacian: standard iso-neutral operator (Madec)
-         CALL tra_ldf_iso  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts,  1, Kmm )
+         CALL tra_ldf_iso  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
       CASE ( np_lap_it )                                 ! laplacian: triad iso-neutral operator (griffies)
-         CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs), jpts,  1, Kmm )
+         CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts,  1 )
       CASE ( np_blp , np_blp_i , np_blp_it )             ! bilaplacian: iso-level & iso-neutral operators
-         CALL tra_ldf_blp  ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, ts(:,:,:,:,Kbb), ts(:,:,:,:,Krhs),             jpts, nldf_tra, Kmm )
+         CALL tra_ldf_blp  ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs),             jpts, nldf_tra )
       END SELECT
       !
       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_ldf, ztrdt )
          CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds )
@@ -86,6 +87,6 @@
       ENDIF
       !                                        !* print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ts(:,:,:,jp_tem,Krhs), clinfo1=' ldf  - 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=' ldf  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_ldf')
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90	(revision 10979)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_iso.F90	(revision 10980)
@@ -48,7 +48,7 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_iso( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
-      &                                                   pgui, pgvi,   &
-      &                                       ptb , ptbb, pta , kjpt, kpass, Kmm )
+  SUBROUTINE tra_ldf_iso( kt, Kmm, kit000, cdtype, pahu, pahv,                    &
+      &                                            pgu , pgv    ,   pgui, pgvi,   &
+      &                                       pt , pt2 , pt_rhs , kjpt  , kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_iso  ***
@@ -87,7 +87,7 @@
       !!         difft = 1/(e1e2t*e3t) dk[ zftw ]
       !!      Add this trend to the general trend (ta,sa):
-      !!         pta = pta + difft
-      !!
-      !! ** Action :   Update pta arrays with the before rotated diffusion
+      !!         pt_rhs = pt_rhs + difft
+      !!
+      !! ** Action :   Update pt_rhs arrays with the before rotated diffusion
       !!----------------------------------------------------------------------
       INTEGER                              , INTENT(in   ) ::   kt         ! ocean time-step index
@@ -100,7 +100,7 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -220,6 +220,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -249,8 +249,8 @@
             !
             !                             !== Vertical tracer gradient
-            zdk1t(:,:) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
+            zdk1t(:,:) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * wmask(:,:,jk+1)     ! level jk+1
             !
             IF( jk == 1 ) THEN   ;   zdkt(:,:) = zdk1t(:,:)                          ! surface: zdkt(jk=1)=zdkt(jk=2)
-            ELSE                 ;   zdkt(:,:) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * wmask(:,:,jk)
+            ELSE                 ;   zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk)
             ENDIF
             DO jj = 1 , jpjm1            !==  Horizontal fluxes
@@ -277,9 +277,9 @@
             END DO
             !
-            DO jj = 2 , jpjm1          !== horizontal divergence and add to pta
+            DO jj = 2 , jpjm1          !== horizontal divergence and add to pt_rhs
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
-                     &                                           + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
-                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji,jj,jk) - zftu(ji-1,jj,jk)      &
+                     &                                                 + zftv(ji,jj,jk) - zftv(ji,jj-1,jk)  )   &
+                     &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -327,6 +327,6 @@
                   DO ji = fs_2, fs_jpim1
                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)   &
-                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
-                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )               &
+                        &                            * (  pt(ji,jj,jk-1,jn) -  pt(ji,jj,jk,jn) )
                   END DO
                END DO
@@ -339,17 +339,17 @@
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)    &
-                           &           + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj)   &
-                           &           * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk)                       &
+                           &           + ah_wslp2(ji,jj,jk)  * e1e2t(ji,jj)   &
+                           &           * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)
                      END DO
                   END DO
                END DO 
-            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+            CASE(  2  )                         ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
-                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                      &
-                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
-                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                        ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk)                  &
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
+                           &                            +         akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
                      END DO
                   END DO
@@ -358,9 +358,9 @@
          ENDIF
          !         
-         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
-                     &                                        * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1)  )   &
+                     &                                              * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
                END DO
             END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10979)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_lap_blp.F90	(revision 10980)
@@ -45,7 +45,7 @@
 CONTAINS
 
-   SUBROUTINE tra_ldf_lap( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
-      &                                                    pgui, pgvi,   &
-      &                                        ptb , pta , kjpt, kpass, Kmm ) 
+   SUBROUTINE tra_ldf_lap( kt, Kmm, kit000, cdtype, pahu, pahv  ,               &
+      &                                             pgu , pgv   , pgui, pgvi,   &
+      &                                             pt  , pt_rhs, kjpt, kpass ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_lap  ***
@@ -59,8 +59,8 @@
       !!          difft = 1/(e1e2t*e3t) {  di-1[ pahu e2u*e3u/e1u di(tb) ]
       !!                                 + dj-1[ pahv e1v*e3v/e2v dj(tb) ] }
-      !!      Add this trend to the general tracer trend pta :
-      !!          pta = pta + difft
-      !!
-      !! ** Action  : - Update pta arrays with the before iso-level 
+      !!      Add this trend to the general tracer trend pt_rhs :
+      !!          pt_rhs = pt_rhs + difft
+      !!
+      !! ** Action  : - Update pt_rhs arrays with the before iso-level 
       !!                harmonic mixing trend.
       !!----------------------------------------------------------------------
@@ -74,6 +74,6 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      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 tracer fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend 
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -114,6 +114,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1
-                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) )
-                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
+                  ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) )
+                  ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) )
                END DO
             END DO
@@ -139,7 +139,7 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
-                     &                                   + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
-                     &                                / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)     &
+                     &                                      +    ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )   &
+                     &                                      / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
                END DO
             END DO
@@ -160,7 +160,7 @@
    
 
-   SUBROUTINE tra_ldf_blp( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
-      &                                                    pgui, pgvi,   &
-      &                                                    ptb , pta , kjpt, kldf, Kmm )
+   SUBROUTINE tra_ldf_blp( kt, Kmm, kit000, cdtype, pahu, pahv  ,             &
+      &                                             pgu , pgv   , pgui, pgvi, &
+      &                                             pt  , pt_rhs, kjpt, kldf )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE tra_ldf_blp  ***
@@ -184,6 +184,6 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu, pgv   ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top levels
-      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 ::   ji, jj, jk, jn   ! dummy loop indices
@@ -205,12 +205,12 @@
       zlap(:,:,:,:) = 0._wp
       !
-      SELECT CASE ( kldf )       !==  1st laplacian applied to ptb (output in zlap)  ==!
+      SELECT CASE ( kldf )       !==  1st laplacian applied to pt (output in zlap)  ==!
       !
       CASE ( np_blp    )               ! iso-level bilaplacian
-         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb,      zlap, kjpt, 1, Kmm )
+         CALL tra_ldf_lap  ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt,     zlap, kjpt, 1 )
       CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
-         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1, Kmm )
+         CALL tra_ldf_iso  ( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
-         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1, Kmm )
+         CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, pt, pt, zlap, kjpt, 1 )
       END SELECT
       !
@@ -221,12 +221,12 @@
       ENDIF
       !
-      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pta)  ==!
+      SELECT CASE ( kldf )       !==  2nd laplacian applied to zlap (output in pt_rhs)  ==!
       !
       CASE ( np_blp    )               ! iso-level bilaplacian
-         CALL tra_ldf_lap  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta,      kjpt, 2, Kmm )
+         CALL tra_ldf_lap  ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt_rhs,         kjpt, 2 )
       CASE ( np_blp_i  )               ! rotated   bilaplacian : standard operator (Madec)
-         CALL tra_ldf_iso  ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2, Kmm )
+         CALL tra_ldf_iso  ( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt    , pt_rhs, kjpt, 2 )
       CASE ( np_blp_it )               ! rotated  bilaplacian : triad operator (griffies)
-         CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2, Kmm )
+         CALL tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pt    , pt_rhs, kjpt, 2 )
       END SELECT
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10979)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90	(revision 10980)
@@ -48,7 +48,7 @@
 CONTAINS
 
-  SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv ,   &
-      &                                                     pgui, pgvi,   &
-      &                                         ptb , ptbb, pta , kjpt, kpass, Kmm )
+  SUBROUTINE tra_ldf_triad( kt, Kmm, kit000, cdtype, pahu, pahv,               &
+      &                                              pgu , pgv  , pgui, pgvi , &
+      &                                         pt , pt2, pt_rhs, kjpt, kpass )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_ldf_triad  ***
@@ -66,5 +66,5 @@
       !!      see documentation for the desciption
       !!
-      !! ** Action :   pta   updated with the before rotated diffusion
+      !! ** Action :   pt_rhs   updated with the before rotated diffusion
       !!               ah_wslp2 ....
       !!               akz   stabilizing vertical diffusivity coefficient (used in trazdf_imp)
@@ -79,7 +79,7 @@
       REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu , pgv  ! tracer gradient at pstep levels
       REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgui, pgvi ! tracer gradient at top   levels
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb        ! tracer (kpass=1) or laplacian of tracer (kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptbb       ! tracer (only used in kpass=2)
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta        ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt         ! tracer (kpass=1) or laplacian of tracer (kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt2        ! tracer (only used in kpass=2)
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs     ! tracer trend
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -230,6 +230,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zdit(ji,jj,jk) = ( ptb(ji+1,jj  ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk)
-                  zdjt(ji,jj,jk) = ( ptb(ji  ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
+                  zdit(ji,jj,jk) = ( pt(ji+1,jj  ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)
+                  zdjt(ji,jj,jk) = ( pt(ji  ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
                END DO
             END DO
@@ -258,9 +258,9 @@
          DO jk = 1, jpkm1
             !                    !==  Vertical tracer gradient at level jk and jk+1
-            zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
+            zdkt3d(:,:,1) = ( pt(:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)
             !
             !                    ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1)
             IF( jk == 1 ) THEN   ;   zdkt3d(:,:,0) = zdkt3d(:,:,1)
-            ELSE                 ;   zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)
+            ELSE                 ;   zdkt3d(:,:,0) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)
             ENDIF
             !
@@ -363,5 +363,5 @@
             DO jj = 2 , jpjm1
                DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  zftu(ji-1,jj,jk) - zftu(ji,jj,jk)       &
                      &                                           + zftv(ji,jj-1,jk) - zftv(ji,jj,jk)   )   &
                      &                                        / (  e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm)  )
@@ -378,5 +378,5 @@
                      ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)   &
                         &                            * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) )             &
-                        &                            * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                        &                            * (  pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                   END DO
                END DO
@@ -389,15 +389,15 @@
                      DO ji = fs_2, fs_jpim1
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)             &
-                           &                            * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) )
+                           &                            * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )
                      END DO
                   END DO
                END DO 
-            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb  and ptbb gradients, resp.
+            CASE(  2  )                            ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt  and pt2 gradients, resp.
                DO jk = 2, jpkm1 
                   DO jj = 1, jpjm1
                      DO ji = fs_2, fs_jpim1
                         ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk)                      &
-                           &                            * (  ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) )   &
-                           &                               + akz     (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) )   )
+                           &                            * (  ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) )   &
+                           &                               + akz     (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) )   )
                      END DO
                   END DO
@@ -406,9 +406,9 @@
          ENDIF
          !
-         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pta  ==!
+         DO jk = 1, jpkm1                 !==  Divergence of vertical fluxes added to pt_rhs  ==!
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1  ! vector opt.
-                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
-                     &                                        / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
+                  pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * (  ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk)  )   &
+                     &                                              / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) )
                END DO
             END DO
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 10979)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10980)
@@ -256,5 +256,5 @@
       IF( lrst_oce .AND. ln_zdfosm ) &
            &             CALL osm_rst( kstp,      Nnn, 'WRITE' )         ! write OSMOSIS outputs + wn (so must do here) to restarts
-                         CALL tra_ldf( kstp, Nbb, Nnn, Nrhs )            ! lateral mixing
+                         CALL tra_ldf( kstp, Nbb, Nnn      , ts, Nrhs )  ! lateral mixing
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90	(revision 10979)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcldf.F90	(revision 10980)
@@ -94,15 +94,15 @@
       !
       CASE ( np_lap   )                                                                                    ! iso-level laplacian
-         CALL tra_ldf_lap  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
-           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs),                 jptra, 1 , Kmm )
+         CALL tra_ldf_lap  ( kt, Kmm, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,            &
+           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs),                  jptra, 1 )
       CASE ( np_lap_i )                                                                                    ! laplacian : standard iso-neutral operator (Madec)
-         CALL tra_ldf_iso  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
-           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, 1 , Kmm )
+         CALL tra_ldf_iso  ( kt, Kmm, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,            &
+           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, 1 )
       CASE ( np_lap_it )                                                                                   ! laplacian : triad iso-neutral operator (griffies)
-         CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
-           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, 1 , Kmm )
+         CALL tra_ldf_triad( kt, Kmm, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,            &
+           &                     tr(:,:,:,:,Kbb), tr(:,:,:,:,Kbb), tr(:,:,:,:,Krhs), jptra, 1 )
       CASE ( np_blp , np_blp_i , np_blp_it )                                                               ! bilaplacian: all operator (iso-level, -neutral)
-         CALL tra_ldf_blp  ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,               &
-           &                     tr(:,:,:,:,Kbb) , tr(:,:,:,:,Krhs),                 jptra, nldf_trc, Kmm )
+         CALL tra_ldf_blp  ( kt, Kmm, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi,            &
+           &                     tr(:,:,:,:,Kbb) , tr(:,:,:,:,Krhs),                 jptra, nldf_trc )
       END SELECT
       !
