Index: trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 363)
+++ trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 367)
@@ -12,5 +12,9 @@
    USE dom_oce         ! ocean space and time domain
    USE in_out_manager  ! I/O manager
+   USE obc_oce         ! ocean open boundary conditions
    USE obcdyn          ! open boundary condition for momentum (obc_dyn routine)
+   USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
+   USE obcvol          ! ocean open boundary condition (obc_vol routines)
+   USE dynspg_oce      ! type of surface pressure gradient
    USE lbclnk          ! lateral boundary condition (or mpp link)
    USE prtctl          ! Print control
@@ -105,4 +109,24 @@
       ! Update (ua,va) along open boundaries (only in the rigid-lid case)
       CALL obc_dyn( kt )
+
+      IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN
+         !Flather boundary condition :
+         !        - Update sea surface height on each open boundary
+         !                 sshn (= after ssh) for explicit case
+         !                 sshn_b (= after ssha_b) for time-splitting case
+         !        - Correct the barotropic velocities
+         CALL obc_dyn_bt( kt )
+
+         !Boundary conditions on sshn ( after ssh)
+         CALL lbc_lnk( sshn, 'T', 1. )
+
+         IF(ln_ctl) THEN         ! print sum trends (used for debugging)
+            CALL prt_ctl(tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask)
+         ENDIF
+
+         IF ( ln_vol_cst ) CALL obc_vol( kt )
+
+      ENDIF
+
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
Index: trunk/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 363)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 367)
@@ -12,4 +12,6 @@
    USE oce            ! ocean dynamics and tracers variables
    USE dom_oce        ! ocean space and time domain variables
+   USE obc_oce        ! ocean open boundary conditions
+   USE dynspg_oce     ! surface pressure gradient variables
    USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
    USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
@@ -27,26 +29,4 @@
    !! * Accessibility
    PUBLIC dyn_spg         ! routine called by step module
-
-   !! * Public variables
-#if   defined key_dynspg_exp   ||  defined key_esopa
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .TRUE.  !: Explicit free surface flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .FALSE. !: Explicit free surface flag
-#endif
-#if   defined key_dynspg_ts   ||  defined key_esopa
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .TRUE.  !: Free surface with time splitting flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .FALSE. !: Free surface with time splitting flag
-#endif
-#if   defined key_dynspg_flt  ||  defined key_esopa
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .TRUE.  !: Filtered free surface cst volume flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .FALSE. !: Filtered free surface cst volume flag
-#endif
-#if   defined key_dynspg_rl
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl  = .TRUE.  !: Rigid-lid flag
-#else
-   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl  = .FALSE. !: Rigid-lid flag
-#endif
 
    !! * module variables
@@ -208,4 +188,16 @@
       ENDIF
 
+#if key_obc
+      ! Conservation of ocean volume (key_dynspg_flt)
+      ! ---------------------------------------------
+      IF( lk_dynspg_flt ) ln_vol_cst = .true.
+
+      ! Application of Flather's algorithm at open boundaries
+      ! -----------------------------------------------------
+      IF( lk_dynspg_flt ) ln_obc_fla = .false.
+      IF( lk_dynspg_exp ) ln_obc_fla = .true.
+      IF( lk_dynspg_ts  ) ln_obc_fla = .true.
+#endif
+
    END SUBROUTINE dyn_spg_ctl
 
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 363)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90	(revision 367)
@@ -20,4 +20,5 @@
    USE ocesbc          ! ocean surface boundary condition
    USE obc_oce         ! Lateral open boundary condition
+   USE obcdta          ! open boundary condition data     (obc_dta_bt routine)
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -101,6 +102,8 @@
       ENDIF
 
-      ! 0. Local constant initialization
-      ! --------------------------------
+      ! 0. Initialization
+      ! -----------------
+      ! read or estimate sea surface height and vertically integrated velocities
+      IF( lk_obc )   CALL obc_dta_bt( kt, 0 )
       z2dt = 2. * rdt                                       ! time step: leap-frog
       IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt       ! time step: Euler if restart from rest
@@ -179,6 +182,6 @@
       ENDIF
 
