Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 4354)
@@ -120,10 +120,4 @@
    REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary
 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh                  !: 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur                  !: 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr                  !: Pointers for barotropic fields 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   pub2d, pun2d, pua2d   !: 
-   REAL(wp), POINTER, DIMENSION(:,:)           ::   pvb2d, pvn2d, pva2d   !: 
-
    !!----------------------------------------------------------------------
    !! open boundary data variables
@@ -134,7 +128,8 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy)
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy)
+!$AGRIF_DO_NOT_TREAT
    TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process)
    TYPE(OBC_DATA) , DIMENSION(jp_bdy), TARGET      ::   dta_bdy           !: bdy external data (local process)
-
+!$AGRIF_END_DO_NOT_TREAT
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -153,5 +148,10 @@
       ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),     &  
          &      STAT=bdy_oce_alloc )
-         !
+      !
+      ! Initialize masks 
+      bdytmask(:,:) = 1._wp
+      bdyumask(:,:) = 1._wp
+      bdyvmask(:,:) = 1._wp
+      ! 
       IF( lk_mpp             )   CALL mpp_sum ( bdy_oce_alloc )
       IF( bdy_oce_alloc /= 0 )   CALL ctl_warn('bdy_oce_alloc: failed to allocate arrays.')
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 4354)
@@ -20,5 +20,4 @@
    !!    bdy_dta_init   : initialise arrays etc for reading of external data
    !!----------------------------------------------------------------------
-   USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
    USE oce             ! ocean dynamics and tracers
@@ -51,7 +50,7 @@
    LOGICAL,           DIMENSION(jp_bdy) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions
                                                                ! =F => baroclinic velocities in 3D boundary conditions
-
+!$AGRIF_DO_NOT_TREAT
    TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read)
-
+!$AGRIF_END_DO_NOT_TREAT
    TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap
 
@@ -103,14 +102,4 @@
          ! Calculate depth-mean currents
          !-----------------------------
-         CALL wrk_alloc(jpi,jpj,pun2d,pvn2d) 
-
-         pun2d(:,:) = 0.e0
-         pvn2d(:,:) = 0.e0
-         DO ik = 1, jpkm1   !! Vertically integrated momentum trends
-             pun2d(:,:) = pun2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik)
-             pvn2d(:,:) = pvn2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik)
-         END DO
-         pun2d(:,:) = pun2d(:,:) * hur(:,:)
-         pvn2d(:,:) = pvn2d(:,:) * hvr(:,:)
          
          DO ib_bdy = 1, nb_bdy
@@ -135,5 +124,5 @@
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%u2d(ib) = pun2d(ii,ij) * umask(ii,ij,1)         
+                     dta_bdy(ib_bdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
                   END DO 
                END IF
@@ -143,5 +132,5 @@
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                      ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                     dta_bdy(ib_bdy)%v2d(ib) = pvn2d(ii,ij) * vmask(ii,ij,1)         
+                     dta_bdy(ib_bdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
                   END DO 
                END IF
@@ -156,5 +145,5 @@
                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pun2d(ii,ij) ) * umask(ii,ij,ik)         
+                        dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
                      END DO
                   END DO 
@@ -166,5 +155,5 @@
                         ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                         ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                        dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pvn2d(ii,ij) ) * vmask(ii,ij,ik)         
+                        dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
                         END DO
                   END DO 
@@ -262,5 +251,4 @@
          ENDDO ! ib_bdy
 
-         CALL wrk_dealloc(jpi,jpj,pun2d,pvn2d) 
 
       ENDIF ! kt .eq. nit000
@@ -919,2 +907,4 @@
    !!==============================================================================
 END MODULE bdydta
+
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 4354)
@@ -60,5 +60,6 @@
       LOGICAL               :: ll_dyn2d, ll_dyn3d, ll_orlanski
       !!
-      REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1     ! inverse depth at u and v points
+      REAL(wp), POINTER, DIMENSION(:,:) :: pua2d, pva2d     ! after barotropic velocities
+      REAL(wp), POINTER, DIMENSION(:,:) :: phura, phvra     ! after inverse depth at u and v points
 
       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn')
@@ -81,9 +82,5 @@
       !-------------------------------------------------------
 
-      pssh => sshn
-      phur => hur
-      phvr => hvr
-      CALL wrk_alloc(jpi,jpj,pua2d,pva2d) 
-      IF ( ll_orlanski ) CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 
+      CALL wrk_alloc(jpi,jpj,pua2d,pva2d,phura,phvra) 
 
       !-------------------------------------------------------
