Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dtauvd.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dtauvd.F90	(revision 10998)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dtauvd.F90	(revision 11001)
@@ -117,5 +117,5 @@
 
 
-   SUBROUTINE dta_uvd( kt, puvd )
+   SUBROUTINE dta_uvd( kt, Kmm, puvd )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE dta_uvd  ***
@@ -133,4 +133,5 @@
       !!----------------------------------------------------------------------
       INTEGER                           , INTENT(in   ) ::   kt     ! ocean time-step
+      INTEGER                           , INTENT(in   ) ::   Kmm    ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk,2), INTENT(  out) ::   puvd   ! U & V current data
       !
@@ -160,5 +161,5 @@
             DO ji = 1, jpi                ! determines the interpolated U & V current profiles at each (i,j) point
                DO jk = 1, jpk
-                  zl = gdept_n(ji,jj,jk)
+                  zl = gdept(ji,jj,jk,Kmm)
                   IF    ( zl < gdept_1d(1  ) ) THEN          ! extrapolate above the first level of data
                      zup(jk) =  puvd(ji,jj,1    ,1)
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyncor_c1d.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyncor_c1d.F90	(revision 10998)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyncor_c1d.F90	(revision 11001)
@@ -56,5 +56,5 @@
 
 
-   SUBROUTINE dyn_cor_c1d( kt )
+   SUBROUTINE dyn_cor_c1d( kt, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE dyn_cor_c1d  ***
@@ -63,5 +63,7 @@
       !!               the general trend of the momentum equation in 1D case.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER                             , INTENT(in   ) ::   kt        ! ocean time-step index
+      INTEGER                             , INTENT(in   ) ::   Kmm, Krhs ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv  ! ocean velocities and RHS of momentum equation
       !!
       INTEGER ::   ji, jj, jk   ! dummy loop indices
@@ -78,6 +80,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * (vn(ji,jj,jk) + vsd(ji,jj,jk))
-                  va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * (un(ji,jj,jk) + usd(ji,jj,jk))
+                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ff_t(ji,jj) * (pvv(ji,jj,jk,Kmm) + vsd(ji,jj,jk))
+                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ff_t(ji,jj) * (puu(ji,jj,jk,Kmm) + usd(ji,jj,jk))
                END DO
             END DO
@@ -87,6 +89,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * vn(ji,jj,jk)
-                  va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * un(ji,jj,jk)
+                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ff_t(ji,jj) * pvv(ji,jj,jk,Kmm)
+                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ff_t(ji,jj) * puu(ji,jj,jk,Kmm)
                END DO
             END DO
@@ -95,6 +97,6 @@
       
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cor  - Ua: ', mask1=umask,  &
-         &                       tab3d_2=va, clinfo2=' Va: '       , mask2=vmask )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cor  - Ua: ', mask1=umask,  &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=' Va: '       , mask2=vmask )
       !
    END SUBROUTINE dyn_cor_c1d
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyndmp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyndmp.F90	(revision 10998)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/dyndmp.F90	(revision 11001)
@@ -130,5 +130,5 @@
 
 
-   SUBROUTINE dyn_dmp( kt )
+   SUBROUTINE dyn_dmp( kt, Kbb, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE dyn_dmp  ***
@@ -140,13 +140,15 @@
       !! ** Method  :   Compute Newtonian damping towards u_dta and v_dta 
       !!      and add to the general momentum trends:
-      !!                     ua = ua + resto_uv * (u_dta - ub)
-      !!                     va = va + resto_uv * (v_dta - vb)
+      !!                     puu(Krhs) = puu(Krhs) + resto_uv * (u_dta - puu(Kbb))
+      !!                     pvv(Krhs) = pvv(Krhs) + resto_uv * (v_dta - pvv(Kbb))
       !!      The trend is computed either throughout the water column
       !!      (nn_zdmp=0), where the vertical mixing is weak (nn_zdmp=1) or
       !!      below the well mixed layer (nn_zdmp=2)
       !!
-      !! ** Action  : - (ua,va)   momentum trends updated with the damping trend
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !! ** Action  : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs))   momentum trends updated with the damping trend
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time-step index
+      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm, Krhs ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv       ! ocean velocities and RHS of momentum equation
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -159,5 +161,5 @@
       !
       !                           !==   read and interpolate U & V current data at kt   ==!
