Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/diawri.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/diawri.F90	(revision 11770)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/diawri.F90	(revision 11771)
@@ -57,6 +57,6 @@
    USE lib_mpp         ! MPP library
    USE timing          ! preformance summary
-   USE diurnal_bulk    ! diurnal warm layer
-   USE cool_skin       ! Cool skin
+   USE diu_bulk        ! diurnal warm layer
+   USE diu_coolskin    ! Cool skin
 
    IMPLICIT NONE
@@ -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
@@ -115,6 +116,4 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   bu, bv   ! volume of u- and v-boxes
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   r1_bt    ! inverse of t-box volume
       !!----------------------------------------------------------------------
       ! 
@@ -123,5 +122,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
@@ -132,27 +131,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)
+         CALL iom_put( "ssh" , (ssh(:,:,Kmm)+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
@@ -160,11 +159,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
@@ -177,8 +176,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) 
                !
@@ -189,11 +188,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
@@ -201,11 +200,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
@@ -213,10 +212,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 )  
@@ -236,7 +235,7 @@
             DO jj = 2, jpjm1                                    ! sal gradient
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zztmp  = tsn(ji,jj,jk,jp_sal)
-                  zztmpx = ( tsn(ji+1,jj,jk,jp_sal) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,jk,jp_sal) ) * r1_e1u(ji-1,jj)
-                  zztmpy = ( tsn(ji,jj+1,jk,jp_sal) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,jk,jp_sal) ) * r1_e2v(ji,jj-1)
+                  zztmp  = ts(ji,jj,jk,jp_sal,Kmm)
+                  zztmpx = ( ts(ji+1,jj,jk,jp_sal,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,jk,jp_sal,Kmm) ) * r1_e1u(ji-1,jj)
+                  zztmpy = ( ts(ji,jj+1,jk,jp_sal,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,jk,jp_sal,Kmm) ) * r1_e2v(ji,jj-1)
                   z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
                      &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)
@@ -253,7 +252,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)
@@ -272,5 +271,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
@@ -284,5 +283,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
@@ -296,5 +295,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) * 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) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -306,9 +305,11 @@
          z3d(:,:,jpk) = 0._wp 
          DO jk = 1, jpkm1
-            DO jj = 2, jpj
-               DO ji = 2, jpi
-                  zztmpx = 0.5 * ( un(ji-1,jj  ,jk) + un(ji,jj,jk) )
-                  zztmpy = 0.5 * ( vn(ji  ,jj-1,jk) + vn(ji,jj,jk) )
-                  z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy )
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  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
@@ -326,9 +327,9 @@
             DO jj = 2, jpj
                DO ji = 2, jpi
-                  z3d(ji,jj,jk) = 0.25_wp * ( un(ji  ,jj,jk) * un(ji  ,jj,jk) * e1e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)  &
-                     &                      + un(ji-1,jj,jk) * un(ji-1,jj,jk) * e1e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)  &
-                     &                      + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1e2v(ji,jj  ) * e3v_n(ji,jj  ,jk)  &
-                     &                      + vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1e2v(ji,jj-1) * e3v_n(ji,jj-1,jk)  )  &
-                     &                    * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
+                  z3d(ji,jj,jk) = 0.25_wp * ( uu(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)  &
+                     &                      + uu(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)  &
+                     &                      + vv(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)  &
+                     &                      + vv(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)  )  &
+                     &                    * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -342,5 +343,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * z3d(ji,jj,jk) * tmask(ji,jj,jk)
+                  z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -350,5 +351,5 @@
       ENDIF
       !
-      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
+      CALL iom_put( "hdiv", hdiv )                  ! Horizontal divergence
 
       IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN
@@ -358,6 +359,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  z3d(ji,jj,jk) = (  e2v(ji+1,jj  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
-                     &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * r1_e1e2f(ji,jj)
+                  z3d(ji,jj,jk) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm)    &
+                     &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm)  ) * r1_e1e2f(ji,jj)
                END DO
             END DO
@@ -378,6 +379,6 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  ze3  = (  e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk)  )
+                  ze3  = (  e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)   &
+                     &    + e3t(ji,jj  ,jk,Kmm)*tmask(ji,jj  ,jk) + e3t(ji+1,jj  ,jk,Kmm)*tmask(ji+1,jj  ,jk)  )
                   IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
                   ELSE                      ;   ze3 = 0._wp
@@ -397,5 +398,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
@@ -409,5 +410,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
@@ -422,5 +423,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
@@ -434,5 +435,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
@@ -444,5 +445,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
@@ -457,5 +458,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
@@ -470,5 +471,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
@@ -482,5 +483,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
@@ -493,7 +494,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')
@@ -521,5 +522,5 @@
 
    
-   SUBROUTINE dia_wri( kt )
+   SUBROUTINE dia_wri( kt, Kmm )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dia_wri  ***
@@ -534,4 +535,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
@@ -551,5 +553,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
@@ -688,9 +690,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
@@ -709,8 +711,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 )
@@ -797,5 +799,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
@@ -810,5 +812,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
@@ -823,5 +825,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
@@ -861,21 +863,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
@@ -884,7 +886,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
@@ -922,5 +924,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
@@ -928,5 +930,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
@@ -941,11 +943,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.
@@ -974,5 +976,5 @@
 #endif
 
-   SUBROUTINE dia_wri_state( cdfile_name )
+   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE dia_wri_state  ***
@@ -987,4 +989,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
       !!
@@ -1003,10 +1006,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
@@ -1024,6 +1027,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/tests/CANAL/MY_SRC/domvvl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/domvvl.F90	(revision 11770)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/domvvl.F90	(revision 11771)
@@ -8,4 +8,5 @@
    !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
    !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
+   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
    !!----------------------------------------------------------------------
 
@@ -13,5 +14,5 @@
    !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
    !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
-   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
+   !!   dom_vvl_sf_update   : Swap vertical scale factors and update the vertical grid
    !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
    !!   dom_vvl_rst      : read/write restart file
@@ -37,5 +38,5 @@
    PUBLIC  dom_vvl_init       ! called by domain.F90
    PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
-   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
+   PUBLIC  dom_vvl_sf_update  ! called by step.F90
    PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
 
@@ -93,5 +94,5 @@
 
 
-   SUBROUTINE dom_vvl_init
+   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_init  ***
@@ -104,10 +105,10 @@
       !!
       !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
-      !!              - Regrid: e3(u/v)_n
-      !!                        e3(u/v)_b       
-      !!                        e3w_n           
-      !!                        e3(u/v)w_b      
-      !!                        e3(u/v)w_n      
-      !!                        gdept_n, gdepw_n and gde3w_n
+      !!              - Regrid: e3[u/v](:,:,:,Kmm)
+      !!                        e3[u/v](:,:,:,Kmm)       
+      !!                        e3w(:,:,:,Kmm)           
+      !!                        e3[u/v]w_b
+      !!                        e3[u/v]w_n      
+      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
       !!              - h(t/u/v)_0
       !!              - frq_rst_e3t and frq_rst_hdv
@@ -115,4 +116,6 @@
       !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
+      !
       INTEGER ::   ji, jj, jk
       INTEGER ::   ii0, ii1, ij0, ij1
@@ -130,33 +133,33 @@
       !
       !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
-      CALL dom_vvl_rst( nit000, 'READ' )
-      e3t_a(:,:,jpk) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
+      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )
+      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
       !
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
       !                                ! Horizontal interpolation of e3t
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )    ! from T to U
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )    ! from T to V 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )    ! from U to F
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V 
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
       !                                ! Vertical interpolation of e3t,u,v 
