Index: NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90	(revision 14769)
+++ NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90	(revision 14789)
@@ -17,4 +17,5 @@
    !!                 !           (A. Coward) suppression of hel, wdj and rot options
    !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
+   !!            4.2  !  2020-12  (M. Bell, A. Young) hpg_djc: revised djc scheme
    !!----------------------------------------------------------------------
 
@@ -72,5 +73,8 @@
    INTEGER, PARAMETER ::   np_isf    =  5   ! s-coordinate similar to sco modify for isf
    !
-   INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
+   INTEGER, PUBLIC  ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
+   !
+   LOGICAL          ::   ln_hpg_djc_vnh, ln_hpg_djc_vnv                 ! flag to specify hpg_djc boundary condition type
+   REAL(wp), PUBLIC ::   aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt !: coefficients for hpg_djc hor and vert boundary conditions
 
    !! * Substitutions
@@ -150,5 +154,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
@@ -156,5 +159,6 @@
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
-         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf
+         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf,     &
+         &                 ln_hpg_djc_vnh, ln_hpg_djc_vnv
       !!----------------------------------------------------------------------
       !
@@ -179,19 +183,15 @@
       ENDIF
       !
-      IF( ln_hpg_djc )   &
-         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method',   &
-         &                 '   currently disabled (bugs under investigation).'        ,   &
-         &                 '   Please select either  ln_hpg_sco or  ln_hpg_prj instead' )
-         !
-      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )          &
-         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   &
-         &                 '   the standard jacobian formulation hpg_sco    or '    ,   &
-         &                 '   the pressure jacobian formulation hpg_prj'             )
-         !
+      IF( .NOT.ln_linssh .AND. (ln_hpg_zco.OR.ln_hpg_zps) )   &
+         &   CALL ctl_stop( 'dyn_hpg_init : non-linear free surface incompatible with hpg_zco or hpg_zps' )
+      !
+      IF( (.NOT.ln_hpg_isf .AND. ln_isfcav) .OR. (ln_hpg_isf .AND. .NOT.ln_isfcav) )                  &
+         &   CALL ctl_stop( 'dyn_hpg_init : ln_hpg_isf=T requires ln_isfcav=T and vice versa' )  
+      !
+#if defined key_qco
       IF( ln_hpg_isf ) THEN
-         IF( .NOT. ln_isfcav )   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
-       ELSE
-         IF(       ln_isfcav )   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
+         CALL ctl_stop( 'dyn_hpg_init : key_qco and ln_hpg_isf not yet compatible' )
       ENDIF
+#endif
       !
       !                               ! Set nhpg from ln_hpg_... flags & consistency check
@@ -220,4 +220,25 @@
       ENDIF
       !                          
+      IF ( ln_hpg_djc ) THEN
+         IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition
+           IF(lwp) WRITE(numout,*) '           horizontal bc: von Neumann '
+           aco_bc_hor = 6.0_wp/5.0_wp
+           bco_bc_hor = 7.0_wp/15.0_wp
+         ELSE ! Linear extrapolation
+           IF(lwp) WRITE(numout,*) '           horizontal bc: linear extrapolation'
+           aco_bc_hor = 3.0_wp/2.0_wp
+           bco_bc_hor = 1.0_wp/2.0_wp
+         END IF
+         IF (ln_hpg_djc_vnv) THEN ! Von Neumann boundary condition
+           IF(lwp) WRITE(numout,*) '           vertical bc: von Neumann '
+           aco_bc_vrt = 6.0_wp/5.0_wp
+           bco_bc_vrt = 7.0_wp/15.0_wp
+         ELSE ! Linear extrapolation
+           IF(lwp) WRITE(numout,*) '           vertical bc: linear extrapolation'
+           aco_bc_vrt = 3.0_wp/2.0_wp
+           bco_bc_vrt = 1.0_wp/2.0_wp
+         END IF
+      END IF
+      !
    END SUBROUTINE dyn_hpg_init
 
@@ -245,5 +266,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 +274,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
       !
@@ -302,6 +318,7 @@
       INTEGER  ::   iku, ikv                         ! temporary integers
       REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
-      REAL(wp), DIMENSION(jpi,jpj) :: zgtsu, zgtsv, zgru, zgrv
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj
+      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zgtsu, zgtsv
+      REAL(wp), DIMENSION(jpi,jpj)     :: zgru, zgrv
       !!----------------------------------------------------------------------
       !
@@ -390,7 +407,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
@@ -402,11 +419,8 @@
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OCE original scheme used'
       ENDIF
       !
       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
@@ -448,22 +462,19 @@
           END IF
         END_2D