-      CALL dta_uvd( kt, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
+      CALL dta_uvd( kt, Kmm, zuv_dta ) !!! NOTE: This subroutine must be altered for use outside
                                   !!!       the C1D context (use of U,V grid variables)
       !
@@ -168,8 +170,8 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
-                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-                  va(ji,jj,jk) = va(ji,jj,jk) + zva
+                  zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
+                  zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
+                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
+                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
                   utrdmp(ji,jj,jk) = zua           ! save the trends
                   vtrdmp(ji,jj,jk) = zva      
@@ -183,12 +185,12 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   IF( avt(ji,jj,jk) <= avt_c ) THEN
-                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
-                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
+                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
+                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
                   ELSE
                      zua = 0._wp
                      zva = 0._wp  
                   ENDIF
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-                  va(ji,jj,jk) = va(ji,jj,jk) + zva
+                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
+                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
                   utrdmp(ji,jj,jk) = zua           ! save the trends
                   vtrdmp(ji,jj,jk) = zva
@@ -201,13 +203,13 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  IF( gdept_n(ji,jj,jk) >= hmlp (ji,jj) ) THEN
-                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - ub(ji,jj,jk) )
-                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - vb(ji,jj,jk) )
+                  IF( gdept(ji,jj,jk,Kmm) >= hmlp (ji,jj) ) THEN
+                     zua = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,1) - puu(ji,jj,jk,Kbb) )
+                     zva = resto_uv(ji,jj,jk) * ( zuv_dta(ji,jj,jk,2) - pvv(ji,jj,jk,Kbb) )
                   ELSE
                      zua = 0._wp
                      zva = 0._wp  
                   ENDIF
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
-                  va(ji,jj,jk) = va(ji,jj,jk) + zva
+                  puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zua
+                  pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zva
                   utrdmp(ji,jj,jk) = zua           ! save the trends
                   vtrdmp(ji,jj,jk) = zva
@@ -219,6 +221,6 @@
       !
       !                           ! Control print
-      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=ua(:,:,:), clinfo1=' dmp  - Ua: ', mask1=umask,   &
-         &                           tab3d_2=va(:,:,:), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF( ln_ctl   )   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' dmp  - Ua: ', mask1=umask,   &
+         &                           tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90	(revision 10998)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90	(revision 11001)
@@ -69,13 +69,13 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Ocean physics update                (ua, va, ta, sa used as workspace)
+      ! Ocean physics update        
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL eos_rab( tsb, rab_b )   ! before local thermal/haline expension ratio at T-points
-                         CALL eos_rab( tsn, rab_n )   ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
-                         CALL bn2( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
+                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )  ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )  ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2( ts(:,:,:,:,Nnn), rab_n, rn2 , Nnn ) ! now    Brunt-Vaisala frequency
       
       !  VERTICAL PHYSICS
-                         CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
+                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs  )    ! vertical physics update (bfr, avt, avs, avm + MLD)
 
       IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )  ! after ssh (includes call to div_hor)
@@ -84,5 +84,5 @@
       IF(.NOT.ln_linssh )   CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )  ! now cross-level velocity 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
