Index: NEMO/branches/UKMO/NEMO_4.0-HEAD_DJC_backwards_port/src/OCE/DYN/dynhpg.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0-HEAD_DJC_backwards_port/src/OCE/DYN/dynhpg.F90	(revision 15034)
+++ NEMO/branches/UKMO/NEMO_4.0-HEAD_DJC_backwards_port/src/OCE/DYN/dynhpg.F90	(revision 15035)
@@ -71,4 +71,6 @@
    !
    INTEGER, PUBLIC ::   nhpg         !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
+   LOGICAL,  PUBLIC ::   ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt
+   REAL(wp), PUBLIC ::   aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt
 
    !! * Substitutions
@@ -147,5 +149,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_vN_hor, ln_hpg_djc_vN_vrt
       !!----------------------------------------------------------------------
       !
@@ -172,10 +175,5 @@
       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) )          &
+      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf.OR.ln_hpg_djc) )          &
          &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ',   &
          &                 '   the standard jacobian formulation hpg_sco    or '    ,   &
@@ -213,4 +211,18 @@
       ENDIF
       !                          
+      IF (ln_hpg_djc_vN_hor) THEN ! Von Neumann boundary condition
+        aco_bc_hor = 6.0_wp/5.0_wp
+        bco_bc_hor = 7.0_wp/15.0_wp
+      ELSE ! 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_vN_vrt) THEN ! Von Neumann boundary condition
+        aco_bc_vrt = 6.0_wp/5.0_wp
+        bco_bc_vrt = 7.0_wp/15.0_wp
+      ELSE ! Linear extrapolation
+        aco_bc_vrt = 3.0_wp/2.0_wp
+        bco_bc_vrt = 1.0_wp/2.0_wp
+      END IF
       IF ( .NOT. ln_isfcav ) THEN     !--- no ice shelf load
          riceload(:,:) = 0._wp
@@ -673,12 +685,19 @@
       !!
       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) ::   rhd_opt
+ 
+      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
       !!----------------------------------------------------------------------
@@ -735,154 +754,265 @@
       ! 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
+
+      IF( .NOT. ln_linssh ) THEN
+         rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option
+      ELSE
+         rhd_opt(:,:,:) = rhd(:,:,:)
+      END IF
       !----------------------------------------------------------------------------------------
-      !  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
+!!bug gm   Not a true bug, but... zdzz=e3w  for zdzx, zdzy verify what it is really
+!! zdzz, zdzx and zdzy changed to heights rather than depths; lower bounds of jj and ji changed from 2 to 1 
 
       DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi   
