Index: branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1399)
+++ branches/dev_004_VVL/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 1405)
@@ -1,6 +1,3 @@
 MODULE dynspg_ts
-   !!======================================================================
-   !!                   ***  MODULE  dynspg_ts  ***
-   !! Ocean dynamics:  surface pressure gradient trend
    !!======================================================================
    !! History :   9.0  !  04-12  (L. Bessieres, G. Madec)  Original code
@@ -85,9 +82,6 @@
       !!          (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b 
       !!          and zvX_b (X specifying after, now or before).
-      !!      -3- Update of sea surface height from time averaged barotropic 
-      !!          variables.
-      !!        - apply lateral boundary conditions on sshn.
-      !!      -4- The new general trend becomes :
-      !!          ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H
+      !!      -3- The new general trend becomes :
+      !!          ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H
       !!
       !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
@@ -107,17 +101,14 @@
          zcu, zcv, zwx, zwy, zhdiv,         &  ! temporary arrays
          zua, zva, zub, zvb,                &  !     "        "
-         zssha_b, zua_b, zva_b,             &  !     "        "
-         zub_e, zvb_e,                      &  !     "        "
-         zun_e, zvn_e                          !     "        "
+         zua_e, zva_e, zssha_e,             &  !     "        "
+         zub_e, zvb_e, zsshb_e,             &  !     "        "
+         zu_sum, zv_sum, zssh_sum
       !! Variable volume
       REAL(wp), DIMENSION(jpi,jpj)     ::   &
          zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfse3un_e, zfse3vn_e  ! 3D workspace
       !!----------------------------------------------------------------------
 
       ! Arrays initialization
       ! ---------------------
-      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0
-      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0
       zhdiv(:,:) = 0.e0
 
@@ -133,10 +124,9 @@
          !                               ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b
 
-         ssha_e(:,:) = sshn(:,:)
-         ua_e(:,:)   = un_b(:,:)
-         va_e(:,:)   = vn_b(:,:)
-         hu_e(:,:)   = hu(:,:)
-         hv_e(:,:)   = hv(:,:)
-
+         zssha_e(:,:) = sshn(:,:)
+         zua_e  (:,:) = un_e(:,:)
+         zva_e  (:,:) = vn_e(:,:)
+         hu_e   (:,:) = hu  (:,:)
+         hv_e   (:,:) = hv  (:,:)
          IF( ln_dynvor_een ) THEN
             ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0
@@ -170,5 +160,5 @@
                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
          END DO
 
@@ -193,5 +183,5 @@
       zfact2 = 0.5 * 0.5
       zraur  = 1. / rauw                                    ! 1 / volumic mass of pure water
-      
+
       ! -----------------------------------------------------------------------------
       !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
@@ -215,6 +205,6 @@
                zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk)
                !                                                           ! Vertically integrated transports (before)
-               zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk)
-               zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk)
+               zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
+               zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
                !                                                           ! Planetary vorticity transport fluxes (now)
                zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk)
@@ -228,6 +218,6 @@
             zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
             !                                                           ! Vertically integrated transports (before)
-            zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk)
-            zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk)
+            zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
+            zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
             !                                                           ! Planetary vorticity (now)
             zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)
@@ -304,24 +294,32 @@
       !----------------
       ! Number of iteration of the barotropic loop
-      icycle = 3  * nn_baro / 2
+      icycle = 2  * nn_baro + 1
 
       ! variables for the barotropic equations