-      CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )  ! from T to W
-      CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )  ! from U to UW
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )  ! from V to UW
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
+      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
+      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
       ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
-      e3t_a(:,:,:) = e3t_n(:,:,:)
-      e3u_a(:,:,:) = e3u_n(:,:,:)
-      e3v_a(:,:,:) = e3v_n(:,:,:)
+      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
+      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
+      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
       !
       !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)       ! reference to the ocean surface (used for MLD and light penetration)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)  ! reference to a common level z=0 for hpg
-      gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1)
-      gdepw_b(:,:,1) = 0.0_wp
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
+      gdepw(:,:,1,Kmm) = 0.0_wp
+      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
+      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
+      gdepw(:,:,1,Kbb) = 0.0_wp
       DO jk = 2, jpk                               ! vertical sum
          DO jj = 1,jpj
@@ -165,13 +168,13 @@
                !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
                !                             ! 0.5 where jk = mikt     
-!!gm ???????   BUG ?  gdept_n as well as gde3w_n  does not include the thickness of ISF ??
+!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
                zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk)) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
-               gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1)
-               gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  &
-                  &                + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk)) 
+               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
+                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
+               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
+               gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
+               gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
+                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
             END DO
          END DO
@@ -179,22 +182,22 @@
       !
       !                    !==  thickness of the water column  !!   (ocean portion only)
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
-      hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
-      hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1)
-      hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
-      hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1)
+      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
+      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
+      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
+      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
+      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
-         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk)
-         hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
-         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk)
+         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
+         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
+         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
+         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
       END DO
       !
       !                    !==  inverse of water column thickness   ==!   (u- and v- points)
-      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
-      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
-      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
+      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
+      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
 
       !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
@@ -266,5 +269,5 @@
 
 
-   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
+   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE dom_vvl_sf_nxt  ***
@@ -288,6 +291,7 @@
       !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )           ::   kt      ! time step
-      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
+      INTEGER, INTENT( in )           ::   kt             ! time step
+      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
+      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
       !
       INTEGER                ::   ji, jj, jk            ! dummy loop indices
@@ -321,8 +325,8 @@
       !                                                ! --------------------------------------------- !
       !
-      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
       DO jk = 1, jpkm1
-         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
-         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
+         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
       END DO
       !
@@ -337,6 +341,6 @@
          zht(:,:)   = 0._wp
          DO jk = 1, jpkm1
-            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
-            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
+            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
          END DO
          zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
@@ -348,5 +352,5 @@
                DO jk = 1, jpkm1
                   hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
-                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
+                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
                END DO
             ENDIF
@@ -361,9 +365,9 @@
          IF( ln_vvl_ztilde ) THEN     ! z_tilde case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
             END DO
          ELSE                         ! layer case
             DO jk = 1, jpkm1
-               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
+               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
             END DO
          ENDIF
@@ -476,7 +480,7 @@
             zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
-         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
+         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
          DO jk = 1, jpkm1
-            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
+            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
          END DO
 
@@ -486,5 +490,5 @@
       !                                           ! ---baroclinic part--------- !
          DO jk = 1, jpkm1
-            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
+            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
          END DO
       ENDIF
@@ -501,37 +505,37 @@
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
          !
          zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
-            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
-         END DO
-         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
-         !
-         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
          CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
       END IF
 
@@ -540,6 +544,6 @@
       ! *********************************** !
 
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
-      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
 
       ! *********************************** !
@@ -547,14 +551,14 @@
       ! *********************************** !
 
-      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
-      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
+      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
+      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
       DO jk = 2, jpkm1
-         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
-         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
+         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
+         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
       END DO
       !                                        ! Inverse of the local depth
 !!gm BUG ?  don't understand the use of umask_i here .....
-      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
-      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
+      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
+      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
       !
       IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
@@ -563,23 +567,23 @@
 
 
-   SUBROUTINE dom_vvl_sf_swp( kt )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE dom_vvl_sf_swp  ***
+   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE dom_vvl_sf_update  ***
       !!                   
-      !! ** Purpose :  compute time filter and swap of scale factors 
+      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors 
       !!               compute all depths and related variables for next time step
       !!               write outputs and restart file
       !!
-      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
+      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
       !!               - reconstruct scale factor at other grid points (interpolate)
       !!               - recompute depths and water height fields
       !!
-      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
+      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
       !!               - Recompute:
       !!                    e3(u/v)_b       
-      !!                    e3w_n           
+      !!                    e3w(:,:,:,Kmm)           
       !!                    e3(u/v)w_b      
       !!                    e3(u/v)w_n      
-      !!                    gdept_n, gdepw_n  and gde3w_n
+      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
       !!                    h(u/v) and h(u/v)r
       !!
@@ -587,5 +591,6 @@
       !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! time step
+      INTEGER, INTENT( in ) ::   kt              ! time step
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -595,10 +600,10 @@
       IF( ln_linssh )   RETURN      ! No calculation in linear free surface
       !
-      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
+      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
       !
       IF( kt == nit000 )   THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
+         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
       ENDIF
       !
@@ -615,10 +620,10 @@
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
       ENDIF
-      gdept_b(:,:,:) = gdept_n(:,:,:)
-      gdepw_b(:,:,:) = gdepw_n(:,:,:)
-
-      e3t_n(:,:,:) = e3t_a(:,:,:)
-      e3u_n(:,:,:) = e3u_a(:,:,:)
-      e3v_n(:,:,:) = e3v_a(:,:,:)
+      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
+      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
+
+      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
+      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)
+      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)
 
       ! Compute all missing vertical scale factor and depths
@@ -626,21 +631,21 @@
       ! Horizontal scale factor interpolations
       ! --------------------------------------
-      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
-      ! - JC - hu_b, hv_b, hur_b, hvr_b also
+      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
+      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
       
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
       
       ! Vertical scale factor interpolations
-      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
-      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
+      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
+      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
+      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
 
       ! t- and w- points depth (set the isf depth as it is in the initial step)
-      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-      gdepw_n(:,:,1) = 0.0_wp
-      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
+      gdepw(:,:,1,Kmm) = 0.0_wp
+      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
       DO jk = 2, jpk
          DO jj = 1,jpj
@@ -649,8 +654,8 @@
                                                                  ! 1 for jk = mikt
                zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
-               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
-               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
-                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
-               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
+               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
+               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
+                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
+               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
             END DO
          END DO
@@ -659,19 +664,17 @@
       ! Local depth and Inverse of the local depth of the water
       ! -------------------------------------------------------
-      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
-      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
-      !
-      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
+      !
+      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
       DO jk = 2, jpkm1
-         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
+         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
       END DO
 
       ! write restart file
       ! ==================
-      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
-      !
-      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
-      !
-   END SUBROUTINE dom_vvl_sf_swp
+      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
+      !
+      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
+      !
+   END SUBROUTINE dom_vvl_sf_update
 
 
@@ -783,5 +786,5 @@
 
 
-   SUBROUTINE dom_vvl_rst( kt, cdrw )
+   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE dom_vvl_rst  ***
@@ -795,6 +798,7 @@
       !!                they are set to 0.
       !!----------------------------------------------------------------------
-      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
-      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
+      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
+      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
       !
       INTEGER ::   ji, jj, jk
@@ -806,5 +810,5 @@
          IF( ln_rstart ) THEN                   !* Read the restart file
             CALL rst_read_open                  !  open the restart file if necessary
-            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
+            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
             !
             id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
@@ -817,39 +821,39 @@
             !                             ! --------- !
             IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
                ! needed to restart if land processor not computed 
