Index: NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynhpg.F90	(revision 13360)
+++ NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/OCE/DYN/dynhpg.F90	(revision 13364)
@@ -150,5 +150,4 @@
       !!
       INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
-      REAL(wp) ::   znad
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zts_top, zrhd   ! hypothesys on isf density
       REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::  zrhdtop_isf    ! density at bottom of ISF
@@ -245,5 +244,5 @@
       INTEGER  ::   ji, jj, jk       ! dummy loop indices
       REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
+      REAL(wp), DIMENSION(jpi,jpj) ::  zhpi, zhpj
       !!----------------------------------------------------------------------
       !
@@ -253,33 +252,28 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
       ENDIF
-
-      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
-
-      ! Surface value
-      DO_2D( 0, 0, 0, 0 )
+      !
+      zcoef0 = - grav * 0.5_wp            ! Local constant initialization
+      !
+      DO_2D( 0, 0, 0, 0 )                 ! Surface value
          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
-         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)
+         !                                   ! hydrostatic pressure gradient
+         zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
+         zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
+         !                                   ! add to the general momentum trend
+         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj)
+         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj)
       END_2D
-
-      !
-      ! interior value (2=<jk=<jpkm1)
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
+      !
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! interior value (2=<jk=<jpkm1)
          zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
-         ! hydrostatic pressure gradient
-         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
-            &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
-            &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
-
-         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
-            &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
-            &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
-         ! add to the general momentum trend
-         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)
+         !                                   ! hydrostatic pressure gradient
+         zhpi(ji,jj) = zhpi(ji,jj) + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )  &
+            &                                  - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
+
+         zhpj(ji,jj) = zhpj(ji,jj) + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )  &
+            &                                  - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
+         !                                   ! add to the general momentum trend
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj)
       END_3D
       !
@@ -390,7 +384,7 @@
       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
-      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
-      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
+      INTEGER  ::   ji, jj, jk, jii, jjj           ! dummy loop indices
+      REAL(wp) ::   zcoef0, zuap, zvap, ztmp       ! local scalars
+      LOGICAL  ::   ll_tmp1, ll_tmp2               ! local logical variables
       REAL(wp), DIMENSION(jpi,jpj,jpk)      ::   zhpi, zhpj
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
@@ -406,7 +400,4 @@
       !
       zcoef0 = - grav * 0.5_wp
-      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
-      ELSE                    ;   znad = 1._wp         ! Variable volume: density
-      ENDIF
       !
       IF( ln_wd_il ) THEN
@@ -450,20 +441,17 @@
         CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
-
-      ! Surface value
-      DO_2D( 0, 0, 0, 0 )
-         ! hydrostatic pressure gradient along s-surfaces
-         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 )   &
+      !
+      DO_2D( 0, 0, 0, 0 )              ! Surface value
+         !                                   ! hydrostatic pressure gradient along s-surfaces
+         zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj)                      &
+            &          * (  e3w(ji+1,jj  ,1,Kmm) * rhd(ji+1,jj  ,1)  &
+            &             - e3w(ji  ,jj  ,1,Kmm) * rhd(ji  ,jj  ,1)  )
+         zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj)                      &
+            &          * (  e3w(ji  ,jj+1,1,Kmm) * rhd(ji  ,jj+1,1)  &
+            &             - e3w(ji  ,jj  ,1,Kmm) * rhd(ji  ,jj  ,1)  )
+         !                                   ! s-coordinate pressure gradient correction
+         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) )   &
             &           * ( 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 )   &
+         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) )   &
             &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
          !
@@ -474,23 +462,21 @@
             zvap = zvap * zcpy(ji,jj)
          ENDIF
-         !
-         ! add to the general momentum trend
+         !                                   ! add to the general momentum trend
          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_2D
-
-      ! interior value (2=<jk=<jpkm1)
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         ! hydrostatic pressure gradient along s-surfaces
-         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
-            &           * (  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(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 )   &
+      !
+      DO_3D( 0, 0, 0, 0, 2, jpkm1 )    ! interior value (2=<jk=<jpkm1)
+         !                                   ! hydrostatic pressure gradient along s-surfaces
+         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)                         &
+            &           * (  e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )  &
+            &              - e3w(ji  ,jj,jk,Kmm) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  )
+         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)                         &
+            &           * (  e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )  &
+            &              - e3w(ji,jj  ,jk,Kmm) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  )
+         !                                   ! s-coordinate pressure gradient correction
+         zuap = -zcoef0 * ( rhd  (ji+1,jj  ,jk) + rhd  (ji,jj,jk) ) &
             &           * ( 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 )   &
+         zvap = -zcoef0 * ( rhd  (ji  ,jj+1,jk) + rhd  (ji,jj,jk) ) &
             &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj)
          !
@@ -535,6 +521,8 @@
       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
-      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
+      INTEGER  ::   ikt ,  ikti1,  iktj1   ! local integer
+      REAL(wp) ::   ze3w, ze3wi1, ze3wj1   ! local scalars
+      REAL(wp) ::   zcoef0, zuap, zvap     !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpk ) ::  zhpi, zhpj
       REAL(wp), DIMENSION(jpi,jpj,jpts) ::  zts_top
@@ -543,6 +531,4 @@
       !
       zcoef0 = - grav * 0.5_wp   ! Local constant initialization
-      !
-      znad=1._wp                 ! To use density and not density anomaly
       !
       !                          ! iniitialised to 0. zhpi zhpi 
@@ -560,54 +546,49 @@
       CALL eos( zts_top, risfdep, zrhdtop_oce )
 
-!==================================================================================     
-!===== Compute surface value ===================================================== 
-!==================================================================================
+      !                     !===========================!
+      !                     !=====  surface value  =====!
+      !                     !===========================!
       DO_2D( 0, 0, 0, 0 )
