Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3490)
@@ -28,4 +28,5 @@
       INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap
       REAL   , POINTER, DIMENSION(:,:)   ::  nbw
+      REAL   , POINTER, DIMENSION(:,:)   ::  nbd
       REAL   , POINTER, DIMENSION(:)     ::  flagu
       REAL   , POINTER, DIMENSION(:)     ::  flagv
@@ -73,9 +74,7 @@
    INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta      	    !: = 0 use the initial state as bdy dta ; 
                                                             !: = 1 read it in a NetCDF file
-   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp
-   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp
-   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp
-   LOGICAL                    ::   ln_tra_dmp_tot
-   LOGICAL                    ::   ln_dyn3d_dmp_tot
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_tra_dmp               !: =T Tracer damping
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_dyn3d_dmp             !: =T Baroclinic velocity damping
+   REAL,    DIMENSION(jp_bdy) ::   rn_time_dmp              !: Damping time scale in days
 
 #if defined key_lim2
@@ -107,5 +106,6 @@
    INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions
                                                                           !: =1 => some data to be read in from data files
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays
+   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)
    TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process)
    TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process)
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3490)
@@ -192,5 +192,5 @@
             IF( PRESENT(jit) ) THEN
                ! Update barotropic boundary conditions only
-               ! jit is optional argument for fld_read and tide_update
+               ! jit is optional argument for fld_read and bdytide_update
                IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
                   IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
@@ -200,47 +200,47 @@
                   ENDIF
                   IF (nn_tra(ib_bdy).ne.4) THEN
-                   IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN
-                     ! For the runoff case, no need to update the forcing (already done in the baroclinic part)
-                     jend = nb_bdy_fld(ib_bdy)
-                     IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2
-                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset )
-                     IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2
-                     ! If full velocities in boundary data then split into barotropic and baroclinic data
-                     IF( ln_full_vel_array(ib_bdy) .and.                                             &
-                    &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN
-                        igrd = 2                      ! zonal velocity
-                        dta_bdy(ib_bdy)%u2d(:) = 0.0
-                        DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                           ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                           ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                           DO ik = 1, jpkm1
-                              dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
-                    &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                     IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN
+                        ! For the runoff case, no need to update the forcing (already done in the baroclinic part)
+                        jend = nb_bdy_fld(ib_bdy)
+                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2
+                        CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset )
+                        IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2
+                        ! If full velocities in boundary data then split into barotropic and baroclinic data
+                        IF( ln_full_vel_array(ib_bdy) .and.                                             &
+                       &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN
+                           igrd = 2                      ! zonal velocity
+                           dta_bdy(ib_bdy)%u2d(:) = 0.0
+                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                              DO ik = 1, jpkm1
+                                 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
+                       &                          + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                              END DO
+                              dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
+                              DO ik = 1, jpkm1
+                                 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)
+                              END DO
                            END DO
-                           dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
-                           DO ik = 1, jpkm1
-                              dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib)
+                           igrd = 3                      ! meridional velocity
+                           dta_bdy(ib_bdy)%v2d(:) = 0.0
+                           DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                              ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                              DO ik = 1, jpkm1
+                                 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
+                       &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
+                              END DO
+                              dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
+                              DO ik = 1, jpkm1
+                                 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)
+                              END DO
                            END DO
-                        END DO
-                        igrd = 3                      ! meridional velocity
-                        dta_bdy(ib_bdy)%v2d(:) = 0.0
-                        DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
-                           ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
-                           ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
-                           DO ik = 1, jpkm1
-                              dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
-                    &                       + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
-                           END DO
-                           dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
-                           DO ik = 1, jpkm1
-                              dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib)
-                           END DO
-                        END DO
-                     ENDIF                    
-                   ENDIF
-                   IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                      CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   & 
-                      &                 jit=jit, time_offset=time_offset )
-                   ENDIF
+                        ENDIF                    
+                     ENDIF
+                     IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                        CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy),   & 
+                          &                 jit=jit, time_offset=time_offset )
+                     ENDIF
                   ENDIF
                ENDIF
@@ -306,5 +306,5 @@
                   ENDIF
                   IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                     CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )
+                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )
                   ENDIF
                ENDIF
@@ -378,18 +378,12 @@
       ! Set nn_dta
       DO ib_bdy = 1, nb_bdy
-         IF (nn_dyn2d_dta(ib_bdy).eq.1.OR.nn_dyn2d_dta(ib_bdy).eq.3) nn_dta(ib_bdy) = 1
-         IF (nn_dyn3d_dta(ib_bdy).ge.1) nn_dta(ib_bdy) = 1
-         IF (nn_tra_dta  (ib_bdy).ge.1) nn_dta(ib_bdy) = 1
-#if defined key_lim2
-         IF (nn_ice_lim2_dta(ib_bdy).eq.1) nn_dta(ib_bdy) = 1
-#endif
-         IF (lwp) THEN
-            WRITE(numout,*) 'Bdy package number ',ib_bdy
-            IF (nn_dta(ib_bdy).eq.1) THEN
-               WRITE(numout,*) 'We use external data'
-            ELSE
-               WRITE(numout,*) 'We do not use external data'
-            ENDIF
-         ENDIF
+         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       &
+                               ,nn_dyn3d_dta(ib_bdy)       &
+                               ,nn_tra_dta(ib_bdy)         &
+#if defined key_lim2
+                               ,nn_ice_lim2_dta(ib_bdy)    &
+#endif
+                              )
+         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1
       END DO
 
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 3490)
@@ -14,4 +14,5 @@
    !!----------------------------------------------------------------------
    USE timing          ! Timing
+   USE wrk_nemo        ! Memory Allocation
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
@@ -27,4 +28,6 @@
    PUBLIC   bdy_dyn3d_dmp ! routine called by step
 
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -222,4 +225,30 @@
       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp')
       !
+      !-------------------------------------------------------
+      ! Remove barotropic part from before velocity
+      !-------------------------------------------------------
+      CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
+
+      pu2d(:,:) = 0.e0
+      pv2d(:,:) = 0.e0
+
+      DO jk = 1, jpkm1
+#if defined key_vvl
+         pu2d(:,:) = pu2d(:,:) + fse3u_b(:,:,jk)* ub(:,:,jk)   *umask(:,:,jk) 
+         pv2d(:,:) = pv2d(:,:) + fse3v_b(:,:,jk)* vb(:,:,jk)   *vmask(:,:,jk)
+#else
+         pu2d(:,:) = pu2d(:,:) + fse3u_0(:,:,jk) * ub(:,:,jk)  * umask(:,:,jk)
+         pv2d(:,:) = pv2d(:,:) + fse3v_0(:,:,jk) * vb(:,:,jk)  * vmask(:,:,jk)
+#endif
+      END DO
+
+      IF( lk_vvl ) THEN
+         pu2d(:,:) = pu2d(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) )
+         pv2d(:,:) = pv2d(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) )
+      ELSE
+         pu2d(:,:) = pv2d(:,:) * hur(:,:)
+         pv2d(:,:) = pu2d(:,:) * hvr(:,:)
+      ENDIF
+
       DO ib_bdy=1, nb_bdy
          IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN
@@ -228,7 +257,7 @@
                ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
                ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
-               zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday )
+               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
                DO jk = 1, jpkm1