-               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
                WHERE ( tmask(:,:,:) == 0.0_wp ) 
-                  e3t_n(:,:,:) = e3t_0(:,:,:)
-                  e3t_b(:,:,:) = e3t_0(:,:,:)
+                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
+                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
                END WHERE
                IF( neuler == 0 ) THEN
-                  e3t_b(:,:,:) = e3t_n(:,:,:)
+                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                ENDIF
             ELSE IF( id1 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
-               e3t_n(:,:,:) = e3t_b(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
+               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
                neuler = 0
             ELSE IF( id2 > 0 ) THEN
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
                IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
-               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                neuler = 0
             ELSE
-               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
+               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
                IF(lwp) write(numout,*) 'Compute scale factor from sshn'
                IF(lwp) write(numout,*) 'neuler is forced to 0'
                DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
                       &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
                       &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
                END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
                neuler = 0
             ENDIF
@@ -888,21 +892,19 @@
                IF( cn_cfg == 'wad' ) THEN
                   ! Wetting and drying test case
-                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
-                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-                  sshn (:,:)     = sshb(:,:)
-                  un   (:,:,:)   = ub  (:,:,:)
-                  vn   (:,:,:)   = vb  (:,:,:)
+                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
+                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
+                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
                ELSE
                   ! if not test case
-                  sshn(:,:) = -ssh_ref
-                  sshb(:,:) = -ssh_ref
+                  ssh(:,:,Kmm) = -ssh_ref
+                  ssh(:,:,Kbb) = -ssh_ref
 
                   DO jj = 1, jpj
                      DO ji = 1, jpi
                         IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
-
-                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
-                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
+                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
                         ENDIF
                      ENDDO
@@ -912,9 +914,9 @@
                ! Adjust vertical metrics for all wad
                DO jk = 1, jpk
-                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
                     &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
                     &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
                END DO
-               e3t_b(:,:,:) = e3t_n(:,:,:)
+               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
 
                DO ji = 1, jpi
@@ -928,17 +930,15 @@
             ELSE
                !
-               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t_b and e3t_n
-               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )  
+               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm)
+               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
                ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn)
                !
                DO jk=1,jpk
-                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
-                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
-                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points
+                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
+                    &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
+                    &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points
                END DO
-               e3t_n(:,:,:) = e3t_b(:,:,:)
-               sshn(:,:) = sshb(:,:)   ! needed later for gde3w
-!!$                e3t_n(:,:,:)=e3t_0(:,:,:)
-!!$                e3t_b(:,:,:)=e3t_0(:,:,:)
+               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
+               ssh(:,:  ,Kmm) = ssh(:,:  ,Kbb)   ! needed later for gde3w
                !
             END IF           ! end of ll_wd edits
@@ -958,6 +958,6 @@
          !                                           ! all cases !
          !                                           ! --------- !
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
          !                                           ! ----------------------- !
          IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/sbcmod.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/sbcmod.F90	(revision 11771)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/sbcmod.F90	(revision 11771)