@@ -97,24 +94,24 @@
       
       IF (lk_vvl) THEN
-         phur(:,:) = 0.
-         phvr(:,:) = 0.
+         phura(:,:) = 0.
+         phvra(:,:) = 0.
          DO jk = 1, jpkm1
-            phur(:,:) = phur(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
-            phvr(:,:) = phvr(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
+            phura(:,:) = phura(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
+            phvra(:,:) = phvra(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
             pua2d(:,:) = pua2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
             pva2d(:,:) = pva2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
          END DO
-         phur(:,:) = umask(:,:,1) / ( phur(:,:) + 1. - umask(:,:,1) )
-         phvr(:,:) = vmask(:,:,1) / ( phvr(:,:) + 1. - vmask(:,:,1) )
-         pua2d(:,:) = pua2d(:,:) * phur(:,:)
-         pva2d(:,:) = pva2d(:,:) * phvr(:,:)
+         phura(:,:) = umask(:,:,1) / ( phura(:,:) + 1. - umask(:,:,1) )
+         phvra(:,:) = vmask(:,:,1) / ( phvra(:,:) + 1. - vmask(:,:,1) )
       ELSE
+         phura(:,:) = hur(:,:)
+         phvra(:,:) = hvr(:,:)
          DO jk = 1, jpkm1
             pua2d(:,:) = pua2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
             pva2d(:,:) = pva2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
          END DO
-         pua2d(:,:) = pua2d(:,:) * phur(:,:)
-         pva2d(:,:) = pva2d(:,:) * phvr(:,:)
       ENDIF
+      pua2d(:,:) = pua2d(:,:) * phura(:,:)
+      pva2d(:,:) = pva2d(:,:) * phvra(:,:)
 
       DO jk = 1 , jpkm1
@@ -126,32 +123,7 @@
 
       IF ( ll_orlanski ) THEN          
-         pub2d(:,:) = 0.e0
-         pvb2d(:,:) = 0.e0
-
-         IF (lk_vvl) THEN
-            phur1(:,:) = 0.
-            phvr1(:,:) = 0.
-            DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-               phur1(:,:) = phur1(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
-               phvr1(:,:) = phvr1(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
-               pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) * ub(:,:,jk)
-               pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk)
-            END DO
-            phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) )
-            phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) )
-            pub2d(:,:) = pub2d(:,:) * phur1(:,:)
-            pvb2d(:,:) = pvb2d(:,:) * phvr1(:,:)
-         ELSE
-            DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-               pub2d(:,:) = pub2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ub(:,:,jk)
-               pvb2d(:,:) = pvb2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * vb(:,:,jk)
-            END DO
-            pub2d(:,:) = pub2d(:,:) * phur(:,:)
-            pvb2d(:,:) = pvb2d(:,:) * phvr(:,:)
-         ENDIF
-
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = (ub(:,:,jk) - pub2d(:,:)) * umask(:,:,jk)
-            vb(:,:,jk) = (vb(:,:,jk) - pvb2d(:,:)) * vmask(:,:,jk)
+            ub(:,:,jk) = (ub(:,:,jk) - ub_b(:,:)) * umask(:,:,jk)
+            vb(:,:,jk) = (vb(:,:,jk) - vb_b(:,:)) * vmask(:,:,jk)
          END DO
       END IF
@@ -162,5 +134,5 @@
       !-------------------------------------------------------
 
-      IF( ll_dyn2d ) CALL bdy_dyn2d( kt )
+      IF( ll_dyn2d ) CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, phura, phvra, ssha )
 
       IF( ll_dyn3d ) CALL bdy_dyn3d( kt )
@@ -177,11 +149,10 @@
       IF ( ll_orlanski ) THEN
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = ( ub(:,:,jk) + pub2d(:,:) ) * umask(:,:,jk)
-            vb(:,:,jk) = ( vb(:,:,jk) + pvb2d(:,:) ) * vmask(:,:,jk)
+            ub(:,:,jk) = ( ub(:,:,jk) + ub_b(:,:) ) * umask(:,:,jk)
+            vb(:,:,jk) = ( vb(:,:,jk) + vb_b(:,:) ) * vmask(:,:,jk)
          END DO
       END IF
 
