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 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90	(revision 10965)
@@ -14,4 +14,5 @@
    !!----------------------------------------------------------------------
    USE step_oce        ! time stepping definition modules 
+   USE step, ONLY : Nbb, Nnn, Naa, Nrhs ! time level indices
 #if defined key_top
    USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
@@ -85,6 +86,6 @@
       ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL dia_wri( kstp )       ! ocean model: outputs
-      IF( lk_diahth  )   CALL dia_hth( kstp )       ! Thermocline depth (20°C)
+                         CALL dia_wri( kstp, Nnn )  ! ocean model: outputs
+      IF( lk_diahth  )   CALL dia_hth( kstp, Nnn )  ! Thermocline depth (20°C)
 
 
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/dia25h.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/dia25h.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/dia25h.F90	(revision 10965)
@@ -39,5 +39,5 @@
 CONTAINS
 
-   SUBROUTINE dia_25h_init 
+   SUBROUTINE dia_25h_init( Kbb )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE dia_25h_init  ***
@@ -47,4 +47,6 @@
       !! ** Method : Read namelist
       !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb       ! Time level index
+      !
       INTEGER ::   ios                 ! Local integer output status for namelist read
       INTEGER ::   ierror              ! Local integer for memory allocation
@@ -95,9 +97,9 @@
       ! ------------------------- !
       cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem)
-      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal)
-      sshn_25h(:,:)   = sshb(:,:)
-      un_25h  (:,:,:) = ub  (:,:,:)
-      vn_25h  (:,:,:) = vb  (:,:,:)
+      tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kbb)
+      sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kbb)
+      sshn_25h(:,:)   = ssh(:,:,Kbb)
+      un_25h  (:,:,:) = uu  (:,:,:,Kbb)
+      vn_25h  (:,:,:) = vv  (:,:,:,Kbb)
       avt_25h (:,:,:) = avt (:,:,:)
       avm_25h (:,:,:) = avm (:,:,:)
@@ -116,5 +118,5 @@
 
 
-   SUBROUTINE dia_25h( kt )  
+   SUBROUTINE dia_25h( kt, Kmm )  
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_25h  ***
@@ -125,4 +127,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !!
       INTEGER ::   ji, jj, jk
@@ -150,5 +153,5 @@
       ! wn_25h could not be initialised in dia_25h_init, so we do it here instead
       IF( kt == nn_it000 ) THEN
-         wn_25h(:,:,:) = wn(:,:,:)
+         wn_25h(:,:,:) = ww(:,:,:)
       ENDIF
 
@@ -161,10 +164,10 @@
          ENDIF
 
-         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem)
-         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal)
-         sshn_25h(:,:)       = sshn_25h(:,:)   + sshn(:,:)
-         un_25h  (:,:,:)     = un_25h  (:,:,:) + un  (:,:,:)
-         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vn  (:,:,:)
-         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + wn  (:,:,:)
+         tn_25h  (:,:,:)     = tn_25h  (:,:,:) + ts (:,:,:,jp_tem,Kmm)
+         sn_25h  (:,:,:)     = sn_25h  (:,:,:) + ts (:,:,:,jp_sal,Kmm)
+         sshn_25h(:,:)       = sshn_25h(:,:)   + ssh(:,:,Kmm)
+         un_25h  (:,:,:)     = un_25h  (:,:,:) + uu  (:,:,:,Kmm)
+         vn_25h  (:,:,:)     = vn_25h  (:,:,:) + vv  (:,:,:,Kmm)
+         wn_25h  (:,:,:)     = wn_25h  (:,:,:) + ww  (:,:,:)
          avt_25h (:,:,:)     = avt_25h (:,:,:) + avt (:,:,:)
          avm_25h (:,:,:)     = avm_25h (:,:,:) + avm (:,:,:)
@@ -245,10 +248,10 @@
          !
          ! After the write reset the values to cnt=1 and sum values equal current value 
-         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem)
-         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal)
-         sshn_25h(:,:)   = sshn(:,:)
-         un_25h  (:,:,:) = un  (:,:,:)
-         vn_25h  (:,:,:) = vn  (:,:,:)
-         wn_25h  (:,:,:) = wn  (:,:,:)
+         tn_25h  (:,:,:) = ts (:,:,:,jp_tem,Kmm)
+         sn_25h  (:,:,:) = ts (:,:,:,jp_sal,Kmm)
+         sshn_25h(:,:)   = ssh(:,:,Kmm)
+         un_25h  (:,:,:) = uu  (:,:,:,Kmm)
+         vn_25h  (:,:,:) = vv  (:,:,:,Kmm)
+         wn_25h  (:,:,:) = ww  (:,:,:)
          avt_25h (:,:,:) = avt (:,:,:)
          avm_25h (:,:,:) = avm (:,:,:)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaar5.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaar5.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaar5.F90	(revision 10965)
@@ -62,5 +62,5 @@
 
 
-   SUBROUTINE dia_ar5( kt )
+   SUBROUTINE dia_ar5( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ar5  ***
@@ -70,4 +70,5 @@
       !
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
       !
       INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
@@ -89,5 +90,5 @@
          ALLOCATE( zrhd(jpi,jpj,jpk) , zrhop(jpi,jpj,jpk) )
          ALLOCATE( ztsn(jpi,jpj,jpk,jpts) )
-         zarea_ssh(:,:) = area(:,:) * sshn(:,:)
+         zarea_ssh(:,:) = area(:,:) * ssh(:,:,Kmm)
       ENDIF
       !
@@ -100,5 +101,5 @@
          CALL iom_put( 'voltot', zvol               )
          CALL iom_put( 'sshtot', zvolssh / area_tot )
-         CALL iom_put( 'sshdyn', sshn(:,:) - (zvolssh / area_tot) )
+         CALL iom_put( 'sshdyn', ssh(:,:,Kmm) - (zvolssh / area_tot) )
          !
       ENDIF
@@ -106,11 +107,11 @@
       IF( iom_use( 'botpres' ) .OR. iom_use( 'sshthster' )  .OR. iom_use( 'sshsteric' )  ) THEN    
          !                     
-         ztsn(:,:,:,jp_tem) = tsn(:,:,:,jp_tem)                    ! thermosteric ssh
+         ztsn(:,:,:,jp_tem) = ts(:,:,:,jp_tem,Kmm)                    ! thermosteric ssh
          ztsn(:,:,:,jp_sal) = sn0(:,:,:)
-         CALL eos( ztsn, zrhd, gdept_n(:,:,:) )                       ! now in situ density using initial salinity
+         CALL eos( ztsn, zrhd, gdept(:,:,:,Kmm) )                       ! now in situ density using initial salinity
          !
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
          DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)
          END DO
          IF( ln_linssh ) THEN
@@ -118,9 +119,9 @@
                DO ji = 1, jpi
                   DO jj = 1, jpj
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
             END IF
 !!gm
@@ -135,5 +136,5 @@
       
          !                                         ! steric sea surface height
-         CALL eos( tsn, zrhd, zrhop, gdept_n(:,:,:) )                 ! now in situ and potential density
+         CALL eos( ts(:,:,:,:,Kmm), zrhd, zrhop, gdept(:,:,:,Kmm) )                 ! now in situ and potential density
          zrhop(:,:,jpk) = 0._wp
          CALL iom_put( 'rhop', zrhop )
@@ -141,5 +142,5 @@
          zbotpres(:,:) = 0._wp                        ! no atmospheric surface pressure, levitating sea-ice
          DO jk = 1, jpkm1
-            zbotpres(:,:) = zbotpres(:,:) + e3t_n(:,:,jk) * zrhd(:,:,jk)
+            zbotpres(:,:) = zbotpres(:,:) + e3t(:,:,jk,Kmm) * zrhd(:,:,jk)
          END DO
          IF( ln_linssh ) THEN
@@ -147,9 +148,9 @@
                DO ji = 1,jpi
                   DO jj = 1,jpj
-                     zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+                     zbotpres(ji,jj) = zbotpres(ji,jj) + ssh(ji,jj,Kmm) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
                   END DO
                END DO
             ELSE
-               zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+               zbotpres(:,:) = zbotpres(:,:) + ssh(:,:,Kmm) * zrhd(:,:,1)
             END IF
          END IF
@@ -162,5 +163,5 @@
          !                                         ! ocean bottom pressure
          zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
-         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
+         zbotpres(:,:) = zztmp * ( zbotpres(:,:) + ssh(:,:,Kmm) + thick0(:,:) )
          CALL iom_put( 'botpres', zbotpres )
          !
@@ -174,7 +175,7 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zztmp = area(ji,jj) * e3t_n(ji,jj,jk)
-                  ztemp = ztemp + zztmp * tsn(ji,jj,jk,jp_tem)
-                  zsal  = zsal  + zztmp * tsn(ji,jj,jk,jp_sal)
+                  zztmp = area(ji,jj) * e3t(ji,jj,jk,Kmm)
+                  ztemp = ztemp + zztmp * ts(ji,jj,jk,jp_tem,Kmm)
+                  zsal  = zsal  + zztmp * ts(ji,jj,jk,jp_sal,Kmm)
                END DO
             END DO
@@ -184,11 +185,11 @@
                DO ji = 1, jpi
                   DO jj = 1, jpj
-                     ztemp = ztemp + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_tem) 
-                     zsal  = zsal  + zarea_ssh(ji,jj) * tsn(ji,jj,mikt(ji,jj),jp_sal) 
+                     ztemp = ztemp + zarea_ssh(ji,jj) * ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) 
+                     zsal  = zsal  + zarea_ssh(ji,jj) * ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) 
                   END DO
                END DO
             ELSE