+      ! diagnostics and outputs       
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL dia_wri( kstp, Nnn )  ! ocean model: outputs
@@ -94,21 +94,21 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                        CALL trc_stp( kstp )       ! time-stepping
+                        CALL trc_stp( kstp, Nbb, Nnn, Nrhs, Naa  )   ! time-stepping
 #endif
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! Active tracers                              (ua, va used as workspace)
+      ! Active tracers                              (uu(:,:,:,Nrhs), vv(:,:,:,Nrhs) used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                        tsa(:,:,:,:) = 0._wp       ! set tracer trends to zero
+                        ts(:,:,:,:,Nrhs) = 0._wp       ! set tracer trends to zero
 
-                        CALL tra_sbc( kstp )       ! surface boundary condition
-      IF( ln_traqsr )   CALL tra_qsr( kstp )       ! penetrative solar radiation qsr
-      IF( ln_tradmp )   CALL tra_dmp( kstp )       ! internal damping trends- tracers
-      IF(.NOT.ln_linssh)CALL tra_adv( kstp )       ! horizontal & vertical advection
-      IF( ln_zdfosm  )  CALL tra_osm( kstp )       ! OSMOSIS non-local tracer fluxes
-                        CALL tra_zdf( kstp )       ! vertical mixing
-                        CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) )   ! now potential density for zdfmxl
-      IF( ln_zdfnpc )   CALL tra_npc( kstp )       ! applied non penetrative convective adjustment on (t,s)
-                        CALL tra_nxt( kstp )       ! tracer fields at next time step
+                        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_tradmp )   CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs  )  ! internal damping trends- tracers
+      IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection
+      IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes
+                        CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   ) ! vertical mixing
+                        CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) )  ! now potential density for zdfmxl
+      IF( ln_zdfnpc )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa   ) ! applied non penetrative convective adjustment on (t,s)
+                        CALL tra_nxt( kstp, Nbb, Nnn, Nrhs,     Naa   ) ! tracer fields at next time step
 
 
@@ -117,22 +117,22 @@
       ! Dynamics                                    (ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                        ua(:,:,:) = 0._wp          ! set dynamics trends to zero
-                        va(:,:,:) = 0._wp
+                        uu(:,:,:,Nrhs) = 0._wp          ! set dynamics trends to zero
+                        vv(:,:,:,Nrhs) = 0._wp
 
-      IF( ln_dyndmp )   CALL dyn_dmp    ( kstp )   ! internal damping trends- momentum
-                        CALL dyn_cor_c1d( kstp )   ! vorticity term including Coriolis
-      IF( ln_zdfosm  )  CALL dyn_osm    ( kstp )   ! OSMOSIS non-local velocity fluxes
-                        CALL dyn_zdf    ( kstp )   ! vertical diffusion
-                        CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step
-      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )   ! swap of sea surface height
+      IF( ln_dyndmp )   CALL dyn_dmp    ( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! internal damping trends- momentum
+                        CALL dyn_cor_c1d( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity term including Coriolis
+      IF( ln_zdfosm  )  CALL dyn_osm    ( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes
+                        CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
+                        CALL dyn_nxt    ( kstp, Nbb, Nnn, Naa )                 ! lateral velocity at next time step
+      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp, Nbb, Nnn, Naa )                 ! swap of sea surface height
 
-      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )! swap of vertical scale factors
+      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa )              ! swap of vertical scale factors
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control and restarts
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                             CALL stp_ctl( kstp, indic )
-      IF( kstp == nit000 )   CALL iom_close( numror )      ! close input  ocean restart file
-      IF( lrst_oce       )   CALL rst_write( kstp )        ! write output ocean restart file
+                             CALL stp_ctl( kstp, Nnn, indic )
+      IF( kstp == nit000 )   CALL iom_close( numror )          ! close input  ocean restart file
+      IF( lrst_oce       )   CALL rst_write( kstp, Nbb, Nnn )  ! write output ocean restart file
       !
 #if defined key_iomput
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 10998)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 11001)
@@ -125,5 +125,5 @@
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-      IF( ln_sto_eos ) CALL sto_pts( tsn  )          ! Random T/S fluctuations
+      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -131,8 +131,8 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  THERMODYNAMICS
-                         CALL eos_rab( tsb, rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
-                         CALL eos_rab( tsn, rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2    ( tsb, rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
-                         CALL bn2    ( tsn, rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
+                         CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn )       ! now    local thermal/haline expension ratio at T-points
+                         CALL bn2    ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
+                         CALL bn2    ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn  ) ! now    Brunt-Vaisala frequency
 
       !  VERTICAL PHYSICS
@@ -142,12 +142,12 @@
       !
       IF( l_ldfslp ) THEN                             ! slope of lateral mixing
-                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
+                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
 
          IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
-            &            CALL zps_hde    ( kstp, Nnn, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
             &                                          rhd, gru , grv    )       ! of t, s, rd at the last ocean level
 
          IF( ln_zps .AND.       ln_isfcav)                                                &
-            &            CALL zps_hde_isf( kstp, Nnn, jpts, tsb, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+            &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
             &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
          IF( ln_traldf_triad ) THEN 
@@ -169,5 +169,5 @@
                             CALL wzv           ( kstp, Nbb, Nnn, ww,  Naa )    ! now cross-level velocity 
       IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
-                            CALL eos    ( tsn, rhd, rhop, gdept_n(:,:,:) )  ! now in situ density for hpg computation
+                            CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
                             
 !!jc: fs simplification
@@ -176,13 +176,13 @@
 !!                                         with previous versions using split-explicit free surface          
             IF( ln_zps .AND. .NOT. ln_isfcav )                                    &
-               &            CALL zps_hde    ( kstp, Nnn, jpts, tsn, gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
+               &            CALL zps_hde    ( kstp, Nnn, jpts, ts(:,:,:,:,Nnn), gtsu, gtsv,   &  ! Partial steps: before horizontal gradient
                &                                          rhd, gru , grv     )       ! of t, s, rd at the last ocean level
             IF( ln_zps .AND.       ln_isfcav )                                               &
-               &            CALL zps_hde_isf( kstp, Nnn, jpts, tsn, gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
+               &            CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nnn), gtsu, gtsv, gtui, gtvi,  &  ! Partial steps for top cell (ISF)
                &                                          rhd, gru , grv , grui, grvi   )       ! of t, s, rd at the first ocean level
 !!jc: fs simplification
                             
-                         ua(:,:,:) = 0._wp            ! set dynamics trends to zero
-                         va(:,:,:) = 0._wp
+                         uu(:,:,:,Nrhs) = 0._wp            ! set dynamics trends to zero
+                         vv(:,:,:,Nrhs) = 0._wp
 
       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
@@ -200,5 +200,5 @@
                          CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
 
-                                                      ! With split-explicit free surface, since now transports have been updated and ssha as well
+                                                      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
       IF( ln_dynspg_ts ) THEN                         ! vertical scale factors and vertical velocity need to be updated
                             CALL div_hor       ( kstp, Nbb, Nnn )    ! Horizontal divergence  (2nd call in time-split case)
@@ -238,29 +238,29 @@
       ! Active tracers                              
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         tsa(:,:,:,:) = 0._wp         ! set tracer trends to zero
+                         ts(:,:,:,:,Nrhs) = 0._wp         ! set tracer trends to zero
 
       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, 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
+         & ln_trainc )   CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs )  ! apply tracer assimilation increment
+                         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())  & 
                &         CALL Agrif_Sponge_tra        ! tracers sponge
 #endif