-      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d) 
-      IF ( ll_orlanski ) CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 
+      CALL wrk_dealloc(jpi,jpj,pua2d,pva2d,phura,phvra) 
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn')
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 4354)
@@ -41,5 +41,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn2d( kt )
+   SUBROUTINE bdy_dyn2d( kt, pua2d, pva2d, pub2d, pvb2d, phur, phvr, pssh  )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn2d  ***
@@ -49,4 +49,8 @@
       !!----------------------------------------------------------------------
       INTEGER,                      INTENT(in) ::   kt   ! Main time step counter
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pub2d, pvb2d
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: phur, phvr
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pssh
       !!
       INTEGER                                  ::   ib_bdy ! Loop counter
@@ -58,11 +62,13 @@
             CYCLE
          CASE('frs')
-            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d )
          CASE('flather')
-            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, pua2d, pva2d, pssh, phur, phvr )
          CASE('orlanski')
-            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
+            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
+                                     & pua2d, pva2d, pub2d, pvb2d, ll_npo=.false.)
          CASE('orlanski_npo')
-            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
+            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, &
+                                     & pua2d, pva2d, pub2d, pvb2d, ll_npo=.true. )
          CASE DEFAULT
             CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' )
@@ -72,5 +78,5 @@
    END SUBROUTINE bdy_dyn2d
 
-   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy )
+   SUBROUTINE bdy_dyn2d_frs( idx, dta, ib_bdy, pua2d, pva2d )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn2d_frs  ***
@@ -86,4 +92,5 @@
       TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
       INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 
       !!
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -118,5 +125,5 @@
 
 
-   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy )
+   SUBROUTINE bdy_dyn2d_fla( idx, dta, ib_bdy, pua2d, pva2d, pssh, phur, phvr )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn2d_fla  ***
@@ -140,4 +147,6 @@
       TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
       INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh, phur, phvr 
 
       INTEGER  ::   jb, igrd                         ! dummy loop indices
@@ -212,5 +221,5 @@
 
 
-   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo )
+   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, pua2d, pva2d, pub2d, pvb2d, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  ***
@@ -226,4 +235,6 @@
       TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
       INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub2d, pvb2d 
       LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version
 
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 4354)
@@ -15,5 +15,4 @@
    !!----------------------------------------------------------------------
    USE timing          ! Timing
-   USE wrk_nemo        ! Memory Allocation
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
@@ -266,5 +265,4 @@
       REAL(wp) ::   zwgt           ! boundary weight
       INTEGER  ::  ib_bdy          ! loop index
-      REAL(wp), POINTER, DIMENSION(:,:) :: phur1, phvr1     ! inverse depth at u and v points
       !!----------------------------------------------------------------------
       !
@@ -272,34 +270,4 @@
       !
       !-------------------------------------------------------
-      ! Remove barotropic part from before velocity
-      !-------------------------------------------------------
-      CALL wrk_alloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 
-
-      pub2d(:,:) = 0.e0
-      pvb2d(:,:) = 0.e0
-
-      phur1(:,:) = 0.
-      phvr1(:,:) = 0.
-      DO jk = 1, jpkm1
-#if defined key_vvl
-         phur1(:,:) = phur1(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
-         phvr1(:,:) = phvr1(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
-         pub2d(:,:) = pub2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
-         pvb2d(:,:) = pvb2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
-#else
-         pub2d(:,:) = pub2d(:,:) + fse3u(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
-         pvb2d(:,:) = pvb2d(:,:) + fse3v(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
-#endif
-      END DO
-
-      IF( lk_vvl ) THEN
-         phur1(:,:) = umask(:,:,1) / ( phur1(:,:) + 1. - umask(:,:,1) )
-         phvr1(:,:) = vmask(:,:,1) / ( phvr1(:,:) + 1. - vmask(:,:,1) )
-         pub2d(:,:) = pub2d(:,:) * umask(:,:,1) * phur1(:,:)
-         pvb2d(:,:) = pvb2d(:,:) * vmask(:,:,1) * phvr1(:,:)
-      ELSE
-         pub2d(:,:) = pvb2d(:,:) * hur(:,:)
-         pvb2d(:,:) = pub2d(:,:) * hvr(:,:)
-      ENDIF
 
       DO ib_bdy=1, nb_bdy
@@ -312,5 +280,5 @@
                DO jk = 1, jpkm1
                   ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
-                                   ub(ii,ij,jk) + pub2d(ii,ij)) ) * umask(ii,ij,jk)
+                                   ub(ii,ij,jk) + ub_b(ii,ij)) ) * umask(ii,ij,jk)
                END DO
             END DO