-        CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
+        CALL lbc_lnk( '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 +485,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 +544,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 +554,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 +569,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)
@@ -630,12 +634,18 @@
       !!
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      INTEGER  ::   iktb, iktt          ! jk indices at tracer points for top and bottom points 
       REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
-      REAL(wp) ::   z1_10, cffu, cffx   !    "         "
-      REAL(wp) ::   z1_12, cffv, cffy   !    "         "
+      REAL(wp) ::   z_grav_10, z1_12
+      REAL(wp) ::   cffu, cffx          !    "         "
+      REAL(wp) ::   cffv, cffy          !    "         "
       LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dzx, dzy, dzz, dzu, dzv, dzw
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   drhox, drhoy, drhoz, drhou, drhov, drhow
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   rho_i, rho_j, rho_k
+ 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdzx, zdzy, zdzz                          ! Primitive grid differences ('delta_xyz')
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdz_i, zdz_j, zdz_k                       ! Harmonic average of primitive grid differences ('d_xyz')
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrhox, zdrhoy, zdrhoz                    ! Primitive rho differences ('delta_rho')
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdrho_i, zdrho_j, zdrho_k                 ! Harmonic average of primitive rho differences ('d_rho')
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z_rho_i, z_rho_j, z_rho_k                 ! Face intergrals
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j    ! temporary arrays 
       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
       !!----------------------------------------------------------------------
@@ -679,5 +689,5 @@
           END IF
         END_2D
-        CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
+        CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       END IF
 
@@ -690,147 +700,219 @@
       ! Local constant initialization
       zcoef0 = - grav * 0.5_wp
-      z1_10  = 1._wp / 10._wp
-      z1_12  = 1._wp / 12._wp
+      z_grav_10  = grav / 10._wp
+      z1_12  = 1.0_wp / 12._wp
 
       !----------------------------------------------------------------------------------------
-      !  compute and store in provisional arrays elementary vertical and horizontal differences
+      !  1. compute and store elementary vertical differences in provisional arrays 
       !----------------------------------------------------------------------------------------
 
-!!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really
-
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (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(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(ji  ,jj+1,jk) - gde3w(ji,jj,jk  )
+!!bug gm   Not a true bug, but... zdzz=e3w  for zdzx, zdzy verify what it is really
+
+      DO_3D( 1, 1, 1, 1, 2, jpkm1 )  
+         zdrhoz(ji,jj,jk) =   rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1)
+         zdzz  (ji,jj,jk) = - gde3w(ji  ,jj  ,jk) + gde3w(ji,jj,jk-1)
       END_3D
 
       !-------------------------------------------------------------------------
-      ! compute harmonic averages using eq. 5.18
+      ! 2. compute harmonic averages for vertical differences using eq. 5.18
       !-------------------------------------------------------------------------
       zep = 1.e-15
 
-!!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
-!!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
-
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1)
-
-         cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  )
-         cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  )
-
-         cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  )
-         cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  )
-
+!! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 
+      zdrho_k(:,:,:) = 0._wp
+      zdz_k  (:,:,:) = 0._wp
+
+      DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 
+         cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1)
          IF( cffw > zep) THEN
