Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90	(revision 10928)
@@ -51,5 +51,5 @@
 CONTAINS
 
-   SUBROUTINE istate_init( Kbb, Kmm )
+   SUBROUTINE istate_init( Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_init  ***
@@ -57,5 +57,6 @@
       !! ** Purpose :   Initialization of the dynamics and tracer fields.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   Kbb, Kmm   ! ocean time level indices
+      INTEGER, INTENT( in )  ::  Kbb, Kmm, Kaa   ! ocean time level indices
+      !
       INTEGER ::   ji, jj, jk   ! dummy loop indices
 !!gm see comment further down
@@ -77,9 +78,9 @@
       rhd  (:,:,:  ) = 0._wp   ;   rhop (:,:,:  ) = 0._wp      ! set one for all to 0 at level jpk
       rn2b (:,:,:  ) = 0._wp   ;   rn2  (:,:,:  ) = 0._wp      ! set one for all to 0 at levels 1 and jpk
-      tsa  (:,:,:,:) = 0._wp                                   ! set one for all to 0 at level jpk
+      ts  (:,:,:,:,Kaa) = 0._wp                                   ! set one for all to 0 at level jpk
       rab_b(:,:,:,:) = 0._wp   ;   rab_n(:,:,:,:) = 0._wp      ! set one for all to 0 at level jpk
 #if defined key_agrif
-      ua   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization
-      va   (:,:,:  ) = 0._wp   ! used in agrif_oce_sponge at initialization    
+      uu   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization
+      vv   (:,:,:  ,Kaa) = 0._wp   ! used in agrif_oce_sponge at initialization    
 #endif
 
@@ -98,9 +99,9 @@
          !
          IF( ln_tsd_init ) THEN               
-            CALL dta_tsd( nit000, tsb )       ! read 3D T and S data at nit000
+            CALL dta_tsd( nit000, ts(:,:,:,:,Kbb) )       ! read 3D T and S data at nit000
             !
-            sshb(:,:)   = 0._wp               ! set the ocean at rest
+            ssh(:,:,Kbb)   = 0._wp               ! set the ocean at rest
             IF( ll_wd ) THEN
-               sshb(:,:) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
+               ssh(:,:,Kbb) =  -ssh_ref  ! Added in 30 here for bathy that adds 30 as Iterative test CEOD 
                !
                ! Apply minimum wetdepth criterion
@@ -108,26 +109,26 @@
                DO jj = 1,jpj
                   DO ji = 1,jpi
-                     IF( ht_0(ji,jj) + sshb(ji,jj)  < rn_wdmin1 ) THEN
-                        sshb(ji,jj) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
+                     IF( ht_0(ji,jj) + ssh(ji,jj,Kbb)  < rn_wdmin1 ) THEN
+                        ssh(ji,jj,Kbb) = tmask(ji,jj,1)*( rn_wdmin1 - (ht_0(ji,jj)) )
                      ENDIF
                   END DO
                END DO 
             ENDIF 
-            ub  (:,:,:) = 0._wp
-            vb  (:,:,:) = 0._wp  
+            uu  (:,:,:,Kbb) = 0._wp
+            vv  (:,:,:,Kbb) = 0._wp  
             !
          ELSE                                 ! user defined initial T and S
-            CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )         
+            CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )         
          ENDIF