@@ -323,5 +291,5 @@
                DO jk = 1, jpkm1
                   va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
-                                   vb(ii,ij,jk) + pvb2d(ii,ij)) ) * vmask(ii,ij,jk)
+                                   vb(ii,ij,jk) + vb_b(ii,ij)) ) * vmask(ii,ij,jk)
                END DO
             END DO
@@ -329,6 +297,4 @@
       ENDDO
       !
-      CALL wrk_dealloc(jpi,jpj,pub2d,pvb2d,phur1,phvr1) 
-      !
       CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
       !
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 4354)
@@ -110,6 +110,4 @@
 
       IF( nn_timing == 1 ) CALL timing_start('bdy_init')
-
-      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
 
       IF(lwp) WRITE(numout,*)
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 4354)
@@ -51,5 +51,7 @@
    END TYPE TIDES_DATA
 
+!$AGRIF_DO_NOT_TREAT
    TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
+!$AGRIF_END_DO_NOT_TREAT
    TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
 
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90	(revision 4354)
@@ -71,5 +71,6 @@
       !!----------------------------------------------------------------------
       ! - ML - needed for initialization of e3t_b
-      INTEGER  ::  jk     ! dummy loop indice
+      INTEGER  ::  ji,jj,jk     ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zhur_b, zhvr_b ! U & Inverse of before depths
       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::  zuvd    ! U & V data workspace
       !!----------------------------------------------------------------------
@@ -155,4 +156,54 @@
       ENDIF
       !
+      ! 
+      ! Initialize "now" and "before" barotropic velocities:
+      ! Do it whatever the free surface method, these arrays
+      ! being eventually used
+      !
+      IF (lk_vvl) THEN 
+         CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b  )
+         zhur_b(:,:) = 0._wp
+         zhvr_b(:,:) = 0._wp
+         DO jk = 1, jpk
+            zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+            zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) )
+         zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) )
+      ENDIF
+      !
+      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
+      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
+      !
+      DO jk = 1, jpkm1
+#if defined key_vectopt_loop
+         DO jj = 1, 1         !Vector opt. => forced unrolling
+            DO ji = 1, jpij
+#else 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+#endif                  
+               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               !
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      un_b(:,:) = un_b(:,:) * hur(:,:)
+      vn_b(:,:) = vn_b(:,:) * hvr(:,:)
+      !
+      IF( lk_vvl ) THEN
+         ub_b(:,:) = ub_b(:,:) * zhur_b(:,:)
+         vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:)
+      ELSE
+         ub_b(:,:) = ub_b(:,:) * hur(:,:)
+         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
+      ENDIF
+      !
+      IF (lk_vvl) CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b  ) 
+      !
       IF( nn_timing == 1 )  CALL timing_stop('istate_init')
       !
@@ -537,2 +588,3 @@
    !!=====================================================================
 END MODULE istate
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 4354)
@@ -102,5 +102,5 @@
       REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zec   ! local scalars
       REAL(wp) ::   zve3a, zve3n, zve3b, zvf        !   -      -
-      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva, zhura, zhvra
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zua, zva, zhura, zhvra, zhurb, zhvrb
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  ze3u_f, ze3v_f 
       !!----------------------------------------------------------------------
@@ -339,5 +339,54 @@
          !
       ENDIF ! neuler =/0