-      sshb_e (:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
-      sshn_e (:,:) = sshn_b(:,:)
-      zub_e  (:,:) = un_b  (:,:)       ! barotropic transports issued from the barotropic equations (before and now)
-      zvb_e  (:,:) = vn_b  (:,:)
-      zun_e  (:,:) = un_b  (:,:)
-      zvn_e  (:,:) = vn_b  (:,:)
-      zssha_b(:,:) = 0.e0
-      zua_b  (:,:) = 0.e0
-      zva_b  (:,:) = 0.e0
-      hu_e   (:,:) = hu   (:,:)     ! (barotropic) ocean depth at u-point
-      hv_e   (:,:) = hv   (:,:)     ! (barotropic) ocean depth at v-point
-      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
-         zfse3un_e(:,:,:) = fse3u(:,:,:)   ! (barotropic) scale factors  at u-point
-         zfse3un_e(:,:,:) = fse3v(:,:,:)   ! (barotropic) scale factors  at v-point
+      zu_sum  (:,:) = 0.e0
+      zv_sum  (:,:) = 0.e0
+      zssh_sum(:,:) = 0.e0
+      hu_e    (:,:) = hu    (:,:)      ! (barotropic) ocean depth at u-point
+      hv_e    (:,:) = hv    (:,:)      ! (barotropic) ocean depth at v-point
+      zsshb_e (:,:) = sshn_e(:,:)      ! (barotropic) sea surface height (before and now)
+      ! vertical sum
+      un_e  (:,:) = 0.e0
+      vn_e  (:,:) = 0.e0
+      IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
+         DO jk = 1, jpkm1
+            DO ji = 1, jpij
+               un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
+               vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
+            END DO
+         END DO
+      ELSE                             ! No  vector opt.
+         DO jk = 1, jpkm1
+            un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)
+            vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+         END DO
       ENDIF
+      zub_e (:,:) = un_e(:,:)
+      zvb_e (:,:) = vn_e(:,:)
+
 
       ! set ssh corrections to 0
@@ -352,8 +350,8 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * zun_e(ji  ,jj)              &
-                  &            -e2u(ji-1,jj  ) * zun_e(ji-1,jj)              &
-                  &            +e1v(ji  ,jj  ) * zvn_e(ji  ,jj)              &
-                  &            -e1v(ji  ,jj-1) * zvn_e(ji  ,jj-1) )          &
+               zhdiv(ji,jj) = ( e2u(ji  ,jj  ) * un_e(ji  ,jj)              &
+                  &            -e2u(ji-1,jj  ) * un_e(ji-1,jj)              &
+                  &            +e1v(ji  ,jj  ) * vn_e(ji  ,jj)              &
+                  &            -e1v(ji  ,jj-1) * vn_e(ji  ,jj-1) )          &
                   &           / (e1t(ji,jj)*e2t(ji,jj))
             END DO
@@ -370,9 +368,9 @@
 
 #if defined key_bdy
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
-               END DO
-            END DO
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
+            END DO
+         END DO
 #endif
 
@@ -381,6 +379,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ssha_e(ji,jj) = ( sshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
-            &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
+               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
+                  &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
             END DO
          END DO
@@ -388,6 +386,6 @@
          ! evolution of the barotropic transport ( following the vorticity scheme used)
          ! ----------------------------------------------------------------------------
-         zwx(:,:) = e2u(:,:) * zun_e(:,:)
-         zwy(:,:) = e1v(:,:) * zvn_e(:,:)
+         zwx(:,:) = e2u(:,:) * un_e(:,:)
+         zwy(:,:) = e1v(:,:) * vn_e(:,:)
 
          IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
@@ -396,8 +394,8 @@
                   ! surface pressure gradient
                   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)
+                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) &
+                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(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(ji,jj) / e2v(ji,jj)
                   ELSE
                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
@@ -412,6 +410,6 @@
                   zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
                   ! after transports
-                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
@@ -422,8 +420,8 @@
                   ! surface pressure gradient
                   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)
+                     zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  ) &
+                        &            - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(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(ji,jj) / e2v(ji,jj)
                   ELSE
                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
@@ -432,12 +430,12 @@
                   ! enstrophy conserving formulation for planetary vorticity term
                   zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
-                                 + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
+                     + zwy(ji  ,jj  ) + zwy(ji+1,jj  ) ) / e1u(ji,jj)
                   zx1 =-zfact1 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1)   &