-         ikt    = mikt(ji,jj)
-         iktp1i = mikt(ji+1,jj)
-         iktp1j = mikt(ji,jj+1)
-         ! 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(ji+1,jj,iktp1i,Kmm)                                    &
-            &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
-            &                                  - 0.5_wp * e3w(ji,jj,ikt,Kmm)                                         &
-            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-            &                                  + ( risfload(ji+1,jj) - risfload(ji,jj))                            ) 
-         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(ji,jj,ikt,Kmm)                                         & 
-            &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
-            &                                  + ( risfload(ji,jj+1) - risfload(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 )   &
+         ikt   = mikt(ji  ,jj  )   ;   ze3w   = e3w(ji  ,jj  ,ikt  ,Kmm)
+         ikti1 = mikt(ji+1,jj  )   ;   ze3wi1 = e3w(ji+1,jj  ,ikti1,Kmm)
+         iktj1 = mikt(ji  ,jj+1)   ;   ze3wj1 = e3w(ji  ,jj+1,iktj1,Kmm)
+         !                          ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
+         !                          ! we assume ISF is in isostatic equilibrium
+         zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * (   risfload(ji+1,jj) - risfload(ji,jj)  &
+            &                                       + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhdtop_oce(ji+1,jj) )     &
+            &                                                  - ze3w   * ( rhd(ji  ,jj,ikt  ) + zrhdtop_oce(ji  ,jj) ) )   )
+         zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * (   risfload(ji,jj+1) - risfload(ji,jj)  &
+            &                                       + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhdtop_oce(ji,jj+1) )      &
+            &                                                  - ze3w   * ( rhd(ji,jj  ,ikt  ) + zrhdtop_oce(ji,jj  ) ) )   ) 
+         !                          ! s-coordinate pressure gradient correction (=0 if z coordinate)
+         zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) )   &
             &           * ( 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 )   &
+         zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) )   &
             &           * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
-         ! add to the general momentum trend
+         !                          ! add to the general momentum trend
          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_2D
-!==================================================================================     
-!===== Compute interior value ===================================================== 
-!==================================================================================
-      ! interior value (2=<jk=<jpkm1)
+      !   
+      !                     !=============================!
+      !                     !=====  interior values  =====!
+      !                     !=============================!
       DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         ! hydrostatic pressure gradient along s-surfaces
+         ze3w   = e3w(ji  ,jj  ,jk,Kmm)
+         ze3wi1 = e3w(ji+1,jj  ,jk,Kmm)
+         ze3wj1 = e3w(ji  ,jj+1,jk,Kmm)
+         !                          ! hydrostatic pressure gradient along s-surfaces
          zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-            &           * (  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)   )
+            &           * (  ze3wi1 * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) * wmask(ji+1,jj,jk)   &
+            &              - ze3w   * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) ) * wmask(ji  ,jj,jk)   )
          zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-            &           * (  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 )   &
+            &           * (  ze3wj1 * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) * wmask(ji,jj+1,jk)   &
+            &              - ze3w   * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) ) * wmask(ji,jj  ,jk)   )
+         !                          ! s-coordinate pressure gradient correction
+         zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) )   &
             &           * ( 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 )   &
+         zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) )   &
             &           * ( gde3w(ji  ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj)
-         ! add to the general momentum trend
+         !                          ! add to the general momentum trend
          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)
@@ -634,4 +615,5 @@
       REAL(wp) ::   z1_12, cffv, cffy   !    "         "
       LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zpgu, zpgv   ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw
@@ -826,5 +808,13 @@
       END_3D
       CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp )
-
+      !
+      ! ---------------
+      !  Surface pressure gradient to be removed
+      ! ---------------
+      DO_2D( 0, 0, 0, 0 )
+         zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
+         zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
+      END_2D
+      !
       ! ---------------
       !  Surface value
@@ -838,6 +828,6 @@
          ENDIF
          ! add to the general momentum trend
-         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)
+         puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) - zpgu(ji,jj)
+         pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) - zpgv(ji,jj)
       END_2D
 
@@ -858,6 +848,6 @@
          ENDIF
          ! add to the general momentum trend
-         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)
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) - zpgu(ji,jj)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) - zpgv(ji,jj)
       END_3D
       !
@@ -892,7 +882,8 @@
       REAL(wp) :: zrhdt1
       REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zpgu, zpgv   ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zsshu_n, zsshv_n
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
-      REAL(wp), DIMENSION(jpi,jpj)   ::   zsshu_n, zsshv_n
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
@@ -907,6 +898,14 @@
       zcoef0 = - grav
       znad = 1._wp
-      IF( ln_linssh )   znad = 0._wp
-
+      IF( ln_linssh )   znad = 1._wp
+      !
+      ! ---------------
+      !  Surface pressure gradient to be removed
+      ! ---------------
+      DO_2D( 0, 0, 0, 0 )
+         zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
+         zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
+      END_2D
+      !
       IF( ln_wd_il ) THEN
          ALLOCATE( zcpx(jpi,jpj) , zcpy(jpi,jpj) )
@@ -974,5 +973,5 @@
       ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
       DO_2D( 1, 1, 1, 1 )
-         zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm) * znad
+         zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm)
       END_2D
 
@@ -1025,6 +1024,6 @@
 
       DO_2D( 0, 0, 0, 0 )
-       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)
+       zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 
+       zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )
       END_2D
 
@@ -1108,5 +1107,5 @@
             zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
          ENDIF
-         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
+         puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 
       ENDIF
 
@@ -1168,5 +1167,5 @@
          ENDIF
 
-         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
+         pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)
       ENDIF
          !