-         tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
-         sshn (:,:)     = sshb(:,:)   
-         un   (:,:,:)   = ub  (:,:,:)
-         vn   (:,:,:)   = vb  (:,:,:)
-         hdivn(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
-         CALL div_hor( 0, Kmm )               ! compute interior hdivn value  
-!!gm                                    hdivn(:,:,:) = 0._wp
+         ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
+         ssh (:,:,Kmm)     = ssh(:,:,Kbb)   
+         uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
+         vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
+         hdiv(:,:,jpk) = 0._wp               ! bottom divergence set one for 0 to zero at jpk level
+         CALL div_hor( 0, Kmm )              ! compute interior hdiv value  
+!!gm                                    hdiv(:,:,:) = 0._wp
 
 !!gm POTENTIAL BUG :
-!!gm  ISSUE :  if sshb /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
+!!gm  ISSUE :  if ssh(:,:,Kbb) /= 0  then, in non linear free surface, the e3._n, e3._b should be recomputed
 !!             as well as gdept and gdepw....   !!!!! 
 !!      ===>>>>   probably a call to domvvl initialisation here....
@@ -139,14 +140,14 @@
 !            ALLOCATE( zuvd(jpi,jpj,jpk,2) )
 !            CALL dta_uvd( nit000, zuvd )
-!            ub(:,:,:) = zuvd(:,:,:,1) ;  un(:,:,:) = ub(:,:,:)
-!            vb(:,:,:) = zuvd(:,:,:,2) ;  vn(:,:,:) = vb(:,:,:)
+!            uu(:,:,:,Kbb) = zuvd(:,:,:,1) ;  uu(:,:,:,Kmm) = uu(:,:,:,Kbb)
+!            vv(:,:,:,Kbb) = zuvd(:,:,:,2) ;  vv(:,:,:,Kmm) = vv(:,:,:,Kbb)
 !            DEALLOCATE( zuvd )
 !         ENDIF
          !
 !!gm This is to be changed !!!!
-!         ! - ML - sshn could be modified by istate_eel, so that initialization of e3t_b is done here
+!         ! - ML - ssh(:,:,Kmm) could be modified by istate_eel, so that initialization of e3t(:,:,:,Kbb) is done here
 !         IF( .NOT.ln_linssh ) THEN
 !            DO jk = 1, jpk
-!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
+!               e3t(:,:,jk,Kbb) = e3t(:,:,jk,Kmm)
 !            END DO
 !         ENDIF
@@ -158,25 +159,25 @@
       ! Do it whatever the free surface method, these arrays being eventually used
       !
-      un_b(:,:) = 0._wp   ;   vn_b(:,:) = 0._wp
-      ub_b(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp
+      uu_b(:,:,Kmm) = 0._wp   ;   vv_b(:,:,Kmm) = 0._wp
+      uu_b(:,:,Kbb) = 0._wp   ;   vv_b(:,:,Kbb) = 0._wp
       !
-!!gm  the use of umsak & vmask is not necessary below as un, vn, ub, vb are always masked
+!!gm  the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
       DO jk = 1, jpkm1
          DO jj = 1, jpj
             DO ji = 1, jpi
-               un_b(ji,jj) = un_b(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
-               vn_b(ji,jj) = vn_b(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
+               vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
                !
-               ub_b(ji,jj) = ub_b(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
+               uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
+               vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
             END DO
          END DO
       END DO
       !
-      un_b(:,:) = un_b(:,:) * r1_hu_n(:,:)
-      vn_b(:,:) = vn_b(:,:) * r1_hv_n(:,:)
+      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu_n(:,:)
+      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv_n(:,:)
       !
-      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
-      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
+      uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu_b(:,:)
+      vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv_b(:,:)
       !
    END SUBROUTINE istate_init
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/divhor.F90	(revision 10928)
@@ -55,11 +55,11 @@
       !!
       !! ** Method  :   the now divergence is computed as :
-      !!         hdivn = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
+      !!         hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
       !!      and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 
       !!
-      !! ** Action  : - update hdivn, the now horizontal divergence
+      !! ** Action  : - update hdiv, the now horizontal divergence
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kmm  ! ocean time-level index
+      INTEGER, INTENT(in) ::   Kmm  ! ocean time level index
       !
       INTEGER  ::   ji, jj, jk    ! dummy loop indices
@@ -78,9 +78,9 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               hdivn(ji,jj,jk) = (  e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) * un(ji  ,jj,jk)      &
-                  &               - e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) * un(ji-1,jj,jk)      &
-                  &               + e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vn(ji,jj  ,jk)      &
-                  &               - e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vn(ji,jj-1,jk)  )   &
-                  &            * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
+               hdiv(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * uu(ji  ,jj,jk,Kmm)      &
+                  &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm)      &
+                  &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vv(ji,jj  ,jk,Kmm)      &
+                  &               - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm)  )   &
+                  &            * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
             END DO  
          END DO  
@@ -88,22 +88,22 @@
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN
-         IF( nbondi == -1 .OR. nbondi == 2 )   hdivn(   2   ,  :   ,:) = 0._wp      ! west
-         IF( nbondi ==  1 .OR. nbondi == 2 )   hdivn( nlci-1,  :   ,:) = 0._wp      ! east
-         IF( nbondj == -1 .OR. nbondj == 2 )   hdivn(   :   ,  2   ,:) = 0._wp      ! south
-         IF( nbondj ==  1 .OR. nbondj == 2 )   hdivn(   :   ,nlcj-1,:) = 0._wp      ! north
+         IF( nbondi == -1 .OR. nbondi == 2 )   hdiv(   2   ,  :   ,:) = 0._wp      ! west
+         IF( nbondi ==  1 .OR. nbondi == 2 )   hdiv( nlci-1,  :   ,:) = 0._wp      ! east
+         IF( nbondj == -1 .OR. nbondj == 2 )   hdiv(   :   ,  2   ,:) = 0._wp      ! south
+         IF( nbondj ==  1 .OR. nbondj == 2 )   hdiv(   :   ,nlcj-1,:) = 0._wp      ! north
       ENDIF
 #endif
       !
-      IF( ln_rnf )   CALL sbc_rnf_div( hdivn, Kmm )                     !==  runoffs    ==!   (update hdivn field)
+      IF( ln_rnf )   CALL sbc_rnf_div( hdiv, Kmm )                     !==  runoffs    ==!   (update hdiv field)
       !
 #if defined key_asminc 
-      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, hdivn )   !==  SSH assimilation  ==!   (update hdivn field)
+      IF( ln_sshinc .AND. ln_asmiau )   CALL ssh_asm_div( kt, hdiv )   !==  SSH assimilation  ==!   (update hdiv field)
       ! 
 #endif
-      IF( ln_isf )   CALL sbc_isf_div( hdivn, Kmm )                     !==  ice shelf  ==!   (update hdivn field)
+      IF( ln_isf )   CALL sbc_isf_div( hdiv, Kmm )                     !==  ice shelf  ==!   (update hdiv field)
       !
-      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdivn )             !==  ice sheet  ==!   (update hdivn field)
+      IF( ln_iscpl .AND. ln_hsb )   CALL iscpl_div( hdiv ) !==  ice sheet  ==!   (update hdiv field)
       !
-      CALL lbc_lnk( 'divhor', hdivn, 'T', 1. )   !   (no sign change)
+      CALL lbc_lnk( 'divhor', hdiv, 'T', 1. )   !   (no sign change)
       !
       IF( ln_timing )   CALL timing_stop('div_hor')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90	(revision 10928)
@@ -141,6 +141,6 @@
       ENDIF
       !                                   ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cen2 adv - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
    END SUBROUTINE dyn_adv_cen2
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90	(revision 10928)
@@ -234,6 +234,6 @@
       ENDIF
       !                                         ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=           ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
    END SUBROUTINE dyn_adv_ubs
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90	(revision 10928)
@@ -81,5 +81,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_hpg( kt )
+   SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_hpg  ***
@@ -88,8 +88,11 @@
       !!              using the scheme defined in the namelist
       !!
-      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
       !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
+      !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
       !!----------------------------------------------------------------------
@@ -97,28 +100,28 @@
       IF( ln_timing )   CALL timing_start('dyn_hpg')
       !