-      !Boundary conditions on sshn
-      CALL lbc_lnk( sshn, 'T', 1. )
+      ! Boundary conditions on sshn
+      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
  
       IF(ln_ctl) THEN         ! print sum trends (used for debugging)
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 367)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 367)
@@ -0,0 +1,54 @@
+MODULE dynspg_oce
+   !!----------------------------------------------------------------------
+   !!                       ***  MODULE dynspg_oce  ***
+   !!       
+   !! ** Purpose :   Define in memory all the ocean space domain variables
+   !!----------------------------------------------------------------------
+   !! Modules used
+   USE par_oce          ! ocean parameters
+
+   IMPLICIT NONE
+   PUBLIC           
+   !!----------------------------------------------------------------------
+   !!  OPA 9.0 , LOCEAN-IPSL (2005) 
+   !! $Header$ 
+   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !!----------------------------------------------------------------------
+
+   !! Surface pressure gradient logicals
+   !! ----------------------------------
+#if   defined key_dynspg_exp   ||  defined key_esopa
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .TRUE.  !: Explicit free surface flag
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_exp = .FALSE. !: Explicit free surface flag
+#endif
+#if   defined key_dynspg_ts   ||  defined key_esopa
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .TRUE.  !: Free surface with time splitting flag
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_ts  = .FALSE. !: Free surface with time splitting flag
+#endif
+#if   defined key_dynspg_flt  ||  defined key_esopa
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .TRUE.  !: Filtered free surface cst volume flag
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_flt = .FALSE. !: Filtered free surface cst volume flag
+#endif
+#if   defined key_dynspg_rl
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl  = .TRUE.  !: Rigid-lid flag
+#else
+   LOGICAL, PUBLIC, PARAMETER ::   lk_dynspg_rl  = .FALSE. !: Rigid-lid flag
+#endif
+
+#if   defined key_dynspg_ts   ||  defined key_esopa
+   !! Time splitting variables
+   !! ------------------------
+      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop
+         sshn_b, sshb_b,               &  ! sea surface heigth (now, before)
+         un_b  , vn_b                     ! vertically integrated horizontal velocities (now)
+      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop
+         sshn_e, ssha_e,               &  ! sea surface heigth (now,after)
+         ua_e  , va_e                     ! vertically integrated horizontal velocities (after)
+#endif
+
+   !!----------------------------------------------------------------------
+
+END MODULE dynspg_oce
Index: trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 363)
+++ trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 367)
@@ -17,4 +17,6 @@
    USE phycst          ! physical constants
    USE ocesbc          ! ocean surface boundary condition
+   USE obcdta          ! open boundary condition data     
+   USE obcfla          ! Flather open boundary condition  
    USE dynvor          ! vorticity term
    USE obc_oce         ! Lateral open boundary condition
@@ -22,4 +24,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
+   USE dynspg_oce      ! surface pressure gradient variables
    USE in_out_manager  ! I/O manager
 
@@ -29,9 +32,4 @@
    !! * Accessibility
    PUBLIC dyn_spg_ts  ! routine called by step.F90
-
-   !! * Module variables
-      REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop
-         sshn_b, sshb_b,               &  ! sea surface heigth (now, before)
-         un_b  , vn_b                     ! vertically integrated horizontal velocities (now)
 
    !! * Substitutions
@@ -61,6 +59,6 @@
       !!          surface gradient and the Coriolis force are updated within
       !!          the barotropic integration.
-      !!      -2- Barotropic loop : updates of sea surface height (zssha_e) and 
-      !!          barotropic transports (zua_e and zva_e) through barotropic 
+      !!      -2- Barotropic loop : updates of sea surface height (ssha_e) and 
+      !!          barotropic transports (ua_e and va_e) through barotropic 
       !!          momentum and continuity integration. Barotropic former 
       !!          variables are time averaging over the full barotropic cycle
@@ -97,6 +95,5 @@
          zssha_b, zua_b, zva_b,             &  !     "        "
          zsshb_e, zub_e, zvb_e,             &  !     "        "
