Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3895)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3896)
@@ -24,4 +24,5 @@
    USE phycst           ! physical constants
    USE dom_oce          ! ocean domain
+   USE domvvl           ! ocean vertical scale factors
    USE dom_ice_2        ! LIM-2: ice domain
    USE ice_2            ! LIM-2: ice variables
@@ -59,4 +60,5 @@
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/LIM2 4.0 , UCL - NEMO Consortium (2011)
@@ -444,4 +446,6 @@
       !!-------------------------------------------------------------------
       !
+      INTEGER :: jk           ! local integer
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'lim_sbc_init_2 : LIM-2 sea-ice - surface boundary condition'
@@ -475,4 +479,35 @@
          sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
          sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+         do jk = 1,jpkm1                     ! adjust initial vertical scale factors
+          fse3t_n(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshn(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+          fse3t_b(:,:,jk) = e3t_0(:,:,jk)*( 1._wp + sshb(:,:)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) )
+         end do
+         fse3t_a(:,:,:) = fse3t_b(:,:,:)
+         ! Reconstruction of all vertical scale factors at now and before time steps
+         ! =============================================================================
+         ! Horizontal scale factor interpolations
+         ! --------------------------------------
+         CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
+         CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
+         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
+         ! Vertical scale factor interpolations
+         ! ------------------------------------
+         CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
+         CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
+         CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
+         CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
+         CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
+         ! t- and w- points depth
+         ! ----------------------
+         fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+         fsdepw_n(:,:,1) = 0.0_wp
+         fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
+         DO jk = 2, jpk
+            fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
+            fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
+            fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
+         END DO
       ENDIF
       !
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 3895)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 3896)
@@ -85,4 +85,6 @@
          IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
          IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
+         un_td = 0.0_wp
+         vn_td = 0.0_wp
       ENDIF
       IF( ln_vvl_ztilde ) THEN
@@ -158,6 +160,6 @@
       ! t- and w- points depth
       ! ----------------------
-      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
-      fsdepw_n(:,:,1) = 0.e0
+      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.0_wp
       fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
       DO jk = 2, jpk
@@ -168,7 +170,7 @@
       ! Reference water column height at t-, u- and v- point
       ! ----------------------------------------------------
-      ht_0(:,:) = 0.e0
-      hu_0(:,:) = 0.e0
-      hv_0(:,:) = 0.e0
+      ht_0(:,:) = 0.0_wp
+      hu_0(:,:) = 0.0_wp
+      hv_0(:,:) = 0.0_wp
       DO jk = 1, jpk
          ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
@@ -181,10 +183,10 @@
       IF( ln_vvl_ztilde ) THEN
          ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
-         frq_rst_e3t(:,:) = 2.e0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
-         frq_rst_hdv(:,:) = 2.e0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
+         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
+         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
          IF( ln_vvl_ztilde_as_zstar ) THEN
             ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
-            frq_rst_e3t(:,:) = 0.e0_wp 
-            frq_rst_hdv(:,:) = 1.e0_wp / rdt
+            frq_rst_e3t(:,:) = 0.0_wp 
+            frq_rst_hdv(:,:) = 1.0_wp / rdt
          ENDIF
       ENDIF
@@ -304,6 +306,6 @@
          ! 3 - Thickness diffusion term
          ! ----------------------------
-         zwu(:,:) = 0.e0
-         zwv(:,:) = 0.e0
+         zwu(:,:) = 0.0_wp
+         zwv(:,:) = 0.0_wp
          ! a - first derivative: diffusive fluxes
          DO jk = 1, jpkm1
@@ -346,5 +348,5 @@
             z2dt =  rdt
          ELSE
-            z2dt = 2.e0 * rdt
+            z2dt = 2.0_wp * rdt
          ENDIF
          CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
@@ -353,5 +355,5 @@
          ! Maximum deformation control
          ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         ze3t(:,:,jpk) = 0.e0
+         ze3t(:,:,jpk) = 0.0_wp
          DO jk = 1, jpkm1
             ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
@@ -421,5 +423,5 @@
          END IF
          !
-         zht(:,:) = 0.e0
+         zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
             zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
@@ -427,7 +429,7 @@
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
-         !
-         zht(:,:) = 0.e0
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
+         !
+         zht(:,:) = 0.0_wp
          DO jk = 1, jpkm1
             zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
@@ -435,6 +437,25 @@
          z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
          IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
-         IF( lwp    ) WRITE(numout, *) 'MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
-         !
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
+         !
+         zht(:,:) = 0.0_wp
+         DO jk = 1, jpkm1
+            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
+         END DO
+         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
+         !
+         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
+         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
+         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
       END IF
 
@@ -502,5 +523,5 @@
             tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
          ELSE
