Index: branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 8111)
+++ branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 8112)
@@ -84,4 +84,5 @@
          <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" />
          <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" />
+         <field id="rho"          long_name="in sity density"                                                           unit="kg/m3"  grid_ref="grid_T_3D" />
 
          <!-- Energy - horizontal divergence -->
@@ -415,4 +416,5 @@
          <field id="uocet"        long_name="ocean transport along i-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_U_3D" />
          <field id="uoces"        long_name="ocean transport along i-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_U_3D" />
+         <field id="hpgi"         long_name="horizontal pressure gradient along i-axis"                                                          unit="N/kg"        grid_ref="grid_U_3D" />
 
          <!-- variables available with MLE -->
@@ -459,4 +461,5 @@
          <field id="vocet"        long_name="ocean transport along j-axis times temperature (CRS)"                                               unit="degree_C*m/s"   grid_ref="grid_V_3D" />
          <field id="voces"        long_name="ocean transport along j-axis times salinity (CRS)"                                                  unit="0.001*m/s"   grid_ref="grid_V_3D" />
+         <field id="hpgj"         long_name="horizontal pressure gradient along j-axis"                                                          unit="N/kg"        grid_ref="grid_V_3D" />
 
          <!-- variables available with MLE -->
@@ -492,4 +495,6 @@
         <field id="woce"         long_name="ocean vertical velocity"              standard_name="upward_sea_water_velocity"   unit="m/s"  />
         <field id="wocetr_eff"   long_name="effective ocean vertical transport"                                               unit="m3/s" />
+        <field id="hpgdepth"     long_name="depth from HPG module" unit="m"  />
+        <field id="pressure"     long_name="in situ pressure"      unit="N/m2"  />
 
         <!-- woce_eiv: available with key_traldf_eiv and key_diaeiv -->
Index: branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 8111)
+++ branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90	(revision 8112)
@@ -1133,11 +1133,11 @@
          e3vw_0(:,:,jk) = e3w_1d(jk)
       END DO
-      DO jk = 1,jpk                         ! Computed as the minimum of neighbooring scale factors
+      DO jk = 1,jpk                         ! Computed as the average of neighbouring scale factors
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
-               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
-               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
-               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
+               e3u_0 (ji,jj,jk) = 0.5 * ( e3t_0(ji,jj,jk) + e3t_0(ji+1,jj,jk) )
+               e3v_0 (ji,jj,jk) = 0.5 * ( e3t_0(ji,jj,jk) + e3t_0(ji,jj+1,jk) )
+               e3uw_0(ji,jj,jk) = 0.5 * ( e3w_0(ji,jj,jk) + e3w_0(ji+1,jj,jk) )
+               e3vw_0(ji,jj,jk) = 0.5 * ( e3w_0(ji,jj,jk) + e3w_0(ji,jj+1,jk) )
             END DO
          END DO
@@ -1173,8 +1173,8 @@
          e3f_0(:,:,jk) = e3t_1d(jk)
       END DO
-      DO jk = 1, jpk                        ! Computed as the minimum of neighbooring V-scale factors
+      DO jk = 1, jpk                        ! Computed as the average of neighbouring V-scale factors
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
+               e3f_0(ji,jj,jk) = 0.5 * ( e3v_0(ji,jj,jk) + e3v_0(ji+1,jj,jk) )
             END DO
          END DO
Index: branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8111)
+++ branches/UKMO/dev_r5518_GO6_package_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8112)
@@ -44,4 +44,5 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+   USE iom
 
    IMPLICIT NONE
@@ -379,9 +380,13 @@
       !!
       INTEGER  ::   ji, jj, jk                 ! dummy loop indices
-      REAL(wp) ::   zcoef0, zuap, zvap, znad   ! temporary scalars
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
+      REAL(wp) ::   zcoef0, zuap, zvap, znad, zvnum   ! temporary scalars
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zunum, zudenom    ! temporary numerators for u and v velocities (for debugging) 
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpb, zpt   ! pressure at top and bottom of cells
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zzb, zzt   ! geopotential height (gz) at top and bottom of cells
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zhpi, zhpj, zp3d, zz3d
+      !!----------------------------------------------------------------------
+      !
+      CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zp3d, zz3d )
+      CALL wrk_alloc( jpi,jpj, zunum, zudenom, zpb, zpt, zzb, zzt )
       !
       IF( kt == nit000 ) THEN
@@ -398,47 +403,72 @@
       ENDIF
 