@@ -0,0 +1,590 @@
+MODULE sbcmod
+   !!======================================================================
+   !!                       ***  MODULE  sbcmod  ***
+   !! Surface module :  provide to the ocean its surface boundary condition
+   !!======================================================================
+   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
+   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
+   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
+   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
+   !!             -   ! 2010-10  (J. Chanut, C. Bricaud, G. Madec)  add the surface pressure forcing
+   !!            3.4  ! 2011-11  (C. Harris) CICE added as an option
+   !!            3.5  ! 2012-11  (A. Coward, G. Madec) Rethink of heat, mass and salt surface fluxes
+   !!            3.6  ! 2014-11  (P. Mathiot, C. Harris) add ice shelves melting
+   !!            4.0  ! 2016-06  (L. Brodeau) new general bulk formulation
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   sbc_init      : read namsbc namelist
+   !!   sbc           : surface ocean momentum, heat and freshwater boundary conditions
+   !!   sbc_final     : Finalize CICE ice model (if used)
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers
+   USE dom_oce        ! ocean space and time domain
+   USE phycst         ! physical constants
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE trc_oce        ! shared ocean-passive tracers variables
+   USE sbc_ice        ! Surface boundary condition: ice fields
+   USE sbcdcy         ! surface boundary condition: diurnal cycle
+   USE sbcssm         ! surface boundary condition: sea-surface mean variables
+   USE sbcflx         ! surface boundary condition: flux formulation
+   USE sbcblk         ! surface boundary condition: bulk formulation
+   USE sbcice_if      ! surface boundary condition: ice-if sea-ice model
+#if defined key_si3
+   USE icestp         ! surface boundary condition: SI3 sea-ice model
+#endif
+   USE sbcice_cice    ! surface boundary condition: CICE sea-ice model
+   USE sbcisf         ! surface boundary condition: ice-shelf
+   USE sbccpl         ! surface boundary condition: coupled formulation
+   USE cpl_oasis3     ! OASIS routines for coupling
+   USE sbcssr         ! surface boundary condition: sea surface restoring
+   USE sbcrnf         ! surface boundary condition: runoffs
+   USE sbcapr         ! surface boundary condition: atmo pressure 
+   USE sbcisf         ! surface boundary condition: ice shelf
+   USE sbcfwb         ! surface boundary condition: freshwater budget
+   USE icbstp         ! Icebergs
+   USE icb_oce  , ONLY : ln_passive_mode      ! iceberg interaction mode
+   USE traqsr         ! active tracers: light penetration
+   USE sbcwave        ! Wave module
+   USE bdy_oce   , ONLY: ln_bdy
+   USE usrdef_sbc     ! user defined: surface boundary condition
+   USE closea         ! closed sea
+   !
+   USE prtctl         ! Print control                    (prt_ctl routine)
+   USE iom            ! IOM library
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE timing         ! Timing
+   USE wet_dry
+   USE diu_bulk, ONLY:   ln_diurnal_only   ! diurnal SST diagnostic
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc        ! routine called by step.F90
+   PUBLIC   sbc_init   ! routine called by opa.F90
+
+   INTEGER ::   nsbc   ! type of surface boundary condition (deduced from namsbc informations)
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: sbcmod.F90 11480 2019-08-29 09:23:25Z davestorkey $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_init( Kbb, Kmm, Kaa )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_init ***
+      !!
+      !! ** Purpose :   Initialisation of the ocean surface boundary computation
+      !!
+      !! ** Method  :   Read the namsbc namelist and set derived parameters
+      !!                Call init routines for all other SBC modules that have one
+      !!
+      !! ** Action  : - read namsbc parameters
+      !!              - nsbc: type of sbc
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa         ! ocean time level indices
+      INTEGER ::   ios, icpt                         ! local integer
+      LOGICAL ::   ll_purecpl, ll_opa, ll_not_nemo   ! local logical
+      !!
+      NAMELIST/namsbc/ nn_fsbc  ,                                                    &
+         &             ln_usr   , ln_flx   , ln_blk       ,                          &
+         &             ln_cpl   , ln_mixcpl, nn_components,                          &
+         &             nn_ice   , ln_ice_embd,                                       &
+         &             ln_traqsr, ln_dm2dc ,                                         &
+         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     &
+         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor   ,     &
+         &             ln_tauw  , nn_lsm, nn_sdrift
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_init : surface boundary condition setting'
+         WRITE(numout,*) '~~~~~~~~ '
+      ENDIF
+      !
+      !                       !**  read Surface Module namelist
+      REWIND( numnam_ref )          !* Namelist namsbc in reference namelist : Surface boundary
+      READ  ( numnam_ref, namsbc, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc in reference namelist', lwp )
+      REWIND( numnam_cfg )          !* Namelist namsbc in configuration namelist : Parameters of the run
+      READ  ( numnam_cfg, namsbc, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc in configuration namelist', lwp )
+      IF(lwm) WRITE( numond, namsbc )
+      !
+#if defined key_mpp_mpi
+      ncom_fsbc = nn_fsbc    ! make nn_fsbc available for lib_mpp
+#endif
+      !                             !* overwrite namelist parameter using CPP key information
+#if defined key_agrif
+      IF( Agrif_Root() ) THEN                ! AGRIF zoom (cf r1242: possibility to run without ice in fine grid)
+         IF( lk_si3  )   nn_ice      = 2
+         IF( lk_cice )   nn_ice      = 3
+      ENDIF
+#else
+      IF( lk_si3  )   nn_ice      = 2
+      IF( lk_cice )   nn_ice      = 3
+#endif
+      !
+      IF(lwp) THEN                  !* Control print
+         WRITE(numout,*) '   Namelist namsbc (partly overwritten with CPP key setting)'
+         WRITE(numout,*) '      frequency update of sbc (and ice)             nn_fsbc       = ', nn_fsbc
+         WRITE(numout,*) '      Type of air-sea fluxes : '
+         WRITE(numout,*) '         user defined formulation                   ln_usr        = ', ln_usr
+         WRITE(numout,*) '         flux         formulation                   ln_flx        = ', ln_flx
+         WRITE(numout,*) '         bulk         formulation                   ln_blk        = ', ln_blk
+         WRITE(numout,*) '      Type of coupling (Ocean/Ice/Atmosphere) : '
+         WRITE(numout,*) '         ocean-atmosphere coupled formulation       ln_cpl        = ', ln_cpl
+         WRITE(numout,*) '         mixed forced-coupled     formulation       ln_mixcpl     = ', ln_mixcpl
+!!gm  lk_oasis is controlled by key_oasis3  ===>>>  It shoud be removed from the namelist 
+         WRITE(numout,*) '         OASIS coupling (with atm or sas)           lk_oasis      = ', lk_oasis
+         WRITE(numout,*) '         components of your executable              nn_components = ', nn_components
+         WRITE(numout,*) '      Sea-ice : '
+         WRITE(numout,*) '         ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice
+         WRITE(numout,*) '         ice embedded into ocean                    ln_ice_embd   = ', ln_ice_embd
+         WRITE(numout,*) '      Misc. options of sbc : '
+         WRITE(numout,*) '         Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr
+         WRITE(numout,*) '            daily mean to diurnal cycle qsr            ln_dm2dc   = ', ln_dm2dc
+         WRITE(numout,*) '         Sea Surface Restoring on SST and/or SSS    ln_ssr        = ', ln_ssr
+         WRITE(numout,*) '         FreshWater Budget control  (=0/1/2)        nn_fwb        = ', nn_fwb
+         WRITE(numout,*) '         Patm gradient added in ocean & ice Eqs.    ln_apr_dyn    = ', ln_apr_dyn
+         WRITE(numout,*) '         runoff / runoff mouths                     ln_rnf        = ', ln_rnf
+         WRITE(numout,*) '         iceshelf formulation                       ln_isf        = ', ln_isf
+         WRITE(numout,*) '         nb of iterations if land-sea-mask applied  nn_lsm        = ', nn_lsm
+         WRITE(numout,*) '         surface wave                               ln_wave       = ', ln_wave
+         WRITE(numout,*) '               Stokes drift corr. to vert. velocity ln_sdw        = ', ln_sdw
+         WRITE(numout,*) '                  vertical parametrization          nn_sdrift     = ', nn_sdrift
+         WRITE(numout,*) '               wave modified ocean stress           ln_tauwoc     = ', ln_tauwoc
+         WRITE(numout,*) '               wave modified ocean stress component ln_tauw       = ', ln_tauw
+         WRITE(numout,*) '               Stokes coriolis term                 ln_stcor      = ', ln_stcor
+         WRITE(numout,*) '               neutral drag coefficient (CORE,NCAR) ln_cdgw       = ', ln_cdgw
+      ENDIF
+      !
+      IF( .NOT.ln_wave ) THEN
+         ln_sdw = .false. ; ln_cdgw = .false. ; ln_tauwoc = .false. ; ln_tauw = .false. ; ln_stcor = .false.
+      ENDIF 
+      IF( ln_sdw ) THEN
+         IF( .NOT.(nn_sdrift==jp_breivik_2014 .OR. nn_sdrift==jp_li_2017 .OR. nn_sdrift==jp_peakfr) ) &
+            CALL ctl_stop( 'The chosen nn_sdrift for Stokes drift vertical velocity must be 0, 1, or 2' )
+      ENDIF
+      ll_st_bv2014  = ( nn_sdrift==jp_breivik_2014 )
+      ll_st_li2017  = ( nn_sdrift==jp_li_2017 )
+      ll_st_bv_li   = ( ll_st_bv2014 .OR. ll_st_li2017 )
+      ll_st_peakfr  = ( nn_sdrift==jp_peakfr )
+      IF( ln_tauwoc .AND. ln_tauw ) &
+         CALL ctl_stop( 'More than one method for modifying the ocean stress has been selected ', &
+                                  '(ln_tauwoc=.true. and ln_tauw=.true.)' )
+      IF( ln_tauwoc ) &
+         CALL ctl_warn( 'You are subtracting the wave stress to the ocean (ln_tauwoc=.true.)' )
+      IF( ln_tauw ) &
+         CALL ctl_warn( 'The wave modified ocean stress components are used (ln_tauw=.true.) ', &
+                              'This will override any other specification of the ocean stress' )
+      !
+      IF( .NOT.ln_usr ) THEN     ! the model calendar needs some specificities (except in user defined case)
+         IF( MOD( rday , rdt ) /= 0. )   CALL ctl_stop( 'the time step must devide the number of second of in a day' )
+         IF( MOD( rday , 2.  ) /= 0. )   CALL ctl_stop( 'the number of second of in a day must be an even number'    )
+         IF( MOD( rdt  , 2.  ) /= 0. )   CALL ctl_stop( 'the time step (in second) must be an even number'           )
+      ENDIF
+      !                       !**  check option consistency
+      !
+      IF(lwp) WRITE(numout,*)       !* Single / Multi - executable (NEMO / OPA+SAS) 
+      SELECT CASE( nn_components )
+      CASE( jp_iam_nemo )
+         IF(lwp) WRITE(numout,*) '   ==>>>   NEMO configured as a single executable (i.e. including both OPA and Surface module)'
+      CASE( jp_iam_opa  )
+         IF(lwp) WRITE(numout,*) '   ==>>>   Multi executable configuration. Here, OPA component'
+         IF( .NOT.lk_oasis )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
+         IF( ln_cpl        )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_cpl = T in OPA'   )
+         IF( ln_mixcpl     )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      CASE( jp_iam_sas  )
+         IF(lwp) WRITE(numout,*) '   ==>>>   Multi executable configuration. Here, SAS component'
+         IF( .NOT.lk_oasis )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but key_oasis3 disabled' )
+         IF( ln_mixcpl     )   CALL ctl_stop( 'sbc_init : OPA-SAS coupled via OASIS, but ln_mixcpl = T in OPA' )
+      CASE DEFAULT
+         CALL ctl_stop( 'sbc_init : unsupported value for nn_components' )
+      END SELECT
+      !                             !* coupled options
+      IF( ln_cpl ) THEN
+         IF( .NOT. lk_oasis )   CALL ctl_stop( 'sbc_init : coupled mode with an atmosphere model (ln_cpl=T)',   &
+            &                                  '           required to defined key_oasis3' )
+      ENDIF
+      IF( ln_mixcpl ) THEN
+         IF( .NOT. lk_oasis )   CALL ctl_stop( 'sbc_init : mixed forced-coupled mode (ln_mixcpl=T) ',   &
+            &                                  '           required to defined key_oasis3' )
+         IF( .NOT.ln_cpl    )   CALL ctl_stop( 'sbc_init : mixed forced-coupled mode (ln_mixcpl=T) requires ln_cpl = T' )
+         IF( nn_components /= jp_iam_nemo )    &
+            &                   CALL ctl_stop( 'sbc_init : the mixed forced-coupled mode (ln_mixcpl=T) ',   &
+            &                                   '          not yet working with sas-opa coupling via oasis' )
+      ENDIF
+      !                             !* sea-ice
+      SELECT CASE( nn_ice )
+      CASE( 0 )                        !- no ice in the domain
+      CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)  
+      CASE( 2 )                        !- SI3  ice model
+      CASE( 3 )                        !- CICE ice model
+         IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl = T' )
+         IF( lk_agrif                    )   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
+      CASE DEFAULT                     !- not supported
+      END SELECT
+      IF( ln_diurnal .AND. .NOT. ln_blk  )   CALL ctl_stop( "sbc_init: diurnal flux processing only implemented for bulk forcing" )
+      !
+      !                       !**  allocate and set required variables
+      !
+      !                             !* allocate sbc arrays
+      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_oce arrays' )
+#if ! defined key_si3 && ! defined key_cice
+      IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop( 'sbc_init : unable to allocate sbc_ice arrays' )
+#endif
+      !
+      IF( .NOT.ln_isf ) THEN        !* No ice-shelf in the domain : allocate and set to zero
+         IF( sbc_isf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_isf arrays' )
+         fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp
+         fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp
+      END IF
+      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero
+         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case
+      ENDIF
+      !
+      sfx   (:,:) = 0._wp           !* salt flux due to freezing/melting
+      fmmflx(:,:) = 0._wp           !* freezing minus melting flux
+
+      taum(:,:) = 0._wp             !* wind stress module (needed in GLS in case of reduced restart)
+
+      !                          ! Choice of the Surface Boudary Condition (set nsbc)
+      IF( ln_dm2dc ) THEN           !* daily mean to diurnal cycle
+         nday_qsr = -1   ! allow initialization at the 1st call
+         IF( .NOT.( ln_flx .OR. ln_blk ) .AND. nn_components /= jp_iam_opa )   &
+            &   CALL ctl_stop( 'qsr diurnal cycle from daily values requires a flux or bulk formulation' )
+      ENDIF
+      !                             !* Choice of the Surface Boudary Condition
+      !                             (set nsbc)
+      !
+      ll_purecpl  = ln_cpl .AND. .NOT.ln_mixcpl
+      ll_opa      = nn_components == jp_iam_opa
+      ll_not_nemo = nn_components /= jp_iam_nemo
+      icpt = 0
+      !
+      IF( ln_usr          ) THEN   ;   nsbc = jp_usr     ; icpt = icpt + 1   ;   ENDIF       ! user defined         formulation
+      IF( ln_flx          ) THEN   ;   nsbc = jp_flx     ; icpt = icpt + 1   ;   ENDIF       ! flux                 formulation
+      IF( ln_blk          ) THEN   ;   nsbc = jp_blk     ; icpt = icpt + 1   ;   ENDIF       ! bulk                 formulation
+      IF( ll_purecpl      ) THEN   ;   nsbc = jp_purecpl ; icpt = icpt + 1   ;   ENDIF       ! Pure Coupled         formulation
+      IF( ll_opa          ) THEN   ;   nsbc = jp_none    ; icpt = icpt + 1   ;   ENDIF       ! opa coupling via SAS module
+      !
+      IF( icpt /= 1 )    CALL ctl_stop( 'sbc_init : choose ONE and only ONE sbc option' )
+      !
+      IF(lwp) THEN                     !- print the choice of surface flux formulation
+         WRITE(numout,*)
+         SELECT CASE( nsbc )
+         CASE( jp_usr     )   ;   WRITE(numout,*) '   ==>>>   user defined forcing formulation'
+         CASE( jp_flx     )   ;   WRITE(numout,*) '   ==>>>   flux formulation'
+         CASE( jp_blk     )   ;   WRITE(numout,*) '   ==>>>   bulk formulation'
+         CASE( jp_purecpl )   ;   WRITE(numout,*) '   ==>>>   pure coupled formulation'
+!!gm abusive use of jp_none ??   ===>>> need to be check and changed by adding a jp_sas parameter
+         CASE( jp_none    )   ;   WRITE(numout,*) '   ==>>>   OPA coupled to SAS via oasis'
+            IF( ln_mixcpl )       WRITE(numout,*) '               + forced-coupled mixed formulation'
+         END SELECT
+         IF( ll_not_nemo  )       WRITE(numout,*) '               + OASIS coupled SAS'
+      ENDIF
+      !
+      !                             !* OASIS initialization
+      !
+      IF( lk_oasis )   CALL sbc_cpl_init( nn_ice )   ! Must be done before: (1) first time step
+      !                                              !                      (2) the use of nn_fsbc
+      !     nn_fsbc initialization if OPA-SAS coupling via OASIS
+      !     SAS time-step has to be declared in OASIS (mandatory) -> nn_fsbc has to be modified accordingly
+      IF( nn_components /= jp_iam_nemo ) THEN
+         IF( nn_components == jp_iam_opa )   nn_fsbc = cpl_freq('O_SFLX') / NINT(rdt)
+         IF( nn_components == jp_iam_sas )   nn_fsbc = cpl_freq('I_SFLX') / NINT(rdt)
+         !
+         IF(lwp)THEN
+            WRITE(numout,*)
+            WRITE(numout,*)"   OPA-SAS coupled via OASIS : nn_fsbc re-defined from OASIS namcouple ", nn_fsbc
+            WRITE(numout,*)
+         ENDIF
+      ENDIF
+      !
+      !                             !* check consistency between model timeline and nn_fsbc
+      IF( MOD( nitend - nit000 + 1, nn_fsbc) /= 0 .OR.   &
+          MOD( nstock             , nn_fsbc) /= 0 ) THEN
+         WRITE(ctmp1,*) 'sbc_init : experiment length (', nitend - nit000 + 1, ') or nstock (', nstock,   &
+            &           ' is NOT a multiple of nn_fsbc (', nn_fsbc, ')'
+         CALL ctl_stop( ctmp1, 'Impossible to properly do model restart' )
+      ENDIF
+      !
+      IF( MOD( rday, REAL(nn_fsbc, wp) * rdt ) /= 0 )   &
+         &  CALL ctl_warn( 'sbc_init : nn_fsbc is NOT a multiple of the number of time steps in a day' )
+      !
+      IF( ln_dm2dc .AND. NINT(rday) / ( nn_fsbc * NINT(rdt) ) < 8  )   &
+         &   CALL ctl_warn( 'sbc_init : diurnal cycle for qsr: the sampling of the diurnal cycle is too small...' )
+      !
+   
+      !                       !**  associated modules : initialization
+      !
+                          CALL sbc_ssm_init ( Kbb, Kmm ) ! Sea-surface mean fields initialization
+      !
+      IF( ln_blk      )   CALL sbc_blk_init              ! bulk formulae initialization
+
+      IF( ln_ssr      )   CALL sbc_ssr_init              ! Sea-Surface Restoring initialization
+      !
+      IF( ln_isf      )   CALL sbc_isf_init( Kmm )       ! Compute iceshelves
+      !
+                          CALL sbc_rnf_init( Kmm )       ! Runof initialization
+      !
+      IF( ln_apr_dyn )    CALL sbc_apr_init              ! Atmo Pressure Forcing initialization
+      !
+#if defined key_si3
+      IF( lk_agrif .AND. nn_ice == 0 ) THEN            ! allocate ice arrays in case agrif + ice-model + no-ice in child grid
+                          IF( sbc_ice_alloc() /= 0 )   CALL ctl_stop('STOP', 'sbc_ice_alloc : unable to allocate arrays' )
+      ELSEIF( nn_ice == 2 ) THEN
+                          CALL ice_init( Kbb, Kmm, Kaa )         ! ICE initialization
+      ENDIF
+#endif
+      IF( nn_ice == 3 )   CALL cice_sbc_init( nsbc, Kbb, Kmm )   ! CICE initialization
+      !
+      IF( ln_wave     )   CALL sbc_wave_init                     ! surface wave initialisation
+      !
+      IF( lwxios ) THEN
+         CALL iom_set_rstw_var_active('utau_b')
+         CALL iom_set_rstw_var_active('vtau_b')
+         CALL iom_set_rstw_var_active('qns_b')
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_set_rstw_var_active('qsr_b')
+         CALL iom_set_rstw_var_active('emp_b')
+         CALL iom_set_rstw_var_active('sfx_b')
+      ENDIF
+
+   END SUBROUTINE sbc_init
+
+
+   SUBROUTINE sbc( kt, Kbb, Kmm )
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc  ***
+      !!
+      !! ** Purpose :   provide at each time-step the ocean surface boundary
+      !!                condition (momentum, heat and freshwater fluxes)
+      !!
+      !! ** Method  :   blah blah  to be written ?????????
+      !!                CAUTION : never mask the surface stress field (tke sbc)
+      !!
+      !! ** Action  : - set the ocean surface boundary condition at before and now
+      !!                time step, i.e.
+      !!                utau_b, vtau_b, qns_b, qsr_b, emp_n, sfx_b, qrp_b, erp_b
+      !!                utau  , vtau  , qns  , qsr  , emp  , sfx  , qrp  , erp
+      !!              - updte the ice fraction : fr_i
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      !
+      LOGICAL ::   ll_sas, ll_opa   ! local logical
+      !
+      REAL(wp) ::     zthscl        ! wd  tanh scale
+      REAL(wp), DIMENSION(jpi,jpj) ::  zwdht, zwght  ! wd dep over wd limit, wgt  
+
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('sbc')
+      !
+      !                                            ! ---------------------------------------- !
+      IF( kt /= nit000 ) THEN                      !          Swap of forcing fields          !
+         !                                         ! ---------------------------------------- !
+         utau_b(:,:) = utau(:,:)                         ! Swap the ocean forcing fields
+         vtau_b(:,:) = vtau(:,:)                         ! (except at nit000 where before fields
+         qns_b (:,:) = qns (:,:)                         !  are set at the end of the routine)
+         emp_b (:,:) = emp (:,:)
+         sfx_b (:,:) = sfx (:,:)
+         IF ( ln_rnf ) THEN
+            rnf_b    (:,:  ) = rnf    (:,:  )
+            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:)
+         ENDIF
+         IF( ln_isf )  THEN
+            fwfisf_b  (:,:  ) = fwfisf  (:,:  )               
+            risf_tsc_b(:,:,:) = risf_tsc(:,:,:)              
+         ENDIF
+        !
+      ENDIF
+      !                                            ! ---------------------------------------- !
+      !                                            !        forcing field computation         !
+      !                                            ! ---------------------------------------- !
+      !
+      ll_sas = nn_components == jp_iam_sas               ! component flags
+      ll_opa = nn_components == jp_iam_opa
+      !
+      IF( .NOT.ll_sas )   CALL sbc_ssm ( kt, Kbb, Kmm )  ! mean ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
+      IF( ln_wave     )   CALL sbc_wave( kt, Kmm )       ! surface waves
+
+      !
+      !                                            !==  sbc formulation  ==!
+      !                                                   
+      SELECT CASE( nsbc )                                ! Compute ocean surface boundary condition
+      !                                                  ! (i.e. utau,vtau, qns, qsr, emp, sfx)
+      CASE( jp_usr   )     ;   CALL usrdef_sbc_oce( kt, Kbb, Kmm )                   ! user defined formulation 
+      CASE( jp_flx     )   ;   CALL sbc_flx       ( kt )                             ! flux formulation
+      CASE( jp_blk     )
+         IF( ll_sas    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! OPA-SAS coupling: SAS receiving fields from OPA
+                               CALL sbc_blk       ( kt )                             ! bulk formulation for the ocean
+                               !
+      CASE( jp_purecpl )   ;   CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! pure coupled formulation
+      CASE( jp_none    )
+         IF( ll_opa    )       CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! OPA-SAS coupling: OPA receiving fields from SAS
+      END SELECT
+      !
+      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice, Kbb, Kmm )  ! forced-coupled mixed formulation after forcing
+      !
+      IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )              ! Wind stress provided by waves 
+      !
+      !                                            !==  Misc. Options  ==!
+      !
+      SELECT CASE( nn_ice )                                       ! Update heat and freshwater fluxes over sea-ice areas
+      CASE(  1 )   ;         CALL sbc_ice_if   ( kt, Kbb, Kmm )   ! Ice-cover climatology ("Ice-if" model)
+#if defined key_si3
+      CASE(  2 )   ;         CALL ice_stp  ( kt, Kbb, Kmm, nsbc ) ! SI3 ice model
+#endif
+      CASE(  3 )   ;         CALL sbc_ice_cice ( kt, nsbc )       ! CICE ice model
+      END SELECT
+
+      IF( ln_icebergs    )   THEN
+                                     CALL icb_stp( kt )           ! compute icebergs
+         ! icebergs may advect into haloes during the icb step and alter emp.
+         ! A lbc_lnk is necessary here to ensure restartability (#2113)
+         IF( .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) ! ensure restartability with icebergs
+      ENDIF
+
+      IF( ln_isf         )   CALL sbc_isf( kt, Kmm )                   ! compute iceshelves
+
+      IF( ln_rnf         )   CALL sbc_rnf( kt )                        ! add runoffs to fresh water fluxes
+
+      IF( ln_ssr         )   CALL sbc_ssr( kt )                        ! add SST/SSS damping term
+
+      IF( nn_fwb    /= 0 )   CALL sbc_fwb( kt, nn_fwb, nn_fsbc, Kmm )  ! control the freshwater budget
+
+      ! Special treatment of freshwater fluxes over closed seas in the model domain
+      ! Should not be run if ln_diurnal_only
+      IF( l_sbc_clo .AND. (.NOT. ln_diurnal_only) )   CALL sbc_clo( kt )   
+
+!!$!RBbug do not understand why see ticket 667
+!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why.
+!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. )
+      IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now
+         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999 
+         zwdht(:,:) = ssh(:,:,Kmm) + ht_0(:,:) - rn_wdmin1   ! do this calc of water
+                                                     ! depth above wd limit once
+         WHERE( zwdht(:,:) <= 0.0 )
+            taum(:,:) = 0.0
+            utau(:,:) = 0.0
+            vtau(:,:) = 0.0
+            qns (:,:) = 0.0
+            qsr (:,:) = 0.0
+            emp (:,:) = min(emp(:,:),0.0) !can allow puddles to grow but not shrink
+            sfx (:,:) = 0.0
+         END WHERE
+         zwght(:,:) = tanh(zthscl*zwdht(:,:))
+         WHERE( zwdht(:,:) > 0.0  .and. zwdht(:,:) < rn_wd_sbcdep ) !  5 m hard limit here is arbitrary
+            qsr  (:,:) =  qsr(:,:)  * zwght(:,:)
+            qns  (:,:) =  qns(:,:)  * zwght(:,:)
+            taum (:,:) =  taum(:,:) * zwght(:,:)
+            utau (:,:) =  utau(:,:) * zwght(:,:)
+            vtau (:,:) =  vtau(:,:) * zwght(:,:)
+            sfx  (:,:) =  sfx(:,:)  * zwght(:,:)
+            emp  (:,:) =  emp(:,:)  * zwght(:,:)
+         END WHERE
+      ENDIF
+      !
+      IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
+         !                                             ! ---------------------------------------- !
+         IF( ln_rstart .AND.    &                               !* Restart: read in restart file
+            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file'
+            CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b, ldxios = lrxios )   ! before i-stress  (U-point)
+            CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b, ldxios = lrxios )   ! before j-stress  (V-point)
+            CALL iom_get( numror, jpdom_autoglo, 'qns_b' , qns_b, ldxios = lrxios  )   ! before non solar heat flux (T-point)
+            ! The 3D heat content due to qsr forcing is treated in traqsr
+            ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point)
+            CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b, ldxios = lrxios  )    ! before     freshwater flux (T-point)
+            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6
+            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN
+               CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b, ldxios = lrxios )  ! before salt flux (T-point)
+            ELSE
+               sfx_b (:,:) = sfx(:,:)
+            ENDIF
+         ELSE                                                   !* no restart: set from nit000 values
+            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields set to nit000'
+            utau_b(:,:) = utau(:,:)
+            vtau_b(:,:) = vtau(:,:)
+            qns_b (:,:) = qns (:,:)
+            emp_b (:,:) = emp (:,:)
+            sfx_b (:,:) = sfx (:,:)
+         ENDIF
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+         !                                             ! ---------------------------------------- !
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
+            &                    'at it= ', kt,' date= ', ndastp
+         IF(lwp) WRITE(numout,*) '~~~~'
+         IF( lwxios ) CALL iom_swap(      cwxios_context          )
+         CALL iom_rstput( kt, nitrst, numrow, 'utau_b' , utau, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'vtau_b' , vtau, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'qns_b'  , qns, ldxios = lwxios  )
+         ! The 3D heat content due to qsr forcing is treated in traqsr
+         ! CALL iom_rstput( kt, nitrst, numrow, 'qsr_b'  , qsr  )
+         CALL iom_rstput( kt, nitrst, numrow, 'emp_b'  , emp, ldxios = lwxios  )
+         CALL iom_rstput( kt, nitrst, numrow, 'sfx_b'  , sfx, ldxios = lwxios  )
+         IF( lwxios ) CALL iom_swap(      cxios_context          )
+      ENDIF
+      !                                                ! ---------------------------------------- !
+      !                                                !        Outputs and control print         !
+      !                                                ! ---------------------------------------- !
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
+         CALL iom_put( "empmr"  , emp    - rnf )                ! upward water flux
+         CALL iom_put( "empbmr" , emp_b  - rnf )                ! before upward water flux ( needed to recalculate the time evolution of ssh in offline )
+         CALL iom_put( "saltflx", sfx  )                        ! downward salt flux (includes virtual salt flux beneath ice in linear free surface case)
+         CALL iom_put( "fmmflx", fmmflx  )                      ! Freezing-melting water flux
+         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux
+         CALL iom_put( "qns"   , qns        )                   ! solar heat flux
+         CALL iom_put( "qsr"   ,       qsr  )                   ! solar heat flux
+         IF( nn_ice > 0 .OR. ll_opa )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction
+         CALL iom_put( "taum"  , taum       )                   ! wind stress module
+         CALL iom_put( "wspd"  , wndm       )                   ! wind speed  module over free ocean or leads in presence of sea-ice
+      ENDIF
+      !
+      CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice)
+      CALL iom_put( "vtau", vtau )   ! j-wind stress
+      !
+      IF(ln_ctl) THEN         ! print mean trends (used for debugging)
+         CALL prt_ctl(tab2d_1=fr_i                , clinfo1=' fr_i    - : ' , mask1=tmask )
+         CALL prt_ctl(tab2d_1=(emp-rnf + fwfisf)  , clinfo1=' emp-rnf - : ' , mask1=tmask )
+         CALL prt_ctl(tab2d_1=(sfx-rnf + fwfisf)  , clinfo1=' sfx-rnf - : ' , mask1=tmask )
+         CALL prt_ctl(tab2d_1=qns                 , clinfo1=' qns      - : ', mask1=tmask )
+         CALL prt_ctl(tab2d_1=qsr                 , clinfo1=' qsr      - : ', mask1=tmask )
+         CALL prt_ctl(tab3d_1=tmask               , clinfo1=' tmask    - : ', mask1=tmask, kdim=jpk )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_tem,Kmm), clinfo1=' sst      - : ', mask1=tmask, kdim=1   )
+         CALL prt_ctl(tab3d_1=ts(:,:,:,jp_sal,Kmm), clinfo1=' sss      - : ', mask1=tmask, kdim=1   )
+         CALL prt_ctl(tab2d_1=utau                , clinfo1=' utau     - : ', mask1=umask,                      &
+            &         tab2d_2=vtau                , clinfo2=' vtau     - : ', mask2=vmask )
+      ENDIF
+
+      IF( kt == nitend )   CALL sbc_final         ! Close down surface module if necessary
+      !
+      IF( ln_timing )   CALL timing_stop('sbc')
+      !
+   END SUBROUTINE sbc
+
+
+   SUBROUTINE sbc_final
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE sbc_final  ***
+      !!
+      !! ** Purpose :   Finalize CICE (if used)
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_ice == 3 )   CALL cice_sbc_final
+      !
+   END SUBROUTINE sbc_final
+
+   !!======================================================================
+END MODULE sbcmod
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/stpctl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/stpctl.F90	(revision 11770)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/stpctl.F90	(revision 11771)
@@ -42,5 +42,5 @@
 CONTAINS
 