-
+      !
+      ! Set "now" and "before" barotropic velocities for next time step:
+      ! JC: Would be more clever to swap variables than to make a full vertical
+      ! integration
+      !
+      IF (lk_vvl) THEN 
+         CALL wrk_alloc( jpi, jpj, zhurb, zhvrb  )
+         zhurb(:,:) = 0._wp
+         zhvrb(:,:) = 0._wp
+         DO jk = 1, jpk
+            zhurb(:,:) = zhurb(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
+            zhvrb(:,:) = zhvrb(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
+         END DO
+         zhurb(:,:) = umask(:,:,1) / ( zhurb(:,:) + 1._wp - umask(:,:,1) )
+         zhvrb(:,:) = vmask(:,:,1) / ( zhvrb(:,:) + 1._wp - vmask(:,:,1) )
+      ENDIF
+      !
+      un_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp
+      ub_b(:,:) = 0._wp ; vb_b(:,:) = 0._wp
+      !
+      DO jk = 1, jpkm1
+#if defined key_vectopt_loop
+         DO jj = 1, 1         !Vector opt. => forced unrolling
+            DO ji = 1, jpij
+#else 
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+#endif                  
+               un_b(ji,jj) = un_b(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)
+               vn_b(ji,jj) = vn_b(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
+               !
+               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
+               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !
+      un_b(:,:) = un_b(:,:) * hur(:,:)
+      vn_b(:,:) = vn_b(:,:) * hvr(:,:)
+      !
+      IF( lk_vvl ) THEN
+         ub_b(:,:) = ub_b(:,:) * zhurb(:,:)
+         vb_b(:,:) = vb_b(:,:) * zhvrb(:,:)
+      ELSE
+         ub_b(:,:) = ub_b(:,:) * hur(:,:)
+         vb_b(:,:) = vb_b(:,:) * hvr(:,:)
+      ENDIF
+      !
+      IF (lk_vvl) CALL wrk_dealloc( jpi, jpj, zhurb, zhvrb  ) 
+      !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
          &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
@@ -352,2 +401,3 @@
    !!=========================================================================
 END MODULE dynnxt
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 4354)
@@ -39,7 +39,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur_e , hvr_e    ! inverse of hu_e and hv_e
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_b           ! before field without time-filter
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_b, va_b     ! after  averaged velocities
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_b, vn_b     ! now    averaged velocities
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b, vb_b     ! before averaged velocities
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   un_adv, vn_adv ! Advection vel. at "now" barocl. step
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub2_b,  vb2_b  ! Advection vel. at "now-0.5" barocl. step
@@ -58,6 +55,4 @@
       ALLOCATE( sshn_e(jpi,jpj) , ua_e(jpi,jpj) , hu_e(jpi,jpj) , hur_e(jpi,jpj) ,      &
          &      ssha_e(jpi,jpj) , va_e(jpi,jpj) , hv_e(jpi,jpj) , hvr_e(jpi,jpj) ,      &
-         &      ub_b(jpi,jpj)   , vb_b(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj)  ,      &
-         &      ua_b(jpi,jpj)   , va_b(jpi,jpj)                                  ,      & 
          &      ub2_b(jpi,jpj)  , vb2_b(jpi,jpj)                                 ,      &
          &      un_adv(jpi,jpj) , vn_adv(jpi,jpj)                                ,      &
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 4354)
@@ -201,5 +201,5 @@
          IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~   free surface with time splitting'
-         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro
+         IF(lwp) WRITE(numout,*)
          !
          IF (neuler==0) ll_init=.TRUE.
@@ -317,11 +317,9 @@
       ! -----------------------------------------------------------------------------
       !      
-      ! Some vertical sums (at now and before time steps) below could be suppressed 
-      ! if one swap barotropic arrays somewhere
-      !
-      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated)
+      !
+      !                                   !* e3*d/dt(Ua) (Vertically integrated)
       !                                   ! --------------------------------------------------
-      zu_frc(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp  ;  un_b(:,:) = 0._wp
-      zv_frc(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp  ;  vn_b(:,:) = 0._wp
+      zu_frc(:,:) = 0._wp
+      zv_frc(:,:) = 0._wp
       !
       DO jk = 1, jpkm1
@@ -332,14 +330,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-#endif
-               !        ! now trend:                                                                   
+#endif                                                                   
                zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
                zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)         
-               !        ! now bt transp:                   
-               un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk)        
-               vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk)
-               !	! before bt transp:
-               ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk)  * umask(ji,jj,jk)
-               vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk)  * vmask(ji,jj,jk)
             END DO
          END DO
@@ -349,11 +340,4 @@
       zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)
       !
-      IF( lk_vvl ) THEN
-          ub_b(:,:) = ub_b(:,:) * zhur_b(:,:)
-          vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:)
-      ELSE
-          ub_b(:,:) = ub_b(:,:) * hur(:,:)
-          vb_b(:,:) = vb_b(:,:) * hvr(:,:)
-      ENDIF
       !
       !                                   !* baroclinic momentum trend (remove the vertical mean trend)
@@ -368,6 +352,6 @@
       !                                   !* barotropic Coriolis trends (vorticity scheme dependent)
       !                                   ! --------------------------------------------------------
-      zwx(:,:) = un_b(:,:) * e2u(:,:)           ! now transport 
-      zwy(:,:) = vn_b(:,:) * e1v(:,:)
+      zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:)        ! now fluxes 
+      zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:)
       !
       IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN      ! energy conserving or mixed scheme