-      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
+      IF( l_trddyn ) THEN                    ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn)
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
-         ztrdu(:,:,:) = ua(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
       ENDIF
       !
       SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
-      CASE ( np_zco )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
-      CASE ( np_zps )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
-      CASE ( np_sco )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
-      CASE ( np_djc )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
-      CASE ( np_prj )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
-      CASE ( np_isf )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf
+      CASE ( np_zco )   ;   CALL hpg_zco    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate
+      CASE ( np_zps )   ;   CALL hpg_zps    ( kt, Kmm, puu, pvv, Krhs )  ! z-coordinate plus partial steps (interpolation)
+      CASE ( np_sco )   ;   CALL hpg_sco    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (standard jacobian formulation)
+      CASE ( np_djc )   ;   CALL hpg_djc    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Density Jacobian with Cubic polynomial)
+      CASE ( np_prj )   ;   CALL hpg_prj    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate (Pressure Jacobian scheme)
+      CASE ( np_isf )   ;   CALL hpg_isf    ( kt, Kmm, puu, pvv, Krhs )  ! s-coordinate similar to sco modify for ice shelf
       END SELECT
       !
       IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )
          DEALLOCATE( ztrdu , ztrdv )
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( ln_timing )   CALL timing_stop('dyn_hpg')
@@ -127,5 +130,5 @@
 
 
-   SUBROUTINE dyn_hpg_init
+   SUBROUTINE dyn_hpg_init( Kmm )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dyn_hpg_init  ***
@@ -137,4 +140,6 @@
       !!      with the type of vertical coordinate used (zco, zps, sco)
       !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: Kmm   ! ocean time level index
+      !
       INTEGER ::   ioptio = 0      ! temporary integer
       INTEGER ::   ios             ! Local integer output status for namelist read
@@ -228,5 +233,5 @@
          !
          DO jk = 1, jpk                   !- compute density of the water displaced by the ice shelf 
-            CALL eos( zts_top(:,:,:), gdept_n(:,:,jk), zrhd(:,:,jk) )
+            CALL eos( zts_top(:,:,:), gdept(:,:,jk,Kmm), zrhd(:,:,jk) )
          END DO
          !
@@ -239,11 +244,11 @@
             DO ji = 1, jpi                      ! divided by 2 later
                ikt = mikt(ji,jj)
-               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))
+               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w(ji,jj,1,Kmm) * (1._wp - tmask(ji,jj,1))
                DO jk = 2, ikt-1
-                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &
+                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w(ji,jj,jk,Kmm) &
                      &                              * (1._wp - tmask(ji,jj,jk))
                END DO
                IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
-                  &                                              * ( risfdep(ji,jj) - gdept_n(ji,jj,ikt-1) )
+                  &                                              * ( risfdep(ji,jj) - gdept(ji,jj,ikt-1,Kmm) )
             END DO
          END DO
@@ -256,5 +261,5 @@
 
 
-   SUBROUTINE hpg_zco( kt )
+   SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE hpg_zco  ***
@@ -266,11 +271,13 @@
       !!      level:    zhpi = grav .....
       !!                zhpj = grav .....
-      !!      add it to the general momentum trend (ua,va).
-      !!            ua = ua - 1/e1u * zhpi
-      !!            va = va - 1/e2v * zhpj
-      !!
-      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
+      !!            puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
+      !!            pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
+      !!
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
@@ -290,11 +297,11 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
+            zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
             ! hydrostatic pressure gradient
             zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
             zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
+            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
+            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
          END DO
       END DO
@@ -305,5 +312,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
+               zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
                ! hydrostatic pressure gradient
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
@@ -315,6 +322,6 @@
                   &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
             END DO
          END DO
@@ -324,5 +331,5 @@
 
 
-   SUBROUTINE hpg_zps( kt )
+   SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE hpg_zps  ***
@@ -330,7 +337,9 @@
       !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
       !!
-      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !! ** Action  : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
       INTEGER  ::   ji, jj, jk                       ! dummy loop indices
@@ -347,5 +356,5 @@
 
       ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level
-!jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv )
+!jc      CALL zps_hde    ( kt, jpts, ts(:,:,:,:,Kmm), gtsu, gtsv, rhd, gru , grv )
 
       ! Local constant initialization
@@ -355,11 +364,11 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
+            zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
             ! hydrostatic pressure gradient
             zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
             zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
+            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
+            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
          END DO
       END DO
@@ -369,5 +378,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
+               zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
                ! hydrostatic pressure gradient
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
@@ -379,6 +388,6 @@
                   &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
             END DO
          END DO
@@ -390,17 +399,17 @@
             iku = mbku(ji,jj)
             ikv = mbkv(ji,jj)
-            zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) )
-            zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) )
+            zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj  ,iku,Kmm) )
+            zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji  ,jj+1,ikv,Kmm) )
             IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
-               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
+               puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku)         ! subtract old value
                zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
                   &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
-               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
+               puu  (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
             ENDIF
             IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
-               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
+               pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv)         ! subtract old value
                zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
                   &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
-               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
+               pvv  (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
             ENDIF
          END DO
@@ -410,5 +419,5 @@
 
 
-   SUBROUTINE hpg_sco( kt )
+   SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE hpg_sco  ***
@@ -422,11 +431,13 @@
       !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
       !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
-      !!      add it to the general momentum trend (ua,va).
-      !!         ua = ua - 1/e1u * zhpi
-      !!         va = va - 1/e2v * zhpj
-      !!
-      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
+      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
+      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
+      !!
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
       INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
@@ -453,10 +464,10 @@
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(  sshn(ji,jj)               ,  sshn(ji+1,jj) ) >                &
+             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)               ,  ssh(ji+1,jj,Kmm) ) >                &
                   &    MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) +  ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
