Index: /branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 3074)
@@ -633,8 +633,16 @@
 $\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true)
 
-Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume 
+$\bullet$ Pressure Jacobian scheme (prj) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_prj}=true)
+
+Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume 
 formulation is activated (\key{vvl}) because in that case, even with a flat bottom, 
 the coordinate surfaces are not horizontal but follow the free surface 
-\citep{Levier2007}. The other pressure gradient options are not yet available.
+\citep{Levier2007}. Only the pressure jacobian scheme (\np{ln\_dynhpg\_prj}=true) is available as an 
+alternative to the default \np{ln\_dynhpg\_sco}=true when \key{vvl} is active.  The pressure Jacobian scheme uses 
+a constrained cubic spline to reconstruct the density profile across the water column. This method
+maintains the monotonicity between the density nodes and is of a higher order than the linear
+interpolation method. The pressure can be calculated by analytical integration of the density profile and
+a pressure Jacobian method is used to solve the horizontal pressure gradient. This method should
+provide a more accurate calculation of the horizontal pressure gradient than the standard scheme.
 
 %--------------------------------------------------------------------------------------------------------------
Index: /branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_hpg
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_hpg	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_hpg	(revision 3074)
@@ -9,4 +9,5 @@
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme)
+   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
    rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
Index: /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 3074)
@@ -535,4 +535,5 @@
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme)
+   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
    rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
Index: /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 3074)
@@ -533,4 +533,5 @@
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme)
+   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
    rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
Index: /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 3074)
@@ -535,4 +535,5 @@
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme)
+   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
    rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
Index: /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 3074)
@@ -536,4 +536,5 @@
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme)
+   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
    rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme)
    ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T)
Index: /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist	(revision 3074)
@@ -535,4 +535,5 @@
    ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial)
    ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme)
+   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme)
    rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme)
    ln_dynhpg_imp = .true.  !  time stepping: semi-implicit time scheme  (T)
Index: /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3073)
+++ /branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3074)
@@ -27,4 +27,5 @@
    !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial)
    !!       hpg_rot  : s-coordinate (ROTated axes scheme)
+   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial)
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -52,8 +53,9 @@
    LOGICAL , PUBLIC ::   ln_hpg_djc    = .FALSE.   !: s-coordinate (Density Jacobian with Cubic polynomial)
    LOGICAL , PUBLIC ::   ln_hpg_rot    = .FALSE.   !: s-coordinate (ROTated axes scheme)
+   LOGICAL , PUBLIC ::   ln_hpg_prj    = .FALSE.   !: s-coordinate (Pressure Jacobian scheme)
    REAL(wp), PUBLIC ::   rn_gamma      = 0._wp     !: weighting coefficient
    LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag
 
-   INTEGER  ::   nhpg  =  0   ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)
+   INTEGER  ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)
 
    !! * Substitutions