-               ztemp = ztemp + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_tem) )
-               zsal  = zsal  + SUM( zarea_ssh(:,:) * tsn(:,:,1,jp_sal) )
+               ztemp = ztemp + SUM( zarea_ssh(:,:) * ts(:,:,1,jp_tem,Kmm) )
+               zsal  = zsal  + SUM( zarea_ssh(:,:) * ts(:,:,1,jp_sal,Kmm) )
             END IF
          ENDIF
@@ -219,8 +220,8 @@
                   DO ji = 1, jpi
                      IF( rn2(ji,jj,jk) > 0._wp ) THEN
-                        zrw =   ( gdepw_n(ji,jj,jk  ) - gdept_n(ji,jj,jk) )   &
-                           &  / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
+                        zrw =   ( gdepw(ji,jj,jk  ,Kmm) - gdept(ji,jj,jk,Kmm) )   &
+                           &  / ( gdept(ji,jj,jk-1,Kmm) - gdept(ji,jj,jk,Kmm) )
 !!gm  this can be reduced to :  (depw-dept) / e3w   (NB idem dans bn2 !)
-!                        zrw =   ( gdept_n(ji,jj,jk) - gdepw_n(ji,jj,jk) ) / e3w_n(ji,jj,jk)
+!                        zrw =   ( gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm) ) / e3w(ji,jj,jk,Kmm)
 !!gm end
                         !
@@ -229,6 +230,6 @@
                         !
                         zpe(ji, jj) = zpe(ji, jj)            &
-                           &        -  grav * (  avt(ji,jj,jk) * zaw * (tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) )  &
-                           &                   - avs(ji,jj,jk) * zbw * (tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal) ) )
+                           &        -  grav * (  avt(ji,jj,jk) * zaw * (ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) )  &
+                           &                   - avs(ji,jj,jk) * zbw * (ts(ji,jj,jk-1,jp_sal,Kmm) - ts(ji,jj,jk,jp_sal,Kmm) ) )
                      ENDIF
                   END DO
@@ -239,5 +240,5 @@
                DO ji = 1, jpi
                   DO jj = 1, jpj
-                     zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w_n(ji, jj, jk)
+                     zpe(ji,jj) = zpe(ji,jj) + avt(ji, jj, jk) * MIN(0._wp,rn2(ji, jj, jk)) * rau0 * e3w(ji, jj, jk,Kmm)
                   END DO
                END DO
@@ -261,5 +262,5 @@
 
 
-   SUBROUTINE dia_ar5_hst( ktra, cptr, pua, pva ) 
+   SUBROUTINE dia_ar5_hst( ktra, cptr, puflx, pvflx ) 
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ar5_htr ***
@@ -270,6 +271,6 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pua   ! 3D input array of advection/diffusion
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: puflx  ! u-flux of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx  ! v-flux of advection/diffusion
       !
       INTEGER    ::  ji, jj, jk
@@ -277,9 +278,9 @@
 
     
-      z2d(:,:) = pua(:,:,1) 
+      z2d(:,:) = puflx(:,:,1) 
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               z2d(ji,jj) = z2d(ji,jj) + pua(ji,jj,jk) 
+               z2d(ji,jj) = z2d(ji,jj) + puflx(ji,jj,jk) 
             END DO
          END DO
@@ -295,9 +296,9 @@
        ENDIF
        !
-       z2d(:,:) = pva(:,:,1) 
+       z2d(:,:) = pvflx(:,:,1) 
        DO jk = 1, jpkm1
           DO jj = 2, jpjm1
              DO ji = fs_2, fs_jpim1   ! vector opt.
-                z2d(ji,jj) = z2d(ji,jj) + pva(ji,jj,jk) 
+                z2d(ji,jj) = z2d(ji,jj) + pvflx(ji,jj,jk) 
              END DO
           END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diacfl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diacfl.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diacfl.F90	(revision 10965)
@@ -45,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE dia_cfl ( kt )
+   SUBROUTINE dia_cfl ( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dia_cfl  ***
@@ -53,4 +53,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       INTEGER                ::   ji, jj, jk                            ! dummy loop indices
@@ -71,7 +72,7 @@
          DO jj = 1, jpj
             DO ji = 1, fs_jpim1   ! vector opt.
-               zCu_cfl(ji,jj,jk) = ABS( un(ji,jj,jk) ) * z2dt / e1u  (ji,jj)      ! for i-direction
-               zCv_cfl(ji,jj,jk) = ABS( vn(ji,jj,jk) ) * z2dt / e2v  (ji,jj)      ! for j-direction
-               zCw_cfl(ji,jj,jk) = ABS( wn(ji,jj,jk) ) * z2dt / e3w_n(ji,jj,jk)   ! for k-direction
+               zCu_cfl(ji,jj,jk) = ABS( uu(ji,jj,jk,Kmm) ) * z2dt / e1u  (ji,jj)      ! for i-direction
+               zCv_cfl(ji,jj,jk) = ABS( vv(ji,jj,jk,Kmm) ) * z2dt / e2v  (ji,jj)      ! for j-direction
+               zCw_cfl(ji,jj,jk) = ABS( ww(ji,jj,jk) ) * z2dt / e3w(ji,jj,jk,Kmm)   ! for k-direction
             END DO
          END DO         
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diadct.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diadct.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diadct.F90	(revision 10965)
@@ -178,5 +178,5 @@
  
  
-  SUBROUTINE dia_dct( kt )
+  SUBROUTINE dia_dct( kt, Kmm )
      !!---------------------------------------------------------------------
      !!               ***  ROUTINE diadct  ***  
@@ -195,5 +195,6 @@
      !!               Reinitialise all relevant arrays to zero 
      !!---------------------------------------------------------------------
-     INTEGER, INTENT(in) ::   kt
+     INTEGER, INTENT(in) ::   kt    ! ocean time step
+     INTEGER, INTENT(in) ::   Kmm   ! time level index
      !
      INTEGER ::   jsec              ! loop on sections
@@ -235,5 +236,5 @@
 
            !Compute transport through section  
-           CALL transport(secs(jsec),lldebug,jsec) 
+           CALL transport(Kmm,secs(jsec),lldebug,jsec) 
 
         ENDDO
@@ -249,5 +250,5 @@
            ! Sum over each class 
            DO jsec=1,nb_sec 
-              CALL dia_dct_sum(secs(jsec),jsec) 
+              CALL dia_dct_sum(Kmm,secs(jsec),jsec) 
            ENDDO 
 
@@ -561,5 +562,5 @@
 
 
-   SUBROUTINE transport(sec,ld_debug,jsec)
+   SUBROUTINE transport(Kmm,sec,ld_debug,jsec)
      !!-------------------------------------------------------------------------------------------
      !!                     ***  ROUTINE transport  ***
@@ -581,4 +582,5 @@
      !!
      !!-------------------------------------------------------------------------------------------
+     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index
      TYPE(SECTION),INTENT(INOUT) :: sec
      LOGICAL      ,INTENT(IN)    :: ld_debug
@@ -676,25 +678,25 @@
             SELECT CASE( sec%direction(jseg) )
                CASE(0,1) 
-                  ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
-                  zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
-                  zrhop = interp(k%I,k%J,jk,'V',rhop) 
-                  zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
-                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1) 
+                  ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) ) 
+                  zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
+                  zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I,k%J+1,Kmm)    ) * vmask(k%I,k%J,1) 
                CASE(2,3) 
-                  ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
-                  zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
-                  zrhop = interp(k%I,k%J,jk,'U',rhop) 
-                  zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
-                  zsshn =  0.5*( sshn(k%I,k%J) + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+                  ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) ) 
+                  zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
+                  zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
+                  zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                  zsshn =  0.5*( ssh(k%I,k%J,Kmm) + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
                END SELECT 
                !
-               zdep= gdept_n(k%I,k%J,jk) 
+               zdep= gdept(k%I,k%J,jk,Kmm) 
   
                SELECT CASE( sec%direction(jseg) )                !compute velocity with the correct direction 
                CASE(0,1)   
                   zumid=0._wp
-                  zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk) 
+                  zvmid=isgnv*vv(k%I,k%J,jk,Kmm)*vmask(k%I,k%J,jk) 
                CASE(2,3) 
-                  zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk) 
+                  zumid=isgnu*uu(k%I,k%J,jk,Kmm)*umask(k%I,k%J,jk) 
                   zvmid=0._wp
                END SELECT 
@@ -702,6 +704,6 @@
                !zTnorm=transport through one cell; 
                !velocity* cell's length * cell's thickness 
-               zTnorm = zumid*e2u(k%I,k%J) * e3u_n(k%I,k%J,jk)     & 
-                  &   + zvmid*e1v(k%I,k%J) * e3v_n(k%I,k%J,jk) 
+               zTnorm = zumid*e2u(k%I,k%J) * e3u(k%I,k%J,jk,Kmm)     & 
+                  &   + zvmid*e1v(k%I,k%J) * e3v(k%I,k%J,jk,Kmm) 
 
 !!gm  THIS is WRONG  no transport due to ssh in linear free surface case !!!!!
@@ -768,5 +770,5 @@
 
 
-  SUBROUTINE dia_dct_sum(sec,jsec) 
+  SUBROUTINE dia_dct_sum(Kmm,sec,jsec) 
      !!------------------------------------------------------------- 
      !! Purpose: Average the transport over nn_dctwri time steps  