-            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) + atfp * ( tilde_e3t_b(:,:,:) - 2.e0 * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
+            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
          ENDIF
          tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
@@ -525,6 +546,6 @@
       ! t- and w- points depth
       ! ----------------------
-      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
-      fsdepw_n(:,:,1) = 0.e0
+      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
+      fsdepw_n(:,:,1) = 0.0_wp
       fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
       DO jk = 2, jpk
@@ -595,7 +616,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
-                     &                    * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
-                     &                        + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
+                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
+                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
+                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
                END DO
             END DO
@@ -613,7 +634,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5 * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
-                     &                    * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                        + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
+                  pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
+                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
+                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
                END DO
             END DO
@@ -631,7 +652,7 @@
             DO jj = 1, jpjm1
                DO ji = 1, fs_jpim1   ! vector opt.
-                  pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
-                     &                    * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
-                     &                        + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
+                  pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
+                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
+                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
                END DO
             END DO
@@ -649,6 +670,6 @@
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
-               &                              +      0.5 * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
+               &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
          END DO
          !               ! -------------------------------------- !
@@ -659,6 +680,6 @@
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
-               &                               +      0.5 * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
+               &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
          END DO
          !               ! -------------------------------------- !
@@ -669,6 +690,6 @@
          ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
          DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
-               &                               +      0.5 * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
+            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
+               &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
          END DO
       END SELECT
@@ -733,6 +754,6 @@
                   CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
                ELSE                            ! one at least array is missing
-                  tilde_e3t_b(:,:,:) = 0.e0
-                  tilde_e3t_n(:,:,:) = 0.e0
+                  tilde_e3t_b(:,:,:) = 0.0_wp
+                  tilde_e3t_n(:,:,:) = 0.0_wp
                ENDIF
                !                          ! ------------ !
@@ -742,5 +763,5 @@
                      CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
                   ELSE                ! array is missing
-                     hdiv_lf(:,:,:) = 0.e0
+                     hdiv_lf(:,:,:) = 0.0_wp
                   ENDIF
                ENDIF
@@ -751,7 +772,7 @@
             fse3t_n(:,:,:) = e3t_0(:,:,:)
             IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
-               tilde_e3t_b(:,:,:) = 0.e0
-               tilde_e3t_n(:,:,:) = 0.e0
-               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.e0
+               tilde_e3t_b(:,:,:) = 0.0_wp
+               tilde_e3t_n(:,:,:) = 0.0_wp
+               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
             END IF
          ENDIF
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 3895)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90	(revision 3896)
@@ -572,4 +572,5 @@
       INTEGER  ::   ierr                                          ! local integer
       REAL(wp) ::   zfac12, zua, zva                              ! local scalars
+      REAL(wp) ::   zmsk, ze3                                     ! local scalars
       !                                                           !  3D workspace 
       REAL(wp), POINTER    , DIMENSION(:,:  )         :: zwx, zwy, zwz
@@ -577,6 +578,5 @@
 #if defined key_vvl
       REAL(wp), POINTER    , DIMENSION(:,:,:)         :: ze3f     !  3D workspace (lk_vvl=T)
-#endif
-#if ! defined key_vvl
+#else
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE   :: ze3f     ! lk_vvl=F, ze3f=1/e3f saved one for all
 #endif
@@ -604,11 +604,13 @@
       ENDIF
 
-      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t)
+      IF( kt == nit000 .OR. lk_vvl ) THEN      ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points)
          DO jk = 1, jpk
             DO jj = 1, jpjm1
                DO ji = 1, jpim1
-                  ze3f(ji,jj,jk) = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
-                     &             + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) ) * 0.25
-                  IF( ze3f(ji,jj,jk) /= 0._wp )   ze3f(ji,jj,jk) = 1._wp / ze3f(ji,jj,jk)
+                  ze3  = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
+                     &   + fse3t(ji,jj  ,jk)*tmask(ji,jj  ,jk) + fse3t(ji+1,jj  ,jk)*tmask(ji+1,jj  ,jk) )
+                  zmsk = (                   tmask(ji,jj+1,jk) +                     tmask(ji+1,jj+1,jk)   &
+                     &                     + tmask(ji,jj  ,jk) +                     tmask(ji+1,jj  ,jk) )
+                  IF( ze3 /= 0._wp )   ze3f(ji,jj,jk) = zmsk / ze3
                END DO
             END DO
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 3895)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90	(revision 3896)
@@ -110,5 +110,5 @@
       ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
       ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
-      z1_rau0 = 0.5 / rau0
+      z1_rau0 = 0.5_wp * r1_rau0
       ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
 
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 3895)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 3896)
@@ -23,5 +23,4 @@
    USE eosbn2          ! equation of state            (eos bn2 routine)
    USE trdmld_oce      ! ocean active mixed layer tracers trends variables
-   USE domvvl          ! variable volume
    USE divcur          ! hor. divergence and curl      (div & cur routines)
 