+                  &    MAX(  ssh(ji,jj,Kmm) +  ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
                   &                                                       > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)              -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (       &
-                  &    MAX(   sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
+             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)              -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (       &
+                  &    MAX(   ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
                   &    MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -464,17 +475,17 @@
                zcpx(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) )
              ELSE
                zcpx(ji,jj) = 0._wp
              END IF
       
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
+             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
+                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
                   &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                &
+             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        &
+                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -482,7 +493,7 @@
                zcpy(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
              ELSE
                zcpy(ji,jj) = 0._wp
@@ -497,13 +508,13 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             ! hydrostatic pressure gradient along s-surfaces
-            zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    &
-               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
-            zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    &
-               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
+            zhpi(ji,jj,1) = zcoef0 * (  e3w(ji+1,jj  ,1,Kmm) * ( znad + rhd(ji+1,jj  ,1) )    &
+               &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
+            zhpj(ji,jj,1) = zcoef0 * (  e3w(ji  ,jj+1,1,Kmm) * ( znad + rhd(ji  ,jj+1,1) )    &
+               &                      - e3w(ji  ,jj  ,1,Kmm) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
             ! s-coordinate pressure gradient correction
             zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
+               &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
             zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
+               &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
             !
             IF( ln_wd_il ) THEN
@@ -515,6 +526,6 @@
             !
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
+            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap
+            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap
          END DO
       END DO
@@ -526,14 +537,14 @@
                ! hydrostatic pressure gradient along s-surfaces
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
-                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
-                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
+                  &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
+                  &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
-                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
-                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
+                  &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
+                  &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
                ! s-coordinate pressure gradient correction
                zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
+                  &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj)
                zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
+                  &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj)
                !
                IF( ln_wd_il ) THEN
@@ -545,6 +556,6 @@
                !
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap
             END DO
          END DO
@@ -556,5 +567,5 @@
 
 
-   SUBROUTINE hpg_isf( kt )
+   SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE hpg_isf  ***
@@ -568,12 +579,14 @@
       !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
       !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
-      !!      add it to the general momentum trend (ua,va).
-      !!         ua = ua - 1/e1u * zhpi
-      !!         va = va - 1/e2v * zhpj
+      !!      add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
+      !!         puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
+      !!         pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
       !!      iceload is added and partial cell case are added to the top and bottom
       !!      
-      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
+      !!----------------------------------------------------------------------
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
       INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
@@ -596,6 +609,6 @@
         DO jj = 1, jpj
           ikt = mikt(ji,jj)
-          zts_top(ji,jj,1) = tsn(ji,jj,ikt,1)
-          zts_top(ji,jj,2) = tsn(ji,jj,ikt,2)
+          zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
+          zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
         END DO
       END DO
@@ -612,22 +625,22 @@
             ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
             ! we assume ISF is in isostatic equilibrium
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    &
+            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w(ji+1,jj,iktp1i,Kmm)                                    &
                &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
-               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &
+               &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
                &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
                &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            ) 
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    &
+            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm)                                    &
                &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
-               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
+               &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         & 
                &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
                &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            ) 
             ! s-coordinate pressure gradient correction (=0 if z coordinate)
             zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
+               &           * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
             zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
-               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
+               &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
-            va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
+            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
+            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
          END DO
       END DO
@@ -641,17 +654,17 @@
                ! hydrostatic pressure gradient along s-surfaces
                zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
-                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
+                  &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
+                  &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
                zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
-                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
+                  &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
+                  &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
                ! s-coordinate pressure gradient correction
                zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)
+                  &           * ( gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj)
                zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)
+                  &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj)
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
             END DO
          END DO
@@ -661,5 +674,5 @@
 
 
-   SUBROUTINE hpg_djc( kt )
+   SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE hpg_djc  ***
@@ -669,5 +682,7 @@
       !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
@@ -687,27 +702,27 @@
         DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
+             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
+                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
                   &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (        &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                &
+             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (        &
+                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
              IF(ll_tmp1) THEN
                zcpx(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji+1,jj) - sshn(ji  ,jj)) )
+               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji+1,jj,Kmm) - ssh(ji  ,jj,Kmm)) )
              ELSE
                zcpx(ji,jj) = 0._wp
              END IF
       
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
+             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
+                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
                   &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (        &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                &
+             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (        &
+                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -715,7 +730,7 @@
                zcpy(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
              ELSE
                zcpy(ji,jj) = 0._wp
@@ -747,9 +762,9 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1)
-               dzz  (ji,jj,jk) = gde3w_n(ji  ,jj  ,jk) - gde3w_n(ji,jj,jk-1)
+               dzz  (ji,jj,jk) = gde3w(ji  ,jj  ,jk) - gde3w(ji,jj,jk-1)
                drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
-               dzx  (ji,jj,jk) = gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk  )
+               dzx  (ji,jj,jk) = gde3w(ji+1,jj  ,jk) - gde3w(ji,jj,jk  )
                drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  )
-               dzy  (ji,jj,jk) = gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk  )
+               dzy  (ji,jj,jk) = gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk  )
             END DO
          END DO
@@ -838,9 +853,9 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) )               &
+            rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) )               &
                &                   * (  rhd(ji,jj,1)                                     &
                &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  &
-               &                              * ( e3w_n  (ji,jj,1) - gde3w_n(ji,jj,1) )  &
-               &                              / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) )  )
+               &                              * ( e3w  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  &
+               &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  )
          END DO
       END DO
@@ -854,8 +869,8 @@
 
                rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   &
-                  &                     * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )                                   &
+                  &                     * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) )                                   &
                   &            - grav * z1_10 * (                                                                   &
                   &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     &
-                  &   * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   &
+                  &   * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   &
                   &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     &
                   &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   &
@@ -863,8 +878,8 @@
 
                rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   &
-                  &                     * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                   &
+                  &                     * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                   &
                   &            - grav* z1_10 * (                                                                    &
                   &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     &
-                  &   * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   &
+                  &   * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   &
                   &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     &
                   &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   &
@@ -872,8 +887,8 @@
 
                rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 &
-                  &                     * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   &
+                  &                     * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                   &
                   &            - grav* z1_10 * (                                                                    &
                   &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     &
-                  &   * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   &
+                  &   * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   &
                   &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     &
                   &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   &
@@ -897,6 +912,6 @@
             ENDIF
             ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
+            puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
+            pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
          END DO
       END DO
@@ -920,6 +935,6 @@
                ENDIF
                ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
             END DO
          END DO
@@ -931,5 +946,5 @@
 
 
-   SUBROUTINE hpg_prj( kt )
+   SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE hpg_prj  ***
@@ -940,8 +955,10 @@
       !!      all vertical coordinate systems
       !!
-      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
+      !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
       !!----------------------------------------------------------------------
       INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
-      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  kt          ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv    ! ocean velocities and RHS of momentum equation
       !!
       INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices
@@ -975,10 +992,10 @@
          DO jj = 2, jpjm1
            DO ji = 2, jpim1 
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji+1,jj) ) >                &
+             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji+1,jj,Kmm) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji+1,jj) + ht_0(ji+1,jj) )  &
+                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) )  &
                   &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji+1,jj) ) > 1.E-12 ) .AND. (         &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji+1,jj) ) >                &