-                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) ) ) * umask(ii,ij,jk)
+                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) + pu2d(ii,ij)) ) * umask(ii,ij,jk)
                END DO
             END DO
@@ -238,11 +267,13 @@
                ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
                ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
-               zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday )
+               zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
                DO jk = 1, jpkm1
-                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) ) ) * vmask(ii,ij,jk)
+                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) + pv2d(ii,ij)) ) * vmask(ii,ij,jk)
                END DO
             END DO
          ENDIF
       ENDDO
+      !
+      CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
       !
       CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3490)
@@ -11,4 +11,5 @@
    !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
    !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
+   !!            3.4  !  2012     (J. Chanut) straight open boundary case update
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -26,4 +27,6 @@
    USE lib_mpp         ! for mpp_sum  
    USE iom             ! I/O
+   USE sbctide, ONLY: lk_tide ! Tidal forcing or not
+   USE phycst, ONLY: rday
 
    IMPLICIT NONE
@@ -32,4 +35,13 @@
    PUBLIC   bdy_init   ! routine called in nemo_init
 
+   INTEGER, PARAMETER          :: jp_nseg = 100
+   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured
+                                               ! open boundary data files
+   ! Straight open boundary segment parameters:
+   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
+   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge
+   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw
+   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn
+   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -53,10 +65,7 @@
       ! namelist variables
       !-------------------
-      INTEGER, PARAMETER          :: jp_nseg = 100
-      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
-      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft
-      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft
-      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft
-      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft
+      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
+      CHARACTER(LEN=1)   ::   ctypebdy
+      INTEGER :: nbdyind, nbdybeg, nbdyend
 
       ! local variables
@@ -66,14 +75,14 @@
       INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
       INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
+      INTEGER  ::   jpbdtau, jpbdtas                       !   -       -
+      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
       INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
       REAL   , POINTER  ::  flagu, flagv                   !    -   -
       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
-      INTEGER, DIMENSION (2)                ::   kdimsz
+      INTEGER, DIMENSION (2)                  ::   kdimsz
       INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays 
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
       INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
-      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask
-      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
-      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid
+      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
       !!
       NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
@@ -86,8 +95,5 @@
          &             ln_vol, nn_volctl, nn_rimwidth
       !!
-      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             &
-                             nbdysegw, jpiwob, jpjwdt, jpjwft,             &
-                             nbdysegn, jpjnob, jpindt, jpinft,             &
-                             nbdysegs, jpjsob, jpisdt, jpisft
+      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
 
       !!----------------------------------------------------------------------
@@ -106,5 +112,5 @@
 
       cgrid= (/'t','u','v'/)
-
+      
       ! -----------------------------------------
       ! Initialise and read namelist parameters
@@ -139,7 +145,4 @@
       ! Check and write out namelist parameters
       ! -----------------------------------------
- 
-      ln_tra_dmp_tot=.false.
-      ln_dyn3d_dmp_tot=.false.
       !                                   ! control prints
       IF(lwp) WRITE(numout,*) '         nambdy'
@@ -164,7 +167,7 @@
         IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
         SELECT CASE( nn_dyn2d(ib_bdy) )                  
-          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
-          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
+          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
         END SELECT
@@ -177,4 +180,7 @@
               CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
            END SELECT
+           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN
+             CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' )
+           ENDIF
         ENDIF
         IF(lwp) WRITE(numout,*)
@@ -182,6 +188,6 @@
         IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
         SELECT CASE( nn_dyn3d(ib_bdy) )                  
-          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
@@ -195,15 +201,18 @@
            END SELECT
         ENDIF
-        IF(lwp) WRITE(numout,*)
 
         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
            IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN
-              IF(lwp) WRITE(numout,*) 'No open boundary condition for the baroclinic velocitues : we force ln_dyn3d_dmp=.false.'
+              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
               ln_dyn3d_dmp(ib_bdy)=.false.
+           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN
+              CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
-              IF(lwp) WRITE(numout,*) '      Damping of the baroclinic velocities along the boundaries'
-              IF(lwp) WRITE(numout,*) '         Time damping :',rn_time_dmp(ib_bdy),' days'
-              ln_dyn3d_dmp_tot=.true.
+              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
+              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
            ENDIF
+        ELSE
+           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
         ENDIF
         IF(lwp) WRITE(numout,*)
@@ -211,6 +220,6 @@
         IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
         SELECT CASE( nn_tra(ib_bdy) )                  
-          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
-          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
@@ -225,15 +234,18 @@
            END SELECT
         ENDIF
-        IF(lwp) WRITE(numout,*)
 
         IF ( ln_tra_dmp(ib_bdy) ) THEN
            IF ( nn_tra(ib_bdy).EQ.0 ) THEN
-              IF(lwp) WRITE(numout,*) 'No open boundary condition for the tracer : we force ln_tra_dmp=.false.'
+              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
               ln_tra_dmp(ib_bdy)=.false.
+           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN
+              CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
-              IF(lwp) WRITE(numout,*) '      Damping of the scalars along the boundaries'
-              IF(lwp) WRITE(numout,*) '         Time dumping :',rn_time_dmp(ib_bdy),' days'
-              ln_tra_dmp_tot=.true.
+              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
+              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
+              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
            ENDIF
+        ELSE
+           IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
         ENDIF
         IF(lwp) WRITE(numout,*)
@@ -256,21 +268,23 @@
 #endif
 
-        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)
+        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
         IF(lwp) WRITE(numout,*)
 
       ENDDO
 
-     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
-       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
-       IF(lwp) WRITE(numout,*)
-       SELECT CASE ( nn_volctl )
-         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
-         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
-         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
-       END SELECT
-       IF(lwp) WRITE(numout,*)
-     ELSE
-       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
-       IF(lwp) WRITE(numout,*)
+     IF (nb_bdy .gt. 0) THEN
+        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
+          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
+          IF(lwp) WRITE(numout,*)
+          SELECT CASE ( nn_volctl )
+            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
+            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
+            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
+          END SELECT
+          IF(lwp) WRITE(numout,*)
+        ELSE
+          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
+          IF(lwp) WRITE(numout,*)
+        ENDIF
      ENDIF
 
@@ -282,83 +296,93 @@
       ! ---------------------------------------------
       REWIND( numnam )                    
-      jpbdta = 1 
+               
       nblendta(:,:) = 0
+      nbdysege = 0
+      nbdysegw = 0
+      nbdysegn = 0
+      nbdysegs = 0
+      icount   = 0 ! count user defined segments
+      ! Dimensions below are used to allocate arrays to read external data
+      jpbdtas = 1 ! Maximum size of boundary data (structured case)
+      jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
+
       DO ib_bdy = 1, nb_bdy
 
          IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
  
+            icount = icount + 1
             ! No REWIND here because may need to read more than one nambdy_index namelist.
             READ  ( numnam, nambdy_index )
 