-   SUBROUTINE stp_ctl( kt, kindic )
+   SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE stp_ctl  ***
@@ -60,4 +60,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER, INTENT(in   ) ::   Kbb, 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
@@ -148,23 +149,23 @@
       !                                   !==  error handling  ==!
       IF( ( ln_ctl .OR. lsomeoce ) .AND. (   &             ! have use mpp_max (because ln_ctl=.T.) or contains some ocean points
-         &  zmax(1) >   50._wp .OR.   &                    ! too large sea surface height ( > 50 m )
-         &  zmax(2) >   20._wp .OR.   &                    ! too large velocity ( > 20 m/s)
-!!$         &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
-!!$         &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 )
-!!$         &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice)
+         &  zmax(1) >   20._wp .OR.   &                    ! too large sea surface height ( > 20 m )
+         &  zmax(2) >   10._wp .OR.   &                    ! too large velocity ( > 10 m/s)
+!!$      &  zmax(3) >=   0._wp .OR.   &                    ! negative or zero sea surface salinity
+!!$      &  zmax(4) >= 100._wp .OR.   &                    ! too large sea surface salinity ( > 100 )
+!!$      &  zmax(4) <    0._wp .OR.   &                    ! too large sea surface salinity (keep this line for sea-ice)
          &  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
          