@@ -787,4 +789,5 @@
      !! 
      !!------------------------------------------------------------- 
+     INTEGER      ,INTENT(IN)    :: Kmm         ! time level index
      TYPE(SECTION),INTENT(INOUT) :: sec 
      INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section 
@@ -848,18 +851,18 @@
               SELECT CASE( sec%direction(jseg) ) 
               CASE(0,1) 
-                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) ) 
-                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) ) 
-                 zrhop = interp(k%I,k%J,jk,'V',rhop) 
-                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0) 
+                 ztn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_tem,Kmm) ) 
+                 zsn   = interp(Kmm,k%I,k%J,jk,'V',ts(:,:,:,jp_sal,Kmm) ) 
+                 zrhop = interp(Kmm,k%I,k%J,jk,'V',rhop) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'V',rhd*rau0+rau0) 
 
               CASE(2,3) 
-                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) ) 
-                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) ) 
-                 zrhop = interp(k%I,k%J,jk,'U',rhop) 
-                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0) 
-                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1)  
+                 ztn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_tem,Kmm) ) 
+                 zsn   = interp(Kmm,k%I,k%J,jk,'U',ts(:,:,:,jp_sal,Kmm) ) 
+                 zrhop = interp(Kmm,k%I,k%J,jk,'U',rhop) 
+                 zrhoi = interp(Kmm,k%I,k%J,jk,'U',rhd*rau0+rau0) 
+                 zsshn =  0.5*( ssh(k%I,k%J,Kmm)    + ssh(k%I+1,k%J,Kmm)    ) * umask(k%I,k%J,1)  
               END SELECT 
  
-              zdep= gdept_n(k%I,k%J,jk) 
+              zdep= gdept(k%I,k%J,jk,Kmm) 
   
               !------------------------------- 
@@ -1104,5 +1107,5 @@
 
 
-   FUNCTION interp(ki, kj, kk, cd_point, ptab)
+   FUNCTION interp(Kmm, ki, kj, kk, cd_point, ptab)
   !!----------------------------------------------------------------------
   !!
@@ -1165,4 +1168,5 @@
   !!----------------------------------------------------------------------
   !*arguments
+  INTEGER, INTENT(IN)                          :: Kmm          ! time level index
   INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
   CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
@@ -1199,7 +1203,7 @@
   IF( ln_sco )THEN   ! s-coordinate case
 
-     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) * 0.5_wp 
-     zdep1 = gdept_n(ii1,ij1,kk) - zdepu
-     zdep2 = gdept_n(ii2,ij2,kk) - zdepu
+     zdepu = ( gdept(ii1,ij1,kk,Kmm) +  gdept(ii2,ij2,kk,Kmm) ) * 0.5_wp 
+     zdep1 = gdept(ii1,ij1,kk,Kmm) - zdepu
+     zdep2 = gdept(ii2,ij2,kk,Kmm) - zdepu
 
      ! weights
@@ -1213,7 +1217,7 @@
   ELSE       ! full step or partial step case 
 
-     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
-     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
-     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
+     ze3t  = e3t(ii2,ij2,kk,Kmm) - e3t(ii1,ij1,kk,Kmm) 
+     zwgt1 = ( e3w(ii2,ij2,kk,Kmm) - e3w(ii1,ij1,kk,Kmm) ) / e3w(ii2,ij2,kk,Kmm)
+     zwgt2 = ( e3w(ii1,ij1,kk,Kmm) - e3w(ii2,ij2,kk,Kmm) ) / e3w(ii1,ij1,kk,Kmm)
 
      IF(kk .NE. 1)THEN
@@ -1253,7 +1257,8 @@
    END SUBROUTINE dia_dct_init
 
-   SUBROUTINE dia_dct( kt )         ! Dummy routine
+   SUBROUTINE dia_dct( kt, Kmm )         ! Dummy routine
       IMPLICIT NONE
       INTEGER, INTENT( in ) :: kt   ! ocean time-step index
+      INTEGER, INTENT( in ) :: Kmm  ! ocean time level index
       WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
    END SUBROUTINE dia_dct
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaharm.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaharm.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaharm.F90	(revision 10965)
@@ -162,5 +162,5 @@
 
 
-   SUBROUTINE dia_harm ( kt )
+   SUBROUTINE dia_harm ( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_harm  ***
@@ -172,4 +172,5 @@
       !!--------------------------------------------------------------------
       INTEGER, INTENT( IN ) :: kt
+      INTEGER, INTENT( IN ) :: Kmm     ! time level index
       !
       INTEGER  :: ji, jj, jh, jc, nhc
@@ -194,7 +195,7 @@
                   DO ji = 1,jpi
                      ! Elevation
-                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj)        
-                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj)
-                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj)
+                     ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*ssh(ji,jj,Kmm)*ssmask (ji,jj)        
+                     ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*uu_b(ji,jj,Kmm)*ssumask(ji,jj)
+                     ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vv_b(ji,jj,Kmm)*ssvmask(ji,jj)
                   END DO
                END DO
@@ -521,6 +522,7 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_diaharm = .FALSE.
 CONTAINS
-   SUBROUTINE dia_harm ( kt )     ! Empty routine
+   SUBROUTINE dia_harm ( kt, Kmm )     ! Empty routine
       INTEGER, INTENT( IN ) :: kt  
+      INTEGER, INTENT( IN ) :: Kmm  
       WRITE(*,*) 'dia_harm: you should not have seen this print'
    END SUBROUTINE dia_harm
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diahsb.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diahsb.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diahsb.F90	(revision 10965)
@@ -58,5 +58,5 @@
 CONTAINS
 
-   SUBROUTINE dia_hsb( kt )
+   SUBROUTINE dia_hsb( kt, Kbb, Kmm )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE dia_hsb  ***
@@ -69,5 +69,6 @@
       !!
       !!---------------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER, INTENT(in) ::   kt         ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
       !
       INTEGER    ::   ji, jj, jk                  ! dummy loop indice
@@ -86,6 +87,6 @@
       IF( ln_timing )   CALL timing_start('dia_hsb')      
       !
-      tsn(:,:,:,1) = tsn(:,:,:,1) * tmask(:,:,:) ; tsb(:,:,:,1) = tsb(:,:,:,1) * tmask(:,:,:) ;
-      tsn(:,:,:,2) = tsn(:,:,:,2) * tmask(:,:,:) ; tsb(:,:,:,2) = tsb(:,:,:,2) * tmask(:,:,:) ;
+      ts(:,:,:,1,Kmm) = ts(:,:,:,1,Kmm) * tmask(:,:,:) ; ts(:,:,:,1,Kbb) = ts(:,:,:,1,Kbb) * tmask(:,:,:) ;
+      ts(:,:,:,2,Kmm) = ts(:,:,:,2,Kmm) * tmask(:,:,:) ; ts(:,:,:,2,Kbb) = ts(:,:,:,2,Kbb) * tmask(:,:,:) ;
       ! ------------------------- !
       ! 1 - Trends due to forcing !
@@ -108,11 +109,11 @@
             DO ji=1,jpi
                DO jj=1,jpj
-                  z2d0(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_tem)
-                  z2d1(ji,jj) = surf(ji,jj) * wn(ji,jj,mikt(ji,jj)) * tsb(ji,jj,mikt(ji,jj),jp_sal)
+                  z2d0(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_tem,Kbb)
+                  z2d1(ji,jj) = surf(ji,jj) * ww(ji,jj,mikt(ji,jj)) * ts(ji,jj,mikt(ji,jj),jp_sal,Kbb)
                END DO
             END DO
          ELSE
-            z2d0(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_tem)
-            z2d1(:,:) = surf(:,:) * wn(:,:,1) * tsb(:,:,1,jp_sal)
+            z2d0(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_tem,Kbb)
+            z2d1(:,:) = surf(:,:) * ww(:,:,1) * ts(:,:,1,jp_sal,Kbb)
          END IF
          z_wn_trd_t = - glob_sum( 'diahsb', z2d0 ) 
@@ -135,5 +136,5 @@
 
       !                    ! volume variation (calculated with ssh)
-      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*sshn(:,:) - surf_ini(:,:)*ssh_ini(:,:) )
+      zdiff_v1 = glob_sum_full( 'diahsb', surf(:,:)*ssh(:,:,Kmm) - surf_ini(:,:)*ssh_ini(:,:) )
 
       !                    ! heat & salt content variation (associated with ssh)
@@ -142,11 +143,11 @@
             DO ji = 1, jpi
                DO jj = 1, jpj
-                  z2d0(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj) - ssh_hc_loc_ini(ji,jj) ) 
-                  z2d1(ji,jj) = surf(ji,jj) * ( tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj) - ssh_sc_loc_ini(ji,jj) ) 
+                  z2d0(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm) - ssh_hc_loc_ini(ji,jj) ) 
+                  z2d1(ji,jj) = surf(ji,jj) * ( ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm) - ssh_sc_loc_ini(ji,jj) ) 
                END DO
             END DO
          ELSE                          ! no under ice-shelf seas
-            z2d0(:,:) = surf(:,:) * ( tsn(:,:,1,jp_tem) * sshn(:,:) - ssh_hc_loc_ini(:,:) ) 
-            z2d1(:,:) = surf(:,:) * ( tsn(:,:,1,jp_sal) * sshn(:,:) - ssh_sc_loc_ini(:,:) ) 
+            z2d0(:,:) = surf(:,:) * ( ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) - ssh_hc_loc_ini(:,:) ) 
+            z2d1(:,:) = surf(:,:) * ( ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm) - ssh_sc_loc_ini(:,:) ) 
          END IF
          z_ssh_hc = glob_sum_full( 'diahsb', z2d0 ) 