-            ! Automatic boundary definition: if nbdysegX = -1
-            ! set boundary to whole side of model domain.
-            IF( nbdysege == -1 ) THEN
-               nbdysege = 1
-               jpieob(1) = jpiglo - 1
-               jpjedt(1) = 2
-               jpjeft(1) = jpjglo - 1
-            ENDIF
-            IF( nbdysegw == -1 ) THEN
-               nbdysegw = 1
-               jpiwob(1) = 2
-               jpjwdt(1) = 2
-               jpjwft(1) = jpjglo - 1
-            ENDIF
-            IF( nbdysegn == -1 ) THEN
-               nbdysegn = 1
-               jpjnob(1) = jpjglo - 1
-               jpindt(1) = 2
-               jpinft(1) = jpiglo - 1
-            ENDIF
-            IF( nbdysegs == -1 ) THEN
-               nbdysegs = 1
-               jpjsob(1) = 2
-               jpisdt(1) = 2
-               jpisft(1) = jpiglo - 1
-            ENDIF
-
-            DO iseg = 1, nbdysege
-               igrd = 1
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
-               igrd = 2
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
-               igrd = 3
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)               
-            ENDDO
-            DO iseg = 1, nbdysegw
-               igrd = 1
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
-               igrd = 2
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
-               igrd = 3
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)               
-            ENDDO
-            DO iseg = 1, nbdysegn
-               igrd = 1
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1               
-               igrd = 2
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)               
-               igrd = 3
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
-            ENDDO
-            DO iseg = 1, nbdysegs
-               igrd = 1
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
-               igrd = 2
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)
-               igrd = 3
-               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
-            ENDDO
-
-            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
+            SELECT CASE ( TRIM(ctypebdy) )
+              CASE( 'N' )
+                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
+                    nbdybeg  = 2
+                    nbdyend  = jpiglo - 1
+                 ENDIF
+                 nbdysegn = nbdysegn + 1
+                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number
+                 jpjnob(nbdysegn) = nbdyind
+                 jpindt(nbdysegn) = nbdybeg
+                 jpinft(nbdysegn) = nbdyend
+                 !
+              CASE( 'S' )
+                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+                    nbdyind  = 2           ! set boundary to whole side of model domain.
+                    nbdybeg  = 2
+                    nbdyend  = jpiglo - 1
+                 ENDIF
+                 nbdysegs = nbdysegs + 1
+                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number
+                 jpjsob(nbdysegs) = nbdyind
+                 jpisdt(nbdysegs) = nbdybeg
+                 jpisft(nbdysegs) = nbdyend
+                 !
+              CASE( 'E' )
+                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
+                    nbdybeg  = 2
+                    nbdyend  = jpjglo - 1
+                 ENDIF
+                 nbdysege = nbdysege + 1 
+                 npckge(nbdysege) = ib_bdy ! Save bdy package number
+                 jpieob(nbdysege) = nbdyind
+                 jpjedt(nbdysege) = nbdybeg
+                 jpjeft(nbdysege) = nbdyend
+                 !
+              CASE( 'W' )
+                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
+                    nbdyind  = 2           ! set boundary to whole side of model domain.
+                    nbdybeg  = 2
+                    nbdyend  = jpjglo - 1
+                 ENDIF
+                 nbdysegw = nbdysegw + 1
+                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number
+                 jpiwob(nbdysegw) = nbdyind
+                 jpjwdt(nbdysegw) = nbdybeg
+                 jpjwft(nbdysegw) = nbdyend
+                 !
+              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
+            END SELECT
+
+            ! For simplicity we assume that in case of straight bdy, arrays have the same length
+            ! (even if it is true that last tangential velocity points
+            ! are useless). This simplifies a little bit boundary data format (and agrees with format
+            ! used so far in obc package)
+
+            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
+            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
+            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
+            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
 
          ELSE            ! Read size of arrays in boundary coordinates file.
-
-
             CALL iom_open( cn_coords_file(ib_bdy), inum )
             DO igrd = 1, jpbgrd
                id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
                nblendta(igrd,ib_bdy) = kdimsz(1)
-            ENDDO
+               jpbdtau = MAX(jpbdtau, kdimsz(1))
+            ENDDO
+            CALL iom_close( inum )
 
          ENDIF 
@@ -366,178 +390,29 @@
       ENDDO ! ib_bdy
 
-      jpbdta = MAXVAL(nblendta(:,:))
-
-      ! Allocate arrays
-      !---------------
-      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
-         &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
-
-      ALLOCATE( dta_global(jpbdta, 1, jpk) )
+      IF (nb_bdy>0) THEN
+         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
+
+         ! Allocate arrays
+         !---------------
+         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
+            &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
+
+         ALLOCATE( dta_global(jpbdtau, 1, jpk) )
+         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
+         ! 
+      ENDIF
+
+      ! Now look for crossings in user (namelist) defined open boundary segments:
+      !--------------------------------------------------------------------------
+      IF ( icount>0 ) CALL bdy_ctl_seg
 
       ! Calculate global boundary index arrays or read in from file
-      !------------------------------------------------------------
-      REWIND( numnam )                    
+      !------------------------------------------------------------               
+      ! 1. Read global index arrays from boundary coordinates file.
       DO ib_bdy = 1, nb_bdy
 
-         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters
-
-            ! No REWIND here because may need to read more than one nambdy_index namelist.
-            READ  ( numnam, nambdy_index )
-
-            ! Automatic boundary definition: if nbdysegX = -1
-            ! set boundary to whole side of model domain.
-            IF( nbdysege == -1 ) THEN
-               nbdysege = 1
-               jpieob(1) = jpiglo - 1
-               jpjedt(1) = 2
-               jpjeft(1) = jpjglo - 1
-            ENDIF
-            IF( nbdysegw == -1 ) THEN
-               nbdysegw = 1
-               jpiwob(1) = 2
-               jpjwdt(1) = 2
-               jpjwft(1) = jpjglo - 1
-            ENDIF
-            IF( nbdysegn == -1 ) THEN
-               nbdysegn = 1
-               jpjnob(1) = jpjglo - 1
-               jpindt(1) = 2
-               jpinft(1) = jpiglo - 1
-            ENDIF
-            IF( nbdysegs == -1 ) THEN
-               nbdysegs = 1
-               jpjsob(1) = 2
-               jpisdt(1) = 2
-               jpisft(1) = jpiglo - 1
-            ENDIF
-
-            ! ------------ T points -------------
-            igrd = 1  
-            icount = 0
-            DO ir = 1, nn_rimwidth(ib_bdy)
-               ! east
-               DO iseg = 1, nbdysege
-                  DO ij = jpjedt(iseg), jpjeft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
-                     nbjdta(icount, igrd, ib_bdy) = ij
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! west
-               DO iseg = 1, nbdysegw
-                  DO ij = jpjwdt(iseg), jpjwft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
-                     nbjdta(icount, igrd, ib_bdy) = ij
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! north
-               DO iseg = 1, nbdysegn
-                  DO ii = jpindt(iseg), jpinft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = ii
-                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! south
-               DO iseg = 1, nbdysegs
-                  DO ii = jpisdt(iseg), jpisft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = ii
-                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-            ENDDO
-
-            ! ------------ U points -------------
-            igrd = 2  
-            icount = 0
-            DO ir = 1, nn_rimwidth(ib_bdy)
-               ! east
-               DO iseg = 1, nbdysege
-                  DO ij = jpjedt(iseg), jpjeft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
-                     nbjdta(icount, igrd, ib_bdy) = ij
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! west
-               DO iseg = 1, nbdysegw
-                  DO ij = jpjwdt(iseg), jpjwft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
-                     nbjdta(icount, igrd, ib_bdy) = ij
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! north
-               DO iseg = 1, nbdysegn
-                  DO ii = jpindt(iseg), jpinft(iseg) - 1
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = ii
-                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! south
-               DO iseg = 1, nbdysegs
-                  DO ii = jpisdt(iseg), jpisft(iseg) - 1
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = ii
-                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-            ENDDO
-
-            ! ------------ V points -------------
-            igrd = 3  
-            icount = 0
-            DO ir = 1, nn_rimwidth(ib_bdy)
-               ! east
-               DO iseg = 1, nbdysege
-                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
-                     nbjdta(icount, igrd, ib_bdy) = ij
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! west
-               DO iseg = 1, nbdysegw
-                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
-                     nbjdta(icount, igrd, ib_bdy) = ij
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! north
-               DO iseg = 1, nbdysegn
-                  DO ii = jpindt(iseg), jpinft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = ii
-                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-               ! south
-               DO iseg = 1, nbdysegs
-                  DO ii = jpisdt(iseg), jpisft(iseg)
-                     icount = icount + 1
-                     nbidta(icount, igrd, ib_bdy) = ii
-                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
-                     nbrdta(icount, igrd, ib_bdy) = ir
-                  ENDDO
-               ENDDO
-            ENDDO
-
-         ELSE            ! Read global index arrays from boundary coordinates file.
-
+         IF( ln_coords_file(ib_bdy) ) THEN
+
+            CALL iom_open( cn_coords_file(ib_bdy), inum )
             DO igrd = 1, jpbgrd
                CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
