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 10740)
+++ 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
                !