+             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. (         &
+                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji+1,jj,Kmm) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -986,7 +1003,7 @@
                zcpx(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
-               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
+               ! no worries about  ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm) = 0, it won't happen ! here
+               zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji+1,jj,Kmm) -  ssh(ji  ,jj,Kmm)) )
               
                 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
@@ -995,10 +1012,10 @@
              END IF
       
-             ll_tmp1 = MIN(  sshn(ji,jj)              ,  sshn(ji,jj+1) ) >                &
+             ll_tmp1 = MIN(  ssh(ji,jj,Kmm)              ,  ssh(ji,jj+1,Kmm) ) >                &
                   &    MAX( -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) .AND.            &
-                  &    MAX(  sshn(ji,jj) + ht_0(ji,jj),  sshn(ji,jj+1) + ht_0(ji,jj+1) )  &
+                  &    MAX(  ssh(ji,jj,Kmm) + ht_0(ji,jj),  ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) )  &
                   &                                                      > rn_wdmin1 + rn_wdmin2
-             ll_tmp2 = ( ABS( sshn(ji,jj)             -  sshn(ji,jj+1) ) > 1.E-12 ) .AND. (      &
-                  &    MAX(   sshn(ji,jj)             ,  sshn(ji,jj+1) ) >                &
+             ll_tmp2 = ( ABS( ssh(ji,jj,Kmm)             -  ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. (      &
+                  &    MAX(   ssh(ji,jj,Kmm)             ,  ssh(ji,jj+1,Kmm) ) >                &
                   &    MAX(  -ht_0(ji,jj)             , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
 
@@ -1006,7 +1023,7 @@
                zcpy(ji,jj) = 1.0_wp
              ELSE IF(ll_tmp2) THEN
-               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
-               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) &
-                           &    / (sshn(ji,jj+1) - sshn(ji,jj  )) )
+               ! no worries about  ssh(ji,jj+1,Kmm) -  ssh(ji,jj  ,Kmm) = 0, it won't happen ! here
+               zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
+                           &    / (ssh(ji,jj+1,Kmm) - ssh(ji,jj  ,Kmm)) )
                 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp)
 
@@ -1031,6 +1048,6 @@
           ELSEIF( jk < jpkm1 ) THEN
              DO jkk = jk+1, jpk
-                zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk  ), gde3w_n(ji,jj,jkk-1),   &
-                   &                      gde3w_n(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2))
+                zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk  ), gde3w(ji,jj,jkk-1),   &
+                   &                      gde3w(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2))
              END DO
           ENDIF
@@ -1041,5 +1058,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad
+            zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad
          END DO
       END DO
@@ -1048,5 +1065,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk)
+               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm)
             END DO
          END DO
@@ -1065,8 +1082,8 @@
         DO ji = 2, jpi
           zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
-             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1)
+             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)
 
           ! assuming linear profile across the top half surface layer
-          zhpi(ji,jj,1) =  0.5_wp * e3w_n(ji,jj,1) * zrhdt1
+          zhpi(ji,jj,1) =  0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1
         END DO
       END DO
@@ -1090,12 +1107,12 @@
         DO ji = 2, jpim1
 !!gm BUG ?    if it is ssh at u- & v-point then it should be:
-!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &
+!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * &
 !                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
-!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &
+!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * &
 !                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
 !!gm not this:
-          zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &
+          zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
                          & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
-          zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &
+          zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
                          & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
         END DO
@@ -1106,6 +1123,6 @@
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
-          zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 
-          zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad)
+          zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 
+          zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad)
         END DO
       END DO
@@ -1114,6 +1131,6 @@
         DO jj = 2, jpjm1
           DO ji = 2, jpim1
-            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk)
+            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)
+            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)
           END DO
         END DO
@@ -1123,6 +1140,6 @@
         DO jj = 2, jpjm1
           DO ji = 2, jpim1
-            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk)
+            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)
+            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)
           END DO
         END DO
@@ -1176,5 +1193,5 @@
                DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
                  IF( jk1 == 1 ) THEN
-                   zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
+                   zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad)
                    zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
                                                      bsp(jid,jj,1),   csp(jid,jj,1), &
@@ -1196,5 +1213,5 @@
                IF( .NOT.ln_linssh ) THEN
                  zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
-                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
+                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) )
                 ELSE
                  zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
@@ -1204,5 +1221,5 @@
                   zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
                ENDIF
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
             ENDIF
 
@@ -1234,5 +1251,5 @@
                DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
                  IF( jk1 == 1 ) THEN
-                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
+                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)
                    zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
                                                      bsp(ji,jjd,1),   csp(ji,jjd,1), &
@@ -1255,5 +1272,5 @@
                IF( .NOT.ln_linssh ) THEN
                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
-                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )
+                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )
                ELSE
                   zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
@@ -1264,5 +1281,5 @@
                ENDIF
 
-               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
             ENDIF
                !
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90	(revision 10928)
@@ -187,6 +187,6 @@
       ENDIF
       !
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' keg  - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( ln_timing )   CALL timing_stop('dyn_keg')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90	(revision 10928)
@@ -43,5 +43,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_ldf( kt )
+   SUBROUTINE dyn_ldf( kt, Kbb, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_ldf  ***
@@ -49,5 +49,7 @@
       !! ** Purpose :   compute the lateral ocean dynamics physics.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
       !
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
@@ -58,25 +60,28 @@
       IF( l_trddyn )   THEN                      ! temporary save of momentum trends
          ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
-         ztrdu(:,:,:) = ua(:,:,:) 
-         ztrdv(:,:,:) = va(:,:,:) 
+         ztrdu(:,:,:) = puu(:,:,:,Krhs) 
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs) 
       ENDIF
 
       SELECT CASE ( nldf_dyn )                   ! compute lateral mixing trend and add it to the general trend
       !
