Index: branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 2873)
+++ branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3068)
@@ -1050,7 +1050,7 @@
       !!
       !! The local varialbes for the correction term
-      INTEGER  :: jke, jkw, jkn, jks, jk1
+      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, pn, ps, dze, dzw, dzn, dzs
+      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
@@ -1081,5 +1081,4 @@
       zhpi(:,:,:) = 0.
       
-
       !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here
 
@@ -1118,11 +1117,9 @@
       END DO
 
-      ! Constrained cubic spline interpolation
-
+      !                 ! Construct the vertical density profile with the 
+      !                 !Constrained cubic spline interpolation
       CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)      
 
-
-
-      ! Calculate the pressure at T(ji,jj,1)
+      ! Calculate the hydrostatic pressure at T(ji,jj,1)
       DO jj = 2, jpj
         DO ji = 2, jpi 
@@ -1130,14 +1127,11 @@
                                          bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) &
                                * 0.5_wp * fsde3w(ji,jj,1)
-          rhdt1 = max(rhdt1, 0.0_wp)
-          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1
+          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
 
-
-
-!      print*, 'max&min fse3w=',maxval(fse3w), minval(fse3w)
-!      print*, 'max&min rhd===',maxval(rhd), minval(rhd)
-            
       ! Calculate the pressure at T(ji,jj,2:jpkm1)
       DO jk = 2, jpkm1                                  
@@ -1156,6 +1150,6 @@
       DO jj = 2, jpjm1     
         DO ji = 2, jpim1  
-          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad)
-          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad)
+          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
@@ -1164,147 +1158,92 @@
         DO jj = 2, jpjm1     
           DO ji = 2, jpim1  
-            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk)
+            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  
-            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)
-            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)
-          END DO
-        END DO
-      END DO
-
-      DO jk = 1, jpkm1                                  
-        DO jj = 2, jpjm1     
-          DO ji = 2, jpim1  
-               
-            pe = 0.; pw = 0.; dze = 0.; dzw = 0.
-            pn = 0.; ps = 0.; dzn = 0.; dzs = 0.
-            bhitwe = 0;      bhitns = 0    
+            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,jk) > zuijk) THEN
-              DO jk1 = jk+1, mbathy(ji+1,jj)
-                IF(-fsde3w(ji+1,jj,jk1) > zuijk) THEN
-                  pe = pe + & 
-                       integ2(fsde3w(ji+1,jj,jk1-1),fsde3w(ji+1,jj,jk1),&
-                              asp(ji+1,jj,jk1-1),bsp(ji+1,jj,jk1-1),&
-                              csp(ji+1,jj,jk1-1),dsp(ji+1,jj,jk1-1))
-                  IF(jk1 == mbathy(ji+1,jj)) bhitwe = 1 
-                 ELSE
-                  zze = -fsde3w(ji+1,jj,jk1-1) - zuijk
-                  pe = pe + & 
-                  integ2(fsde3w(ji+1,jj,jk1-1),-zuijk,&
-                         asp(ji+1,jj,jk1-1),bsp(ji+1,jj,jk1-1),&
-                         csp(ji+1,jj,jk1-1),dsp(ji+1,jj,jk1-1))
+            !!!!!     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
-                ENDIF
+                END IF
               END DO
-
-              IF(jk == mbathy(ji+1,jj))  bhitwe = 1 
-
-              IF(bhitwe == 1) zuijk = -fsde3w(ji+1,jj,mbathy(ji+1,jj))  
-
-              DO jk1 = jk-1, 1, -1
-                IF(-fsde3w(ji,jj,jk1) < zuijk) THEN
-                  pw = pw + &
-                           integ2(fsde3w(ji,jj,jk1),fsde3w(ji,jj,jk1+1),&
-                                   asp(ji,jj,jk1),bsp(ji,jj,jk1),&
-                                   csp(ji,jj,jk1),dsp(ji,jj,jk1))
-                  IF(jk1 == 1) THEN
-                    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)) &
-                                         * (zuijk + fsde3w(ji,jj,1))
-                    rhdt1 = max(rhdt1,0.0_wp)
-                    rrw = rhdt1
-                    zzw = zuijk + fsde3w(ji,jj,1)
-                    pw = pw + min(0.5 * (rhd(ji,jj,1) + rrw) * zzw, zhpi(ji,jj,1))
-                  END IF
-                 ELSE
-                  zzw = zuijk + fsde3w(ji,jj,jk1+1)
-                  pw = pw + &
-                           integ2(-zuijk,fsde3w(ji,jj,jk1+1),&
-                                   asp(ji,jj,jk1),bsp(ji,jj,jk1),&
-                                   csp(ji,jj,jk1),dsp(ji,jj,jk1))
+            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
-                ENDIF
+                END IF
               END DO