@@ -155,13 +156,13 @@
       !
       DO jk = 1, jpkm1           ! volume variation (calculated with scale factors)
-         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk) - surf_ini(:,:)*e3t_ini(:,:,jk) ) * tmask(:,:,jk)
+         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm) - surf_ini(:,:)*e3t_ini(:,:,jk) ) * tmask(:,:,jk)
       END DO
       zdiff_v2 = glob_sum_full( 'diahsb', zwrk(:,:,:) )
       DO jk = 1, jpkm1           ! heat content variation
-         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_tem) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
+         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_tem,Kmm) - surf_ini(:,:)*hc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
       END DO
       zdiff_hc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
       DO jk = 1, jpkm1           ! salt content variation
-         zwrk(:,:,jk) = ( surf(:,:)*e3t_n(:,:,jk)*tsn(:,:,jk,jp_sal) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
+         zwrk(:,:,jk) = ( surf(:,:)*e3t(:,:,jk,Kmm)*ts(:,:,jk,jp_sal,Kmm) - surf_ini(:,:)*sc_loc_ini(:,:,jk) ) * tmask(:,:,jk)
       END DO
       zdiff_sc = glob_sum_full( 'diahsb', zwrk(:,:,:) )
@@ -185,5 +186,5 @@
       ! ----------------------- !
       DO jk = 1, jpkm1           ! total ocean volume (calculated with scale factors)
-         zwrk(:,:,jk) = surf(:,:) * e3t_n(:,:,jk) * tmask(:,:,jk)
+         zwrk(:,:,jk) = surf(:,:) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
       zvol_tot = glob_sum_full( 'diahsb', zwrk(:,:,:) )
@@ -191,5 +192,5 @@
 !!gm to be added ?
 !      IF( ln_linssh ) THEN            ! fixed volume, add the ssh contribution
-!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * sshn(:,:) )
+!        zvol_tot = zvol_tot + glob_sum( 'diahsb', surf(:,:) * ssh(:,:,Kmm) )
 !      ENDIF
 !!gm end
@@ -233,5 +234,5 @@
       ENDIF
       !
-      IF( lrst_oce )   CALL dia_hsb_rst( kt, 'WRITE' )
+      IF( lrst_oce )   CALL dia_hsb_rst( kt, Kmm, 'WRITE' )
       !
       IF( ln_timing )   CALL timing_stop('dia_hsb')
@@ -240,5 +241,5 @@
 
 
-   SUBROUTINE dia_hsb_rst( kt, cdrw )
+   SUBROUTINE dia_hsb_rst( kt, Kmm, cdrw )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE dia_hsb_rst  ***
@@ -249,4 +250,5 @@
       !!----------------------------------------------------------------------
       INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kmm    ! ocean time level index
       CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
       !
@@ -281,10 +283,10 @@
             IF(lwp) WRITE(numout,*)
             surf_ini(:,:) = e1e2t(:,:) * tmask_i(:,:)         ! initial ocean surface
-            ssh_ini(:,:) = sshn(:,:)                          ! initial ssh
+            ssh_ini(:,:) = ssh(:,:,Kmm)                          ! initial ssh
             DO jk = 1, jpk
               ! if ice sheet/oceqn coupling, need to mask ini variables here (mask could change at the next NEMO instance).
-               e3t_ini   (:,:,jk) = e3t_n(:,:,jk)                      * tmask(:,:,jk)  ! initial vertical scale factors
-               hc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_tem) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial heat content
-               sc_loc_ini(:,:,jk) = tsn(:,:,jk,jp_sal) * e3t_n(:,:,jk) * tmask(:,:,jk)  ! initial salt content
+               e3t_ini   (:,:,jk) = e3t(:,:,jk,Kmm)                      * tmask(:,:,jk)  ! initial vertical scale factors
+               hc_loc_ini(:,:,jk) = ts(:,:,jk,jp_tem,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial heat content
+               sc_loc_ini(:,:,jk) = ts(:,:,jk,jp_sal,Kmm) * e3t(:,:,jk,Kmm) * tmask(:,:,jk)  ! initial salt content
             END DO
             frc_v = 0._wp                                           ! volume       trend due to forcing
@@ -295,11 +297,11 @@
                   DO ji = 1, jpi
                      DO jj = 1, jpj
-                        ssh_hc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_tem) * sshn(ji,jj)   ! initial heat content in ssh
-                        ssh_sc_loc_ini(ji,jj) = tsn(ji,jj,mikt(ji,jj),jp_sal) * sshn(ji,jj)   ! initial salt content in ssh
+                        ssh_hc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_tem,Kmm) * ssh(ji,jj,Kmm)   ! initial heat content in ssh
+                        ssh_sc_loc_ini(ji,jj) = ts(ji,jj,mikt(ji,jj),jp_sal,Kmm) * ssh(ji,jj,Kmm)   ! initial salt content in ssh
                      END DO
                    END DO
                 ELSE
-                  ssh_hc_loc_ini(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:)   ! initial heat content in ssh
-                  ssh_sc_loc_ini(:,:) = tsn(:,:,1,jp_sal) * sshn(:,:)   ! initial salt content in ssh
+                  ssh_hc_loc_ini(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm)   ! initial heat content in ssh
+                  ssh_sc_loc_ini(:,:) = ts(:,:,1,jp_sal,Kmm) * ssh(:,:,Kmm)   ! initial salt content in ssh
                END IF
                frc_wn_t = 0._wp                                       ! initial heat content misfit due to free surface
@@ -338,5 +340,5 @@
 
 
-   SUBROUTINE dia_hsb_init
+   SUBROUTINE dia_hsb_init( Kmm )
       !!---------------------------------------------------------------------------
       !!                  ***  ROUTINE dia_hsb  ***
@@ -350,4 +352,6 @@
       !!             - Compute coefficients for conversion
       !!---------------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm ! time level index
+      !
       INTEGER ::   ierror, ios   ! local integer
       !!
@@ -417,5 +421,5 @@
       ! 4 - initial conservation variables !
       ! ---------------------------------- !
-      CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
+      CALL dia_hsb_rst( nit000, Kmm, 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE dia_hsb_init
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diahth.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diahth.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diahth.F90	(revision 10965)
@@ -60,5 +60,5 @@
 
 
-   SUBROUTINE dia_hth( kt )
+   SUBROUTINE dia_hth( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_hth  ***
@@ -81,4 +81,5 @@
       !!-------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
       !!
       INTEGER                          ::   ji, jj, jk            ! dummy loop arguments
@@ -139,5 +140,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
+            zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
             hth     (ji,jj) = zztmp
             zabs2   (ji,jj) = zztmp
@@ -150,5 +151,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
+               zztmp = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 
                zrho0_3(ji,jj) = zztmp
                zrho0_1(ji,jj) = zztmp
@@ -162,11 +163,11 @@
          DO ji = 1, jpi
             IF( tmask(ji,jj,nla10) == 1. ) THEN
-               zu  =  1779.50 + 11.250 * tsn(ji,jj,nla10,jp_tem) - 3.80   * tsn(ji,jj,nla10,jp_sal)                             &
-                  &                                              - 0.0745 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)   &
-                  &                                              - 0.0100 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_sal)
-               zv  =  5891.00 + 38.000 * tsn(ji,jj,nla10,jp_tem) + 3.00   * tsn(ji,jj,nla10,jp_sal)                             &
-                  &                                              - 0.3750 * tsn(ji,jj,nla10,jp_tem) * tsn(ji,jj,nla10,jp_tem)
-               zut =    11.25 -  0.149 * tsn(ji,jj,nla10,jp_tem) - 0.01   * tsn(ji,jj,nla10,jp_sal)
-               zvt =    38.00 -  0.750 * tsn(ji,jj,nla10,jp_tem)
+               zu  =  1779.50 + 11.250 * ts(ji,jj,nla10,jp_tem,Kmm) - 3.80   * ts(ji,jj,nla10,jp_sal,Kmm)                             &
+                  &                                              - 0.0745 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)   &
+                  &                                              - 0.0100 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_sal,Kmm)
+               zv  =  5891.00 + 38.000 * ts(ji,jj,nla10,jp_tem,Kmm) + 3.00   * ts(ji,jj,nla10,jp_sal,Kmm)                             &
+                  &                                              - 0.3750 * ts(ji,jj,nla10,jp_tem,Kmm) * ts(ji,jj,nla10,jp_tem,Kmm)
+               zut =    11.25 -  0.149 * ts(ji,jj,nla10,jp_tem,Kmm) - 0.01   * ts(ji,jj,nla10,jp_sal,Kmm)
+               zvt =    38.00 -  0.750 * ts(ji,jj,nla10,jp_tem,Kmm)
                zw  = (zu + 0.698*zv) * (zu + 0.698*zv)
                zdelr(ji,jj) = ztem2 * (1000.*(zut*zv - zvt*zu)/zw)
@@ -187,6 +188,6 @@
             DO ji = 1, jpi
                !