-            drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   &
-               &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
+            zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) )
+         ENDIF
+         zdz_k(ji,jj,jk) = 2._wp *   zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1)   &
+            &                  / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) )
+      END_3D
+
+      !----------------------------------------------------------------------------------
+      ! 3. apply boundary conditions at top and bottom using 5.36-5.37
+      !----------------------------------------------------------------------------------
+
+! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition
+      zdrho_k(:,:,1) = aco_bc_vrt * ( rhd    (:,:,2) - rhd    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2)
+      zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2)
+
+      DO_2D( 1, 1, 1, 1 )
+         IF ( mbkt(ji,jj)>1 ) THEN
+            iktb = mbkt(ji,jj)
+            zdrho_k(ji,jj,iktb) = aco_bc_vrt * (     rhd(ji,jj,iktb) -     rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)
+            zdz_k  (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k  (ji,jj,iktb-1) 
+         END IF
+      END_2D
+
+      !--------------------------------------------------------------
+      ! 4. Compute side face integrals
+      !-------------------------------------------------------------
+
+!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights  
+!! rho_k stores grav * FX / rho_0  
+
+      !--------------------------------------------------------------
+      ! 4. a) Upper half of top-most grid box, compute and store
+      !-------------------------------------------------------------
+! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1,Kmm)
+      DO_2D( 0, 1, 0, 1)
+         z_rho_k(ji,jj,1) =  grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) )                        & 
+            &                     * (  rhd(ji,jj,1)                                        &
+            &                     + 0.5_wp * (   rhd    (ji,jj,2) - rhd    (ji,jj,1) ) &
+            &                              * (   ssh   (ji,jj,Kmm) + gde3w(ji,jj,1) )          &
+            &                              / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) )  )
+      END_2D
+
+      !--------------------------------------------------------------
+      ! 4. b) Interior faces, compute and store
+      !-------------------------------------------------------------
+
+      DO_3D( 0, 1, 0, 1, 2, jpkm1 )
+         z_rho_k(ji,jj,jk) = zcoef0 * (   rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   &
+            &                       * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) )                                               &
+            &                       + z_grav_10 * (                                                                           &
+            &     (   zdrho_k  (ji,jj,jk) - zdrho_k  (ji,jj,jk-1) )                                                           &
+            &   * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k  (ji,jj,jk) + zdz_k  (ji,jj,jk-1) ) )             &
+            &   - ( zdz_k    (ji,jj,jk) - zdz_k    (ji,jj,jk-1) )                                                             &
+            &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   &
+            &                             )
+      END_3D
+
+      !----------------------------------------------------------------------------------------
+      !  5. compute and store elementary horizontal differences in provisional arrays 
+      !----------------------------------------------------------------------------------------
+
+      DO_3D( 1, 0, 1, 0, 1, jpkm1 )
+         zdrhox(ji,jj,jk) =   rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
+         zdzx  (ji,jj,jk) = - gde3w(ji+1,jj  ,jk) + gde3w(ji,jj,jk  )
+         zdrhoy(ji,jj,jk) =   rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  )
+         zdzy  (ji,jj,jk) = - gde3w(ji  ,jj+1,jk) + gde3w(ji,jj,jk  )
+      END_3D
+
+      CALL lbc_lnk( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 
+
+      !-------------------------------------------------------------------------
+      ! 6. compute harmonic averages using eq. 5.18
+      !-------------------------------------------------------------------------
+
+      DO_3D( 0, 1, 0, 1, 1, jpkm1 )
+         cffu = 2._wp * zdrhox(ji-1,jj  ,jk) * zdrhox(ji,jj,jk  )
+         IF( cffu > zep ) THEN
+            zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) )
          ELSE
-            drhow(ji,jj,jk) = 0._wp
-         ENDIF
-
-         dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   &
-            &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
-
-         IF( cffu > zep ) THEN
-            drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   &
-               &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
+            zdrho_i(ji,jj,jk ) = 0._wp
+         ENDIF
+
+         cffx = 2._wp * zdzx  (ji-1,jj  ,jk) * zdzx  (ji,jj,jk  )
+         IF( cffx > zep ) THEN
+            zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) )
          ELSE
-            drhou(ji,jj,jk ) = 0._wp
-         ENDIF
-
-         IF( cffx > zep ) THEN
-            dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   &
-               &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
+            zdz_i(ji,jj,jk) = 0._wp
+         ENDIF
+
+         cffv = 2._wp * zdrhoy(ji  ,jj-1,jk) * zdrhoy(ji,jj,jk  )
+         IF( cffv > zep ) THEN
+            zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) )
          ELSE
-            dzu(ji,jj,jk) = 0._wp
-         ENDIF
-
-         IF( cffv > zep ) THEN
-            drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   &
-               &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
+            zdrho_j(ji,jj,jk) = 0._wp
+         ENDIF
+
+         cffy = 2._wp * zdzy  (ji  ,jj-1,jk) * zdzy  (ji,jj,jk  )
+         IF( cffy > zep ) THEN
+            zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) )
          ELSE
-            drhov(ji,jj,jk) = 0._wp
-         ENDIF
-
-         IF( cffy > zep ) THEN
-            dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   &
-               &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
-         ELSE
-            dzv(ji,jj,jk) = 0._wp
-         ENDIF
-
-      END_3D
+            zdz_j(ji,jj,jk) = 0._wp
+         ENDIF
+      END_3D
+      
+!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point      
 
       !----------------------------------------------------------------------------------
-      ! apply boundary conditions at top and bottom using 5.36-5.37
+      ! 6B. apply boundary conditions at side boundaries using 5.36-5.37
       !----------------------------------------------------------------------------------