-         zsshn_e, zun_e, zvn_e,             &  !     "        "
-         zssha_e, zua_e, zva_e                 !     "        "
+         zun_e, zvn_e                          !     "        "
       REAL(wp), DIMENSION(jpi,jpj),SAVE ::  &
          ztnw, ztne, ztsw, ztse
@@ -105,6 +102,6 @@
       ! Arrays initialization
       ! ---------------------
-      zua_b(:,:) = 0.e0   ;   zub_e(:,:) = 0.e0   ;   zun_e(:,:) = 0.e0   ;   zua_e(:,:) = 0.e0
-      zva_b(:,:) = 0.e0   ;   zvb_e(:,:) = 0.e0   ;   zvn_e(:,:) = 0.e0   ;   zva_e(:,:) = 0.e0
+      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
 
@@ -138,7 +135,7 @@
             ENDIF
          ENDIF
-         zssha_e(:,:) = sshn(:,:)
-         zua_e  (:,:) = un_b(:,:)
-         zva_e  (:,:) = vn_b(:,:)
+         ssha_e(:,:) = sshn(:,:)
+         ua_e(:,:)   = un_b(:,:)
+         va_e(:,:)   = vn_b(:,:)
 
          IF( ln_dynvor_een ) THEN
@@ -278,12 +275,21 @@
       ! variables for the barotropic equations
       zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
-      zsshn_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(:,:) = sshn(:,:)        ! time averaged variables over all sub-timesteps
-      zua_b(:,:) = un_b(:,:)   
-      zva_b(:,:) = vn_b(:,:)
+      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(:,:) = sshn  (:,:)       ! time averaged variables over all sub-timesteps
+      zua_b  (:,:) = un_b  (:,:)   
+      zva_b  (:,:) = vn_b  (:,:)
+
+      ! set ssh corrections to 0
+      ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
+#if defined key_obc
+      IF( lp_obc_east  )   sshfoe_b(:,:) = 0.e0
+      IF( lp_obc_west  )   sshfow_b(:,:) = 0.e0
+      IF( lp_obc_south )   sshfos_b(:,:) = 0.e0
+      IF( lp_obc_north )   sshfon_b(:,:) = 0.e0
+#endif
 
       ! Barotropic integration over 2 baroclinic time steps
@@ -296,4 +302,7 @@
          z2dt_e = 2. * rdtbt
          IF ( jit == 1 )   z2dt_e = rdtbt
+
+         ! Time interpolation of open boundary condition data
+         IF( lk_obc )   CALL obc_dta_bt( kt, jit )
 
          ! Horizontal divergence of barotropic transports
@@ -312,8 +321,8 @@
          ! 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
-         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
-         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
-         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
-         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
+         IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0      ! east
+         IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0      ! west
+         IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0      ! north
+         IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0      ! south
 #endif
 
@@ -322,5 +331,5 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
+               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
             &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
             END DO
@@ -336,6 +345,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  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)
                   ! energy conserving formulation for planetary vorticity term
                   zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
@@ -346,6 +355,6 @@
                   zcvbt =-zfact2 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 )
                   ! after transports
-                  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)
+                  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)
                END DO
             END DO
@@ -355,6 +364,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  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)
                   ! enstrophy conserving formulation for planetary vorticity term
                   zy1 = zfact1 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1)   &
@@ -365,6 +374,6 @@
                   zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) )
                   ! after transports
-                  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)
+                  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)
                END DO
             END DO
@@ -375,6 +384,6 @@
                DO ji = fs_2, fs_jpim1   ! vector opt.
                   ! surface pressure gradient
-                  zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)
-                  zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)
+                  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)
                   ! energy/enstrophy conserving formulation for planetary vorticity term
                   zcubt = + zfac25 / e1u(ji,jj) * (  ztne(ji,jj  ) * zwy(ji  ,jj  ) + ztnw(ji+1,jj) * zwy(ji+1,jj  )   &
@@ -383,37 +392,43 @@
                      &                             + ztnw(ji,jj  ) * zwx(ji-1,jj  ) + ztne(ji,jj  ) * zwx(ji  ,jj  ) )
                   ! after transports
-                  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)
+                  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)
                END DO
             END DO
          ENDIF
 