+               zdrhoz(ji,jj,jk) =   rhd_opt    (ji  ,jj  ,jk) - rhd_opt    (ji,jj,jk-1)
+               zdzz  (ji,jj,jk) = - gde3w_n(ji  ,jj  ,jk) + gde3w_n(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+
+      !-------------------------------------------------------------------------
+      ! 2. compute harmonic averages for vertical differences using eq. 5.18
+      !-------------------------------------------------------------------------
+      zep = 1.e-15
+
+!!bug  gm  zdrhoz not defined at level 1 and used (jk-1 with jk=2) ; this issue has now been addressed
+!!bug  gm  idem for zdrhox, zdrhoy et ji=jpi and jj=jpj; idem
+
+!! zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 
+!! DO loops broken up so that zdrho_k and zdz_k are calculated only for jk = 2 to jpk - 2 
+      zdrho_k(:,:,:) = 0._wp
+      zdz_k  (:,:,:) = 0._wp
+
+      DO jk = 2, jpk - 2
+         DO jj = 1, jpj
+            DO ji = 1, jpi   
+               cffw = 2._wp * zdrhoz(ji  ,jj  ,jk) * zdrhoz(ji,jj,jk+1)
+               IF( cffw > zep) THEN
+                  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 DO
+         END DO
+      END DO
+
+      !----------------------------------------------------------------------------------
+      ! 3. apply boundary conditions at top and bottom using 5.36-5.37
+      !----------------------------------------------------------------------------------
+
+! 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_opt    (:,:,2) - rhd_opt    (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2)
+      zdz_k  (:,:,1) = aco_bc_vrt * (-gde3w_n(:,:,2) + gde3w_n(:,:,1) ) - bco_bc_vrt * zdz_k  (:,:,2)
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( mbkt(ji,jj)>1 ) THEN
+               iktb = mbkt(ji,jj)
+               zdrho_k(ji,jj,iktb) = aco_bc_vrt * (     rhd_opt(ji,jj,iktb) -     rhd_opt(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)
+               zdz_k  (ji,jj,iktb) = aco_bc_vrt * (-gde3w_n(ji,jj,iktb) + gde3w_n(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k  (ji,jj,iktb-1) 
+            END IF
+         END DO
+      END DO
+
+      !--------------------------------------------------------------
+      ! 4. Compute side face integrals
+      !-------------------------------------------------------------
+
+!! sshn replaces e3w_n ; gde3w_n 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
+      !-------------------------------------------------------------
+      ! Concerns that zdrho_k might be oddly defined (just -1.5rho) for single celled columns are resolved by the fact that z_rho_k is defined explicity for the surface layer
+      DO jj = 2, jpj
+         DO ji = 2, jpi 
+            z_rho_k(ji,jj,1) =  grav * ( sshn(ji,jj) + gde3w_n(ji,jj,1) )                  & ! *** AY sshn included in djc but not in sco
+               &                     * (  rhd_opt(ji,jj,1)                                     &
+               &                     + 0.5_wp * (   rhd_opt    (ji,jj,2) - rhd_opt    (ji,jj,1) )  &
+               &                              * (   sshn   (ji,jj  ) + gde3w_n(ji,jj,1) )  &
+               &                              / ( - gde3w_n(ji,jj,2) + gde3w_n(ji,jj,1) )  )
+         END DO
+      END DO
+
+      !--------------------------------------------------------------
+      ! 4. b) Interior faces, compute and store
+      !-------------------------------------------------------------
+
+      DO jk = 2, jpkm1
+         DO jj = 2, jpj
+            DO ji = 2, jpi 
+               z_rho_k(ji,jj,jk) = zcoef0 * (   rhd_opt    (ji,jj,jk) + rhd_opt    (ji,jj,jk-1) )                                   &
+                  &                       * ( - gde3w_n(ji,jj,jk) + gde3w_n(ji,jj,jk-1) )                                   &
+                  &                       + z_grav_10 * (                                                                   &
+                  &     (   zdrho_k  (ji,jj,jk) - zdrho_k  (ji,jj,jk-1) )                                                   &
+                  &   * ( - gde3w_n(ji,jj,jk) + gde3w_n(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_opt    (ji,jj,jk) - rhd_opt    (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) )   &
+                  &                             )
+            END DO
+         END DO
+      END DO
+
+      !----------------------------------------------------------------------------------------
+      !  5. compute and store elementary horizontal differences in provisional arrays 
+      !----------------------------------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zdrhox(ji,jj,jk) =   rhd_opt    (ji+1,jj  ,jk) - rhd_opt    (ji,jj,jk  )
+               zdzx  (ji,jj,jk) = - gde3w_n(ji+1,jj  ,jk) + gde3w_n(ji,jj,jk  )
+               zdrhoy(ji,jj,jk) =   rhd_opt    (ji  ,jj+1,jk) - rhd_opt    (ji,jj,jk  )
+               zdzy  (ji,jj,jk) = - gde3w_n(ji  ,jj+1,jk) + gde3w_n(ji,jj,jk  )
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 
+
+      !-------------------------------------------------------------------------
+      ! 6. compute harmonic averages using eq. 5.18
+      !-------------------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 2, jpj
+            DO ji = fs_2, jpi   ! vector opt.
+
+               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
+                  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
+                  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
+                  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
+                  zdz_j(ji,jj,jk) = 0._wp
+               ENDIF
+
+            END DO
+         END DO
+      END DO
+      
+!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point      
+
+      !----------------------------------------------------------------------------------
+      ! 6B. apply boundary conditions at side boundaries using 5.36-5.37
+      !----------------------------------------------------------------------------------
+
+      DO jk = 1, jpkm1
+         DO jj = 2, jpj
+            DO ji = 2, jpi   ! vector opt.
+               zz_drho_i(ji,jj) = zdrho_i(ji,jj,jk)
+               zz_dz_i  (ji,jj) = zdz_i  (ji,jj,jk)
+               zz_drho_j(ji,jj) = zdrho_j(ji,jj,jk)
+               zz_dz_j  (ji,jj) = zdz_j  (ji,jj,jk)
+               ! 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. tmask(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_opt    (ji+1,jj,jk) - rhd_opt    (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 
+                     zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w_n(ji+1,jj,jk) + gde3w_n(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_opt    (ji,jj,jk) - rhd_opt    (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)  
+                     zz_dz_i  (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj,jk) + gde3w_n(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. tmask(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_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)
+                     zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj+1,jk) + gde3w_n(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_opt    (ji,jj,jk) - rhd_opt    (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 
+                     zz_dz_j  (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj,jk) + gde3w_n(ji,jj-1,jk) ) - bco_bc_hor * zdz_j  (ji,jj-1,jk)
+                  END IF
+               END IF
+
+            END DO
+         END DO
+         
+         DO jj = 2, jpj
+            DO ji = 2, jpi   ! vector opt.
+              zdrho_i(ji,jj,jk) = zz_drho_i(ji,jj)
+              zdz_i  (ji,jj,jk) = zz_dz_i  (ji,jj)
+              zdrho_j(ji,jj,jk) = zz_drho_j(ji,jj)
+              zdz_j  (ji,jj,jk) = zz_dz_j  (ji,jj)
+            END DO
+         END DO
+
+      END DO
+
+      !--------------------------------------------------------------
+      ! 7. Calculate integrals on side faces  
+      !-------------------------------------------------------------
+
+      DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             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)
-               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  )
-               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  )
-            END DO
-         END DO
-      END DO
-
-      !-------------------------------------------------------------------------
-      ! compute harmonic averages 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 jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               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  )
-
-               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) )
-               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) )
-               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) )
-               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) )
-               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 DO
-         END DO
-      END DO
-
-      !----------------------------------------------------------------------------------
-      ! apply boundary conditions at top and bottom 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)
-
+
+! 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_opt    (ji+1,jj,jk) + rhd_opt    (ji,jj,jk) )                                     &
+                   &                    * ( gde3w_n(ji+1,jj,jk) - gde3w_n(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_n(ji+1,jj,jk) + gde3w_n(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_opt    (ji+1,jj,jk) - rhd_opt    (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_opt    (ji,jj+1,jk) + rhd_opt    (ji,jj,jk) )                                     &
+                   &                    * ( gde3w_n(ji,jj+1,jk) - gde3w_n(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_n(ji,jj+1,jk) + gde3w_n(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_opt    (ji,jj+1,jk) - rhd_opt    (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) )  &
+                   &                                                 )
+               END IF
+
+            END DO
+         END DO
+      END DO
 
       !--------------------------------------------------------------
-      ! Upper half of top-most grid box, compute and store
+      ! 8. Integrate in the vertical   
       !-------------------------------------------------------------
-
-!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   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 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) )               &
-               &                   * (  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) )  )
-         END DO
-      END DO
-
-!!bug gm    : here also, simplification is possible
-!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
-
-      DO jk = 2, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-
-               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) )                                   &
-                  &            - 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) ) )   &
-                  &   - ( 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_n(ji+1,jj,jk) - gde3w_n(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) ) )   &
-                  &   - ( 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_n(ji,jj+1,jk) - gde3w_n(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) ) )   &
-                  &   - ( 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 DO
-         END DO
-      END DO
-      CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1., rho_i, 'U', 1., rho_j, 'V', 1. )
 
       ! ---------------
@@ -891,6 +1021,6 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, fs_jpim1   ! vector opt.
-            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)
@@ -911,9 +1041,9 @@
                ! 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)
+                  &           + (  ( 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)                                &
-                  &           + (  ( 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)
+                  &           + (  ( 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)