@@ -560,5 +435,4 @@
                IF (ibr_max < nn_rimwidth(ib_bdy))   &
                      CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
-
             END DO
             CALL iom_close( inum )
@@ -566,13 +440,236 @@
          ENDIF 
 
-      ENDDO 
+      ENDDO      
+    
+      ! 2. Now fill indices corresponding to straight open boundary arrays:
+      ! East
+      !-----
+      DO iseg = 1, nbdysege
+         ib_bdy = npckge(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjedt(iseg), jpjeft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjedt(iseg), jpjeft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
+            DO ij = jpjedt(iseg), jpjeft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+      ENDDO
+      !
+      ! West
+      !-----
+      DO iseg = 1, nbdysegw
+         ib_bdy = npckgw(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjwdt(iseg), jpjwft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ij = jpjwdt(iseg), jpjwft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
+            DO ij = jpjwdt(iseg), jpjwft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+               nbjdta(icount, igrd, ib_bdy) = ij
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+      ENDDO
+      !
+      ! North
+      !-----
+      DO iseg = 1, nbdysegn
+         ib_bdy = npckgn(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpindt(iseg), jpinft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+!            DO ii = jpindt(iseg), jpinft(iseg) - 1
+            DO ii = jpindt(iseg), jpinft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpindt(iseg), jpinft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+      ! South
+      !-----
+      DO iseg = 1, nbdysegs
+         ib_bdy = npckgs(iseg)
+         !
+         ! ------------ T points -------------
+         igrd=1
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpisdt(iseg), jpisft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+         !
+         ! ------------ U points -------------
+         igrd=2
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+!            DO ii = jpisdt(iseg), jpisft(iseg) - 1
+            DO ii = jpisdt(iseg), jpisft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
+         ENDDO
+         !
+         ! ------------ V points -------------
+         igrd=3
+         icount=0
+         DO ir = 1, nn_rimwidth(ib_bdy)
+            DO ii = jpisdt(iseg), jpisft(iseg)
+               icount = icount + 1
+               nbidta(icount, igrd, ib_bdy) = ii
+               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
+               nbrdta(icount, igrd, ib_bdy) = ir
+            ENDDO
+         ENDDO
+      ENDDO
+
+      !  Deal with duplicated points
+      !-----------------------------
+      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
+      ! if their distance to the bdy is greater than the other
+      ! If their distance are the same, just keep only one to avoid updating a point twice
+      DO igrd = 1, jpbgrd
+         DO ib_bdy1 = 1, nb_bdy
+            DO ib_bdy2 = 1, nb_bdy
+               IF (ib_bdy1/=ib_bdy2) THEN
+                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
+                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
+                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
+                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
+!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', & 
+!                                                       &              nbidta(ib1, igrd, ib_bdy1),      & 
+!                                                       &              nbjdta(ib2, igrd, ib_bdy2)
+                           ! keep only points with the lowest distance to boundary:
+                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
+                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
+                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
+                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
+                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                           ! Arbitrary choice if distances are the same:
+                           ELSE
+                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
+                           ENDIF
+                        END IF
+                     END DO
+                  END DO
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       ! Work out dimensions of boundary data on each processor
       ! ------------------------------------------------------
-     
-      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
-      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
-      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
-      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
+
+      ! Rather assume that boundary data indices are given on global domain
+      ! TO BE DISCUSSED ?
+!      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
+!      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
+!      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
+!      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1      
+      iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2
+      ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1
+      is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2
+      in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1
 
       DO ib_bdy = 1, nb_bdy
@@ -610,4 +707,5 @@
          ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
@@ -629,6 +727,11 @@
                      !
                      icount = icount  + 1
-                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
-                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
+
+                     ! Rather assume that boundary data indices are given on global domain
+                     ! TO BE DISCUSSED ?
+!                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
+!                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
+                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom
+                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom
                      idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
                      idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
@@ -643,7 +746,17 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
                nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
-!               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
-               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.      ! quadratic
-!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
+               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic
+!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear
+            END DO
+         END DO 
+
+         ! Compute damping coefficients
+         ! ----------------------------
+         DO igrd = 1, jpbgrd
+            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
+               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
+               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
             END DO
          END DO 
@@ -660,56 +773,47 @@
       !          = 0  elsewhere   
  
-      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution
-         zmask(         :                ,:) = 0.e0
-         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0          
-      ELSE IF( ln_mask_file ) THEN
+      IF( ln_mask_file ) THEN
          CALL iom_open( cn_mask_file, inum )
-         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
+         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
          CALL iom_close( inum )
-      ELSE
-         zmask(:,:) = 1.e0
-      ENDIF
-
-      DO ij = 1, nlcj      ! Save mask over local domain      
-         DO ii = 1, nlci
-            bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
+
+         ! Derive mask on U and V grid from mask on T grid
+         bdyumask(:,:) = 0.e0
+         bdyvmask(:,:) = 0.e0
+         DO ij=1, jpjm1
+            DO ii=1, jpim1
+               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
+               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)  
+            END DO
          END DO
-      END DO
-
-      ! Derive mask on U and V grid from mask on T grid
-      bdyumask(:,:) = 0.e0
-      bdyvmask(:,:) = 0.e0
-      DO ij=1, jpjm1
-         DO ii=1, jpim1
-            bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
-            bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1)  
+         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
+
+
+         ! Mask corrections
+         ! ----------------
+         DO ik = 1, jpkm1
+            DO ij = 1, jpj
+               DO ii = 1, jpi
+                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
+                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
+                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
+                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
+               END DO      
+            END DO
          END DO
-      END DO
-      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
-
-
-      ! Mask corrections
-      ! ----------------
-      DO ik = 1, jpkm1
-         DO ij = 1, jpj
-            DO ii = 1, jpi
-               tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
-               umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
-               vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
-               bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
-            END DO      
+
+         DO ik = 1, jpkm1
+            DO ij = 2, jpjm1
+               DO ii = 2, jpim1
+                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
+                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
+               END DO      
+            END DO
          END DO
-      END DO
-
-      DO ik = 1, jpkm1
-         DO ij = 2, jpjm1
-            DO ii = 2, jpim1
-               fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
-                  &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
-            END DO      
-         END DO
-      END DO
-
-      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)             
+
+         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
+
+      ENDIF ! ln_mask_file=.TRUE.
+      
       bdytmask(:,:) = tmask(:,:,1)
 