-      CASE ( np_lap   )    ;   CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 )      ! iso-level    laplacian
-      CASE ( np_lap_i )    ;   CALL dyn_ldf_iso( kt )                         ! rotated      laplacian
-      CASE ( np_blp   )    ;   CALL dyn_ldf_blp( kt, ub, vb, ua, va    )      ! iso-level bi-laplacian
+      CASE ( np_lap   )  
+         CALL dyn_ldf_lap( kt, Kbb, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs), 1 ) ! iso-level    laplacian
+      CASE ( np_lap_i ) 
+         CALL dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs    )                                                   ! rotated      laplacian
+      CASE ( np_blp   )  
+         CALL dyn_ldf_blp( kt, Kbb, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs)    ) ! iso-level bi-laplacian
       !
       END SELECT
 
       IF( l_trddyn ) THEN                        ! save the horizontal diffusive trends for further diagnostics
-         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
-         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
+         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
+         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
          CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt )
          DEALLOCATE ( ztrdu , ztrdv )
       ENDIF
       !                                          ! print sum trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf  - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldf  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( ln_timing )   CALL timing_stop('dyn_ldf')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_iso.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_iso.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_iso.F90	(revision 10928)
@@ -60,5 +60,5 @@
 
 
-   SUBROUTINE dyn_ldf_iso( kt )
+   SUBROUTINE dyn_ldf_iso( kt, Kbb, Kmm, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE dyn_ldf_iso  ***
@@ -81,27 +81,29 @@
       !!      horizontal fluxes associated with the rotated lateral mixing:
       !!      u-component:
-      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ ub ]
-      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
-      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
-      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(ub)) ]
+      !!         ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t  di[ uu ]
+      !!               -  ahmt              e2t * mi-1(uslp) dk[ mi(mk(uu)) ]
+      !!         zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f  dj[ uu ]
+      !!               -  ahmf              e1f * mi(vslp)   dk[ mj(mk(uu)) ]
       !!      v-component:
-      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vb ]
-      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vb)) ]
-      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ ub ]
-      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
+      !!         zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t  di[ vv ]
+      !!               -  ahmf              e2t * mj(uslp)   dk[ mi(mk(vv)) ]
+      !!         zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f  dj[ vv ]
+      !!               -  ahmt              e1f * mj-1(vslp) dk[ mj(mk(vv)) ]
       !!      take the horizontal divergence of the fluxes:
       !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
       !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
-      !!      Add this trend to the general trend (ua,va):
-      !!         ua = ua + diffu
+      !!      Add this trend to the general trend (uu(rhs),vv(rhs)):
+      !!         uu(rhs) = uu(rhs) + diffu
       !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
       !!      should be modified for applications others than orca_r2 (!!bug)
       !!
       !! ** Action :
-      !!       -(ua,va) updated with the before geopotential harmonic mixing trend
+      !!       -(puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the before geopotential harmonic mixing trend
       !!       -(akzu,akzv) to accompt for the diagonal vertical component
       !!                    of the rotated operator in dynzdf module
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  kt               ! ocean time-step index
+      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv         ! ocean velocities and RHS of momentum equation
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -128,8 +130,8 @@
             DO jj = 2, jpjm1
                DO ji = 2, jpim1
-                  uslp (ji,jj,jk) = - ( gdept_b(ji+1,jj,jk) - gdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
-                  vslp (ji,jj,jk) = - ( gdept_b(ji,jj+1,jk) - gdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
-                  wslpi(ji,jj,jk) = - ( gdepw_b(ji+1,jj,jk) - gdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
-                  wslpj(ji,jj,jk) = - ( gdepw_b(ji,jj+1,jk) - gdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
+                  uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
+                  vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kbb) - gdept(ji ,jj ,jk,Kbb) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
+                  wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kbb) - gdepw(ji-1,jj,jk,Kbb) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5
+                  wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kbb) - gdepw(ji,jj-1,jk,Kbb) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5
                END DO
             END DO
@@ -151,6 +153,6 @@
          !                             zdkv(jk=1)=zdkv(jk=2)
 
-         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
-         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
+         zdk1u(:,:) = ( puu(:,:,jk,Kbb) -puu(:,:,jk+1,Kbb) ) * umask(:,:,jk+1)
+         zdk1v(:,:) = ( pvv(:,:,jk,Kbb) -pvv(:,:,jk+1,Kbb) ) * vmask(:,:,jk+1)
 
          IF( jk == 1 ) THEN
@@ -158,6 +160,6 @@
             zdkv(:,:) = zdk1v(:,:)
          ELSE
-            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
-            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
+            zdku(:,:) = ( puu(:,:,jk-1,Kbb) - puu(:,:,jk,Kbb) ) * umask(:,:,jk)
+            zdkv(:,:) = ( pvv(:,:,jk-1,Kbb) - pvv(:,:,jk,Kbb) ) * vmask(:,:,jk)
          ENDIF
 
@@ -171,5 +173,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u_n(ji,jj,jk), e3u_n(ji-1,jj,jk) ) * r1_e1t(ji,jj)
+                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj)
 
                   zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)     &
@@ -178,5 +180,5 @@
                   zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
    
-                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )    &
+                  ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )    &
                      &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)      &
                      &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
