Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 4279)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 4280)
@@ -25,4 +25,6 @@
 
    PUBLIC   dom_rea    ! routine called by inidom.F90
+  !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OFF 3.3 , NEMO Consortium (2010)
@@ -173,8 +175,8 @@
             CALL iom_get( inum4, jpdom_unknown, 'esigw', esigw )
 
-            CALL iom_get( inum4, jpdom_data, 'e3t', e3t ) ! scale factors
-            CALL iom_get( inum4, jpdom_data, 'e3u', e3u )
-            CALL iom_get( inum4, jpdom_data, 'e3v', e3v )
-            CALL iom_get( inum4, jpdom_data, 'e3w', e3w )
+            CALL iom_get( inum4, jpdom_data, 'e3t', fse3t_n(:,:,:) ) ! scale factors
+            CALL iom_get( inum4, jpdom_data, 'e3u', fse3u_n(:,:,:) )
+            CALL iom_get( inum4, jpdom_data, 'e3v', fse3v_n(:,:,:) )
+            CALL iom_get( inum4, jpdom_data, 'e3w', fse3w_n(:,:,:) )
 
             CALL iom_get( inum4, jpdom_unknown, 'gdept_1d', gdept_1d ) ! depth
@@ -190,8 +192,8 @@
             !
             IF( nmsh <= 6 ) THEN                                        ! 3D vertical scale factors
-               CALL iom_get( inum4, jpdom_data, 'e3t', e3t )
-               CALL iom_get( inum4, jpdom_data, 'e3u', e3u )
-               CALL iom_get( inum4, jpdom_data, 'e3v', e3v )
-               CALL iom_get( inum4, jpdom_data, 'e3w', e3w )
+               CALL iom_get( inum4, jpdom_data, 'e3t', fse3t_n(:,:,:) )
+               CALL iom_get( inum4, jpdom_data, 'e3u', fse3u_n(:,:,:) )
+               CALL iom_get( inum4, jpdom_data, 'e3v', fse3v_n(:,:,:) )
+               CALL iom_get( inum4, jpdom_data, 'e3w', fse3w_n(:,:,:) )
             ELSE                                                        ! 2D bottom scale factors
                CALL iom_get( inum4, jpdom_data, 'e3t_ps', e3tp )
@@ -199,14 +201,14 @@
                !                                                        ! deduces the 3D scale factors
                DO jk = 1, jpk