-      drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  )
-      drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  )
-      drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  )
-
-      drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
-      drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
-      drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
-
+
+      DO jk = 1, jpkm1
+         zz_drho_i(:,:) = zdrho_i(:,:,jk)
+         zz_dz_i  (:,:) = zdz_i  (:,:,jk)
+         zz_drho_j(:,:) = zdrho_j(:,:,jk)
+         zz_dz_j  (:,:) = zdz_j  (:,:,jk)
+         DO_2D( 0, 1, 0, 1)
+            ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
+            IF (ji < jpi) THEN
+               IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp)  THEN  
+                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 
+                  zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i  (ji+1,jj,jk)
+               END IF
+            END IF
+            ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
+            IF (ji > 2) THEN
+               IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN
+                  zz_drho_i(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)  
+                  zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i  (ji-1,jj,jk)
+               END IF
+            END IF
+            ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
+            IF (jj < jpj) THEN
+               IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp)  THEN
+                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
+                  zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j  (ji,jj+1,jk)
+               END IF
+            END IF 
+            ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
+            IF (jj > 2) THEN
+               IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 
+                  zz_drho_j(ji,jj) = aco_bc_hor * ( rhd    (ji,jj,jk) - rhd    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 
+                  zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk)
+               END IF
+            END IF
+         END_2D
+         zdrho_i(:,:,jk) = zz_drho_i(:,:)
+         zdz_i  (:,:,jk) = zz_dz_i  (:,:)
+         zdrho_j(:,:,jk) = zz_drho_j(:,:)
+         zdz_j  (:,:,jk) = zz_dz_j  (:,:)
+      END DO
 
       !--------------------------------------------------------------
-      ! Upper half of top-most grid box, compute and store
+      ! 7. Calculate integrals on side faces  
       !-------------------------------------------------------------
 
-!!bug gm   :  e3w-gde3w(:,:,:) = 0.5*e3w  ....  and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) ....   to be verified
-!          true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be
-
-      DO_2D( 0, 0, 0, 0 )
-         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  (ji,jj,1,Kmm) - gde3w(ji,jj,1) )  &
-            &                              / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) )  )
-      END_2D
-
-!!bug gm    : here also, simplification is possible
-!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
-
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-
-         rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (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(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) ) )   &
-            &                             )
-
-         rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (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(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) ) )   &
-            &                            )
-
-         rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (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(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) ) )   &
-            &                            )
-
-      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 )
-
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height)
+         z_rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                       &
+             &                    * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )                                    
+         IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
+            z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * (                                                               &
+             &     (   zdrho_i  (ji+1,jj,jk) - zdrho_i  (ji,jj,jk) )                                                            &
+             &   * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i  (ji+1,jj,jk) + zdz_i  (ji,jj,jk) ) )              &
+             &   - (   zdz_i    (ji+1,jj,jk) - zdz_i    (ji,jj,jk) )                                                            &
+             &   * (   rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) )  &
+             &                                               )
+         END IF
+  
+         z_rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                       &
+             &                    * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )                                  
+         IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
+            z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * (                                                               &
+             &     (   zdrho_j  (ji,jj+1,jk) - zdrho_j  (ji,jj,jk) )                                                            &
+             &   * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j  (ji,jj+1,jk) + zdz_j  (ji,jj,jk) ) )              &
+             &   - (   zdz_j    (ji,jj+1,jk) - zdz_j    (ji,jj,jk) )                                                            &
+             &   * (   rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  &
+             &                                                 )
+         END IF
+      END_3D
+
+      !--------------------------------------------------------------
+      ! 8. Integrate in the vertical   
+      !-------------------------------------------------------------
+      !
       ! ---------------
       !  Surface value
       ! ---------------
       DO_2D( 0, 0, 0, 0 )
-         zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
-         zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
+         zhpi(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj  ,1) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
+         zhpj(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji  ,jj+1,1) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
          IF( ln_wd_il ) THEN
            zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
@@ -847,10 +929,10 @@
       DO_3D( 0, 0, 0, 0, 2, jpkm1 )
          ! hydrostatic pressure gradient along s-surfaces
-         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
-            &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    &
-            &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj)
-         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                &
-            &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
-            &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj)
+         zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                                     &
+            &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk  ) )                     &
+            &              - ( z_rho_i(ji,jj,jk) - z_rho_i(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
+         zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                                     &
+            &           + (  ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk  ) )                     &
+            &               -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
          IF( ln_wd_il ) THEN
            zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
@@ -892,7 +974,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 +990,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) )
@@ -952,5 +1043,5 @@
             ENDIF
          END_2D
-         CALL lbc_lnk_multi( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
+         CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )
       ENDIF
 
@@ -974,5 +1065,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
 
@@ -1022,9 +1113,9 @@
       END_2D
 
-      CALL lbc_lnk_multi ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )
+      CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp )
 
       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 +1199,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 +1259,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
          !
