Index: branches/UKMO/dev_3.6_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/UKMO/dev_3.6_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 7980)
+++ branches/UKMO/dev_3.6_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 8476)
@@ -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
