Index: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 575)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 592)
@@ -34,4 +34,5 @@
    USE iom
    USE restart         ! only for lrst_oce
+   USE domvvl          ! variable volume
 
    IMPLICIT NONE
@@ -102,4 +103,8 @@
          zsshb_e, zub_e, zvb_e,             &  !     "        "
          zun_e, zvn_e                          !     "        "
+      !! Variable volume
+      REAL(wp), DIMENSION(jpi,jpj)     ::   &
+         zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e         ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfse3un_e, zfse3vn_e  ! 3D workspace
       !!----------------------------------------------------------------------
 
@@ -137,4 +142,37 @@
          ENDIF
          !
+      ENDIF
+
+      ! Substract the surface pressure gradient from the momentum
+      ! ---------------------------------------------------------
+      IF( lk_vvl ) THEN    ! Variable volume
+
+         ! 0. Local initialization
+         ! -----------------------
+         zspgu_1(:,:) = 0.e0
+         zspgv_1(:,:) = 0.e0
+
+         ! 1. Surface pressure gradient (now)
+         ! ----------------------------
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  ) &
+                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e1u(ji,jj)
+               zspgv_1(ji,jj) = -grav * ( ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1) &
+                  &                     - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  ) ) / e2v(ji,jj)
+            END DO 
+         END DO
+
+         ! 2. Substract the surface pressure trend from the general trend
+         ! ------------------------------------------------------
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj)
+                  va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj)
+               END DO
+            END DO
+         END DO
+
       ENDIF
 
@@ -269,4 +307,12 @@
       zua_b  (:,:) = un_b  (:,:)   
       zva_b  (:,:) = vn_b  (:,:)
+      IF( lk_vvl ) THEN
+         zsshun_e(:,:)    = sshu (:,:)     ! (barotropic) sea surface height (now) at u-point
+         zsshvn_e(:,:)    = sshv (:,:)     ! (barotropic) sea surface height (now) at v-point
+         hu_e(:,:)        = hu   (:,:)     ! (barotropic) ocean depth at u-point
+         hv_e(:,:)        = hv   (:,:)     ! (barotropic) ocean depth at v-point
+         zfse3un_e(:,:,:) = fse3u(:,:,:)   ! (barotropic) scale factors  at u-point
+         zfse3un_e(:,:,:) = fse3v(:,:,:)   ! (barotropic) scale factors  at v-point
+      ENDIF
 
       ! set ssh corrections to 0
@@ -330,6 +376,13 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  IF( lk_vvl) THEN
+                     zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) &
+                        &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)
+                     zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &
+                        &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)
+                  ELSE
+                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
+                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  ENDIF
                   ! energy conserving formulation for planetary vorticity term
                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
@@ -349,6 +402,13 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  IF( lk_vvl) THEN
+                     zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) &
+                        &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)
+                     zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &
+                        &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)
+                  ELSE
+                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
+                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  ENDIF
                   ! enstrophy conserving formulation for planetary vorticity term
                   zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
@@ -369,6 +429,13 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  IF( lk_vvl) THEN
+                     zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) &
+                        &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)
+                     zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &
+                        &  - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)
+                  ELSE
+                     zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
+                     zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  ENDIF
                   ! energy/enstrophy conserving formulation for planetary vorticity term
                   zcubt = + zfac25 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   &
@@ -403,5 +470,5 @@
          ! Time filter and swap of dynamics arrays
          ! ---------------------------------------
-         IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
+         IF( neuler == 0 .AND. jit == 1 ) THEN   ! Euler (forward) time stepping
             zsshb_e(:,:) = sshn_e(:,:)
             zub_e  (:,:) = zun_e (:,:)
@@ -419,4 +486,45 @@
          ENDIF
 
+         IF( lk_vvl ) THEN ! Variable volume
+
+            ! Sea surface elevation
+            ! ---------------------
+            DO jj = 1, jpjm1
+               DO ji = 1,jpim1
+
+                  ! Sea Surface Height at u-point before
+                  zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )   &
+                     &              * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn_e(ji  ,jj)    &
+                     &                + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
+
+                  ! Sea Surface Height at v-point before
+                  zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )   &
+                     &              * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn_e(ji,jj  )    &
+                     &                + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
+
+               END DO
+            END DO
+
+            ! Boundaries conditions
+            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;    CALL lbc_lnk( zsshvn_e, 'V', 1. )
+
+            ! Scale factors at before and after time step
+            ! -------------------------------------------
+            DO jk = 1, jpkm1
+               zfse3un_e(:,:,jk)  = fsve3u(:,:,jk) * ( 1 + zsshun_e(:,:) * muu(:,:,jk) )
+               zfse3vn_e(:,:,jk)  = fsve3v(:,:,jk) * ( 1 + zsshvn_e(:,:) * muv(:,:,jk) )
+            END DO
+
+            ! Ocean depth at U- and V-points
+            hu_e(:,:) = 0.e0
+            hv_e(:,:) = 0.e0
+
+            DO jk = 1, jpk
+               hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk)
+               hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk)
+            END DO
+
+         ENDIF ! End variable volume case
+
          !                                                 ! ==================== !
       END DO                                               !        end loop      !
@@ -444,13 +552,15 @@
       ! Horizontal divergence of time averaged barotropic transports
       !-------------------------------------------------------------
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
-           &                +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
-           &             / ( e1t(ji,jj) * e2t(ji,jj) )
-         END DO
-      END DO
-
-#if defined key_obc
+      IF( .NOT. lk_vvl ) THEN
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj  ) * un_b(ji-1,jj  )     &
+                  &            +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji  ,jj-1) * vn_b(ji  ,jj-1) )   &
+                  &           / ( e1t(ji,jj) * e2t(ji,jj) )
+            END DO
+         END DO
+      ENDIF
+
+#if defined key_obc && ! defined key_vvl
       ! open boundaries (div must be zero behind the open boundary)
       !  mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
@@ -463,6 +573,12 @@
       ! sea surface height
       !-------------------
-      sshb(:,:) = sshn(:,:)
-      sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      IF( lk_vvl ) THEN
+         sshbb(:,:) = sshb(:,:)
+         sshb (:,:) = sshn(:,:)
+         sshn (:,:) = ssha(:,:)
+      ELSE
+         sshb (:,:) = sshn(:,:)
+         sshn(:,:) = (  sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) )  ) * tmask(:,:,1)
+      ENDIF
 
       ! ... Boundary conditions on sshn
@@ -477,4 +593,5 @@
       ! ------------------------------------------
       sshb_b(:,:) = sshn_b (:,:)
+      IF ( kt == nit000 ) zssha_b(:,:) = sshn(:,:)
       sshn_b(:,:) = zssha_b(:,:)
       un_b  (:,:) = zua_b  (:,:) 