-         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 50 m  or  |U| > 20 m/s  or  NaN encounter in the tests'
+         WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m  or  |U| > 10 m/s  or  S <= 0  or  S >= 100  or  NaN encounter in the tests'
          WRITE(ctmp2,9100) kt,   zmax(1), ih(1) , ih(2)
          WRITE(ctmp3,9200) kt,   zmax(2), iu(1) , iu(2) , iu(3)
@@ -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
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/trazdf.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/trazdf.F90	(revision 11770)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/trazdf.F90	(revision 11771)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE tra_zdf( kt )
+   SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf  ***
@@ -50,5 +50,7 @@
       !! ** Purpose :   compute the vertical ocean tracer physics.
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                                  , INTENT(in)    :: kt                  ! ocean time-step index
+      INTEGER                                  , INTENT(in)    :: Kbb, Kmm, Krhs, Kaa ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts                 ! active tracers and RHS of tracer equation
       !
       INTEGER  ::   jk   ! Dummy loop indices
@@ -70,10 +72,10 @@
       IF( l_trdtra )   THEN                  !* Save ta and sa trends
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
-         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa)
+         ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa)
       ENDIF
       !
       !                                      !* compute lateral mixing trend and add it to the general trend
-      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts ) 
+      CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, Krhs, pts, Kaa, jpts ) 
 
 !!gm WHY here !   and I don't like that !