@@ -800,4 +904,5 @@
       ! Compute total lateral surface for volume correction:
       ! ----------------------------------------------------
+      ! JC: this must be done at each time step with key_vvl
       bdysurftot = 0.e0 
       IF( ln_vol ) THEN  
@@ -806,5 +911,5 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
-               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
                flagu => idx_bdy(ib_bdy)%flagu(ib)
                bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
@@ -819,5 +924,5 @@
             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
-               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
                flagv => idx_bdy(ib_bdy)%flagv(ib)
                bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
@@ -833,9 +938,419 @@
       ! Tidy up
       !--------
-      DEALLOCATE(nbidta, nbjdta, nbrdta)
+      IF (nb_bdy>0) THEN
+         DEALLOCATE(nbidta, nbjdta, nbrdta)
+      ENDIF
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
 
    END SUBROUTINE bdy_init
+
+   SUBROUTINE bdy_ctl_seg
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_ctl_seg  ***
+      !!
+      !! ** Purpose :   Check straight open boundary segments location
+      !!
+      !! ** Method  :   - Look for open boundary corners
+      !!                - Check that segments start or end on land 
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest  
+      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns  
+      REAL(wp), DIMENSION(2) ::   ztestmask
+      !!----------------------------------------------------------------------
+      !
+      IF (lwp) WRITE(numout,*) ' '
+      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
+      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      !
+      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
+      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
+      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
+      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
+      ! 1. Check bounds
+      !----------------
+      DO ib = 1, nbdysegn
+         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
+         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
+            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
+         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
+         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
+      END DO
+      !
+      DO ib = 1, nbdysegs
+         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
+         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
+            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
+         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
+         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
+      END DO
+      !
+      DO ib = 1, nbdysege
+         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
+         IF ((jpieob(ib).ge.jpiglo-1).or.& 
+            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
+         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
+         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
+      END DO
+      !
+      DO ib = 1, nbdysegw
+         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
+         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
+            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
+         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
+         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
+         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
+      ENDDO
+      !
+      !      
+      ! 2. Look for segment crossings
+      !------------------------------ 
+      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
+      !
+      itest = 0 ! corner number
+      !
+      ! flag to detect if start or end of open boundary belongs to a corner
+      ! if not (=0), it must be on land.
+      ! if a corner is detected, save bdy package number for further tests
+      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
+      ! South/West crossings
+      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
+         DO ib1 = 1, nbdysegw        
+            DO ib2 = 1, nbdysegs
+               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
+                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
+                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
+                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
+                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
+                     ! We have a possible South-West corner                      
+!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1) 
+!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
+                     icornw(ib1,1) = npckgs(ib2)
+                     icorns(ib2,1) = npckgw(ib1)
+                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     &                                     jpisft(ib2), jpjwft(ib1)
+                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
+                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
+                     &                                                    ' and South segment: ',npckgs(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  ELSE
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
+                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , &
+                     &                                                    ' and South segment: ',npckgs(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop+1
+                  END IF
+               END IF
+            END DO
+         END DO
+      END IF
+      !
+      ! South/East crossings
+      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
+         DO ib1 = 1, nbdysege
+            DO ib2 = 1, nbdysegs
+               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
+                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
+                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
+                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
+                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
+                     ! We have a possible South-East corner 
+!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1) 
+!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
+                     icorne(ib1,1) = npckgs(ib2)
+                     icorns(ib2,2) = npckge(ib1)
+                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     &                                     jpisdt(ib2), jpjeft(ib1)
+                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
+                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
+                     &                                                    ' and South segment: ',npckgs(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  ELSE
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
+                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
+                     &                                                    ' and South segment: ',npckgs(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  END IF
+               END IF
+            END DO
+         END DO
+      END IF
+      !
+      ! North/West crossings
+      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
+         DO ib1 = 1, nbdysegw        
+            DO ib2 = 1, nbdysegn
+               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
+                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
+                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
+                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
+                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
+                     ! We have a possible North-West corner 
+!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1) 
+!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
+                     icornw(ib1,2) = npckgn(ib2)
+                     icornn(ib2,1) = npckgw(ib1)
+                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     &                                     jpinft(ib2), jpjwdt(ib1)
+                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
+                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &
+                     &                                                    ' and North segment: ',npckgn(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  ELSE
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
+                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), &
+                     &                                                    ' and North segment: ',npckgn(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  END IF
+               END IF
+            END DO
+         END DO
+      END IF
+      !
+      ! North/East crossings
+      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
+         DO ib1 = 1, nbdysege        
+            DO ib2 = 1, nbdysegn
+               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
+                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
+                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
+                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
+                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
+                     ! We have a possible North-East corner 
+!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
+!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
+                     icorne(ib1,2) = npckgn(ib2)
+                     icornn(ib2,2) = npckge(ib1)
+                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     &                                     jpindt(ib2), jpjedt(ib1)
+                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
+                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
+                     &                                                    ' and North segment: ',npckgn(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  ELSE
+                     IF(lwp) WRITE(numout,*)
+                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
+                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
+                     &                                                    ' and North segment: ',npckgn(ib2)
+                     IF(lwp) WRITE(numout,*)
+                     nstop = nstop + 1
+                  END IF
+               END IF
+            END DO
+         END DO
+      END IF
+      !
+      ! 3. Check if segment extremities are on land
+      !-------------------------------------------- 
+      !
+      ! West segments
+      DO ib = 1, nbdysegw
+         ! get mask at boundary extremities:
+         ztestmask(1:2)=0.
+         DO ji = 1, jpi
+            DO jj = 1, jpj             
+              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
+               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
+              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
+               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
+
+         IF (ztestmask(1)==1) THEN 
+            IF (icornw(ib,1)==0) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
+               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                  
+               IF(lwp) WRITE(numout,*)
+               nstop = nstop + 1
+            ELSE
+               ! This is a corner
+               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
+               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
+               itest=itest+1
+            ENDIF
+         ENDIF
+         IF (ztestmask(2)==1) THEN
+            IF (icornw(ib,2)==0) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
+               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                  
+               IF(lwp) WRITE(numout,*)
+               nstop = nstop + 1
+            ELSE
+               ! This is a corner
+               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
+               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
+               itest=itest+1
+            ENDIF
+         ENDIF
+      END DO
+      !
+      ! East segments
+      DO ib = 1, nbdysege
+         ! get mask at boundary extremities:
+         ztestmask(1:2)=0.
+         DO ji = 1, jpi
+            DO jj = 1, jpj             
+              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
+               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
+              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
+               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
+
+         IF (ztestmask(1)==1) THEN
+            IF (icorne(ib,1)==0) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
+               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                  
+               IF(lwp) WRITE(numout,*)
+               nstop = nstop + 1 
+            ELSE
+               ! This is a corner
+               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
+               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
+               itest=itest+1
+            ENDIF
+         ENDIF
+         IF (ztestmask(2)==1) THEN
+            IF (icorne(ib,2)==0) THEN
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
+               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                  
+               IF(lwp) WRITE(numout,*)
+               nstop = nstop + 1
+            ELSE
+               ! This is a corner
+               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
+               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
+               itest=itest+1
+            ENDIF
+         ENDIF
+      END DO
+      !
+      ! South segments
+      DO ib = 1, nbdysegs
+         ! get mask at boundary extremities:
+         ztestmask(1:2)=0.
+         DO ji = 1, jpi
+            DO jj = 1, jpj             
+              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
+               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
+              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
+               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
+
+         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
+            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                  
+            IF(lwp) WRITE(numout,*)
+            nstop = nstop + 1
+         ENDIF
+         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
+            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                  
+            IF(lwp) WRITE(numout,*)
+            nstop = nstop + 1
+         ENDIF
+      END DO
+      !
+      ! North segments
+      DO ib = 1, nbdysegn
+         ! get mask at boundary extremities:
+         ztestmask(1:2)=0.
+         DO ji = 1, jpi
+            DO jj = 1, jpj             
+              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
+               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
+              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
+               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1)  
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
+
+         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
+            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                  
+            IF(lwp) WRITE(numout,*)
+            nstop = nstop + 1
+         ENDIF
+         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
+            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                  
+            IF(lwp) WRITE(numout,*)
+            nstop = nstop + 1
+         ENDIF
+      END DO
+      !
+      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
+      !
+      ! Other tests TBD: 
+      ! segments completly on land
+      ! optimized open boundary array length according to landmask
+      ! Nudging layers that overlap with interior domain
+      !
+   END SUBROUTINE bdy_ctl_seg
+
+   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_ctl_corn  ***
+      !!
+      !! ** Purpose :   Check numerical schemes consistency between
+      !!                segments having a common corner
+      !!
+      !! ** Method  :   
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  ::   ib1, ib2
+      INTEGER :: itest
+      !!----------------------------------------------------------------------
+      itest = 0
+
+      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1
+      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1
+      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1
+      !
+      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
+      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
+      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
+      !
+      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1   
+      !
+      IF ( itest>0 ) THEN
+         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
+         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                  
+         IF(lwp) WRITE(numout,*)
+         nstop = nstop + 1
+      ENDIF
+      !
+   END SUBROUTINE bdy_ctl_corn
 
 #else
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3490)
@@ -8,4 +8,5 @@
    !!            3.0  !  2008-04  (NEMO team)  add in the reference version
    !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
+   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -27,6 +28,7 @@
    USE bdy_oce         ! ocean open boundary conditions
    USE daymod          ! calendar
+   USE wrk_nemo        ! Memory allocation
    USE tideini
-   USE tide_mod
+!   USE tide_mod       ! Useless ??
    USE fldread, ONLY: fld_map
 
@@ -34,6 +36,6 @@
    PRIVATE
 
-   PUBLIC   bdytide_init  ! routine called in nemo_init
-   PUBLIC   tide_update   ! routine called in bdydyn
+   PUBLIC   bdytide_init     ! routine called in bdy_init
+   PUBLIC   bdytide_update   ! routine called in bdy_dta
 
    TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
@@ -47,5 +49,5 @@
 
    TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
-   
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -66,22 +68,31 @@
       !!-------------------
       CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
-      LOGICAL                                   ::   ln_conjug = .FALSE. !: F/T : tidal data in complex/complex conjugate
+      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
+      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
       !!
       INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
+      INTEGER                                   ::   ii, ij              !: dummy loop indices
       INTEGER                                   ::   inum, igrd
-      INTEGER, DIMENSION(3)                     ::   ilen0               !: length of boundary data (from OBC arrays)
+      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
+      INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts
       CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file 
       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
+      REAL(wp), POINTER, DIMENSION(:,:)         ::   ztr, zti            !:  "     "    "   "   "   "        "      " 
       !!
       TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
       !!
-      NAMELIST/nambdy_tide/filtide, ln_conjug
+      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
       !!----------------------------------------------------------------------
 
       IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
 
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+      IF (nb_bdy>0) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+
+      ln_bdytide_2ddta = .FALSE.
+      ln_bdytide_conj  = .FALSE.
 
       REWIND(numnam)
@@ -90,4 +101,6 @@
 
             td => tides(ib_bdy)
+            nblen => idx_bdy(ib_bdy)%nblen
+            nblenrim => idx_bdy(ib_bdy)%nblenrim
 
             ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
@@ -97,65 +110,145 @@
             READ  ( numnam, nambdy_tide )
             !                                               ! Parameter control and print
+            IF(lwp) WRITE(numout,*) '  '
             IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
-            IF(lwp) WRITE(numout,*) '             tidal components specified ', nb_harmo
-            IF(lwp) WRITE(numout,*) '                ', Wave(ntide(1:nb_harmo))%cname_tide
-            IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
-            IF(lwp) WRITE(numout,*) '                ', omega_tide(1:nb_harmo)
+            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
+            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
+            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
+            IF(lwp) THEN 
+                    WRITE(numout,*) '             Tidal cpt name    -     Phase speed (deg/hr)'            
+               DO itide = 1, nb_harmo
+                  WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)   
+               END DO
+            ENDIF 
+            IF(lwp) WRITE(numout,*) ' '
 
             ! Allocate space for tidal harmonics data - get size from OBC data arrays
             ! -----------------------------------------------------------------------
 
-            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 
+            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
+            ! relaxation area      
+            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+               ilen0(:)=nblen(:)
+            ELSE
+               ilen0(:)=nblenrim(:)
+            ENDIF
+
             ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
             ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
 
-            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 
             ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
             ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
 
-            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 
             ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
             ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
 
-            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) )
-
-            ! Open files and read in tidal forcing data
-            ! -----------------------------------------
-
-            DO itide = 1, nb_harmo
-               !                                                              ! SSH fields
-               clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
-               CALL iom_open( clfile, inum )
-               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
-               td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
-               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
-               td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
+            td%ssh0(:,:,:) = 0.e0
+            td%ssh(:,:,:) = 0.e0
+            td%u0(:,:,:) = 0.e0
+            td%u(:,:,:) = 0.e0
+            td%v0(:,:,:) = 0.e0
+            td%v(:,:,:) = 0.e0
+
+            IF (ln_bdytide_2ddta) THEN
+               ! It is assumed that each data file contains all complex harmonic amplitudes
+               ! given on the data domain (ie global, jpidta x jpjdta)
+               !
+               CALL wrk_alloc( jpi, jpj, zti, ztr )
+               !
+               ! SSH fields
+               clfile = TRIM(filtide)//'_grid_T.nc'
+               CALL iom_open (clfile , inum ) 
+               igrd = 1                       ! Everything is at T-points here
+               DO itide = 1, nb_harmo
+                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
+                  DO ib = 1, ilen0(igrd)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     td%ssh0(ib,itide,1) = ztr(ii,ij)
+                     td%ssh0(ib,itide,2) = zti(ii,ij)
+                  END DO
+               END DO 
                CALL iom_close( inum )