-                                 + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
+                     + zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj)
                   zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) )
                   zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
                   ! after transports
-                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
@@ -449,8 +447,8 @@
                   ! surface pressure gradient
                   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)
+                     zspgu = -grav * ( ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  ) &
+                        &            - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  ) ) * hu(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(ji,jj) / e2v(ji,jj)
                   ELSE
                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
@@ -463,6 +461,6 @@
                      &                             + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) )
                   ! after transports
-                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
-                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
+                  zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)
+                  zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)
                END DO
             END DO
@@ -477,35 +475,33 @@
 
          ! ... Boundary conditions on ua_e, va_e, ssha_e
-         CALL lbc_lnk( ua_e  , 'U', -1. )
-         CALL lbc_lnk( va_e  , 'V', -1. )
-         CALL lbc_lnk( ssha_e, 'T',  1. )
+         CALL lbc_lnk( zua_e  , 'U', -1. )
+         CALL lbc_lnk( zva_e  , 'V', -1. )
+         CALL lbc_lnk( zssha_e, 'T',  1. )
 
          ! temporal sum
          !-------------
-         IF( jit >= nn_baro / 2 ) THEN
-            zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:)
-            zua_b  (:,:) = zua_b  (:,:) + ua_e  (:,:)
-            zva_b  (:,:) = zva_b  (:,:) + va_e  (:,:) 
-         ENDIF
+         zu_sum  (:,:) = zu_sum  (:,:) + zua_e  (:,:)
+         zv_sum  (:,:) = zv_sum  (:,:) + zva_e  (:,:) 
+         zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 
 
          ! Time filter and swap of dynamics arrays
          ! ---------------------------------------
          IF( jit == 1 ) THEN   ! Euler (forward) time stepping
-            sshb_e (:,:) = sshn_e(:,:)
-            zub_e  (:,:) = zun_e (:,:)
-            zvb_e  (:,:) = zvn_e (:,:)
-            sshn_e (:,:) = ssha_e(:,:)
-            zun_e  (:,:) = ua_e  (:,:)
-            zvn_e  (:,:) = va_e  (:,:)
+            zsshb_e(:,:) = sshn_e (:,:)
+            zub_e  (:,:) = un_e   (:,:)
+            zvb_e  (:,:) = vn_e   (:,:)
+            sshn_e (:,:) = zssha_e(:,:)
+            un_e   (:,:) = zua_e  (:,:)
+            vn_e   (:,:) = zva_e  (:,:)
          ELSE                                        ! Asselin filtering
-            sshb_e (:,:) = atfp * ( sshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
-            zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:)
-            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:)
-            sshn_e (:,:) = ssha_e(:,:)
-            zun_e  (:,:) = ua_e  (:,:)
-            zvn_e  (:,:) = va_e  (:,:)
+            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:)
+            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * un_e  (:,:)
+            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * vn_e  (:,:)
+            sshn_e (:,:) = zssha_e(:,:)
+            un_e   (:,:) = zua_e  (:,:)
+            vn_e   (:,:) = zva_e  (:,:)
          ENDIF
 
-         IF( lk_vvl ) THEN ! Variable volume
+         IF( lk_vvl ) THEN ! Variable volume   ! needed for BDY ???
 
             ! Sea surface elevation
@@ -528,27 +524,23 @@
 
             ! Boundaries conditions
-            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;    CALL lbc_lnk( zsshvn_e, 'V', 1. )
+            CALL lbc_lnk( zsshun_e, 'U', 1. )    ;  CALL lbc_lnk( zsshvn_e, 'V', 1. )
 
             ! 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
-
+            ! -------------------------------------------
+            hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:)
+            hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:)
+
+            !
+         ENDIF
          !                                                 ! ==================== !
       END DO                                               !        end loop      !
       !                                                    ! ==================== !
 
-
       ! Time average of after barotropic variables
