Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4296)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4297)
@@ -38,5 +38,5 @@
    USE cpl_oasis3, ONLY : lk_cpl
    USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget
-   USE oce,        ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
+   USE oce,        ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
    USE dom_ice,    ONLY : tms
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -434,45 +434,38 @@
             sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
             sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-            !
-            ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
-            !       which were previously set in domvvl
-            IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
-                     zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-                     zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-                     zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-                     sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                     sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                     sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
-                     sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
-                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
-                  END DO
-               END DO
-               CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
-               CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-               DO jj = 1, jpjm1
-                  DO ji = 1, jpim1      ! NO Vector Opt.
-                     sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
-                          &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                          &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                          &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
-                  END DO
-               END DO
-               CALL lbc_lnk( sshf_n, 'F', 1. )
-            ENDIF
+            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
       ENDIF ! .NOT. ln_rstart
-      !
-!!?      IF( .NOT. ln_rstart ) THEN           ! delete the initial ssh below sea-ice area
-!!?         !
-!!?         zarea     = glob_sum( e1e2t(:,:) )           ! interior global domain surface
-!!?         zsum      = glob_sum( e1e2t(:,:) * ( snwice_mass(:,:) ) ) / zarea * r1_rau0
-!!?         sshn(:,:) = sshn(:,:) - zsum 
-!!?         sshb(:,:) = sshb(:,:) - zsum
-!!?      ENDIF
       !
 