-               !                                                              ! U fields
-               clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
-               CALL iom_open( clfile, inum )
-               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
-               td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
-               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
-               td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
+               !
+               ! U fields
+               clfile = TRIM(filtide)//'_grid_U.nc'
+               CALL iom_open (clfile , inum ) 
+               igrd = 2                       ! Everything is at U-points here
+               DO itide = 1, nb_harmo
+                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
+                  DO ib = 1, ilen0(igrd)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     td%u0(ib,itide,1) = ztr(ii,ij)
+                     td%u0(ib,itide,2) = zti(ii,ij)
+                  END DO
+               END DO
                CALL iom_close( inum )
-               !                                                              ! V fields
-               clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
-               CALL iom_open( clfile, inum )
-               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
-               td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
-               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
-               td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
+               !
+               ! V fields
+               clfile = TRIM(filtide)//'_grid_V.nc'
+               CALL iom_open (clfile , inum ) 
+               igrd = 3                       ! Everything is at V-points here
+               DO itide = 1, nb_harmo
+                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
+                  DO ib = 1, ilen0(igrd)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     td%v0(ib,itide,1) = ztr(ii,ij)
+                     td%v0(ib,itide,2) = zti(ii,ij)
+                  END DO
+               END DO  
                CALL iom_close( inum )