@@ -412,6 +396,4 @@
       ENDIF 
       !
-      un_b (:,:) = un_b(:,:) * hur(:,:)         ! Revert now transport to barotropic velocities
-      vn_b (:,:) = vn_b(:,:) * hvr(:,:)  
       !                                   !* Right-Hand-Side of the barotropic momentum equation
       !                                   ! ----------------------------------------------------
@@ -511,9 +493,9 @@
       !                                             ! ==================== !  
       ! Initialize barotropic variables:    
-      IF (ln_bt_fw) THEN                  ! FORWARD integration:  start from NOW fields                             
+      IF (ln_bt_fw) THEN                  ! FORWARD integration: start from NOW fields                             
          sshn_e (:,:) = sshn (:,:)            
          zun_e  (:,:) = un_b (:,:)            
          zvn_e  (:,:) = vn_b (:,:)
-      ELSE                                ! CENTERED integration: start from BEFORE fields
+      ELSE                                ! CENTRED integration: start from BEFORE fields
          sshn_e (:,:) = sshb (:,:)
          zun_e  (:,:) = ub_b (:,:)         
@@ -807,14 +789,6 @@
 
 #if defined key_bdy  
-
-         pssh => ssha_e
-         phur => hur_e
-         phvr => hvr_e
-         pua2d => ua_e
-         pva2d => va_e
-         pub2d => zun_e
-         pvb2d => zvn_e
-                                      
-         IF( lk_bdy )   CALL bdy_dyn2d( kt )               ! open boundaries
+                                                           ! open boundaries
+         IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )
 #endif
 #if defined key_agrif
@@ -1118,5 +1092,5 @@
          IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true.  => Forward integration of barotropic variables '
       ELSE
-         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables '
+         IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
       ENDIF
       !
@@ -1168,2 +1142,3 @@
    !!======================================================================
 END MODULE dynspg_ts
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90	(revision 4354)
@@ -26,5 +26,5 @@
    USE timing          ! Timing
    USE dynadv          ! dynamics: vector invariant versus flux form
-   USE dynspg_oce, ONLY: lk_dynspg_ts, ua_b, va_b
+   USE dynspg_oce, ONLY: lk_dynspg_ts
    USE dynspg_ts
 
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 4354)
@@ -618,4 +618,7 @@
       USE diadct    , ONLY: diadct_alloc 
 #endif 
+#if defined key_bdy
+      USE bdy_oce   , ONLY: bdy_oce_alloc
+#endif
       !
       INTEGER :: ierr
@@ -634,4 +637,7 @@
       ierr = ierr + diadct_alloc    ()          ! 
 #endif 
+#if defined key_bdy
+      ierr = ierr + bdy_oce_alloc   ()          ! bdy masks (incl. initialization)
+#endif
       !
       IF( lk_mpp    )   CALL mpp_sum( ierr )
@@ -819,2 +825,4 @@
    !!======================================================================
 END MODULE nemogcm
+
+
Index: /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 4353)
+++ /branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 4354)
@@ -35,6 +35,8 @@
 
    !! free surface                                      !  before  ! now    ! after  !
-   !! ------------                                      !  fields  ! fields ! trends !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m]
+   !! ------------                                      !  fields  ! fields ! fields !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ub_b   ,  un_b  ,  ua_b  !: Barotropic velocities at u-point [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vb_b   ,  vn_b  ,  va_b  !: Barotropic velocities at v-point [m/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   ,  sshn  ,  ssha  !: sea surface height at t-point [m]
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient
@@ -83,5 +85,7 @@
       ALLOCATE(rhd (jpi,jpj,jpk) ,                                         &
          &     rhop(jpi,jpj,jpk) ,                                         &
-         &     sshb  (jpi,jpj)   , sshn  (jpi,jpj) , ssha  (jpi,jpj) ,     &
+         &     sshb(jpi,jpj)     , sshn(jpi,jpj)   , ssha(jpi,jpj)   ,     &
+         &     ub_b(jpi,jpj)     , un_b(jpi,jpj)   , ua_b(jpi,jpj)   ,     &
+         &     vb_b(jpi,jpj)     , vn_b(jpi,jpj)   , va_b(jpi,jpj)   ,     &
          &     spgu  (jpi,jpj)   , spgv(jpi,jpj)   ,                       &
          &     gtsu(jpi,jpj,jpts), gtsv(jpi,jpj,jpts),                     &