-               zzdep = gdepw_n(ji,jj,jk)
-               zztmp = ( tsn(ji,jj,jk-1,jp_tem) - tsn(ji,jj,jk,jp_tem) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
+               zzdep = gdepw(ji,jj,jk,Kmm)
+               zztmp = ( ts(ji,jj,jk-1,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm) ) / zzdep * tmask(ji,jj,jk)   ! vertical gradient of temperature (dT/dz)
                zzdep = zzdep * tmask(ji,jj,1)
 
@@ -223,7 +224,7 @@
             DO ji = 1, jpi
                !
-               zzdep = gdepw_n(ji,jj,jk) * tmask(ji,jj,1)
+               zzdep = gdepw(ji,jj,jk,Kmm) * tmask(ji,jj,1)
                !
-               zztmp = tsn(ji,jj,nla10,jp_tem) - tsn(ji,jj,jk,jp_tem)  ! - delta T(10m)
+               zztmp = ts(ji,jj,nla10,jp_tem,Kmm) - ts(ji,jj,jk,jp_tem,Kmm)  ! - delta T(10m)
                IF( ABS(zztmp) > ztem2 )      zabs2   (ji,jj) = zzdep   ! abs > 0.2
                IF(     zztmp  > ztem2 )      ztm2    (ji,jj) = zzdep   ! > 0.2
@@ -257,5 +258,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zztmp = tsn(ji,jj,jk,jp_tem)
+               zztmp = ts(ji,jj,jk,jp_tem,Kmm)
                IF( zztmp >= 20. )   ik20(ji,jj) = jk
                IF( zztmp >= 28. )   ik28(ji,jj) = jk
@@ -270,12 +271,12 @@
          DO ji = 1, jpi
             !
-            zzdep = gdepw_n(ji,jj,mbkt(ji,jj)+1)       ! depth of the oean bottom
+            zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)       ! depth of the oean bottom
             !
             iid = ik20(ji,jj)
             IF( iid /= 1 ) THEN 
-               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
-                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
-                  &  * ( 20.*tmask(ji,jj,iid+1) - tsn(ji,jj,iid,jp_tem)                       )   &
-                  &  / ( tsn(ji,jj,iid+1,jp_tem) - tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
+               zztmp =      gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation
+                  &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   &
+                  &  * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm)                       )   &
+                  &  / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) )
                hd20(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1)       ! bound by the ocean depth
             ELSE 
@@ -285,8 +286,8 @@
             iid = ik28(ji,jj)
             IF( iid /= 1 ) THEN 
-               zztmp =      gdept_n(ji,jj,iid  )   &                     ! linear interpolation
-                  &  + (    gdept_n(ji,jj,iid+1) - gdept_n(ji,jj,iid)                       )   &
-                  &  * ( 28.*tmask(ji,jj,iid+1) -    tsn(ji,jj,iid,jp_tem)                       )   &
-                  &  / (  tsn(ji,jj,iid+1,jp_tem) -    tsn(ji,jj,iid,jp_tem) + (1.-tmask(ji,jj,1)) )
+               zztmp =      gdept(ji,jj,iid  ,Kmm)   &                     ! linear interpolation
+                  &  + (    gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm)                       )   &
+                  &  * ( 28.*tmask(ji,jj,iid+1) -    ts(ji,jj,iid,jp_tem,Kmm)                       )   &
+                  &  / (  ts(ji,jj,iid+1,jp_tem,Kmm) -    ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) )
                hd28(ji,jj) = MIN( zztmp , zzdep ) * tmask(ji,jj,1)      ! bound by the ocean depth
             ELSE 
@@ -311,11 +312,11 @@
       END DO
       ! surface boundary condition
-      IF( ln_linssh ) THEN   ;   zthick(:,:) = sshn(:,:)   ;   htc3(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)  
+      IF( ln_linssh ) THEN   ;   zthick(:,:) = ssh(:,:,Kmm)   ;   htc3(:,:) = ts(:,:,1,jp_tem,Kmm) * ssh(:,:,Kmm) * tmask(:,:,1)  
       ELSE                   ;   zthick(:,:) = 0._wp       ;   htc3(:,:) = 0._wp                                   
       ENDIF
       ! integration down to ilevel
       DO jk = 1, ilevel
-         zthick(:,:) = zthick(:,:) + e3t_n(:,:,jk)
-         htc3  (:,:) = htc3  (:,:) + e3t_n(:,:,jk) * tsn(:,:,jk,jp_tem) * tmask(:,:,jk)
+         zthick(:,:) = zthick(:,:) + e3t(:,:,jk,Kmm)
+         htc3  (:,:) = htc3  (:,:) + e3t(:,:,jk,Kmm) * ts(:,:,jk,jp_tem,Kmm) * tmask(:,:,jk)
       END DO
       ! deepest layer
@@ -323,6 +324,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            htc3(ji,jj) = htc3(ji,jj) + tsn(ji,jj,ilevel+1,jp_tem)                  &
-               &                      * MIN( e3t_n(ji,jj,ilevel+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
+            htc3(ji,jj) = htc3(ji,jj) + ts(ji,jj,ilevel+1,jp_tem,Kmm)                  &
+               &                      * MIN( e3t(ji,jj,ilevel+1,Kmm), zthick(ji,jj) ) * tmask(ji,jj,ilevel+1)
          END DO
       END DO
@@ -342,7 +343,8 @@
    LOGICAL , PUBLIC, PARAMETER ::   lk_diahth = .FALSE.  !: thermocline-20d depths flag
 CONTAINS
-   SUBROUTINE dia_hth( kt )         ! Empty routine
+   SUBROUTINE dia_hth( kt, Kmm )         ! Empty routine
       IMPLICIT NONE
       INTEGER, INTENT( in ) :: kt
+      INTEGER, INTENT( in ) :: Kmm
       WRITE(*,*) 'dia_hth: You should not have seen this print! error?', kt
    END SUBROUTINE dia_hth
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaptr.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaptr.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diaptr.F90	(revision 10965)
@@ -71,8 +71,9 @@
 CONTAINS
 
-   SUBROUTINE dia_ptr( pvtr )
+   SUBROUTINE dia_ptr( Kmm, pvtr )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dia_ptr  ***
       !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in)           ::   Kmm    ! time level index
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pvtr   ! j-effective transport
       !
@@ -90,5 +91,5 @@
       REAL(wp), DIMENSION(jpj,jpk,nptr) ::   sjk  , r1_sjk ! i-mean i-k-surface and its inverse
       REAL(wp), DIMENSION(jpj,jpk,nptr) ::   v_msf, sn_jk  , tn_jk ! i-mean T and S, j-Stream-Function
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvn   ! 3D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zvv   ! 3D workspace
 
 
@@ -126,13 +127,13 @@
             zmask(:,:,:) = 0._wp
             zts(:,:,:,:) = 0._wp
-            zvn(:,:,:) = 0._wp
+            zvv(:,:,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpjm1
                   DO ji = 1, jpi
-                     zvfc = e1v(ji,jj) * e3v_n(ji,jj,jk)
+                     zvfc = e1v(ji,jj) * e3v(ji,jj,jk,Kmm)
                      zmask(ji,jj,jk)      = vmask(ji,jj,jk)      * zvfc
-                     zts(ji,jj,jk,jp_tem) = (tsn(ji,jj,jk,jp_tem)+tsn(ji,jj+1,jk,jp_tem)) * 0.5 * zvfc  !Tracers averaged onto V grid
-                     zts(ji,jj,jk,jp_sal) = (tsn(ji,jj,jk,jp_sal)+tsn(ji,jj+1,jk,jp_sal)) * 0.5 * zvfc
-                     zvn(ji,jj,jk)        = vn(ji,jj,jk)         * zvfc
+                     zts(ji,jj,jk,jp_tem) = (ts(ji,jj,jk,jp_tem,Kmm)+ts(ji,jj+1,jk,jp_tem,Kmm)) * 0.5 * zvfc  !Tracers averaged onto V grid
+                     zts(ji,jj,jk,jp_sal) = (ts(ji,jj,jk,jp_sal,Kmm)+ts(ji,jj+1,jk,jp_sal,Kmm)) * 0.5 * zvfc
+                     zvv(ji,jj,jk)        = vv(ji,jj,jk,Kmm)         * zvfc
                   ENDDO
                ENDDO
@@ -147,5 +148,5 @@
              tn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_tem) ) * r1_sjk(:,:,1)
              sn_jk(:,:,1) = ptr_sjk( zts(:,:,:,jp_sal) ) * r1_sjk(:,:,1)
-             v_msf(:,:,1) = ptr_sjk( zvn(:,:,:) )
+             v_msf(:,:,1) = ptr_sjk( zvv(:,:,:) )
 
              htr_ove(:,1) = SUM( v_msf(:,:,1)*tn_jk(:,:,1) ,2 )
@@ -173,5 +174,5 @@
                     tn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_tem), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
                     sn_jk(:,:,jn) = ptr_sjk( zts(:,:,:,jp_sal), btmsk(:,:,jn) ) * r1_sjk(:,:,jn)
-                    v_msf(:,:,jn) = ptr_sjk( zvn(:,:,:), btmsk(:,:,jn) ) 
+                    v_msf(:,:,jn) = ptr_sjk( zvv(:,:,:), btmsk(:,:,jn) ) 
                     htr_ove(:,jn) = SUM( v_msf(:,:,jn)*tn_jk(:,:,jn) ,2 )
                     str_ove(:,jn) = SUM( v_msf(:,:,jn)*sn_jk(:,:,jn) ,2 )
@@ -198,5 +199,5 @@
              WHERE( sjk(:,1,1) /= 0._wp )   r1_sjk(:,1,1) = 1._wp / sjk(:,1,1)
             