-               IF ( ln_conjug ) THEN
-                  IF(lwp) WRITE(numout,*) '       The tidal input data are written in complex conjugate'
-                  td%ssh0(:,itide,2) = - td%ssh0(:,itide,2)
-                  td%u0  (:,itide,2) = - td%u0  (:,itide,2)
-                  td%v0  (:,itide,2) = - td%v0  (:,itide,2)
-               ENDIF
-               !
-            END DO ! end loop on tidal components
+               !
+               CALL wrk_dealloc( jpi, jpj, ztr, zti ) 
+               !
+            ELSE            
+               !
+               ! Read tidal data only on bdy segments
+               ! 
+               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
+
+               ! Open files and read in tidal forcing data
+               ! -----------------------------------------
+
+               DO itide = 1, nb_harmo
+                  !                                                              ! SSH fields
+                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
+                  CALL iom_open( clfile, inum )
+                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
+                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
+                  CALL iom_close( inum )
+                  !                                                              ! U fields
+                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
+                  CALL iom_open( clfile, inum )
+                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
+                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
+                  CALL iom_close( inum )
+                  !                                                              ! V fields
+                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
+                  CALL iom_open( clfile, inum )
+                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
+                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
+                  CALL iom_close( inum )
+                  !
+               END DO ! end loop on tidal components
+               !
+               DEALLOCATE( dta_read )
+            ENDIF ! ln_bdytide_2ddta=.true.
             !
-            DEALLOCATE( dta_read )
+            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
+               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
+               td%u0  (:,:,2) = - td%u0  (:,:,2)
+               td%v0  (:,:,2) = - td%v0  (:,:,2)
+            ENDIF
             !
          ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
@@ -167,7 +260,7 @@
    END SUBROUTINE bdytide_init
 
-   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE tide_update  ***
+   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdytide_update  ***
       !!                
       !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
@@ -186,11 +279,16 @@
                                                         ! etc.
       !!
+      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
       INTEGER                          :: itide, igrd, ib   ! dummy loop indices
       INTEGER                          :: time_add          ! time offset in units of timesteps
-      REAL(wp)                         :: z_arg, z_sarg, zflag      
+      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp      
       REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
       !!----------------------------------------------------------------------
 
-      IF( nn_timing == 1 ) CALL timing_start('tide_update')
+      IF( nn_timing == 1 ) CALL timing_start('bdytide_update')
+
+      ilen0(1) =  SIZE(td%ssh(:,1,1))
+      ilen0(2) =  SIZE(td%u(:,1,1))
+      ilen0(3) =  SIZE(td%v(:,1,1))
 
       zflag=1
@@ -205,6 +303,6 @@
         IF(lwp) THEN
            WRITE(numout,*)
-           WRITE(numout,*) 'bdy_tide : (re)Initialization of the tidal bdy forcing at kt=',kt
-           WRITE(numout,*) '~~~~~~~ '
+           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
+           WRITE(numout,*) '~~~~~~~~~~~~~~ '
         ENDIF
         !
@@ -225,4 +323,8 @@
       ENDIF
 
+      ! Linear ramp on tidal component at open boundaries 
+      zramp = 1.
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.)
+
       DO itide = 1, nb_harmo
          z_sarg = z_arg * omega_tide(itide)
@@ -233,20 +335,20 @@
       DO itide = 1, nb_harmo
          igrd=1                              ! SSH on tracer grid
-         DO ib = 1, idx%nblenrim(igrd)
-            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)
+         DO ib = 1, ilen0(igrd)
+            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
          END DO
          igrd=2                              ! U grid
-         DO ib=1, idx%nblenrim(igrd)
-            dta%u2d(ib) = dta%u2d(ib) + td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)
+         DO ib = 1, ilen0(igrd)
+            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
          END DO
          igrd=3                              ! V grid
-         DO ib=1, idx%nblenrim(igrd)
-            dta%v2d(ib) = dta%v2d(ib) + td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)
+         DO ib = 1, ilen0(igrd) 
+            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
          END DO
       END DO
       !
-      IF( nn_timing == 1 ) CALL timing_stop('tide_update')
+      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update')
       !
-   END SUBROUTINE tide_update
+   END SUBROUTINE bdytide_update
 
    SUBROUTINE tide_init_elevation( idx, td )
@@ -257,21 +359,25 @@
       TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data
       !! * Local declarations
+      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays)
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
 
-      igrd=1                                 ! SSH on tracer grid.
-
-      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
+      igrd=1   
+                              ! SSH on tracer grid.
+   
+      ilen0(1) =  SIZE(td%ssh0(:,1,1))
+
+      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
 
       DO itide = 1, nb_harmo
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
             phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
          END DO
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1 , ilen0(igrd)
             mod_tide(ib)=mod_tide(ib)*ftide(itide)
             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
          ENDDO
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1 , ilen0(igrd)
             td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
             td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
@@ -290,21 +396,25 @@
       TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
       !! * Local declarations
+      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
       REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
       INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
 
+      ilen0(2) =  SIZE(td%u0(:,1,1))
+      ilen0(3) =  SIZE(td%v0(:,1,1))
+
       igrd=2                                 ! U grid.
 
-      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
+      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
 
       DO itide = 1, nb_harmo
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
             phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
          END DO
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             mod_tide(ib)=mod_tide(ib)*ftide(itide)
             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
          ENDDO
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
             td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
@@ -314,18 +424,18 @@
       DEALLOCATE(mod_tide,phi_tide)
 
-      igrd=3                                 ! U grid.
-
-      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
+      igrd=3                                 ! V grid.
+
+      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
 
       DO itide = 1, nb_harmo
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
             phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
          END DO
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             mod_tide(ib)=mod_tide(ib)*ftide(itide)
             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
          ENDDO