@@ -81,24 +83,24 @@
       ! JMM avoid negative salinities near river outlet ! Ugly fix
       ! JMM : restore negative salinities to small salinities:
-!!$      WHERE( tsa(:,:,:,jp_sal) < 0._wp )   tsa(:,:,:,jp_sal) = 0.1_wp
+!!$   WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp )   pts(:,:,:,jp_sal,Kaa) = 0.1_wp
 !!gm
 
       IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
          DO jk = 1, jpkm1
-            ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &
-               &          / (e3t_n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)
-            ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &
-              &           / (e3t_n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)
+            ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) &
+               &          / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk)
+            ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) &
+              &           / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk)
          END DO
 !!gm this should be moved in trdtra.F90 and done on all trends
          CALL lbc_lnk_multi( 'trazdf', ztrdt, 'T', 1. , ztrds, 'T', 1. )
 !!gm
-         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
-         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt )
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds )
          DEALLOCATE( ztrdt , ztrds )
       ENDIF
       !                                          ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
-         &                       tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kaa), clinfo1=' zdf  - Ta: ', mask1=tmask,               &
+         &                       tab3d_2=pts(:,:,:,jp_sal,Kaa), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
       !
       IF( ln_timing )   CALL timing_stop('tra_zdf')
@@ -107,5 +109,5 @@
 
  
-   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt ) 
+   SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt ) 
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE tra_zdf_imp  ***
@@ -125,13 +127,13 @@
       !!      If iso-neutral mixing, add to avt the contribution due to lateral mixing.
       !!