-            vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,1))
+            vsum = ptr_sj( zvv(:,:,:), btmsk(:,:,1))
             tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,1) )
             tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,1) )
@@ -220,5 +221,5 @@
                     r1_sjk(:,1,jn) = 0._wp
                     WHERE( sjk(:,1,jn) /= 0._wp )   r1_sjk(:,1,jn) = 1._wp / sjk(:,1,jn)
-                    vsum = ptr_sj( zvn(:,:,:), btmsk(:,:,jn))
+                    vsum = ptr_sj( zvv(:,:,:), btmsk(:,:,jn))
                     tssum(:,jp_tem) = ptr_sj( zts(:,:,:,jp_tem), btmsk(:,:,jn) )
                     tssum(:,jp_sal) = ptr_sj( zts(:,:,:,jp_sal), btmsk(:,:,jn) )
@@ -247,8 +248,8 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     zsfc = e1t(ji,jj) * e3t_n(ji,jj,jk)
+                     zsfc = e1t(ji,jj) * e3t(ji,jj,jk,Kmm)
                      zmask(ji,jj,jk)      = tmask(ji,jj,jk)      * zsfc
-                     zts(ji,jj,jk,jp_tem) = tsn(ji,jj,jk,jp_tem) * zsfc
-                     zts(ji,jj,jk,jp_sal) = tsn(ji,jj,jk,jp_sal) * zsfc
+                     zts(ji,jj,jk,jp_tem) = ts(ji,jj,jk,jp_tem,Kmm) * zsfc
+                     zts(ji,jj,jk,jp_sal) = ts(ji,jj,jk,jp_sal,Kmm) * zsfc
                   END DO
                END DO
@@ -459,5 +460,5 @@
 
 
-   SUBROUTINE dia_ptr_hst( ktra, cptr, pva ) 
+   SUBROUTINE dia_ptr_hst( ktra, cptr, pvflx ) 
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE dia_ptr_hst ***
@@ -468,18 +469,18 @@
       INTEGER                         , INTENT(in )  :: ktra  ! tracer index
       CHARACTER(len=3)                , INTENT(in)   :: cptr  ! transport type  'adv'/'ldf'/'eiv'
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pva   ! 3D input array of advection/diffusion
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)   :: pvflx   ! 3D input array of advection/diffusion
       INTEGER                                        :: jn    !
 
       IF( cptr == 'adv' ) THEN
-         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_tem )  htr_adv(:,1) = ptr_sj( pvflx )
+         IF( ktra == jp_sal )  str_adv(:,1) = ptr_sj( pvflx )
       ENDIF
       IF( cptr == 'ldf' ) THEN
-         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_tem )  htr_ldf(:,1) = ptr_sj( pvflx )
+         IF( ktra == jp_sal )  str_ldf(:,1) = ptr_sj( pvflx )
       ENDIF
       IF( cptr == 'eiv' ) THEN
-         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pva(:,:,:) )
-         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pva(:,:,:) )
+         IF( ktra == jp_tem )  htr_eiv(:,1) = ptr_sj( pvflx )
+         IF( ktra == jp_sal )  str_eiv(:,1) = ptr_sj( pvflx )
       ENDIF
       !
@@ -489,10 +490,10 @@
              IF( ktra == jp_tem ) THEN 
                 DO jn = 2, nptr
-                   htr_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                   htr_adv(:,jn) = ptr_sj( pvflx, btmsk(:,:,jn) )
                 END DO
              ENDIF
              IF( ktra == jp_sal ) THEN 
                 DO jn = 2, nptr
-                   str_adv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                   str_adv(:,jn) = ptr_sj( pvflx, btmsk(:,:,jn) )
                 END DO
              ENDIF
@@ -501,10 +502,10 @@
              IF( ktra == jp_tem ) THEN 
                 DO jn = 2, nptr
-                    htr_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                    htr_ldf(:,jn) = ptr_sj( pvflx, btmsk(:,:,jn) )
                  END DO
              ENDIF
              IF( ktra == jp_sal ) THEN 
                 DO jn = 2, nptr