@@ -100,4 +102,5 @@
       CASE (  5 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
       CASE (  6 )   ;   CALL hpg_rot    ( kt )      ! s-coordinate (ROTated axes scheme)
+      CASE (  7 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
       END SELECT
       !
@@ -129,5 +132,6 @@
       !!
       NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,    &
-         &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma  , ln_dynhpg_imp
+         &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj,    &
+         &                 rn_gamma  , ln_dynhpg_imp
       !!----------------------------------------------------------------------
       !
@@ -147,10 +151,13 @@
          WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc
          WRITE(numout,*) '      s-coord. (ROTated axes scheme)                    ln_hpg_rot    = ', ln_hpg_rot
+         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj
          WRITE(numout,*) '      weighting coeff. (wdj scheme)                     rn_gamma      = ', rn_gamma
          WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp
       ENDIF
       !
-      IF( lk_vvl .AND. .NOT. ln_hpg_sco )   &
-         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco')
+      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   &
+         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&
+                           & the standard jacobian formulation hpg_sco or &
+                           & the pressure jacobian formulation hpg_prj')
       !
       !                               ! Set nhpg from ln_hpg_... flags
@@ -162,4 +169,5 @@
       IF( ln_hpg_djc )   nhpg = 5
       IF( ln_hpg_rot )   nhpg = 6
+      IF( ln_hpg_prj )   nhpg = 7
       !
       !                               ! Consitency check
@@ -172,4 +180,5 @@
       IF( ln_hpg_djc )   ioptio = ioptio + 1
       IF( ln_hpg_rot )   ioptio = ioptio + 1
+      IF( ln_hpg_prj )   ioptio = ioptio + 1
       IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
       !
@@ -1005,4 +1014,543 @@
    END SUBROUTINE hpg_rot
 
+   
+   SUBROUTINE hpg_prj( kt )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE hpg_prj  ***
+      !!
+      !! ** Method  :   s-coordinate case.
+      !!      Reformulate the horizontal hydrostatical pressure gradient
+      !!      term using Pressure Jacobian. A new correction term 
+      !!      is developed to eliminate the sigma-coordinate error.
+      !!
+      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
+      !!             - Save the trend (l_trddyn=T)
+      !!
+      !!----------------------------------------------------------------------
+
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE oce     , ONLY:   fsde3w_tmp => ta       ! (ta,sa) used as 3D workspace
+      USE oce     , ONLY:      rhd_tmp => sa 
+      USE wrk_nemo, ONLY:         zhpi => wrk_3d_3 
+      USE wrk_nemo, ONLY:           zu => wrk_3d_4 
+      USE wrk_nemo, ONLY:           zv => wrk_3d_5
+      USE wrk_nemo, ONLY:          fsp => wrk_3d_6 
+      USE wrk_nemo, ONLY:          xsp => wrk_3d_7
+      USE wrk_nemo, ONLY:          asp => wrk_3d_8
+      USE wrk_nemo, ONLY:          bsp => wrk_3d_9
+      USE wrk_nemo, ONLY:          csp => wrk_3d_10
+      USE wrk_nemo, ONLY:          dsp => wrk_3d_11
+      !!
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
+      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk, jkk                ! dummy loop indices
+      REAL(wp) ::   zcoef0, znad               ! temporary scalars
+      !!
+      !! The local varialbes for the correction term
+      INTEGER  :: jke, jkw, jkn, jks, jk1, jis, jid, jjs, jjd
+      REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs
+      REAL(wp) :: zuijk, zvijk, pe, pw, pwes, pwed, pn, ps, pnss, pnsd, deps
+      REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2
+      REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2
+      INTEGER  :: bhitwe, bhitns
+      !!----------------------------------------------------------------------
+
+      IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN
+         CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable')   ;   RETURN
+      ENDIF
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian'
+      ENDIF
+
+      !!----------------------------------------------------------------------
+      ! Local constant initialization
+      zcoef0 = - grav 
+      znad = 0.0_wp
+      IF(lk_vvl) znad = 1._wp
+
+      ! Save fsde3w and rhd
+      fsde3w_tmp(:,:,:) = fsde3w(:,:,:) 
+      rhd_tmp(:,:,:)    = rhd(:,:,:)   
+
+      ! Clean 3-D work arraies
+      zhpi(:,:,:) = 0.
+      
+      !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here
+
+      ! Preparing vertical density profile for hybrid-sco coordinate
+      DO jj = 1, jpj
+        DO ji = 1, jpi   
+          jk = mbathy(ji,jj)
+          IF(jk <= 0) THEN; rhd(ji,jj,:) = 0._wp
+            ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)
+            ELSE IF(jk < jpkm1) THEN
+            DO jkk = jk+1, jpk
+              rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),&
+                                       fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))
+            END DO 
+          END IF
+        END DO
+      END DO
+
+      DO jj = 1, jpj
+        DO ji = 1, jpi
+          fsde3w(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1)
+          fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad
+          DO jk = 2, jpk
+            fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+
+      DO jk = 1, jpkm1
+        DO jj = 1, jpj
+          DO ji = 1, jpi
+            fsp(ji,jj,jk) = rhd(ji,jj,jk)
+            xsp(ji,jj,jk) = fsde3w(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+
+      !                 ! Construct the vertical density profile with the 
+      !                 !Constrained cubic spline interpolation
+      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)      
+
+      ! Calculate the hydrostatic pressure at T(ji,jj,1)
+      DO jj = 2, jpj
+        DO ji = 2, jpi 
+          rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), &
+                                         bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) &
+                               * 0.5_wp * fsde3w(ji,jj,1)
+          rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+
+     !                  ! assuming linear profile across the top half surface layer
+          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1  
+        END DO
+      END DO
+
+      ! Calculate the pressure at T(ji,jj,2:jpkm1)
+      DO jk = 2, jpkm1                                  
+        DO jj = 2, jpj     
+          DO ji = 2, jpi
+            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &
+                             integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),&
+                                    asp(ji,jj,jk-1),bsp(ji,jj,jk-1),&
+                                    csp(ji,jj,jk-1),dsp(ji,jj,jk-1))
+          END DO
+        END DO
+      END DO
+
+      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
+
+      DO jj = 2, jpjm1     
+        DO ji = 2, jpim1  
+          zu(ji,jj,1) = - ( 0.5_wp * fse3uw(ji,jj,1) - sshu_n(ji,jj) * znad)
+          zv(ji,jj,1) = - ( 0.5_wp * fse3vw(ji,jj,1) - sshv_n(ji,jj) * znad)
+        END DO
+      END DO
+
+      DO jk = 2, jpkm1                                  
+        DO jj = 2, jpjm1     
+          DO ji = 2, jpim1  
+            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3uw(ji,jj,jk)
+            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3vw(ji,jj,jk)
+          END DO
+        END DO
+      END DO
+               
+      !  Start pressure integration
+
+      DO jk = 1, jpkm1                                  
+        DO jj = 2, jpjm1     
+          DO ji = 2, jpim1  
+            pwes = 0._wp; pwed = 0._wp
+            pnss = 0._wp; pnsd = 0._wp
+            zuijk = zu(ji,jj,jk)
+            zvijk = zv(ji,jj,jk)
+
+            !!!!!     for u equation
+            IF( -fsde3w(ji+1,jj,mbathy(ji+1,jj)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN
+              jis = ji + 1; jid = ji
+            ELSE
+              jis = ji;     jid = ji +1
+            ENDIF
+
+            !     !integrate the pressure on the shallow side
+            jk1 = jk 
+            bhitwe = 0
+            IF( jk1 == mbathy(jis,jj) ) THEN
+              bhitwe = 1
+            ELSE
+              DO WHILE ( -fsde3w(jis,jj,jk1+1) > zuijk )
+                pwes = pwes + & 
+                       integ2(fsde3w(jis,jj,jk1), fsde3w(jis,jj,jk1+1),&
+                              asp(jis,jj,jk1)   , bsp(jis,jj,jk1),     &
+                              csp(jis,jj,jk1)   , dsp(jis,jj,jk1))
+                jk1 = jk1 + 1
+                IF( jk1 == mbathy(jis,jj) ) THEN
+                  bhitwe = 1
+                  EXIT
+                END IF
+              END DO
+            ENDIF
+            
+            IF( bhitwe /= 1 ) THEN
+              pwes = pwes + & 
+                     integ2(fsde3w(jis,jj,jk1), -zuijk,      &
+                            asp(jis,jj,jk1), bsp(jis,jj,jk1),&
+                            csp(jis,jj,jk1), dsp(jis,jj,jk1))
+            ELSE
+               zuijk = -fsde3w(jis,jj,jk1)
+            ENDIF
+
+            !     !integrate the pressure on the deep side
+            jk1 = jk 
+            bhitwe = 0
+            IF( jk1 == 1 ) THEN
+              bhitwe = 1
+            ELSE
+              DO WHILE ( -fsde3w(jid,jj,jk1-1) < zuijk )
+                pwed = pwed + & 
+                       integ2(fsde3w(jid,jj,jk1-1), fsde3w(jid,jj,jk1),&
+                              asp(jid,jj,jk1-1),    bsp(jid,jj,jk1-1), &
+                              csp(jid,jj,jk1-1),    dsp(jid,jj,jk1-1))
+                jk1 = jk1 - 1
+                IF( jk1 == 1 ) THEN
+                  bhitwe = 1
+                  EXIT
+                END IF
+              END DO
+            ENDIF
+            
+            IF( bhitwe /= 1 ) THEN
+              pwed = pwed + & 
+                     integ2(-zuijk,             fsde3w(jid,jj,jk1),&
+                             asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), &
+                             csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1))
+            ELSE
+              deps  = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad)
+              rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1), asp(jid,jj,1), &
+                                              bsp(jid,jj,1),    csp(jid,jj,1), &
+                                              dsp(jid,jj,1) ) * deps
+              rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+              pwed  = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps
+            ENDIF
+
+            IF( jid > jis ) THEN
+              pe = pwed; pw = pwes
+            ELSE
+              pe = pwes; pw = pwed
+            ENDIF
+
+            dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
+            IF( lk_vvl ) THEN
+              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj))) 
+             ELSE
+              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe) 
+            ENDIF
+
+            ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * &
+               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
+  
+            !!!!!     for v equation
+
+            IF( -fsde3w(ji,jj+1,mbathy(ji,jj+1)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN
+              jjs = jj + 1; jjd = jj
+            ELSE
+              jjs = jj    ; jjd = jj + 1
+            ENDIF
+
+            !     !integrate the pressure on the shallow side
+            jk1 = jk 
+            bhitns = 0
+            IF( jk1 == mbathy(ji,jjs) ) THEN
+              bhitns = 1
+            ELSE
+              DO WHILE ( -fsde3w(ji,jjs,jk1+1) > zvijk )
+                pnss = pnss + & 
+                       integ2(fsde3w(ji,jjs,jk1), fsde3w(ji,jjs,jk1+1),&
+                              asp(ji,jjs,jk1),    bsp(ji,jjs,jk1),     &
+                              csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
+                jk1 = jk1 + 1
+                IF( jk1 == mbathy(ji,jjs) ) THEN
+                  bhitns = 1
+                  EXIT
+                END IF
+              END DO
+            ENDIF
+            
+            IF( bhitns /= 1 ) THEN
+              pnss = pnss + & 
+                     integ2(fsde3w(ji,jjs,jk1), -zvijk,          &
+                            asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
+                            csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
+            ELSE
+               zvijk = -fsde3w(ji,jjs,jk1)
+            ENDIF
+
+            !     !integrate the pressure on the deep side
+            jk1 = jk 
+            bhitns = 0
+            IF( jk1 == 1 ) THEN
+              bhitns = 1
+            ELSE
+              DO WHILE ( -fsde3w(ji,jjd,jk1-1) < zvijk )
+                pnsd = pnsd + & 
+                       integ2(fsde3w(ji,jjd,jk1-1), fsde3w(ji,jjd,jk1), &
+                              asp(ji,jjd,jk1-1),    bsp(ji,jjd,jk1-1),  &
+                              csp(ji,jjd,jk1-1),    dsp(ji,jjd,jk1-1) )
+                jk1 = jk1 - 1
+                IF( jk1 == 1 ) THEN
+                  bhitns = 1
+                  EXIT
+                END IF
+              END DO
+            ENDIF
+            
+            IF( bhitns /= 1 ) THEN
+              pnsd = pnsd + & 
+                     integ2(-zvijk,            fsde3w(ji,jjd,jk1), &
+                            asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1),  &
+                            csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
+            ELSE
+              deps  = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad)
+              rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1), asp(ji,jjd,1), &
+                                              bsp(ji,jjd,1),    csp(ji,jjd,1), &
+                                              dsp(ji,jjd,1) ) * deps
+              rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+              pnsd  = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps
+            ENDIF
+
+            IF( jjd > jjs ) THEN
+              pn = pnsd; ps = pnss
+            ELSE
+              pn = pnss; ps = pnsd
+            ENDIF
+
+            dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
+            if( lk_vvl ) then
+                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj))) 
+              else
+                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn ) 
+            end if
+
+            va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)* &
+            &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)
+
+                    
+           END DO
+        END DO
+      END DO
+
+      ! Restore fsde3w and rhd
+      fsde3w(:,:,:) = fsde3w_tmp(:,:,:)
+      rhd(:,:,:)    = rhd_tmp(:,:,:)
+
+      !
+      IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) )   &
+         CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays')
+      !
+   END SUBROUTINE hpg_prj
+
+   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type)
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE cspline  ***
+      !!       
+      !! ** Purpose :   constrained cubic spline interpolation
+      !!          
+      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 
+      !! Reference: K.W. Brodlie, A review of mehtods for curve and function
+      !!                          drawing, 1980
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), DIMENSION(:,:,:), INTENT(in)  :: fsp, xsp           ! value and coordinate
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 
+                                                                    ! the interpoated function
+      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline 
+                                                                    ! 2: Linear
+
+      ! Local Variables      
+      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
+      INTEGER  ::   jpi, jpj, jpkm1
+      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
+      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha
+      REAL(wp) ::   zdf(size(fsp,3))
+      !!----------------------------------------------------------------------
+
+      jpi   = size(fsp,1)
+      jpj   = size(fsp,2)
+      jpkm1 = size(fsp,3) - 1
+
+      ! Clean output arrays
+      asp = 0.0_wp
+      bsp = 0.0_wp
+      csp = 0.0_wp
+      dsp = 0.0_wp
+      
+      
+      Do ji = 1, jpi
+      Do jj = 1, jpj
+
+      If (polynomial_type == 1) Then     !Constrained Cubic Spline
+          Do jk = 2, jpkm1-1
+             zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)  
+             zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
+             zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
+             zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
+  
+             zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
+             
+             If(zdf1 * zdf2 <= 0._wp) Then
+               zdf(jk) = 0._wp
+             Else
+               zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
+             End If
+          End Do
+  
+          zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - &
+                    &  0.5_wp * zdf(2)
+          zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
+                    &           ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - &
+                    & 0.5_wp * zdf(jpkm1 - 1)
+  
+          Do jk = 1, jpkm1 - 1
+             zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 
+             ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
+             ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
+             zddf1  = -2._wp * ztmp1 + ztmp2 
+             ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
+             zddf2  =  2._wp * ztmp1 - ztmp2 
+  
+             dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
+             csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
+             bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 
+                           & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
+                           & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + &
+                           &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + &
+                           &                   xsp(ji,jj,jk)**2 )
+             asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - &
+                           &                 csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - &
+                           &                 dsp(ji,jj,jk) * xsp(ji,jj,jk)**3
+          End Do
+
+       Else If (polynomial_type == 2) Then     !Linear
+
+        Do jk = 1, jpkm1-1
+           zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 
+           ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
+  
+           dsp(ji,jj,jk) = 0._wp
+           csp(ji,jj,jk) = 0._wp
+           bsp(ji,jj,jk) = ztmp1 / zdxtmp
+           asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
+        End Do
+
+       Else
+        CALL ctl_stop( 'invalid polynomial type in cspline' )
+      End If
+
+      End Do
+      End Do
+      
+   End Subroutine cspline
+
+
+   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE interp1  ***
+      !!       
+      !! ** Purpose :   1-d linear interpolation
+      !!          
+      !! ** Method  :  
+      !!                interpolation is straight forward
+      !!                extrapolation is also permitted (no value limit) 
+      !!
+      !! H.Liu, Jan 2009,  POL 
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr   
+      REAL(wp)             ::  f ! result of the interpolation (extrapolation)
+      REAL(wp)             ::  zdeltx
+      !!----------------------------------------------------------------------
+
+      zdeltx = xr - xl
+      IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN
+        f = 0.5_wp * (fl + fr)
+       ELSE
+        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
+      END IF
+      
+   END FUNCTION interp1
+
+   FUNCTION interp2(x, a, b, c, d)  RESULT(f) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE interp1  ***
+      !!       
+      !! ** Purpose :   1-d constrained cubic spline interpolation
+      !!          
+      !! ** Method  :  cubic spline interpolation
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  x, a, b, c, d   
+      REAL(wp)             ::  f ! value from the interpolation
+      !!----------------------------------------------------------------------
+
+      f = a + x* ( b + x * ( c + d * x ) ) 
+
+   END FUNCTION interp2
+
+
+   FUNCTION interp3(x, a, b, c, d)  RESULT(f) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE interp1  ***
+      !!       
+      !! ** Purpose :   deriavtive of a cubic spline function
+      !!          
+      !! ** Method  :  
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  x, a, b, c, d   
+      REAL(wp)             ::  f ! value from the interpolation
+      !!----------------------------------------------------------------------
+
+      f = b + x * ( 2._wp * c + 3._wp * d * x)
+
+   END FUNCTION interp3
+
+   
+   FUNCTION integ2(xl, xr, a, b, c, d)  RESULT(f) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE interp1  ***
+      !!       
+      !! ** Purpose :   1-d constrained cubic spline integration
+      !!          
+      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr 
+      !!
+      !!----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d   
+      REAL(wp)             ::  za1, za2, za3      
+      REAL(wp)             ::  f                   ! integration result
+      !!----------------------------------------------------------------------
+
+      za1 = 0.5_wp * b 
+      za2 = c / 3.0_wp 
+      za3 = 0.25_wp * d 
+
+      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
+         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
+
+   END FUNCTION integ2
+
+
    !!======================================================================
 END MODULE dynhpg