-      zcoef =  1.e0 / ( nn_baro + 1 ) 
-      zssha_b(:,:) = zcoef * zssha_b(:,:) 
-      zua_b  (:,:) = zcoef *  zua_b (:,:) 
-      zva_b  (:,:) = zcoef *  zva_b (:,:) 
+      zcoef =  1.e0 / ( 2 * nn_baro  + 1 ) 
+      un_e  (:,:) = zcoef *  zu_sum  (:,:) 
+      vn_e  (:,:) = zcoef *  zv_sum  (:,:) 
+      sshn_e(:,:) = zcoef *  zssh_sum(:,:) 
+
 #if defined key_obc
       IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
@@ -557,29 +549,17 @@
       IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
 #endif
-     
-
-      ! ---------------------------------------------------------------------------
-      ! Phase 3 : Update sea surface height from time averaged barotropic variables
-      ! ---------------------------------------------------------------------------
-!RB_vvl now done in ssh_wzv and ssh_nxt
 
       ! -----------------------------------------------------------------------------
-      ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)
+      ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step)
       ! -----------------------------------------------------------------------------
 
-      ! Swap on time averaged barotropic variables
-      ! ------------------------------------------
-      sshb_b(:,:) = sshn_b (:,:)
-      IF ( neuler == 0 .AND. kt == nit000 ) zssha_b(:,:) = sshn(:,:)
-      sshn_b(:,:) = zssha_b(:,:)
-      un_b  (:,:) = zua_b  (:,:) 
-      vn_b  (:,:) = zva_b  (:,:) 
-   
+
+
       ! add time averaged barotropic coriolis and surface pressure gradient
       ! terms to the general momentum trend
       ! --------------------------------------------------------------------
       DO jk=1,jpkm1
-         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b
-         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b
+         ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b
+         va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b
       END DO
 
@@ -588,6 +568,4 @@
       IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' )
 
-      ! print sum trends (used for debugging)
-      IF( ln_ctl )     CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
       !
    END SUBROUTINE dyn_spg_ts
@@ -607,35 +585,31 @@
       !
       IF( TRIM(cdrw) == 'READ' ) THEN
-         IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) )   ! free surface issued
-            CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) )   ! from time-splitting loop
-            CALL iom_get( numror, jpdom_autoglo, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
-            CALL iom_get( numror, jpdom_autoglo, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
-            IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:)
+         IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN
+            CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) )   ! free surface and
+            CALL iom_get( numror, jpdom_autoglo, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued
+            CALL iom_get( numror, jpdom_autoglo, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop
          ELSE
-            sshb_b(:,:) = sshb(:,:)
-            sshn_b(:,:) = sshn(:,:)
-            un_b  (:,:) = 0.e0
-            vn_b  (:,:) = 0.e0
+            sshn_e(:,:) = sshn(:,:)
+            un_e  (:,:) = 0.e0
+            vn_e  (:,:) = 0.e0
             ! vertical sum
             IF( lk_vopt_loop ) THEN          ! vector opt., forced unroll
                DO jk = 1, jpkm1
                   DO ji = 1, jpij
-                     un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
-                     vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
+                     un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
+                     vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
                   END DO
                END DO
             ELSE                             ! No  vector opt.
                DO jk = 1, jpkm1
-                  un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
-                  vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
+                  un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)
+                  vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
                END DO
             ENDIF
          ENDIF
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
-         CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) )   ! free surface issued
-         CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) )   ! from barotropic loop
-         CALL iom_rstput( kt, nitrst, numrow, 'un_b'  , un_b  (:,:) )   ! horizontal transports issued
-         CALL iom_rstput( kt, nitrst, numrow, 'vn_b'  , vn_b  (:,:) )   ! from barotropic loop
+         CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) )   ! free surface and
+         CALL iom_rstput( kt, nitrst, numrow, 'un_e'  , un_e  (:,:) )   ! horizontal transports issued
+         CALL iom_rstput( kt, nitrst, numrow, 'vn_e'  , vn_e  (:,:) )   ! from barotropic loop
       ENDIF
       !