-
-              IF(jk == 1) THEN
-                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)) &
-                                     * (zuijk + fsde3w(ji,jj,1))
-                rhdt1 = max(rhdt1,0.0_wp)
-                rrw = rhdt1
-                zzw = zuijk + fsde3w(ji,jj,1)
-                pw = pw + min(0.5 * (rhd(ji,jj,1) + rrw) * zzw, zhpi(ji,jj,1))
-              END IF
-
-              ELSE
-
-              DO jk1 = jk+1, mbathy(ji,jj)
-                IF(-fsde3w(ji,jj,jk1) > zuijk) THEN
-                  pw = pw - &
-                          integ2(fsde3w(ji,jj,jk1-1),fsde3w(ji,jj,jk1),&
-                                  asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),&
-                                  csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1))
-                  IF(jk1 == mbathy(ji,jj)) bhitwe = 1 
-                ELSE
-                  zzw = -fsde3w(ji,jj,jk1-1) - zuijk
-                  pw = pw - &
-                           integ2(fsde3w(ji,jj,jk1-1),-zuijk,&
-                                  asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),&
-                                  csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1))
-                  EXIT
-                ENDIF
-              END DO
-
-              IF(jk == mbathy(ji,jj)) bhitwe = 1 
-
-              IF(bhitwe == 1) zuijk = -fsde3w(ji,jj,mbathy(ji,jj)) 
-
-              DO jk1 = jk-1, 1, -1
-                IF(-fsde3w(ji+1,jj,jk1) < zuijk) THEN
-                  pe = pe - &
-                           integ2(fsde3w(ji+1,jj,jk1),fsde3w(ji+1,jj,jk1+1),&
-                                   asp(ji+1,jj,jk1),bsp(ji+1,jj,jk1),&
-                                   csp(ji+1,jj,jk1),dsp(ji+1,jj,jk1))
-                  IF(jk1 == 1) THEN
-                    rhdt1 = rhd(ji+1,jj,1) - interp3(fsde3w(ji+1,jj,1),asp(ji+1,jj,1), &
-                                           bsp(ji+1,jj,1),csp(ji+1,jj,1),dsp(ji+1,jj,1)) &
-                                         * (zuijk + fsde3w(ji+1,jj,1))
-                    rhdt1 = max(rhdt1,0.0_wp)
-                    rre = rhdt1
-                    zze = zuijk + fsde3w(ji+1,jj,1)
-                    pe = pe - min(0.5 * (rhd(ji+1,jj,1) + rre) * zze, zhpi(ji+1,jj,1))
-                  END IF
-                 ELSE
-                  zze = zuijk + fsde3w(ji+1,jj,jk1+1)
-                  pe = pe - &
-                           integ2(-zuijk,fsde3w(ji+1,jj,jk1+1),&
-                                   asp(ji+1,jj,jk1),bsp(ji+1,jj,jk1),&
-                                   csp(ji+1,jj,jk1),dsp(ji+1,jj,jk1))
-                  EXIT
-                ENDIF
-              END DO
-
-              IF(jk == 1) THEN
-                 rhdt1 = rhd(ji+1,jj,1) - interp3(fsde3w(ji+1,jj,1),asp(ji+1,jj,1), &
-                                          bsp(ji+1,jj,1),csp(ji+1,jj,1),dsp(ji+1,jj,1)) &
-                                        * (zuijk + fsde3w(ji+1,jj,1))
-                rhdt1 = max(rhdt1,0.0_wp)
-                rre = rhdt1
-                zze = zuijk + fsde3w(ji+1,jj,1)
-                pe = pe - min(0.5 * (rhd(ji+1,jj,1) + rre) * zze, zhpi(ji+1,jj,1))
-              END IF
-
             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))
@@ -1317,138 +1256,88 @@
             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,jk) > zvijk) THEN
-              DO jk1 = jk+1, mbathy(ji,jj+1)
-                IF(-fsde3w(ji,jj+1,jk1) > zvijk) THEN
-                  pn = pn + &
-                           integ2(fsde3w(ji,jj+1,jk1-1),fsde3w(ji,jj+1,jk1),&
-                                   asp(ji,jj+1,jk1-1),bsp(ji,jj+1,jk1-1),&
-                                   csp(ji,jj+1,jk1-1),dsp(ji,jj+1,jk1-1))
-                  IF(jk1 == mbathy(ji,jj+1)) bhitns = 1 
-                 ELSE
-                  zzn = -fsde3w(ji,jj+1,jk1-1) - zvijk
-                  pn = pn + &
-                           integ2(fsde3w(ji,jj+1,jk1-1),-zvijk,&
-                                   asp(ji,jj+1,jk1-1),bsp(ji,jj+1,jk1-1),&
-                                   csp(ji,jj+1,jk1-1),dsp(ji,jj+1,jk1-1))
+  
+            !!!!!     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
-                ENDIF
+                END IF
               END DO