-         ! ... Boundary conditions on zua_e, zva_e, zssha_e
-         CALL lbc_lnk( zua_e, 'U', -1. )
-         CALL lbc_lnk( zva_e, 'V', -1. )
-         CALL lbc_lnk( zssha_e, 'T', 1. )
+         ! Flather's boundary condition for the barotropic loop :
+         !         - Update sea surface height on each open boundary
+         !         - Correct the barotropic transports
+         IF( lk_obc )   CALL obc_fla_ts( kt )
+
+
+         ! ... 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. )
 
          ! temporal sum
          !-------------
-         zssha_b(:,:) = zssha_b(:,:) + zssha_e(:,:)
-         zua_b  (:,:) = zua_b  (:,:) + zua_e  (:,:)
-         zva_b  (:,:) = zva_b  (:,:) + zva_e  (:,:) 
+         zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:)
+         zua_b  (:,:) = zua_b  (:,:) + ua_e  (:,:)
+         zva_b  (:,:) = zva_b  (:,:) + va_e  (:,:) 
 
          ! Time filter and swap of dynamics arrays
          ! ---------------------------------------
          IF( neuler == 0 .AND. kt == nit000 ) THEN   ! Euler (forward) time stepping
-            zsshb_e(:,:) = zsshn_e(:,:)
-            zub_e  (:,:) = zun_e  (:,:)
-            zvb_e  (:,:) = zvn_e  (:,:)
-            zsshn_e(:,:) = zssha_e(:,:)
-            zun_e  (:,:) = zua_e  (:,:)
-            zvn_e  (:,:) = zva_e  (:,:)
+            zsshb_e(:,:) = sshn_e(:,:)
+            zub_e  (:,:) = zun_e (:,:)
+            zvb_e  (:,:) = zvn_e (:,:)
+            sshn_e (:,:) = ssha_e(:,:)
+            zun_e  (:,:) = ua_e  (:,:)
+            zvn_e  (:,:) = va_e  (:,:)
          ELSE                                        ! Asselin filtering
-            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * zsshn_e(:,:)
-            zub_e  (:,:) = atfp * ( zub_e  (:,:) + zua_e  (:,:) ) + atfp1 * zun_e  (:,:)
-            zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + zva_e  (:,:) ) + atfp1 * zvn_e  (:,:)
-            zsshn_e(:,:) = zssha_e(:,:)
-            zun_e  (:,:) = zua_e  (:,:)
-            zvn_e  (:,:) = zva_e  (:,:)
+            zsshb_e(:,:) = atfp * ( zsshb_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  (:,:)
          ENDIF
 
@@ -426,6 +441,12 @@
       zcoef =  1.e0 / (  FLOAT( icycle +1 )  )
       zssha_b(:,:) = zcoef * zssha_b(:,:) 
-      zua_b  (:,:) = zcoef * zua_b  (:,:) 
-      zva_b  (:,:) = zcoef * zva_b  (:,:) 
+      zua_b  (:,:) = zcoef *  zua_b (:,:) 
+      zva_b  (:,:) = zcoef *  zva_b (:,:) 
+#if defined key_obc
+         IF( lp_obc_east  )   sshfoe_b(:,:) = zcoef * sshfoe_b(:,:)
+         IF( lp_obc_west  )   sshfow_b(:,:) = zcoef * sshfow_b(:,:)
+         IF( lp_obc_north )   sshfon_b(:,:) = zcoef * sshfon_b(:,:)
+         IF( lp_obc_south )   sshfos_b(:,:) = zcoef * sshfos_b(:,:)
+#endif
      
 
@@ -448,8 +469,8 @@
       ! 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
-      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1) = 0.e0      ! east
-      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1) = 0.e0      ! west
+      IF( lp_obc_east  )   zhdiv(nie0p1:nie1p1,nje0  :nje1)   = 0.e0    ! east
+      IF( lp_obc_west  )   zhdiv(niw0  :niw1  ,njw0  :njw1)   = 0.e0    ! west
       IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
-      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1) = 0.e0      ! south
+      IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
 #endif
 
@@ -460,5 +481,5 @@
 
       ! ... Boundary conditions on sshn
-      CALL lbc_lnk( sshn, 'T', 1. )
+      IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )
 
 