-                  e3t(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
-                  e3u(:,:,jk) = e3t_1d(jk)
-                  e3v(:,:,jk) = e3t_1d(jk)
-                  e3w(:,:,jk) = e3w_1d(jk)
+                  fse3t_n(:,:,jk) = e3t_1d(jk)                                    ! set to the ref. factors
+                  fse3u_n(:,:,jk) = e3t_1d(jk)
+                  fse3v_n(:,:,jk) = e3t_1d(jk)
+                  fse3w_n(:,:,jk) = e3w_1d(jk)
                END DO
                DO jj = 1,jpj                                                  ! adjust the deepest values
                   DO ji = 1,jpi
                      ik = mbkt(ji,jj)
-                     e3t(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
-                     e3w(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(1) * ( 1._wp - tmask(ji,jj,1) )
+                     fse3t_n(ji,jj,ik) = e3tp(ji,jj) * tmask(ji,jj,1) + e3t_1d(1) * ( 1._wp - tmask(ji,jj,1) )
+                     fse3w_n(ji,jj,ik) = e3wp(ji,jj) * tmask(ji,jj,1) + e3w_1d(1) * ( 1._wp - tmask(ji,jj,1) )
                   END DO
                END DO
@@ -214,21 +216,21 @@
                   DO jj = 1, jpjm1
                      DO ji = 1, jpim1
-                        e3u(ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
-                        e3v(ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
+                        fse3u_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji+1,jj,jk) )
+                        fse3v_n(ji,jj,jk) = MIN( fse3t_n(ji,jj,jk), fse3t_n(ji,jj+1,jk) )
                      END DO
                   END DO
                END DO
-               CALL lbc_lnk( e3u , 'U', 1._wp )   ;   CALL lbc_lnk( e3uw, 'U', 1._wp )   ! lateral boundary conditions
-               CALL lbc_lnk( e3v , 'V', 1._wp )   ;   CALL lbc_lnk( e3vw, 'V', 1._wp )
+               CALL lbc_lnk( fse3u_n(:,:,:) , 'U', 1._wp )   ;   CALL lbc_lnk( fse3uw_n(:,:,:), 'U', 1._wp )   ! lateral boundary conditions
+               CALL lbc_lnk( fse3v_n(:,:,:) , 'V', 1._wp )   ;   CALL lbc_lnk( fse3vw_n(:,:,:), 'V', 1._wp )
                !
                DO jk = 1, jpk                        ! set to z-scale factor if zero (i.e. along closed boundaries)
-                  WHERE( e3u(:,:,jk) == 0._wp )   e3u(:,:,jk) = e3t_1d(jk)
-                  WHERE( e3v(:,:,jk) == 0._wp )   e3v(:,:,jk) = e3t_1d(jk)
+                  WHERE( fse3u_n(:,:,jk) == 0._wp )   fse3u_n(:,:,jk) = e3t_1d(jk)
+                  WHERE( fse3v_n(:,:,jk) == 0._wp )   fse3v_n(:,:,jk) = e3t_1d(jk)
                END DO
             END IF
 
             IF( iom_varid( inum4, 'gdept', ldstop = .FALSE. ) > 0 ) THEN   ! 3D depth of t- and w-level
-               CALL iom_get( inum4, jpdom_data, 'gdept', gdept )
-               CALL iom_get( inum4, jpdom_data, 'gdepw', gdepw )
+               CALL iom_get( inum4, jpdom_data, 'gdept', fsdept_n(:,:,:) )
+               CALL iom_get( inum4, jpdom_data, 'gdepw', fsdepw_n(:,:,:) )
             ELSE                                                           ! 2D bottom depth
                CALL iom_get( inum4, jpdom_data, 'hdept', zprt )
@@ -236,6 +238,6 @@
                !
                DO jk = 1, jpk                                              ! deduces the 3D depth
-                  gdept(:,:,jk) = gdept_1d(jk)
-                  gdepw(:,:,jk) = gdepw_1d(jk)
+                  fsdept_n(:,:,jk) = gdept_1d(jk)
+                  fsdepw_n(:,:,jk) = gdepw_1d(jk)
                END DO
                DO jj = 1, jpj
@@ -243,7 +245,7 @@
                      ik = mbkt(ji,jj)
                      IF( ik > 0 ) THEN
-                        gdepw(ji,jj,ik+1) = zprw(ji,jj)
-                        gdept(ji,jj,ik  ) = zprt(ji,jj)
-                        gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
+                        fsdepw_n(ji,jj,ik+1) = zprw(ji,jj)
+                        fsdept_n(ji,jj,ik  ) = zprt(ji,jj)
+                        fsdept_n(ji,jj,ik+1) = fsdept_n(ji,jj,ik) + fse3t_n(ji,jj,ik)
                      ENDIF
                   END DO
@@ -259,10 +261,10 @@
             CALL iom_get( inum4, jpdom_unknown, 'e3w_1d'  , e3w_1d   )
             DO jk = 1, jpk
-               e3t  (:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
-               e3u  (:,:,jk) = e3t_1d(jk)
-               e3v  (:,:,jk) = e3t_1d(jk)
-               e3w  (:,:,jk) = e3w_1d(jk)
-               gdept(:,:,jk) = gdept_1d(jk)
-               gdepw(:,:,jk) = gdepw_1d(jk)
+               fse3t_n(:,:,jk) = e3t_1d(jk)                              ! set to the ref. factors
+               fse3u_n(:,:,jk) = e3t_1d(jk)
+               fse3v_n(:,:,jk) = e3t_1d(jk)
+               fse3w_n(:,:,jk) = e3w_1d(jk)
+               fsdept_n(:,:,jk) = gdept_1d(jk)
+               fsdepw_n(:,:,jk) = gdepw_1d(jk)
             END DO
          ENDIF
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4279)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4280)
@@ -56,5 +56,5 @@
 
    ! Potential namelist parameters below to be read in dyn_spg_ts_init
-   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.TRUE.        !: Forward integration of barotropic sub-stepping
+   LOGICAL,  PUBLIC,  PARAMETER :: ln_bt_fw=.true.        !: Forward integration of barotropic sub-stepping
    LOGICAL,  PRIVATE, PARAMETER :: ln_bt_av=.TRUE.        !: Time averaging of barotropic variables
    LOGICAL,  PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE.  !: Set number of iterations automatically
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 4279)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/PISCES/trcwri_pisces.F90	(revision 4280)
@@ -21,4 +21,5 @@
    PUBLIC trc_wri_pisces 
 
+#  include "top_substitute.h90"
 CONTAINS
 
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/trcdta.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 4279)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 4280)
@@ -194,5 +194,5 @@
                   DO ji = 1, jpi
                      DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                        zl = fsdept_0(ji,jj,jk)