@@ -186,5 +188,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, jpi   ! vector opt.
-                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t_n(ji,jj,jk) * r1_e1t(ji,jj)
+                  zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj)
 
                   zmskt = 1._wp / MAX(   umask(ji-1,jj,jk  ) + umask(ji,jj,jk+1)     &
@@ -193,5 +195,5 @@
                   zcof1 = - zaht_0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
 
-                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
+                  ziut(ji,jj) = (  zabe1 * ( puu(ji,jj,jk,Kbb) - puu(ji-1,jj,jk,Kbb) )   &
                      &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
                      &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
@@ -203,5 +205,5 @@
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f_n(ji,jj,jk) * r1_e2f(ji,jj)
+               zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj)
 
                zmskf = 1._wp / MAX(   umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)     &
@@ -210,5 +212,5 @@
                zcof2 = - zaht_0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
 
-               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
+               zjuf(ji,jj) = (  zabe2 * ( puu(ji,jj+1,jk,Kbb) - puu(ji,jj,jk,Kbb) )   &
                   &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
                   &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
@@ -224,5 +226,5 @@
          DO jj = 2, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f_n(ji,jj,jk) * r1_e1f(ji,jj)
+               zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj)
 
                zmskf = 1._wp / MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)     &
@@ -231,5 +233,5 @@
                zcof1 = - zaht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
 
-               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )    &
+               zivf(ji,jj) = (  zabe1 * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji,jj,jk,Kbb) )    &
                   &           + zcof1 * (  zdkv (ji,jj) + zdk1v(ji+1,jj)      &
                   &                      + zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
@@ -241,5 +243,5 @@
             DO jj = 2, jpj
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v_n(ji,jj,jk), e3v_n(ji,jj-1,jk) ) * r1_e2t(ji,jj)
+                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj)
 
                   zmskt = 1._wp / MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)     &
@@ -248,5 +250,5 @@
                   zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
 
-                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )    &
+                  zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )    &
                      &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)      &
                      &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
@@ -256,5 +258,5 @@
             DO jj = 2, jpj
                DO ji = 1, fs_jpim1   ! vector opt.
-                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t_n(ji,jj,jk) * r1_e2t(ji,jj)
+                  zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj)
 
                   zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
@@ -263,5 +265,5 @@
                   zcof2 = - zaht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
 
-                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
+                  zjvt(ji,jj) = (  zabe2 * ( pvv(ji,jj,jk,Kbb) - pvv(ji,jj-1,jk,Kbb) )   &
                      &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
                      &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
@@ -275,8 +277,8 @@
          DO jj = 2, jpjm1
             DO ji = 2, jpim1          !!gm Question vectop possible??? !!bug
-               ua(ji,jj,jk) = ua(ji,jj,jk) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
-                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
-                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (  ziut(ji+1,jj) - ziut(ji,jj  )    &
+                  &                           + zjuf(ji  ,jj) - zjuf(ji,jj-1)  ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (  zivf(ji,jj  ) - zivf(ji-1,jj)    &
+                  &                           + zjvt(ji,jj+1) - zjvt(ji,jj  )  ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
             END DO
          END DO
@@ -286,6 +288,6 @@
 
       ! print sum trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ldfh - Ua: ', mask1=umask, &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
 
 
@@ -306,11 +308,11 @@
             DO ji = 2, jpi
                ! i-gradient of u at jj
-               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
+               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji-1,jj  ,jk,Kbb) )
                ! j-gradient of u and v at jj
-               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
-               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
+               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( puu(ji,jj+1,jk,Kbb) - puu(ji  ,jj  ,jk,Kbb) )
+               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pvv(ji,jj  ,jk,Kbb) - pvv(ji  ,jj-1,jk,Kbb) )
                ! j-gradient of u and v at jj+1
-               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
-               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
+               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( puu(ji,jj  ,jk,Kbb) - puu(ji  ,jj-1,jk,Kbb) )
+               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pvv(ji,jj+1,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
             END DO
          END DO
@@ -318,5 +320,5 @@
             DO ji = 1, jpim1
                ! i-gradient of v at jj
-               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
+               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pvv(ji+1,jj,jk,Kbb) - pvv(ji  ,jj  ,jk,Kbb) )
             END DO
          END DO
@@ -391,6 +393,6 @@
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
-               ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
+               puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+               pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
             END DO
          END DO
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90	(revision 10928)
@@ -35,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )
+   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
       !!----------------------------------------------------------------------
       !!                     ***  ROUTINE dyn_ldf_lap  ***
@@ -45,10 +45,11 @@
       !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 
       !!
-      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
+      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
       INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s]
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -76,12 +77,12 @@
 !!gm open question here : e3f  at before or now ?    probably now...
 !!gm note that ahmf has already been multiplied by fmask
-               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
-                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
-                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
+               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
+                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
+                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
                !                                      ! ahm * div        (computed from 2 to jpi/jpj)
 !!gm note that ahmt has already been multiplied by tmask
-               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
-                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
-                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
+               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         &
+                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
+                  &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
             END DO  
          END DO  
@@ -89,10 +90,10 @@
          DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 &
-                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
+               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 &
+                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
                   &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
                   !
-               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 &
-                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
+               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                                 &
+                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
                   &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
             END DO
@@ -105,5 +106,5 @@
 
 
-   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
+   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dyn_ldf_blp  ***
@@ -116,9 +117,10 @@
       !!      It is computed by two successive calls to dyn_ldf_lap routine
       !!
-      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
+      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
       !!----------------------------------------------------------------------
       INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
+      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
       !
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
@@ -134,9 +136,9 @@
       zvlap(:,:,:) = 0._wp
       !
-      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
+      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
       !
       CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
       !
-      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
+      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
       !
    END SUBROUTINE dyn_ldf_blp
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg.F90	(revision 10928)
@@ -175,6 +175,6 @@
       ENDIF
       !                                      ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' spg  - Ua: ', mask1=umask, &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' spg  - Ua: ', mask1=umask, &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( ln_timing )   CALL timing_stop('dyn_spg')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90	(revision 10928)
@@ -181,6 +181,6 @@
       !
       !                       ! print sum trends (used for debugging)