-                         CALL tra_adv( kstp, Nbb, Nnn      , ts, Nrhs )  ! hor. + vert. advection	==> RHS
-      IF( ln_zdfosm  )   CALL tra_osm( kstp,      Nnn      , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
+                         CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
+      IF( ln_zdfosm  )   CALL tra_osm    ( kstp, Nnn     , ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
       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      , ts, Nrhs )  ! lateral mixing
+           &             CALL osm_rst    ( kstp, Nnn, 'WRITE' )        ! write OSMOSIS outputs + ww (so must do here) to restarts
+                         CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
       IF( ln_diaptr  )   CALL dia_ptr( Nnn )                 ! Poleward adv/ldf TRansports diagnostics
 !!gm
-                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer	==> after
-      IF( ln_zdfnpc  )   CALL tra_npc( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
+                         CALL tra_zdf    ( kstp, Nbb, Nnn, Nrhs, ts, Naa  )  ! vert. mixing & after tracer	==> after
+      IF( ln_zdfnpc  )   CALL tra_npc    ( kstp,      Nnn, Nrhs, ts, Naa  )  ! update after fields by non-penetrative convection
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -280,5 +280,5 @@
 !!    place.
 !! 
-!!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
+!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
                          CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap
                          CALL dyn_nxt       ( kstp, Nbb, Nnn, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