-      ! Surface value
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            ! hydrostatic pressure gradient along s-surfaces
-            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )   &
-               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
-            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )   &
-               &                                  - fse3w(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) ) )
-            ! s-coordinate pressure gradient correction
-            zuap = -zcoef0 * ( rhd   (ji+1,jj,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
-               &           * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
-            zvap = -zcoef0 * ( rhd   (ji,jj+1,1) + rhd   (ji,jj,1) + 2._wp * znad )   &
-               &           * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
-            ! add to the general momentum trend
-            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
-            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
-         END DO
-      END DO
-
-      ! interior value (2=<jk=<jpkm1)
-      DO jk = 2, jpkm1
+! set the pressure values to zero at the top of the model (the atmospheric surface pressure is included in the barotropic step ?)
+      zpb(:,:) = 0._wp
+      zzb(:,:) = grav * sshn(:,:)     ! surface height had been incorrectly set to zero; corrected on 24 Apr 2017
+
+      if(lwp) write(numout,*) 'hpg_sco: grav is ',grav
+
+      DO jk = 1, jpkm1
+
+         zpt(:,:) = zpb(:,:)
+         zzt(:,:) = zzb(:,:) 
+
+         zp3d(:,:,jk) = zpt(:,:)
+         zz3d(:,:,jk) = zzt(:,:)
+
+         ! integrate the pressure from the top of the cell to the bottom of the cell; a factor rho_0 has been removed from the pressure  
+         ! corrected to calculate for loops from 1 to jpj and 1 to jpi 25 Apr 2017 
+         DO jj = 1, jpj
+            DO ji = 1, jpi   ! vector opt.
+!!$               zpb(ji,jj) = zpt(ji,jj) + grav * fse3w_n(ji,jj,jk) * ( znad + rhd(ji,jj,jk) )    
+               zpb(ji,jj) = zpt(ji,jj) + grav * fse3t_n(ji,jj,jk) * ( znad + rhd(ji,jj,jk) )
+               zzb(ji,jj) = grav * ( sshn(ji,jj) - fsdepw_n(ji,jj,jk+1) )   ! surface height had been set to zero; corrected on 24 Apr 2017
+            END DO
+         END DO
+
+         ! calculate the Jacobian of the pressure forces
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! hydrostatic pressure gradient along s-surfaces
-               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
-                  &           * (  fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
-                  &              - fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
-               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
-                  &           * (  fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
-                  &              - fse3w(ji,jj  ,jk) * ( 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 )   &
-                  &           * ( fsde3w(ji+1,jj  ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
-               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
-                  &           * ( fsde3w(ji  ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
-               ! add to the general momentum trend
-               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
-               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
+               zunum(ji, jj) = ( zpb(ji+1,jj) - zpt(ji  , jj) )*( zzb(ji  , jj) - zzt(ji+1, jj) ) - ( zpb(ji  , jj) - zpt(ji+1, jj) )*( zzb(ji+1, jj) - zzt(ji , jj) )   
+               zudenom(ji, jj) =  zpb(ji+1, jj) - zpt(ji  , jj) + zpb(ji  , jj) - zpt(ji+1, jj) 
+               zhpi(ji, jj, jk) = zunum(ji, jj) * r1_e1u(ji, jj) / zudenom(ji, jj)   
+               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
+
+               zvnum = ( zpb(ji,jj+1) - zpt(ji, jj  ) )*( zzb(ji, jj  ) - zzt(ji, jj+1) ) - ( zpb(ji, jj  ) - zpt(ji, jj+1) )*( zzb(ji, jj+1) - zzt(ji, jj  ) )   
+               zhpj(ji, jj, jk) = zvnum * r1_e2v(ji, jj) / ( zpb(ji, jj+1) - zpt(ji, jj  ) + zpb(ji, jj  ) - zpt(ji, jj+1) )   
+               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
             END DO
          END DO
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
+
+! temporary diagnostic outputs (to be used on just a few timesteps on a small configuration) 
+!!$         DO jj = 1, jpj
+!!$            write(numout, *) 'fse3w_n,   jk, jj = ', jk, jj, (fse3w_n(ji,jj,jk), ji = 1, jpi)
+!!$            write(numout, *) 'rhd,     jk, jj = ', jk, jj, (rhd(ji,jj,jk), ji = 1, jpi)
+!!$            write(numout, *) 'gdepw_n, jk+1, jj = ', jk+1, jj, (gdepw_n(ji,jj,jk+1), ji = 1, jpim1)
+!!$            write(numout, *) 'r1_e1u,  jk, jj = ', jk, jj, (r1_e1u(ji,jj), ji = 1, jpim1)
+!!$            write(numout, *) 'zpt,     jk, jj = ', jk, jj, (zpt(ji,jj), ji = 1, jpi)
+!!$            write(numout, *) 'zpb,     jk, jj = ', jk, jj, (zpb(ji,jj), ji = 1, jpi)
+!!$            write(numout, *) 'zzt,     jk, jj = ', jk, jj, (zzt(ji,jj), ji = 1, jpi)
+!!$            write(numout, *) 'zzb,     jk, jj = ', jk, jj, (zzb(ji,jj), ji = 1, jpi)
+!!$            write(numout, *) 'zunum,   jk, jj = ', jk, jj, (zunum(ji,jj), ji = 1, jpi)
+!!$            write(numout, *) 'zudenom, jk, jj = ', jk, jj, (zudenom(ji,jj), ji = 1, jpi)
+!!$            write(numout, *) 'tmask,   jk, jj = ', jk, jj, (tmask(ji,jj,jk), ji = 1, jpi)
+!!$            write(numout, *) 'umask,   jk, jj = ', jk, jj, (umask(ji,jj,jk), ji = 1, jpi)
+!!$            write(numout, *) 'vmask,   jk, jj = ', jk, jj, (vmask(ji,jj,jk), ji = 1, jpi)
+!!$
+!!$            write(numout, *) 'zhpi,    jk, jj = ', jk, jj, (zhpi(ji,jj,jk), ji = 1, jpi)
+!!$            write(numout, *) 'zhpj,    jk, jj = ', jk, jj, (zhpj(ji,jj,jk), ji = 1, jpi)
+!!$         END DO 
+
+      END DO
+      !
+      CALL iom_put("pressure" , zp3d(:,:,:) )
+      CALL iom_put("hpgdepth" , zz3d(:,:,:) )
+      CALL iom_put("rho" , rhd(:,:,:) )
+      CALL iom_put("hpgi", zhpi(:,:,:) )
+      CALL iom_put("hpgj", zhpj(:,:,:) )
+      !
+      CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zp3d, zz3d )
+      CALL wrk_dealloc( jpi,jpj, zpb, zpt, zzb, zzt )
       !
    END SUBROUTINE hpg_sco