-      !! ** Action  : - pta  becomes the after tracer
-      !!---------------------------------------------------------------------
-      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
-      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
-      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
-      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
-      REAL(wp)                             , INTENT(in   ) ::   p2dt     ! tracer time-step
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! in: tracer trend ; out: after tracer field
+      !! ** Action  : - pt(:,:,:,:,Kaa)  becomes the after tracer
+      !!---------------------------------------------------------------------
+      INTEGER                                  , INTENT(in   ) ::   kt       ! ocean time-step index
+      INTEGER                                  , INTENT(in   ) ::   Kbb, Kmm, Krhs, Kaa  ! ocean time level indices
+      INTEGER                                  , INTENT(in   ) ::   kit000   ! first time step index
+      CHARACTER(len=3)                         , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
+      INTEGER                                  , INTENT(in   ) ::   kjpt     ! number of tracers
+      REAL(wp)                                 , INTENT(in   ) ::   p2dt     ! tracer time-step
+      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::   pt       ! tracers and RHS of tracer equation
       !
       INTEGER  ::  ji, jj, jk, jn   ! dummy loop indices
@@ -181,7 +183,7 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt. (ensure same order of calculation as below if wi=0.)
-                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk  )
-                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
-                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zzwi - zzws   &
+                        zzwi = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk  ,Kmm)
+                        zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws   &
                            &                 + p2dt * ( MAX( wi(ji,jj,jk  ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
                         zwi(ji,jj,jk) = zzwi + p2dt *   MIN( wi(ji,jj,jk  ) , 0._wp )
@@ -194,7 +196,7 @@
                   DO jj = 2, jpjm1
                      DO ji = fs_2, fs_jpim1   ! vector opt.
-                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w_n(ji,jj,jk)
-                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1)
-                        zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)
+                        zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk  ) / e3w(ji,jj,jk,Kmm)
+                        zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+                        zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
                     END DO
                   END DO
@@ -218,5 +220,5 @@
             !   The solution will be in the 4d array pta.
             !   The 3d array zwt is used as a work space array.
-            !   En route to the solution pta is used a to evaluate the rhs and then 
+            !   En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 
             !   used as a work space array: its value is modified.
             !
@@ -238,5 +240,5 @@
          DO jj = 2, jpjm1           !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
             DO ji = fs_2, fs_jpim1
-               pta(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)
+               pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
             END DO
          END DO
@@ -244,6 +246,6 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  zrhs = e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn)   ! zrhs=right hand side
-                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
+                  zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs)   ! zrhs=right hand side
+                  pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
                END DO
             END DO
@@ -252,5 +254,5 @@
          DO jj = 2, jpjm1           !* 3d recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk   (result is the after tracer)
             DO ji = fs_2, fs_jpim1
-               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
+               pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
             END DO
          END DO
@@ -258,5 +260,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
-                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
+                  pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) )   &
                      &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
                END DO
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/usrdef_sbc.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/usrdef_sbc.F90	(revision 11770)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/usrdef_sbc.F90	(revision 11771)
@@ -40,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE usrdef_sbc_oce( kt, Kbb )
+   SUBROUTINE usrdef_sbc_oce( kt, Kmm, Kbb )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE usr_def_sbc  ***
@@ -54,6 +54,6 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time step
-      INTEGER, INTENT(in) ::   Kbb  ! ocean time index
+      INTEGER, INTENT(in) ::   kt        ! ocean time step
+      INTEGER, INTENT(in) ::   Kbb, Kmm  ! ocean time index
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) :: zrhoair = 1.22     ! approximate air density [Kg/m3]
@@ -88,11 +88,11 @@
          
          WHERE( ABS(gphit) <= rn_windszy/2. )
-            zwndrel(:,:) = rn_u10 - rn_uofac * un(:,:,1)
+            zwndrel(:,:) = rn_u10 - rn_uofac * uu(:,:,1,Kmm)
          ELSEWHERE
-            zwndrel(:,:) =        - rn_uofac * un(:,:,1)
+            zwndrel(:,:) =        - rn_uofac * uu(:,:,1,Kmm)
          END WHERE
          utau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:)
 
-         zwndrel(:,:) = - rn_uofac * vn(:,:,1)
+         zwndrel(:,:) = - rn_uofac * vv(:,:,1,Kmm)
          vtau(:,:) = zrhocd * zwndrel(:,:) * zwndrel(:,:)
 