-                   str_ldf(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                   str_ldf(:,jn) = ptr_sj( pvflx, btmsk(:,:,jn) )
                 END DO
              ENDIF
@@ -513,10 +514,10 @@
              IF( ktra == jp_tem ) THEN 
                 DO jn = 2, nptr
-                    htr_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                    htr_eiv(:,jn) = ptr_sj( pvflx, btmsk(:,:,jn) )
                  END DO
              ENDIF
              IF( ktra == jp_sal ) THEN 
                 DO jn = 2, nptr
-                   str_eiv(:,jn) = ptr_sj( pva(:,:,:), btmsk(:,:,jn) )
+                   str_eiv(:,jn) = ptr_sj( pvflx, btmsk(:,:,jn) )
                 END DO
              ENDIF
@@ -554,5 +555,5 @@
 
 
-   FUNCTION ptr_sj_3d( pva, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sj_3d( pvflx, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sj_3d  ***
@@ -560,10 +561,10 @@
       !! ** Purpose :   i-k sum computation of a j-flux array
       !!
-      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
-      !!              pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
-      !!
-      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!----------------------------------------------------------------------
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pva   ! mask flux array at V-point
+      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
+      !!              pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
+      !!
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk)       ::   pvflx   ! mask flux array at V-point
       REAL(wp), INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
       !
@@ -581,5 +582,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! Vector opt.
-                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
+                  p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * tmask_i(ji,jj) * pmsk(ji,jj)
                END DO
             END DO
@@ -589,5 +590,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! Vector opt.
-                  p_fval(jj) = p_fval(jj) + pva(ji,jj,jk) * tmask_i(ji,jj) 
+                  p_fval(jj) = p_fval(jj) + pvflx(ji,jj,jk) * tmask_i(ji,jj) 
                END DO
             END DO
@@ -601,16 +602,16 @@
 
 
-   FUNCTION ptr_sj_2d( pva, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sj_2d( pvflx, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sj_2d  ***
       !!
-      !! ** Purpose :   "zonal" and vertical sum computation of a i-flux array
-      !!
-      !! ** Method  : - i-k sum of pva using the interior 2D vmask (vmask_i).
-      !!      pva is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
-      !!
-      !! ** Action  : - p_fval: i-k-mean poleward flux of pva
-      !!----------------------------------------------------------------------
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pva   ! mask flux array at V-point
+      !! ** Purpose :   "zonal" and vertical sum computation of a j-flux array
+      !!
+      !! ** Method  : - i-k sum of pvflx using the interior 2D vmask (vmask_i).
+      !!      pvflx is supposed to be a masked flux (i.e. * vmask*e1v*e3v)
+      !!
+      !! ** Action  : - p_fval: i-k-mean poleward flux of pvflx
+      !!----------------------------------------------------------------------
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)           ::   pvflx   ! mask flux array at V-point
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj), OPTIONAL ::   pmsk   ! Optional 2D basin mask
       !
@@ -627,5 +628,5 @@
          DO jj = 2, jpjm1
             DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
-               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
+               p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * tmask_i(ji,jj) * pmsk(ji,jj)
             END DO
          END DO
@@ -633,5 +634,5 @@
          DO jj = 2, jpjm1
             DO ji = nldi, nlei   ! No vector optimisation here. Better use a mask ?
-               p_fval(jj) = p_fval(jj) + pva(ji,jj) * tmask_i(ji,jj)
+               p_fval(jj) = p_fval(jj) + pvflx(ji,jj) * tmask_i(ji,jj)
             END DO
          END DO
@@ -644,5 +645,5 @@
 
 
-   FUNCTION ptr_sjk( pta, pmsk )   RESULT ( p_fval )
+   FUNCTION ptr_sjk( pfld, pmsk )   RESULT ( p_fval )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE ptr_sjk  ***
@@ -650,11 +651,11 @@
       !! ** Purpose :   i-sum computation of an array
       !!
-      !! ** Method  : - i-sum of pva using the interior 2D vmask (vmask_i).
-      !!
-      !! ** Action  : - p_fval: i-mean poleward flux of pva
+      !! ** Method  : - i-sum of field using the interior 2D vmask (pmsk).
+      !!
+      !! ** Action  : - p_fval: i-sum of masked field
       !!----------------------------------------------------------------------
       !!
       IMPLICIT none
-      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pta    ! mask flux array at V-point
+      REAL(wp) , INTENT(in), DIMENSION(jpi,jpj,jpk)           ::   pfld   ! input field to be summed
       REAL(wp) , INTENT(in), DIMENSION(jpi,jpj)    , OPTIONAL ::   pmsk   ! Optional 2D basin mask
       !!
@@ -678,5 +679,5 @@
 !!gm here, use of tmask_i  ==> no need of loop over nldi, nlei....
                DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
-                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * pmsk(ji,jj)
+                  p_fval(jj,jk) = p_fval(jj,jk) + pfld(ji,jj,jk) * pmsk(ji,jj)
                END DO
             END DO
@@ -686,5 +687,5 @@
             DO jj = 2, jpjm1
                DO ji =  nldi, nlei   ! No vector optimisation here. Better use a mask ?
-                  p_fval(jj,jk) = p_fval(jj,jk) + pta(ji,jj,jk) * tmask_i(ji,jj)
+                  p_fval(jj,jk) = p_fval(jj,jk) + pfld(ji,jj,jk) * tmask_i(ji,jj)
                END DO
             END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diatmb.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diatmb.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diatmb.F90	(revision 10965)
@@ -94,5 +94,5 @@
 
 
-   SUBROUTINE dia_tmb
+   SUBROUTINE dia_tmb( Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_tmb  ***
@@ -103,14 +103,16 @@
       !!
       !!--------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kmm     ! time level index
+      !
       REAL(wp) ::   zmdi =1.e+20     ! land value
       REAL(wp), DIMENSION(jpi,jpj,3) :: zwtmb    ! workspace 
       !!--------------------------------------------------------------------
       !
-      CALL dia_calctmb( tsn(:,:,:,jp_tem), zwtmb )
+      CALL dia_calctmb( ts(:,:,:,jp_tem,Kmm), zwtmb )
       !ssh already output but here we output it masked
       IF( ll_wd ) THEN
-         CALL iom_put( "sshnmasked", (sshn(:,:)+ssh_ref)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
+         CALL iom_put( "sshnmasked", (ssh(:,:,Kmm)+ssh_ref)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
       ELSE
-         CALL iom_put( "sshnmasked", sshn(:,:)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
+         CALL iom_put( "sshnmasked", ssh(:,:,Kmm)*tmask(:,:,1) + zmdi*(1.0 - tmask(:,:,1)) )
       ENDIF
 
@@ -119,15 +121,15 @@
       CALL iom_put( "bot_temp"  , zwtmb(:,:,3) )    ! tmb Temperature
       !
-      CALL dia_calctmb( tsn(:,:,:,jp_sal), zwtmb )
+      CALL dia_calctmb( ts(:,:,:,jp_sal,Kmm), zwtmb )
       CALL iom_put( "top_sal"   , zwtmb(:,:,1) )    ! tmb Salinity 
       CALL iom_put( "mid_sal"   , zwtmb(:,:,2) )    ! tmb Salinity
       CALL iom_put( "bot_sal"   , zwtmb(:,:,3) )    ! tmb Salinity
       !
-      CALL dia_calctmb( un(:,:,:), zwtmb )
+      CALL dia_calctmb( uu(:,:,:,Kmm), zwtmb )
       CALL iom_put( "top_u"     , zwtmb(:,:,1) )    ! tmb  U Velocity
       CALL iom_put( "mid_u"     , zwtmb(:,:,2) )    ! tmb  U Velocity
       CALL iom_put( "bot_u"     , zwtmb(:,:,3) )    ! tmb  U Velocity
       !
-      CALL dia_calctmb( vn(:,:,:), zwtmb )
+      CALL dia_calctmb( vv(:,:,:,Kmm), zwtmb )
       CALL iom_put( "top_v"     , zwtmb(:,:,1) )    ! tmb  V Velocity
       CALL iom_put( "mid_v"     , zwtmb(:,:,2) )    ! tmb  V Velocity
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DIA/diawri.F90	(revision 10965)
@@ -97,5 +97,5 @@
 
    
-   SUBROUTINE dia_wri( kt )
+   SUBROUTINE dia_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_wri  ***
@@ -107,4 +107,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
       !!
       INTEGER ::   ji, jj, jk       ! dummy loop indices
@@ -120,5 +121,5 @@
       ! Output the initial state and forcings
       IF( ninist == 1 ) THEN                       
-         CALL dia_wri_state( 'output.init' )
+         CALL dia_wri_state( Kmm, 'output.init' )
          ninist = 0
       ENDIF
@@ -129,27 +130,27 @@
       CALL iom_put("e3v_0", e3v_0(:,:,:) )
       !
-      CALL iom_put( "e3t" , e3t_n(:,:,:) )
-      CALL iom_put( "e3u" , e3u_n(:,:,:) )
-      CALL iom_put( "e3v" , e3v_n(:,:,:) )
-      CALL iom_put( "e3w" , e3w_n(:,:,:) )
+      CALL iom_put( "e3t" , e3t(:,:,:,Kmm) )
+      CALL iom_put( "e3u" , e3u(:,:,:,Kmm) )
+      CALL iom_put( "e3v" , e3v(:,:,:,Kmm) )
+      CALL iom_put( "e3w" , e3w(:,:,:,Kmm) )
       IF( iom_use("e3tdef") )   &
-         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
+         CALL iom_put( "e3tdef"  , ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
 
       IF( ll_wd ) THEN
          CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
       ELSE
-         CALL iom_put( "ssh" , sshn )              ! sea surface height
+         CALL iom_put( "ssh" , ssh(:,:,Kmm) )              ! sea surface height
       ENDIF
 
       IF( iom_use("wetdep") )   &                  ! wet depth
-         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
+         CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )
       
-      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
-      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
+      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature
+      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature
       IF ( iom_use("sbt") ) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
                ikbot = mbkt(ji,jj)
-               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
+               z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
             END DO
          END DO
@@ -157,11 +158,11 @@
       ENDIF
       
-      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
-      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
+      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity
+      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity
       IF ( iom_use("sbs") ) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
                ikbot = mbkt(ji,jj)
-               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
+               z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
             END DO
          END DO
@@ -174,8 +175,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
-                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
-                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
-                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
+               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   &
+                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   &
+                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   &
+                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2
                z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
                !
@@ -186,11 +187,11 @@
       ENDIF
          
-      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
-      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
+      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current
+      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current
       IF ( iom_use("sbu") ) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
                ikbot = mbku(ji,jj)
-               z2d(ji,jj) = un(ji,jj,ikbot)
+               z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
             END DO
          END DO
@@ -198,11 +199,11 @@
       ENDIF
       
-      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
-      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
+      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current
+      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current
       IF ( iom_use("sbv") ) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
                ikbot = mbkv(ji,jj)
-               z2d(ji,jj) = vn(ji,jj,ikbot)
+               z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
             END DO
          END DO
@@ -210,10 +211,10 @@
       ENDIF
 
-      CALL iom_put( "woce", wn )                   ! vertical velocity
+      CALL iom_put( "woce", ww )                   ! vertical velocity
       IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
          ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
          z2d(:,:) = rau0 * e1e2t(:,:)
          DO jk = 1, jpk
-            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
+            z3d(:,:,jk) = ww(:,:,jk) * z2d(:,:)
          END DO
          CALL iom_put( "w_masstr" , z3d )  
@@ -231,7 +232,7 @@
          DO jj = 2, jpjm1                                    ! sst gradient
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zztmp  = tsn(ji,jj,1,jp_tem)
-               zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
-               zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
+               zztmp  = ts(ji,jj,1,jp_tem,Kmm)
+               zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
+               zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
                z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
                   &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
@@ -250,5 +251,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -262,5 +263,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -274,9 +275,9 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
-                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
-                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
-                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
-                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
+                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+                  z3d(ji,jj,jk) = zztmp * (  uu(ji-1,jj,jk,Kmm)**2 * e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)   &
+                     &                     + uu(ji  ,jj,jk,Kmm)**2 * e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)   &
+                     &                     + vv(ji,jj-1,jk,Kmm)**2 * e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)   &
+                     &                     + vv(ji,jj  ,jk,Kmm)**2 * e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)   )
                END DO
             END DO
@@ -286,5 +287,5 @@
       ENDIF
       !
-      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
+      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence
       !
       IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
@@ -292,5 +293,5 @@
          z2d(:,:) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
+            z3d(:,:,jk) = rau0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
             z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
          END DO
@@ -304,5 +305,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
                END DO
             END DO
@@ -317,5 +318,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
                END DO
             END DO
@@ -329,5 +330,5 @@
          z3d(:,:,jpk) = 0.e0
          DO jk = 1, jpkm1
-            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
+            z3d(:,:,jk) = rau0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
          END DO
          CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
@@ -339,5 +340,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
                END DO
             END DO
@@ -352,5 +353,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
+                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
                END DO
             END DO
@@ -365,5 +366,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) *  ts(ji,jj,jk,jp_tem,Kmm)
                END DO
             END DO
@@ -377,5 +378,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
                END DO
             END DO
@@ -388,7 +389,7 @@
       !
 
-      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values 
+      IF (ln_diatmb)   CALL dia_tmb( Kmm )            ! tmb values 
           
-      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
+      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )        ! 25h averaging
 
       IF( ln_timing )   CALL timing_stop('dia_wri')
@@ -416,5 +417,5 @@
 
    
-   SUBROUTINE dia_wri( kt )
+   SUBROUTINE dia_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_wri  ***
@@ -429,4 +430,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
       !
       LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
@@ -446,5 +448,5 @@
       !
       IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
-         CALL dia_wri_state( 'output.init' )
+         CALL dia_wri_state( Kmm, 'output.init' )
          ninist = 0
       ENDIF
@@ -583,9 +585,9 @@
             &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
          IF(  .NOT.ln_linssh  ) THEN
-            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
+            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
+            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
-            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
+            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t(:,:,:,Kmm)
             &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
          ENDIF
@@ -604,8 +606,8 @@
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          IF(  ln_linssh  ) THEN
-            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
+            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm)
             &                                                                  , "KgC/m2/s",  &  ! sosst_cd
             &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
+            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm)
             &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
             &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -692,5 +694,5 @@
 
          !                                                                                      !!! nid_U : 3D
-         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
+         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm)
             &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
          IF( ln_wave .AND. ln_sdw) THEN
@@ -705,5 +707,5 @@
 
          !                                                                                      !!! nid_V : 3D
-         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
+         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm)
             &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
          IF( ln_wave .AND. ln_sdw) THEN
@@ -718,5 +720,5 @@
 
          !                                                                                      !!! nid_W : 3D