-         DO ib = 1, idx%nblenrim(igrd)
+         DO ib = 1, ilen0(igrd)
             td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
             td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
@@ -341,11 +451,10 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE bdytide_init                ! Empty routine
+   SUBROUTINE bdytide_init             ! Empty routine
+      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?'
    END SUBROUTINE bdytide_init
-   SUBROUTINE tide_data                ! Empty routine
-   END SUBROUTINE tide_data
-   SUBROUTINE tide_update( kt, jit )   ! Empty routine
-      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit
-   END SUBROUTINE tide_update
+   SUBROUTINE bdytide_update( kt, jit )   ! Empty routine
+      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit
+   END SUBROUTINE bdytide_update
 #endif
 
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3490)
@@ -22,5 +22,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
-   USE phycst
+
 
    IMPLICIT NONE
@@ -28,5 +28,5 @@
 
    PUBLIC bdy_tra      ! routine called in tranxt.F90 
-   PUBLIC bdy_tra_dmp  ! routine called in tranxt.F90 
+   PUBLIC bdy_tra_dmp  ! routine called in step.F90 
 
    !!----------------------------------------------------------------------
@@ -65,4 +65,8 @@
          END SELECT
       ENDDO
+      !
+      ! Boundary points should be updated
+      IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )
+      IF (nb_bdy>0) CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )
 
    END SUBROUTINE bdy_tra
@@ -98,6 +102,4 @@
       END DO 
       !
-      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
-      !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
       !
@@ -133,6 +135,4 @@
          END DO
       END DO
-      !
-      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
       !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
@@ -190,6 +190,4 @@
       END DO
       !
-      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
-      !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
       !
@@ -229,6 +227,4 @@
          END DO
       END DO
-      !
-      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
       !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )
@@ -262,14 +258,14 @@
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-               zwgt = idx_bdy(ib_bdy)%nbw(ib,igrd) / ( rn_time_dmp(ib_bdy) * rday )
+               zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
                DO ik = 1, jpkm1
-                  tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik)
-                  tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik)
+                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik)
+                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik)
+                  tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta
+                  tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa
                END DO
             END DO
          ENDIF
       ENDDO
-      !
-      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )  ;  CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
       !
       IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp')
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 3490)
@@ -26,5 +26,6 @@
        ftide
 
-  LOGICAL, PUBLIC :: ln_tide_pot = .false.
+  LOGICAL, PUBLIC :: ln_tide_pot = .false., ln_tide_ramp = .false.
+  REAL(wp), PUBLIC :: rdttideramp 
   INTEGER, PUBLIC :: nb_harmo
   INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide
@@ -46,5 +47,5 @@
     CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
     !
-    NAMELIST/nam_tide/ln_tide_pot,clname
+    NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname
     !!----------------------------------------------------------------------
 
@@ -54,5 +55,5 @@
           WRITE(numout,*)
           WRITE(numout,*) 'tide_init : Initialization of the tidal components'
-          WRITE(numout,*) '~~~~~~~ '
+          WRITE(numout,*) '~~~~~~~~~ '
        ENDIF
        !
@@ -77,4 +78,10 @@
           WRITE(numout,*) '        Namelist nam_tide'
           WRITE(numout,*) '        nb_harmo    = ', nb_harmo
+          WRITE(numout,*) '        ln_tide_ramp = ', ln_tide_ramp 
+          WRITE(numout,*) '        rdttideramp = ', rdttideramp
+          IF (ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp)) &
+          & CALL ctl_stop('rdttideramp must be lower than run duration')
+          IF (ln_tide_ramp.AND.(rdttideramp<0.)) &
+          & CALL ctl_stop('rdttideramp must be positive')
           CALL flush(numout)
        ENDIF
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90	(revision 3490)
@@ -13,4 +13,5 @@
   USE sbctide
   USE dynspg_oce
+  USE tideini, ONLY: ln_tide_ramp, rdttideramp
 
   IMPLICIT NONE
@@ -33,15 +34,21 @@
     INTEGER, INTENT( in ) ::   kt,kit      ! ocean time-step index
     INTEGER  :: ji,jj,jk
+    REAL (wp) :: zramp
     REAL (wp), DIMENSION(nb_harmo) :: zwt 
-    !...............................................................................
-    ! Potentiel astronomique
     !...............................................................................
 
     pot_astro(:,:)=0.e0
+    zramp = 1.e0
 
     IF (lk_dynspg_ts) THEN
        zwt(:) = omega_tide(:)* ((kt-kt_tide)*rdt + kit*(rdt/REAL(nn_baro,wp)))
+       IF (ln_tide_ramp) THEN
+          zramp = MIN(MAX( ((kt-nit000)*rdt + kit*(rdt/REAL(nn_baro,wp)))/(rdttideramp*rday),0.),1.)
+       ENDIF
     ELSE
        zwt(:) = omega_tide(:)*(kt-kt_tide)*rdt
+       IF (ln_tide_ramp) THEN
+          zramp = MIN(MAX( ((kt-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
+       ENDIF  
     ENDIF
 
@@ -49,5 +56,5 @@
        do ji=1,jpi
           do jj=1,jpj
-             pot_astro(ji,jj)=pot_astro(ji,jj) + (amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))      
+             pot_astro(ji,jj)=pot_astro(ji,jj) + zramp*(amp_pot(ji,jj,jk)*COS(zwt(jk)+phi_pot(ji,jj,jk)))      
           enddo
        enddo
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3490)
@@ -76,4 +76,6 @@
    USE mod_ioclient
 #endif
+   USE sbctide, ONLY: lk_tide
+
 
    IMPLICIT NONE
@@ -318,5 +320,5 @@
                             CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
-                            CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics
+      IF( lk_tide       )   CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics
 
       IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation
Index: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3489)
+++ branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3490)
@@ -95,5 +95,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
-      IF( kstp /= nit000 )   CALL tide_init ( kstp )
+      IF( lk_tide.AND.(kstp /= nit000 ))   CALL tide_init ( kstp )
       IF( lk_tide    )   CALL sbc_tide( kstp )
       IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
@@ -188,6 +188,5 @@
       IF( lk_trabbl      )   CALL tra_bbl    ( kstp )       ! advective (and/or diffusive) bottom boundary layer scheme
       IF( ln_tradmp      )   CALL tra_dmp    ( kstp )       ! internal damping trends
-      IF( lk_bdy .AND.     &
-         & ln_tra_dmp_tot)   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
+      IF( lk_bdy         )   CALL bdy_tra_dmp( kstp )       ! bdy damping trends
                              CALL tra_adv    ( kstp )       ! horizontal & vertical advection
       IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
@@ -222,6 +221,5 @@
          & ln_dyninc       )   CALL dyn_asm_inc( kstp )     ! apply dynamics assimilation increment
       IF( ln_neptsimp )        CALL dyn_nept_cor( kstp )    ! subtract Neptune velocities (simplified)
-      IF( lk_bdy .AND.     &
-         & ln_dyn3d_dmp_tot)   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
+      IF( lk_bdy           )   CALL bdy_dyn3d_dmp(kstp )    ! bdy damping trends
                                CALL dyn_adv( kstp )         ! advection (vector or flux form)
                                CALL dyn_vor( kstp )         ! vorticity term including Coriolis