-      IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' vor  - Ua: ', mask1=umask,               &
-         &                     tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor  - Ua: ', mask1=umask,               &
+         &                     tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( ln_timing )   CALL timing_stop('dyn_vor')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90	(revision 10928)
@@ -116,6 +116,6 @@
       ENDIF
       !                             ! Control print
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zad  - Ua: ', mask1=umask,   &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' zad  - Ua: ', mask1=umask,   &
+         &                       tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
       !
       IF( ln_timing )   CALL timing_stop('dyn_zad')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90	(revision 10928)
@@ -130,6 +130,6 @@
       IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
          DO jk = 1, jpkm1        ! remove barotropic velocities
-            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - ua_b(:,:) ) * umask(:,:,jk)
-            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - va_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk)
          END DO
          DO jj = 2, jpjm1        ! Add bottom/top stress due to barotropic component only
@@ -139,6 +139,6 @@
                ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
                ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
-               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua
-               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va
+               puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
+               pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
             END DO
          END DO
@@ -150,6 +150,6 @@
                   ze3ua =  ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa)
                   ze3va =  ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa)
-                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua
-                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va
+                  puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
+                  pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
                END DO
             END DO
@@ -494,6 +494,6 @@
       ENDIF
       !                                          ! print mean trends (used for debugging)
-      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf  - Ua: ', mask1=umask,               &
-         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
+      IF(ln_ctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf  - Ua: ', mask1=umask,               &
+         &                       tab3d_2=pvv(:,:,:,Kaa), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
          !
       IF( ln_timing )   CALL timing_stop('dyn_zdf')
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90	(revision 10928)
@@ -428,5 +428,5 @@
       !                            
       IF( ln_diurnal_only ) THEN                   ! diurnal only: a subset of the initialisation routines
-         CALL  istate_init( Nbb, Nnn )                ! ocean initial state (Dynamics and tracers)
+         CALL  istate_init( Nbb, Nnn, Naa )           ! ocean initial state (Dynamics and tracers)
          CALL     sbc_init( Nbb, Nnn )                ! Forcings : surface module
          CALL tra_qsr_init                            ! penetrative solar radiation qsr
@@ -440,5 +440,5 @@
       ENDIF
       
-                           CALL  istate_init( Nbb, Nnn ) ! ocean initial state (Dynamics and tracers)
+                           CALL  istate_init( Nbb, Nnn, Naa )    ! ocean initial state (Dynamics and tracers)
 
       !                                      ! external forcing 
@@ -464,10 +464,10 @@
 
       !                                      ! Dynamics
-      IF( lk_c1d       )   CALL dyn_dmp_init      ! internal momentum damping
-                           CALL dyn_adv_init      ! advection (vector or flux form)
-                           CALL dyn_vor_init      ! vorticity term including Coriolis
-                           CALL dyn_ldf_init      ! lateral mixing
-                           CALL dyn_hpg_init      ! horizontal gradient of Hydrostatic pressure
-                           CALL dyn_spg_init      ! surface pressure gradient
+      IF( lk_c1d       )   CALL dyn_dmp_init         ! internal momentum damping
+                           CALL dyn_adv_init         ! advection (vector or flux form)
+                           CALL dyn_vor_init         ! vorticity term including Coriolis
+                           CALL dyn_ldf_init         ! lateral mixing
+                           CALL dyn_hpg_init( Nnn )  ! horizontal gradient of Hydrostatic pressure
+                           CALL dyn_spg_init         ! surface pressure gradient
 
 #if defined key_top
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10928)
@@ -195,8 +195,8 @@
                          CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! advection (VF or FF)	==> RHS
                          CALL dyn_vor( kstp,      Nnn      , uu, vv, Nrhs )  ! vorticity           	==> RHS
-                         CALL dyn_ldf       ( kstp )  ! lateral mixing
+                         CALL dyn_ldf( kstp, Nbb, Nnn      , uu, vv, Nrhs )  ! lateral mixing
       IF( ln_zdfosm  )   CALL dyn_osm( kstp,      Nnn      , uu, vv, Nrhs )  ! OSMOSIS non-local velocity fluxes ==> RHS
-                         CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
-                         CALL dyn_spg       ( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
+                         CALL dyn_hpg( kstp,      Nnn,       uu, vv, Nrhs )  ! horizontal gradient of Hydrostatic pressure
+                         CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa )  ! surface pressure gradient
 
                                                       ! With split-explicit free surface, since now transports have been updated and ssha as well
Index: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10927)
+++ NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90	(revision 10928)
@@ -315,5 +315,5 @@
       IF( ln_ctl       )   CALL prt_ctl_init        ! Print control
 
-                           CALL  istate_init    ! ocean initial state (Dynamics and tracers)
+                           CALL  istate_init( Nnn, Naa )    ! ocean initial state (Dynamics and tracers)
 
                            CALL     sbc_init( Nbb, Nnn )    ! Forcings : surface module
@@ -524,5 +524,5 @@
    END SUBROUTINE nemo_set_cfctl
 
-   SUBROUTINE istate_init
+   SUBROUTINE istate_init( Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE istate_init  ***
@@ -530,11 +530,12 @@
       !! ** Purpose :   Initialization to zero of the dynamics and tracers.
       !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   Kmm, Kaa  ! ocean time level indices
       !
       !     now fields         !     after fields      !
-      un   (:,:,:)   = 0._wp   ;   ua(:,:,:) = 0._wp   !
-      vn   (:,:,:)   = 0._wp   ;   va(:,:,:) = 0._wp   !
-      wn   (:,:,:)   = 0._wp   !                       !
+      uu   (:,:,:,Kmm)   = 0._wp   ;   uu(:,:,:,Kaa) = 0._wp   !
+      vv   (:,:,:,Kmm)   = 0._wp   ;   vv(:,:,:,Kaa) = 0._wp   !
+      ww   (:,:,:)   = 0._wp   !                       !
       hdivn(:,:,:)   = 0._wp   !                       !
-      tsn  (:,:,:,:) = 0._wp   !                       !
+      ts  (:,:,:,:,Kmm) = 0._wp   !                       !
       !
       rhd  (:,:,:) = 0.e0