-         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
+         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww
             &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
          CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
@@ -756,21 +758,21 @@
 
       IF( .NOT.ln_linssh ) THEN
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
+         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content
+         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content
+         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content
+         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content
       ELSE
-         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
-         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
-         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
-         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
+         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
+         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
+         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
+         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
       ENDIF
       IF( .NOT.ln_linssh ) THEN
-         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
-         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
-         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
+         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
+         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness
+         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
          CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
       ENDIF
-      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
+      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
       CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
       CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
@@ -779,7 +781,7 @@
                                                                                   ! in linear free surface case)
       IF( ln_linssh ) THEN
-         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
+         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
          CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
-         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
+         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
          CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
       ENDIF
@@ -817,5 +819,5 @@
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
@@ -823,5 +825,5 @@
          CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
          CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
-         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
+         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
          CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
       ENDIF
@@ -836,11 +838,11 @@
 #endif
 
-      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
+      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
       CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
 
-      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
+      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
       CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
 
-      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
+      CALL histwrite( nid_W, "vovecrtz", it, ww             , ndim_T, ndex_T )    ! vert. current
       CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
       CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
@@ -869,5 +871,5 @@
 #endif
 
-   SUBROUTINE dia_wri_state( cdfile_name )
+   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE dia_wri_state  ***
@@ -882,4 +884,5 @@
       !!      File 'output.abort.nc' is created in case of abnormal job end
       !!----------------------------------------------------------------------
+      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
       CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
       !!
@@ -898,10 +901,10 @@
 #endif
 
-      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
-      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
-      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
-      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
-      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
-      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
+      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
+      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
+      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
+      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
+      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
+      CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww                )    ! now k-velocity
       IF( ALLOCATED(ahtu) ) THEN
          CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
@@ -919,6 +922,6 @@
       CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
       IF(  .NOT.ln_linssh  ) THEN             
-         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth 
-         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness  
+         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth 
+         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness  
       END IF
       IF( ln_wave .AND. ln_sdw ) THEN
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv.F90	(revision 10965)
@@ -136,5 +136,5 @@
       !
 !!gm ???
-      IF( ln_diaptr )   CALL dia_ptr( zvv )                                    ! diagnose the effective MSF 
+      IF( ln_diaptr )   CALL dia_ptr( Kmm, zvv )                               ! diagnose the effective MSF 
 !!gm ???
       !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10965)
@@ -489,9 +489,9 @@
                            CALL dia_ptr_init    ! Poleward TRansports initialization
       IF( lk_diadct    )   CALL dia_dct_init    ! Sections tranports
-                           CALL dia_hsb_init    ! heat content, salt content and volume budgets
+                           CALL dia_hsb_init( Nnn )    ! heat content, salt content and volume budgets
                            CALL     trd_init( Nnn )    ! Mixed-layer/Vorticity/Integral constraints trends
                            CALL dia_obs_init( Nnn )    ! Initialize observational data
                            CALL dia_tmb_init    ! TMB outputs
-                           CALL dia_25h_init    ! 25h mean  outputs
+                           CALL dia_25h_init( Nbb )    ! 25h mean  outputs
       IF( ln_diaobs    )   CALL dia_obs( nit000-1, Nnn )   ! Observation operator for restart
 
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 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10965)
@@ -219,10 +219,10 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_floats  )   CALL flo_stp ( kstp )        ! drifting Floats
-      IF( ln_diacfl  )   CALL dia_cfl ( kstp )        ! Courant number diagnostics
-      IF( lk_diahth  )   CALL dia_hth ( kstp )        ! Thermocline depth (20 degres isotherm depth)
-      IF( lk_diadct  )   CALL dia_dct ( kstp )        ! Transports
-                         CALL dia_ar5 ( kstp )        ! ar5 diag
-      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
-                         CALL dia_wri ( kstp )        ! ocean model: outputs
+      IF( ln_diacfl  )   CALL dia_cfl ( kstp, Nnn )   ! Courant number diagnostics
+      IF( lk_diahth  )   CALL dia_hth ( kstp, Nnn )   ! Thermocline depth (20 degres isotherm depth)
+      IF( lk_diadct  )   CALL dia_dct ( kstp, Nnn )   ! Transports
+                         CALL dia_ar5 ( kstp, Nnn )   ! ar5 diag
+      IF( lk_diaharm )   CALL dia_harm( kstp, Nnn )   ! Tidal harmonic analysis
+                         CALL dia_wri ( kstp, Nnn )   ! ocean model: outputs
       !
       IF( ln_crs     )   CALL crs_fld       ( kstp )  ! ocean model: online field coarsening & output
@@ -259,5 +259,5 @@
 
 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?)
-      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics
+      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
@@ -286,5 +286,5 @@
       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
       !
-      IF( ln_diahsb  )   CALL dia_hsb       ( kstp )  ! - ML - global conservation diagnostics
+      IF( ln_diahsb  )   CALL dia_hsb       ( kstp, Nbb, Nnn )  ! - ML - global conservation diagnostics
 
 !!gm : This does not only concern the dynamics ==>>> add a new title
@@ -309,5 +309,5 @@
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL stp_ctl      ( kstp, indic )
+                         CALL stp_ctl      ( kstp, Nnn, indic )
                          
       IF( kstp == nit000 ) THEN                          ! 1st time step only
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/stpctl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/stpctl.F90	(revision 10964)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/stpctl.F90	(revision 10965)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE stp_ctl( kt, kindic )
+   SUBROUTINE stp_ctl( kt, Kmm, kindic )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE stp_ctl  ***
@@ -60,4 +60,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kmm      ! ocean time level index
       INTEGER, INTENT(inout) ::   kindic   ! error indicator
       !!
@@ -111,13 +112,13 @@
       !                                   !==  test of extrema  ==!
       IF( ll_wd ) THEN
-         zmax(1) = MAXVAL(  ABS( sshn(:,:) + ssh_ref*tmask(:,:,1) )  )        ! ssh max 
+         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) + ssh_ref*tmask(:,:,1) )  )        ! ssh max 
       ELSE
-         zmax(1) = MAXVAL(  ABS( sshn(:,:) )  )                               ! ssh max
+         zmax(1) = MAXVAL(  ABS( ssh(:,:,Kmm) )  )                               ! ssh max
       ENDIF
-      zmax(2) = MAXVAL(  ABS( un(:,:,:) )  )                                  ! velocity max (zonal only)
-      zmax(3) = MAXVAL( -tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
-      zmax(4) = MAXVAL(  tsn(:,:,:,jp_sal) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
-      zmax(5) = MAXVAL( -tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max
-      zmax(6) = MAXVAL(  tsn(:,:,:,jp_tem) , mask = tmask(:,:,:) == 1._wp )   !       temperature max
+      zmax(2) = MAXVAL(  ABS( uu(:,:,:,Kmm) )  )                                  ! velocity max (zonal only)
+      zmax(3) = MAXVAL( -ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus salinity max
+      zmax(4) = MAXVAL(  ts(:,:,:,jp_sal,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       salinity max
+      zmax(5) = MAXVAL( -ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   ! minus temperature max
+      zmax(6) = MAXVAL(  ts(:,:,:,jp_tem,Kmm) , mask = tmask(:,:,:) == 1._wp )   !       temperature max
       zmax(7) = REAL( nstop , wp )                                            ! stop indicator
       IF( ln_zad_Aimp ) THEN
@@ -155,13 +156,13 @@
          &  ISNAN( zmax(1) + zmax(2) + zmax(3) ) ) ) THEN   ! NaN encounter in the tests
          IF( lk_mpp .AND. ln_ctl ) THEN
-            CALL mpp_maxloc( 'stpctl', ABS(sshn)        , ssmask(:,:)  , zzz, ih  )
-            CALL mpp_maxloc( 'stpctl', ABS(un)          , umask (:,:,:), zzz, iu  )
-            CALL mpp_minloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 )
-            CALL mpp_maxloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 )
+            CALL mpp_maxloc( 'stpctl', ABS(ssh(:,:,Kmm))        , ssmask(:,:)  , zzz, ih  )
+            CALL mpp_maxloc( 'stpctl', ABS(uu(:,:,:,Kmm))          , umask (:,:,:), zzz, iu  )
+            CALL mpp_minloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is1 )
+            CALL mpp_maxloc( 'stpctl', ts(:,:,:,jp_sal,Kmm), tmask (:,:,:), zzz, is2 )
          ELSE
-            ih(:)  = MAXLOC( ABS( sshn(:,:)   )                              ) + (/ nimpp - 1, njmpp - 1    /)
-            iu(:)  = MAXLOC( ABS( un  (:,:,:) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /)
-            is1(:) = MINLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
-            is2(:) = MAXLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            ih(:)  = MAXLOC( ABS( ssh(:,:,Kmm)   )                              ) + (/ nimpp - 1, njmpp - 1    /)
+            iu(:)  = MAXLOC( ABS( uu  (:,:,:,Kmm) )                              ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            is1(:) = MINLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
+            is2(:) = MAXLOC( ts(:,:,:,jp_sal,Kmm), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /)
          ENDIF
          
@@ -173,5 +174,5 @@
          WRITE(ctmp6,*) '      ===> output of last computed fields in output.abort.nc file'
          
-         CALL dia_wri_state( 'output.abort' )     ! create an output.abort file
+         CALL dia_wri_state( Kmm, 'output.abort' )     ! create an output.abort file
          
          IF( .NOT. ln_ctl ) THEN
