Index: /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
===================================================================
--- /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3433)
+++ /trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90	(revision 3434)
@@ -678,11 +678,10 @@
       REAL(wp) :: zrhdt1 
       REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
-      INTEGER  :: zbhitwe, zbhitns
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdeptht, zrhh 
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdept, zrhh 
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
       !!----------------------------------------------------------------------
       !
       CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
-      CALL wrk_alloc( jpi,jpj,jpk, zdeptht, zrhh ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 
       !
       IF( kt == nit000 ) THEN
@@ -717,23 +716,15 @@
       END DO
 
-      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)"
-      DO jj = 1, jpj
-        DO ji = 1, jpi
-          zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1)
-          zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad
-          DO jk = 2, jpk
-             zdeptht(ji,jj,jk) = zdeptht(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) = zrhh(ji,jj,jk)
-            xsp(ji,jj,jk) = zdeptht(ji,jj,jk)
-          END DO
-        END DO
-      END DO
+      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
+      DO jj = 1, jpj;   DO ji = 1, jpi
+          zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad
+      END DO        ;   END DO
+
+      DO jk = 2, jpk;   DO jj = 1, jpj;   DO ji = 1, jpi
+          zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)
+      END DO        ;   END DO        ;   END DO
+
+      fsp(:,:,:) = zrhh(:,:,:)
+      xsp(:,:,:) = zdept(:,:,:)
 
       ! Construct the vertical density profile with the 
@@ -745,8 +736,7 @@
       DO jj = 2, jpj
         DO ji = 2, jpi 
-          zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), &
+          zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), &
                                          bsp(ji,jj,1),   csp(ji,jj,1), &
-                                         dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1)
-          zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water
+                                         dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1)
 
           ! assuming linear profile across the top half surface layer
@@ -760,5 +750,5 @@
           DO ji = 2, jpi
             zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                          &
-                             integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),&
+                             integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&
                                     asp(ji,jj,jk-1),    bsp(ji,jj,jk-1), &
                                     csp(ji,jj,jk-1),    dsp(ji,jj,jk-1))
@@ -793,4 +783,16 @@
       END DO
 
+      DO jk = 1, jpkm1
+        DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+            zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
+            zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
+            zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
+            zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
+          END DO
+        END DO
+      END DO
+
+
       DO jk = 1, jpkm1                                  
         DO jj = 2, jpjm1     
@@ -803,5 +805,5 @@
             !!!!!     for u equation
             IF( jk <= mbku(ji,jj) ) THEN
-               IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN
+               IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
                  jis = ji + 1; jid = ji
                ELSE
@@ -811,13 +813,12 @@
                ! integrate the pressure on the shallow side
                jk1 = jk 
-               zbhitwe = 0
-               DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk )
+               DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
                  IF( jk1 == mbku(ji,jj) ) THEN
-                   zbhitwe = 1
+                   zuijk = -zdept(jis,jj,jk1)
                    EXIT
                  ENDIF
-                 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk)
+                 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
                  zpwes = zpwes +                                    & 
-                      integ2(zdeptht(jis,jj,jk1), zdeps,            &
+                      integ_spline(zdept(jis,jj,jk1), zdeps,            &
                              asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
                              csp(jis,jj,jk1),    dsp(jis,jj,jk1))
@@ -825,19 +826,18 @@
                END DO
             
-               IF(zbhitwe == 1) THEN
-                 zuijk = -zdeptht(jis,jj,jk1)
-               ENDIF
-
                ! integrate the pressure on the deep side
                jk1 = jk 
-               zbhitwe = 0
-               DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk )
+               DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
                  IF( jk1 == 1 ) THEN
-                   zbhitwe = 1
+                   zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
+                   zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
+                                                     bsp(jid,jj,1),   csp(jid,jj,1), &
+                                                     dsp(jid,jj,1)) * zdeps
+                   zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
                    EXIT
                  ENDIF
-                 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk)
+                 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
                  zpwed = zpwed +                                        & 
-                        integ2(zdeps,              zdeptht(jid,jj,jk1), &
+                        integ_spline(zdeps,              zdept(jid,jj,jk1), &
                                asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
                                csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
@@ -845,13 +845,4 @@
                END DO
             
-               IF( zbhitwe == 1 ) THEN
-                 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
-                 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), &
-                                                 bsp(jid,jj,1),    csp(jid,jj,1), &
-                                                 dsp(jid,jj,1)) * zdeps
-                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water
-                 zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
-               ENDIF
-
                ! update the momentum trends in u direction
 
@@ -870,5 +861,5 @@
             !!!!!     for v equation
             IF( jk <= mbkv(ji,jj) ) THEN
-               IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN
+               IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
                  jjs = jj + 1; jjd = jj
                ELSE
@@ -878,13 +869,12 @@
                ! integrate the pressure on the shallow side
                jk1 = jk 
-               zbhitns = 0
-               DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk )
+               DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
                  IF( jk1 == mbkv(ji,jj) ) THEN
-                   zbhitns = 1
+                   zvijk = -zdept(ji,jjs,jk1)
                    EXIT
                  ENDIF
-                 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk)
+                 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
                  zpnss = zpnss +                                      & 
-                        integ2(zdeptht(ji,jjs,jk1), zdeps,            &
+                        integ_spline(zdept(ji,jjs,jk1), zdeps,            &
                                asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
                                csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
@@ -892,19 +882,18 @@
                END DO
             
-               IF(zbhitns == 1) THEN
-                 zvijk = -zdeptht(ji,jjs,jk1)
-               ENDIF
-
                ! integrate the pressure on the deep side
                jk1 = jk 
-               zbhitns = 0
-               DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk )
+               DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
                  IF( jk1 == 1 ) THEN
-                   zbhitns = 1
+                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
+                   zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
+                                                     bsp(ji,jjd,1),   csp(ji,jjd,1), &
+                                                     dsp(ji,jjd,1) ) * zdeps
+                   zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
                    EXIT
                  ENDIF
-                 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk)
+                 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
                  zpnsd = zpnsd +                                        & 
-                        integ2(zdeps,              zdeptht(ji,jjd,jk1), &
+                        integ_spline(zdeps,              zdept(ji,jjd,jk1), &
                                asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
                                csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
@@ -912,12 +901,4 @@
                END DO
             
-               IF( zbhitns == 1 ) THEN
-                 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
-                 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), &
-                                                 bsp(ji,jjd,1),    csp(ji,jjd,1), &
-                                                 dsp(ji,jjd,1) ) * zdeps
-                 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0)        ! no lighter than fresh water
-                 zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
-               ENDIF
 
                ! update the momentum trends in v direction
@@ -941,5 +922,5 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 
-      CALL wrk_dealloc( jpi,jpj,jpk, zdeptht, zrhh ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 
       !
    END SUBROUTINE hpg_prj
@@ -1121,5 +1102,5 @@
 
    
-   FUNCTION integ2(xl, xr, a, b, c, d)  RESULT(f) 
+   FUNCTION integ_spline(xl, xr, a, b, c, d)  RESULT(f) 
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE interp1  ***
@@ -1143,7 +1124,8 @@
          & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
 
-   END FUNCTION integ2
+   END FUNCTION integ_spline
 
 
    !!======================================================================
 END MODULE dynhpg
+