-
-              IF(jk == mbathy(ji,jj+1)) bhitns = 1
-
-
-              IF(bhitns == 1) zvijk = -fsde3w(ji,jj+1,mbathy(ji,jj+1)) 
-
-
-              DO jk1 = jk-1, 1, -1
-                IF(-fsde3w(ji,jj,jk1) < zvijk) THEN
-                  ps = ps + &
-                           integ2(fsde3w(ji,jj,jk1),fsde3w(ji,jj,jk1+1),&
-                                  asp(ji,jj,jk1),bsp(ji,jj,jk1),&
-                                  csp(ji,jj,jk1),dsp(ji,jj,jk1))
-                  IF(jk1 == 1) THEN
-                    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)) &
-                                         * (zvijk + fsde3w(ji,jj,1))
-                    rhdt1 = max(rhdt1,0.0_wp)
-                    rrs = rhdt1
-                    zzs = zvijk + fsde3w(ji,jj,1)
-                    ps = ps + min(0.5 * (rhd(ji,jj,1) + rrs ) * zzs, zhpi(ji,jj,1))
-                  END IF
-                ELSE
-                  zzs = zvijk + fsde3w(ji,jj,jk1+1)
-                  ps = ps + &
-                           integ2(-zvijk,fsde3w(ji,jj,jk1+1),&
-                                   asp(ji,jj,jk1),bsp(ji,jj,jk1),&
-                                   csp(ji,jj,jk1),dsp(ji,jj,jk1))
+            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
-                ENDIF
+                END IF
               END DO
-
-              IF(jk == 1) THEN
-                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)) &
-                                     * (zvijk + fsde3w(ji,jj,1))
-                rhdt1 = max(rhdt1,0.0_wp)
-                rrs = rhdt1
-                zzs = zvijk + fsde3w(ji,jj,1)
-                ps = ps + min(0.5 * (rhd(ji,jj,1) + rrs) * zzs, zhpi(ji,jj,1))
-              END IF
-
+            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
-
-              DO jk1 = jk+1, mbathy(ji,jj)
-                IF(-fsde3w(ji,jj,jk1) > zvijk) THEN
-                  ps = ps - &
-                           integ2(fsde3w(ji,jj,jk1-1),fsde3w(ji,jj,jk1),&
-                                   asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),&
-                                   csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1))
-                  IF(jk1 == mbathy(ji,jj)) bhitns = 1
-                 ELSE
-                  zzs = -fsde3w(ji,jj,jk1-1) - zvijk
-                  ps = ps - &
-                           integ2(fsde3w(ji,jj,jk1-1),-zvijk,&
-                                   asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),&
-                                   csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1))
-                  EXIT
-                ENDIF
-              END DO
-
-              IF(jk == mbathy(ji,jj)) bhitns = 1 
-
-
-              IF(bhitns == 1) zvijk = -fsde3w(ji,jj,mbathy(ji,jj))
-
-                DO jk1 = jk-1, 1, -1
-                  IF(-fsde3w(ji,jj+1,jk1) < zvijk) THEN
-                    pn = pn - &
-                             integ2(fsde3w(ji,jj+1,jk1),fsde3w(ji,jj+1,jk1+1),&
-                                    asp(ji,jj+1,jk1),bsp(ji,jj+1,jk1),&
-                                    csp(ji,jj+1,jk1),dsp(ji,jj+1,jk1))
-                    IF(jk1 == 1) THEN
-                      rhdt1 = rhd(ji,jj+1,1) - interp3(fsde3w(ji,jj+1,1),asp(ji,jj+1,1), &
-                                             bsp(ji,jj+1,1),csp(ji,jj+1,1),dsp(ji,jj+1,1)) &
-                                           * (zvijk + fsde3w(ji,jj+1,1))
-                      rhdt1 = max(rhdt1,0.0_wp)
-                      rrn = rhdt1
-                      zzn = zvijk + fsde3w(ji,jj+1,1)
-                      pn = pn - min(0.5 * (rhd(ji,jj+1,1) + rrn) * zzn, zhpi(ji,jj+1,1))
-                    END IF
-
-                   ELSE
-                    zzn = zvijk + fsde3w(ji,jj+1,jk1+1)
-                    pn = pn - &
-                             integ2(-zvijk,fsde3w(ji,jj+1,jk1+1),&
-                                     asp(ji,jj+1,jk1),bsp(ji,jj+1,jk1),&
-                                     csp(ji,jj+1,jk1),dsp(ji,jj+1,jk1))
-                    EXIT
-                  ENDIF
-                END DO
-
-                IF(jk == 1) THEN
-                  rhdt1 = rhd(ji,jj+1,1) - interp3(fsde3w(ji,jj+1,1),asp(ji,jj+1,1), &
-                                          bsp(ji,jj+1,1),csp(ji,jj+1,1),dsp(ji,jj+1,1)) &
-                                        * (zvijk + fsde3w(ji,jj+1,1))
-                  rhdt1 = max(rhdt1,0.0_wp)
-                  rrn = rhdt1
-                  zzn = zvijk + fsde3w(ji,jj+1,1)
-                  pn = pn - min(0.5 * (rhd(ji,jj+1,1) + rrn) * zzn, zhpi(ji,jj+1,1))
-                END IF
-
-              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)
+              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)
 
                     
@@ -1458,5 +1347,4 @@
 
       ! Restore fsde3w and rhd
-
       fsde3w(:,:,:) = fsde3w_tmp(:,:,:)
       rhd(:,:,:)    = rhd_tmp(:,:,:)
@@ -1466,5 +1354,4 @@
          CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays')
       !
-
    END SUBROUTINE hpg_prj
 