+                        zl = fsdept_n(ji,jj,jk)
                         IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
                            ztp(jk) =  ptrc(ji,jj,1    ,jn)
@@ -226,5 +226,5 @@
                         ik = mbkt(ji,jj) 
                         IF( ik > 1 ) THEN
-                           zl = ( gdept_1d(ik) - fsdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                           zl = ( gdept_1d(ik) - fsdept_n(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
                            ptrc(ji,jj,ik,jn) = (1.-zl) * ptrc(ji,jj,ik,jn) + zl * ptrc(ji,jj,ik-1,jn)
                         ENDIF
Index: /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
===================================================================
--- /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 4279)
+++ /branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 4280)
@@ -79,6 +79,6 @@
            !
            sshb_hold  (:,:) = sshn  (:,:)
-           sshu_b_hold(:,:) = sshu_n(:,:)
-           sshv_b_hold(:,:) = sshv_n(:,:)
+!!Z~       sshu_b_hold(:,:) = sshu_n(:,:)
+!!Z~       sshv_b_hold(:,:) = sshv_n(:,:)
            emp_b_hold (:,:) = emp_b (:,:)
            !
@@ -117,6 +117,6 @@
           !
           sshn_tm  (:,:)         = sshn_tm  (:,:)         + sshn  (:,:) 
-          sshu_n_tm(:,:)         = sshu_n_tm(:,:)         + sshu_n(:,:) 
-          sshv_n_tm(:,:)         = sshv_n_tm(:,:)         + sshv_n(:,:) 
+!!Z~      sshu_n_tm(:,:)         = sshu_n_tm(:,:)         + sshu_n(:,:) 
+!!Z~      sshv_n_tm(:,:)         = sshv_n_tm(:,:)         + sshv_n(:,:) 
           rnf_tm   (:,:)         = rnf_tm   (:,:)         + rnf   (:,:) 
           h_rnf_tm (:,:)         = h_rnf_tm (:,:)         + h_rnf (:,:) 
@@ -197,13 +197,13 @@
 # endif
          sshn_temp  (:,:)        = sshn  (:,:)
-         sshu_n_temp(:,:)        = sshu_n(:,:)
-         sshv_n_temp(:,:)        = sshv_n(:,:)
-         sshf_n_temp(:,:)        = sshf_n(:,:)
+!!Z~     sshu_n_temp(:,:)        = sshu_n(:,:)
+!!Z~     sshv_n_temp(:,:)        = sshv_n(:,:)
+!!Z~     sshf_n_temp(:,:)        = sshf_n(:,:)
          sshb_temp  (:,:)        = sshb  (:,:)
-         sshu_b_temp(:,:)        = sshu_b(:,:)
-         sshv_b_temp(:,:)        = sshv_b(:,:)
+!!Z~     sshu_b_temp(:,:)        = sshu_b(:,:)
+!!Z~     sshv_b_temp(:,:)        = sshv_b(:,:)
          ssha_temp  (:,:)        = ssha  (:,:)
-         sshu_a_temp(:,:)        = sshu_a(:,:)
-         sshv_a_temp(:,:)        = sshv_a(:,:)
+!!Z~     sshu_a_temp(:,:)        = sshu_a(:,:)
+!!Z~     sshv_a_temp(:,:)        = sshv_a(:,:)
          rnf_temp   (:,:)        = rnf   (:,:)
          h_rnf_temp (:,:)        = h_rnf (:,:)
@@ -309,6 +309,6 @@
 # endif
          sshn_tm  (:,:)          = sshn_tm    (:,:)       + sshn  (:,:) 
-         sshu_n_tm(:,:)          = sshu_n_tm  (:,:)       + sshu_n(:,:) 
-         sshv_n_tm(:,:)          = sshv_n_tm  (:,:)       + sshv_n(:,:) 
+!!Z~     sshu_n_tm(:,:)          = sshu_n_tm  (:,:)       + sshu_n(:,:) 
+!!Z~     sshv_n_tm(:,:)          = sshv_n_tm  (:,:)       + sshv_n(:,:) 
          rnf_tm   (:,:)          = rnf_tm     (:,:)       + rnf   (:,:) 
          h_rnf_tm (:,:)          = h_rnf_tm   (:,:)       + h_rnf (:,:) 
@@ -321,9 +321,9 @@
          !
          sshn     (:,:)          = sshn_tm    (:,:) * r1_ndttrcp1 
-         sshu_n   (:,:)          = sshu_n_tm  (:,:) * r1_ndttrcp1  
-         sshv_n   (:,:)          = sshv_n_tm  (:,:) * r1_ndttrcp1  
+!!Z~     sshu_n   (:,:)          = sshu_n_tm  (:,:) * r1_ndttrcp1  
+!!Z~     sshv_n   (:,:)          = sshv_n_tm  (:,:) * r1_ndttrcp1  
          sshb     (:,:)          = sshb_hold  (:,:)
-         sshu_b   (:,:)          = sshu_b_hold(:,:)
-         sshv_b   (:,:)          = sshv_b_hold(:,:)
+!!Z~     sshu_b   (:,:)          = sshu_b_hold(:,:)
+!!Z~     sshv_b   (:,:)          = sshv_b_hold(:,:)
          rnf      (:,:)          = rnf_tm     (:,:) * r1_ndttrcp1 
          h_rnf    (:,:)          = h_rnf_tm   (:,:) * r1_ndttrcp1 
@@ -486,13 +486,13 @@
 #endif
          CALL lbc_lnk( sshn  (:,:)         , 'T', 1. ) 
-         CALL lbc_lnk( sshu_n(:,:)         , 'U', 1. ) 
-         CALL lbc_lnk( sshv_n(:,:)         , 'V', 1. ) 
-         CALL lbc_lnk( sshf_n(:,:)         , 'F', 1. ) 
+!!Z~     CALL lbc_lnk( sshu_n(:,:)         , 'U', 1. ) 
+!!Z~     CALL lbc_lnk( sshv_n(:,:)         , 'V', 1. ) 
+!!Z~     CALL lbc_lnk( sshf_n(:,:)         , 'F', 1. ) 
          CALL lbc_lnk( sshb  (:,:)         , 'T', 1. ) 
-         CALL lbc_lnk( sshu_b(:,:)         , 'U', 1. ) 
-         CALL lbc_lnk( sshv_b(:,:)         , 'V', 1. ) 
+!!Z~     CALL lbc_lnk( sshu_b(:,:)         , 'U', 1. ) 
+!!Z~     CALL lbc_lnk( sshv_b(:,:)         , 'V', 1. ) 
          CALL lbc_lnk( ssha  (:,:)         , 'T', 1. ) 
-         CALL lbc_lnk( sshu_a(:,:)         , 'U', 1. ) 
-         CALL lbc_lnk( sshv_a(:,:)         , 'V', 1. ) 
+!!Z~     CALL lbc_lnk( sshu_a(:,:)         , 'U', 1. ) 
+!!Z~     CALL lbc_lnk( sshv_a(:,:)         , 'V', 1. ) 
          CALL lbc_lnk( rnf   (:,:)         , 'T', 1. ) 
          CALL lbc_lnk( h_rnf (:,:)         , 'T', 1. ) 
@@ -592,6 +592,6 @@
 #endif
       sshn_tm  (:,:) = sshn  (:,:) 
-      sshu_n_tm(:,:) = sshu_n(:,:) 
-      sshv_n_tm(:,:) = sshv_n(:,:) 
+!!Z~  sshu_n_tm(:,:) = sshu_n(:,:) 
+!!Z~  sshv_n_tm(:,:) = sshv_n(:,:) 
       rnf_tm   (:,:) = rnf   (:,:) 
       h_rnf_tm (:,:) = h_rnf (:,:) 
@@ -695,11 +695,11 @@
       sshb  (:,:)     =  sshb_temp  (:,:)
       ssha  (:,:)     =  ssha_temp  (:,:)
-      sshu_n(:,:)     =  sshu_n_temp(:,:)
-      sshu_b(:,:)     =  sshu_b_temp(:,:)
-      sshu_a(:,:)     =  sshu_a_temp(:,:)
-      sshv_n(:,:)     =  sshv_n_temp(:,:)
-      sshv_b(:,:)     =  sshv_b_temp(:,:)
-      sshv_a(:,:)     =  sshv_a_temp(:,:)
-      sshf_n(:,:)     =  sshf_n_temp(:,:)
+!!Z~  sshu_n(:,:)     =  sshu_n_temp(:,:)
+!!Z~  sshu_b(:,:)     =  sshu_b_temp(:,:)
+!!Z~  sshu_a(:,:)     =  sshu_a_temp(:,:)
+!!Z~  sshv_n(:,:)     =  sshv_n_temp(:,:)
+!!Z~  sshv_b(:,:)     =  sshv_b_temp(:,:)
+!!Z~  sshv_a(:,:)     =  sshv_a_temp(:,:)
+!!Z~  sshf_n(:,:)     =  sshf_n_temp(:,:)
       rnf   (:,:)     =  rnf_temp   (:,:)
       h_rnf (:,:)     =  h_rnf_temp (:,:)
@@ -816,10 +816,10 @@
       !
       sshb_hold  (:,:) = sshn  (:,:)
-      sshu_b_hold(:,:) = sshu_n(:,:)
-      sshv_b_hold(:,:) = sshv_n(:,:)
+!!Z~  sshu_b_hold(:,:) = sshu_n(:,:)
+!!Z~  sshv_b_hold(:,:) = sshv_n(:,:)
       emp_b_hold (:,:) = emp   (:,:)
       sshn_tm    (:,:) = sshn  (:,:) 
-      sshu_n_tm  (:,:) = sshu_n(:,:) 
-      sshv_n_tm  (:,:) = sshv_n(:,:) 
+!!Z~  sshu_n_tm  (:,:) = sshu_n(:,:) 
+!!Z~  sshv_n_tm  (:,:) = sshv_n(:,:) 
       rnf_tm     (:,:) = rnf   (:,:) 
       h_rnf_tm   (:,:) = h_rnf (:,:) 
@@ -951,6 +951,6 @@
          END DO
          !
-         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
-         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
+!!Z~     hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
+!!Z~     hv(:,:) = hv_0(:,:) + sshv_n(:,:)
          !                                            ! now masked inverse of the ocean depth (at u- and v-points)
          hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
@@ -992,17 +992,17 @@
 
       !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
-      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1      ! NO Vector Opt.
-               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
-                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
-                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
-               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
-                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
-                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
-            END DO
-         END DO
-         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
-      ENDIF
+!!Z~  IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
+!!Z~     DO jj = 1, jpjm1
+!!Z~        DO ji = 1, jpim1      ! NO Vector Opt.
+!!Z~           sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
+!!Z~              &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
+!!Z~              &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
+!!Z~           sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
+!!Z~              &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
+!!Z~              &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
+!!Z~        END DO
+!!Z~     END DO
+!!Z~     CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
+!!Z~  ENDIF
       
 
