Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/CONFIG/AMM12/EXP00/namelist
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/CONFIG/AMM12/EXP00/namelist	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/CONFIG/AMM12/EXP00/namelist	(revision 3583)
@@ -371,17 +371,20 @@
 &nam_tide      !   tide parameters (#ifdef key_tide)
 !-----------------------------------------------------------------------
-   ln_tide_pot   = .true.   !  use tidal potential forcing
-   nb_harmo      =    11    !  number of constituents used
-   clname(1)     =   'M2'   !  name of constituent
-   clname(2)     =   'S2'
-   clname(3)     =   'N2'
-   clname(4)     =   'K1'
-   clname(5)     =   'O1'
-   clname(6)     =   'Q1'
-   clname(7)     =   'M4'
-   clname(8)     =   'K2'
-   clname(9)     =   'P1'
-   clname(10)    =   'Mf'
-   clname(11)    =   'Mm'
+   ln_tide_pot   = .false.   !  use tidal potential forcing
+   clname(1)     =   'Q1'   !  name of constituent
+   clname(2)     =   'O1'
+   clname(3)     =   'P1'
+   clname(4)     =   'S1'
+   clname(5)     =   'K1'
+   clname(6)     =   '2N2'
+   clname(7)     =   'MU2'
+   clname(8)     =   'N2'
+   clname(9)     =   'NU2'
+   clname(10)    =   'M2'
+   clname(11)    =   'L2'
+   clname(12)    =   'T2'
+   clname(13)    =   'S2'
+   clname(14)    =   'K2'
+   clname(15)    =   'M4'
 /
 !-----------------------------------------------------------------------
@@ -427,35 +430,4 @@
 !-----------------------------------------------------------------------
    filtide      = 'bdydta/amm12_bdytide_'         !  file name root of tidal forcing files
-    tide_cpt(1)   ='Q1'  !  names of tidal components used
-    tide_cpt(2)   ='O1'  !  names of tidal components used
-    tide_cpt(3)   ='P1'  !  names of tidal components used
-    tide_cpt(4)   ='S1'  !  names of tidal components used
-    tide_cpt(5)   ='K1'  !  names of tidal components used
-    tide_cpt(6)   ='2N2' !  names of tidal components used
-    tide_cpt(7)   ='MU2' !  names of tidal components used
-    tide_cpt(8)   ='N2'  !  names of tidal components used
-    tide_cpt(9)   ='NU2' !  names of tidal components used
-    tide_cpt(10)   ='M2'  !  names of tidal components used
-    tide_cpt(11)   ='L2'  !  names of tidal components used
-    tide_cpt(12)   ='T2'  !  names of tidal components used
-    tide_cpt(13)   ='S2'  !  names of tidal components used
-    tide_cpt(14)   ='K2'  !  names of tidal components used
-    tide_cpt(15)   ='M4'  !  names of tidal components used
-    tide_speed(1)   = 13.398661 !  phase speeds of tidal components (deg/hour)
-    tide_speed(2)   = 13.943036 !  phase speeds of tidal components (deg/hour)
-    tide_speed(3)   = 14.958932 !  phase speeds of tidal components (deg/hour)
-    tide_speed(4)   = 15.000001 !  phase speeds of tidal components (deg/hour)
-    tide_speed(5)   = 15.041069 !  phase speeds of tidal components (deg/hour)
-    tide_speed(6)   = 27.895355 !  phase speeds of tidal components (deg/hour)
-    tide_speed(7)   = 27.968210 !  phase speeds of tidal components (deg/hour)
-    tide_speed(8)   = 28.439730 !  phase speeds of tidal components (deg/hour)
-    tide_speed(9)   = 28.512585 !  phase speeds of tidal components (deg/hour)
-    tide_speed(10)   = 28.984106 !  phase speeds of tidal components (deg/hour)
-    tide_speed(11)   = 29.528479 !  phase speeds of tidal components (deg/hour)
-    tide_speed(12)   = 29.958935 !  phase speeds of tidal components (deg/hour)
-    tide_speed(13)   = 30.000002 !  phase speeds of tidal components (deg/hour)
-    tide_speed(14)   = 30.082138 !  phase speeds of tidal components (deg/hour)
-    tide_speed(15)   = 57.968212 !  phase speeds of tidal components (deg/hour)
-    ln_tide_date = .true.               !  adjust tidal harmonics for start date of run
 /
 !!======================================================================
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3583)
@@ -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,4 +74,8 @@
    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               !: =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
    INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables 
@@ -101,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_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3583)
@@ -32,4 +32,5 @@
    USE ice_2
 #endif
+   USE sbcapr
 
    IMPLICIT NONE
@@ -108,9 +109,5 @@
 
             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 1
                DO ib = 1, ilen1(igrd)
@@ -134,9 +131,5 @@
 
             IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 2 
                DO ib = 1, ilen1(igrd)
@@ -158,9 +151,5 @@
 
             IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 1                       ! Everything is at T-points here
                DO ib = 1, ilen1(igrd)
@@ -176,9 +165,5 @@
 #if defined key_lim2
             IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(:) = nblen(:)
-               ELSE
-                  ilen1(:) = nblenrim(:)
-               ENDIF
+               ilen1(:) = nblen(:)
                igrd = 1                       ! Everything is at T-points here
                DO ib = 1, ilen1(igrd)
@@ -207,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
@@ -214,71 +199,131 @@
                      dta_bdy(ib_bdy)%v2d(:) = 0.0
                   ENDIF
-                  IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data
-                     jend = jstart + 2
-                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),   &
-                     &              jit=jit, time_offset=time_offset )
-                  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 )
+                  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)
+                              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
+                           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 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
             ELSE
-               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                  dta_bdy(ib_bdy)%ssh(:) = 0.0
-                  dta_bdy(ib_bdy)%u2d(:) = 0.0
-                  dta_bdy(ib_bdy)%v2d(:) = 0.0
+               IF (nn_tra(ib_bdy).eq.4) then      ! runoff condition
+                  jend = nb_bdy_fld(ib_bdy)
+                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
+                  !
+                  igrd = 2                      ! zonal velocity
+                  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)
+                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+                  END DO
+                  !
+                  igrd = 3                      ! meridional velocity
+                  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)
+                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+                  END DO
+               ELSE
+                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     dta_bdy(ib_bdy)%ssh(:) = 0.0
+                     dta_bdy(ib_bdy)%u2d(:) = 0.0
+                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                  ENDIF
+                  IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data
+                     jend = nb_bdy_fld(ib_bdy)
+                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
+                  ENDIF
+                  ! 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
+                     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
+                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. 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), time_offset=time_offset )
+                  ENDIF
                ENDIF
-               IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data
-                  jend = jstart + nb_bdy_fld(ib_bdy) - 1
-                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
-               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 )
-               ENDIF
             ENDIF
             jstart = jend+1
-
-            ! If full velocities in boundary data then split into barotropic and baroclinic data
-            ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same
-            ! time as the dynspg_ts option). 
-
-            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)
+         END IF ! nn_dta(ib_bdy) = 1
+      END DO  ! ib_bdy
+
+      IF ( ln_apr_obc ) THEN
+         DO ib_bdy = 1, nb_bdy
+            IF (nn_tra(ib_bdy).NE.4)THEN
+               igrd = 1                      ! meridional velocity
+               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(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
-
-               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
-
-         END IF ! nn_dta(ib_bdy) = 1
-      END DO  ! ib_bdy
+                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij)
+               ENDDO
+            ENDIF
+         ENDDO
+      ENDIF
 
       IF( nn_timing == 1 ) CALL timing_stop('bdy_dta')
@@ -326,4 +371,9 @@
       IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init')
 
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+      IF(lwp) WRITE(numout,*) ''
+
       ! Set nn_dta
       DO ib_bdy = 1, nb_bdy
@@ -357,4 +407,5 @@
          ENDIF
 #endif               
+         IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy)
       ENDDO            
 
@@ -408,9 +459,4 @@
             ln_full_vel_array(ib_bdy) = ln_full_vel
 
-            IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN
-               CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',&
-            &                  'with dynspg_ts option' )   ;   RETURN  
-            ENDIF             
-
             nblen => idx_bdy(ib_bdy)%nblen
             nblenrim => idx_bdy(ib_bdy)%nblenrim
@@ -420,24 +466,19 @@
             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
 
-               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
+               IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
                   ibdy(jfld) = ib_bdy
                   igrid(jfld) = 1
-                  ilen1(jfld) = nblenrim(igrid(jfld))
+                  ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
                ENDIF
 
                IF( .not. ln_full_vel_array(ib_bdy) ) THEN
-
                   jfld = jfld + 1
                   blf_i(jfld) = bn_u2d
                   ibdy(jfld) = ib_bdy
                   igrid(jfld) = 2
-                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                     ilen1(jfld) = nblen(igrid(jfld))
-                  ELSE
-                     ilen1(jfld) = nblenrim(igrid(jfld))
-                  ENDIF
+                  ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
 
@@ -446,11 +487,6 @@
                   ibdy(jfld) = ib_bdy
                   igrid(jfld) = 3
-                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                     ilen1(jfld) = nblen(igrid(jfld))
-                  ELSE
-                     ilen1(jfld) = nblenrim(igrid(jfld))
-                  ENDIF
+                  ilen1(jfld) = nblen(igrid(jfld))
                   ilen3(jfld) = 1
-
                ENDIF
 
@@ -466,9 +502,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 2
-               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -477,9 +509,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 3
-               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -493,9 +521,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -504,9 +528,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = jpk
 
@@ -521,9 +541,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = 1
 
@@ -532,9 +548,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = 1
 
@@ -543,9 +555,5 @@
                ibdy(jfld) = ib_bdy
                igrid(jfld) = 1
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
+               ilen1(jfld) = nblen(igrid(jfld))
                ilen3(jfld) = 1
 
@@ -566,5 +574,4 @@
       ENDDO ! ib_bdy
 
-
       DO jfld = 1, nb_bdy_fld_sum
          ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) )
@@ -577,5 +584,5 @@
       jstart = 1
       DO ib_bdy = 1, nb_bdy
-         jend = jstart + nb_bdy_fld(ib_bdy) - 1
+         jend = nb_bdy_fld(ib_bdy) 
          CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta',   &
          &              'open boundary conditions', 'nambdy_dta' )
@@ -596,12 +603,13 @@
          IF (nn_dyn2d(ib_bdy) .gt. 0) THEN
             IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN
-               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
-                  ilen0(1:3) = nblen(1:3)
-               ELSE
-                  ilen0(1:3) = nblenrim(1:3)
-               ENDIF
-               ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) )
+               ilen0(1:3) = nblen(1:3)
                ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) )
                ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) )
+               IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN
+                  jfld = jfld + 1
+                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)
+               ELSE
+                  ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) )
+               ENDIF
             ELSE
                IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
@@ -617,9 +625,5 @@
 
          IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN
-            IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
-               ilen0(1:3) = nblen(1:3)
-            ELSE
-               ilen0(1:3) = nblenrim(1:3)
-            ENDIF
+            ilen0(1:3) = nblen(1:3)
             ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) )
             ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) )
@@ -636,9 +640,5 @@
          IF (nn_tra(ib_bdy) .gt. 0) THEN
             IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN
-               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
-                  ilen0(1:3) = nblen(1:3)
-               ELSE
-                  ilen0(1:3) = nblenrim(1:3)
-               ENDIF
+               ilen0(1:3) = nblen(1:3)
                ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) )
                ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) )
@@ -654,9 +654,5 @@
          IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN
             IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN
-               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
-                  ilen0(1:3) = nblen(1:3)
-               ELSE
-                  ilen0(1:3) = nblenrim(1:3)
-               ENDIF
+               ilen0(1:3) = nblen(1:3)
                ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) )
                ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) )
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 3583)
@@ -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
@@ -19,4 +20,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
+   Use phycst
 
    IMPLICIT NONE
@@ -24,5 +26,8 @@
 
    PUBLIC   bdy_dyn3d     ! routine called by bdy_dyn
-
+   PUBLIC   bdy_dyn3d_dmp ! routine called by step
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -55,4 +60,8 @@
          CASE(jp_frs)
             CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(2)
+            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(3)
+            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
          CASE DEFAULT
             CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
@@ -62,4 +71,92 @@
    END SUBROUTINE bdy_dyn3d
 
+   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
+      !!
+      !! ** Purpose : - Apply a specified value for baroclinic velocities
+      !!                at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !!
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe')
+      !
+      igrd = 2                      ! Relaxation of zonal velocity
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            ij   = idx%nbj(jb,igrd)
+            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
+         END DO
+      END DO
+      !
+      igrd = 3                      ! Relaxation of meridional velocity
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            ij   = idx%nbj(jb,igrd)
+            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
+         END DO
+      END DO
+      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe')
+
+   END SUBROUTINE bdy_dyn3d_spe
+
+   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
+      !!
+      !! ** Purpose : - baroclinic velocities = 0. at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !!
+      INTEGER  ::   ib, ik         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd, zcoef   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro')
+      !
+      igrd = 2                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ua(ii,ij,ik) = 0._wp
+         END DO
+      END DO
+
+      igrd = 3                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            va(ii,ij,ik) = 0._wp
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro')
+
+   END SUBROUTINE bdy_dyn3d_zro
+
    SUBROUTINE bdy_dyn3d_frs( idx, dta, kt )
       !!----------------------------------------------------------------------
@@ -111,4 +208,78 @@
    END SUBROUTINE bdy_dyn3d_frs
 
+   SUBROUTINE bdy_dyn3d_dmp( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
+      !!
+      !! ** Purpose : Apply damping for baroclinic velocities at open boundaries.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER                     ::   kt
+      !!
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::  ib_bdy          ! loop index
+      !!----------------------------------------------------------------------
+      !
+      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
+            igrd = 2                      ! Relaxation of zonal velocity
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+               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) + pu2d(ii,ij)) ) * umask(ii,ij,jk)
+               END DO
+            END DO
+            !
+            igrd = 3                      ! Relaxation of meridional velocity
+            DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+               ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+               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) + 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
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp')
+
+   END SUBROUTINE bdy_dyn3d_dmp
 
 #else
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3583)
@@ -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,18 +75,19 @@
       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,             &
          &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
          &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &  
+         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              &
 #if defined key_lim2
          &             nn_ice_lim2, nn_ice_lim2_dta,                       &
@@ -85,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
 
       !!----------------------------------------------------------------------
@@ -105,5 +112,5 @@
 
       cgrid= (/'t','u','v'/)
-
+      
       ! -----------------------------------------
       ! Initialise and read namelist parameters
@@ -121,4 +128,7 @@
       nn_tra(:)         = 0
       nn_tra_dta(:)     = -1  ! uninitialised flag
+      ln_tra_dmp(:)     = .false.
+      ln_dyn3d_dmp(:)   = .false.
+      rn_time_dmp(:)    = 1.
 #if defined key_lim2
       nn_ice_lim2(:)    = 0
@@ -135,5 +145,4 @@
       ! Check and write out namelist parameters
       ! -----------------------------------------
-
       !                                   ! control prints
       IF(lwp) WRITE(numout,*) '         nambdy'
@@ -158,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
@@ -171,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,*)
@@ -176,6 +188,8 @@
         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)'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
         END SELECT
@@ -187,10 +201,28 @@
            END SELECT
         ENDIF
+
+        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
+           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN
+              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,*) '      + 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,*)
 
         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'
+          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
         END SELECT
@@ -201,4 +233,19 @@
               CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
            END SELECT
+        ENDIF
+
+        IF ( ln_tra_dmp(ib_bdy) ) THEN
+           IF ( nn_tra(ib_bdy).EQ.0 ) THEN
+              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,*) '      + 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,*)
@@ -221,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
 
@@ -247,87 +296,93 @@
       ! ---------------------------------------------
       REWIND( numnam )                    
+               
+      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
 
-         jpbdta = 1
          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
-
-            nblendta(:,ib_bdy) = 0
-            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)
-            jpbdta = MAXVAL(nblendta(:,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 )
-            jpbdta = 1
             DO igrd = 1, jpbgrd
                id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
                nblendta(igrd,ib_bdy) = kdimsz(1)
-               jpbdta = MAX(jpbdta, kdimsz(1))
-            ENDDO
+               jpbdtau = MAX(jpbdtau, kdimsz(1))
+            ENDDO
+            CALL iom_close( inum )
 
          ENDIF 
@@ -335,176 +390,29 @@
       ENDDO ! ib_bdy
 
-      ! 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) )
@@ -527,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 )
@@ -533,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
@@ -577,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) )
@@ -596,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
@@ -611,6 +747,16 @@
                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+1-nbr)/FLOAT(nn_rimwidth))**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) = (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 
@@ -627,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)
 
@@ -753,5 +890,5 @@
             END IF
          END DO
- 
+
          IF( icount /= 0 ) THEN
             IF(lwp) WRITE(numout,*)
@@ -767,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  
@@ -800,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_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3583)
@@ -8,5 +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  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
+   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -15,12 +15,6 @@
    !!----------------------------------------------------------------------
    !!   PUBLIC
-   !!      tide_init     : read of namelist and initialisation of tidal harmonics data
+   !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data
    !!      tide_update   : calculation of tidal forcing at each timestep
-   !!   PRIVATE
-   !!      uvset         :\
-   !!      vday          : | Routines to correct tidal harmonics forcing for
-   !!      shpen         : | start time of integration
-   !!      ufset         : |
-   !!      vset          :/
    !!----------------------------------------------------------------------
    USE timing          ! Timing
@@ -34,4 +28,7 @@
    USE bdy_oce         ! ocean open boundary conditions
    USE daymod          ! calendar
+   USE wrk_nemo        ! Memory allocation
+   USE tideini
+!   USE tide_mod       ! Useless ??
    USE fldread, ONLY: fld_map
 
@@ -39,19 +36,18 @@
    PRIVATE
 
-   PUBLIC   tide_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
-      INTEGER                                ::   ncpt       !: Actual number of tidal components
-      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr)
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U
-      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.)
    END TYPE TIDES_DATA
 
-   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents
-
-   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data
-   
+   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -61,7 +57,7 @@
 CONTAINS
 
-   SUBROUTINE tide_init
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE tide_init  ***
+   SUBROUTINE bdytide_init
+      !!----------------------------------------------------------------------
+      !!                    ***  SUBROUTINE bdytide_init  ***
       !!                     
       !! ** Purpose : - Read in namelist for tides and initialise external
@@ -71,28 +67,32 @@
       !! namelist variables
       !!-------------------
-      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date
-      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files
-      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used.
-      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
+      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
+      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
+      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
       !!
-      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components
-      INTEGER                                   ::   ib_bdy, itide, ib  !: dummy loop indices
+      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)
-      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file 
-      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t           
-      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc
-      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data
+      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   
+      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
       !!
-      NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed
-      !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('tide_init')
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
+
+      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)
@@ -101,166 +101,154 @@
 
             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
-            ln_tide_date = .false.
             filtide(:) = ''
-            tide_speed(:) = 0.0
-            tide_cpt(:) = ''
 
             ! Don't REWIND here - may need to read more than one of these namelists.
             READ  ( numnam, nambdy_tide )
-            !                                               ! Count number of components specified
-            td%ncpt = 0
-            DO itide = 1, jptides_max
-              IF( tide_cpt(itide) /= '' ) THEN
-                 td%ncpt = td%ncpt + 1
-              ENDIF
-            END DO
-
-            ! Fill in phase speeds from namelist
-            ALLOCATE( td%speed(td%ncpt) )
-            td%speed = tide_speed(1:td%ncpt)
-
-            ! Find constituents in standard list
-            DO itide = 1, td%ncpt
-               nindx(itide) = 0
-               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
-               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
-               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
-               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
-               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
-               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
-               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
-               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
-               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
-               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
-               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
-               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
-               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
-               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
-               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
-               IF( nindx(itide) == 0  .AND. lwp ) THEN
-                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
-                  CALL ctl_warn( ctmp1 )
-               ENDIF
-            END DO
             !                                               ! Parameter control and print
-            IF( td%ncpt < 1 ) THEN 
-               CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' )
+            IF(lwp) WRITE(numout,*) '  '
+            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
+            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
+            ! -----------------------------------------------------------------------
+
+            ! 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
-               IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
-               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt
-               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt)
-               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
-               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt)
+               ilen0(:)=nblenrim(:)
             ENDIF
 
-
-            ! Allocate space for tidal harmonics data - 
-            ! get size from OBC data arrays
-            ! ---------------------------------------
-
-            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 
-            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) )
-
-            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 
-            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) )
-
-            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 
-            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) )
-
-            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) )
-
-
-            ! Open files and read in tidal forcing data
-            ! -----------------------------------------
-
-            DO itide = 1, td%ncpt
-               !                                                              ! SSH fields
-               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
-               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-               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%ssh(:,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%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1)
+            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
+            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
+
+            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
+            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
+
+            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
+            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
+
+            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(tide_cpt(itide))//'_grid_U.nc'
-               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-               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%u(:,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%u(:,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(tide_cpt(itide))//'_grid_V.nc'
-               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-               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%v(:,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%v(:,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 )
                !
-            END DO ! end loop on tidal components
-
-            IF( ln_tide_date ) THEN      ! correct for date factors
-
-!! used nmonth, nyear and nday from daymod....
-               ! Calculate date corrects for 15 standard consituents
-               ! This is the initialisation step, so nday, nmonth, nyear are the 
-               ! initial date/time of the integration.
-                 print *, nday,nmonth,nyear
-                 nyear  = int(ndate0 / 10000  )                          ! initial year
-                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month
-                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100)
-
-               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
-
-               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
-
-               DO itide = 1, td%ncpt       ! 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 )
                   !
-                  IF( nindx(itide) /= 0 ) THEN
-!!gm use rpi  and rad global variable  
-                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
-                     z_atde=z_ftc(nindx(itide))*cos(z_arg)
-                     z_btde=z_ftc(nindx(itide))*sin(z_arg)
-                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   &
-                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
-                  ELSE
-                     z_atde = 1.0_wp
-                     z_btde = 0.0_wp
-                  ENDIF
-                  !                                         !  elevation         
-                  igrd = 1
-                  DO ib = 1, ilen0(igrd)                
-                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2)
-                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1)
-                     td%ssh(ib,itide,1) = z1t
-                     td%ssh(ib,itide,2) = z2t
-                  END DO
-                  !                                         !  u       
-                  igrd = 2
-                  DO ib = 1, ilen0(igrd)                
-                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2)
-                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1)
-                     td%u(ib,itide,1) = z1t
-                     td%u(ib,itide,2) = z2t
-                  END DO
-                  !                                         !  v       
-                  igrd = 3
-                  DO ib = 1, ilen0(igrd)                
-                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2)
-                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1)
-                     td%v(ib,itide,1) = z1t
-                     td%v(ib,itide,2) = z2t
-                  END DO
-                  !
-               END DO   ! end loop on tidal components
-               !
-            ENDIF ! date correction
+               END DO ! end loop on tidal components
+               !
+               DEALLOCATE( dta_read )
+            ENDIF ! ln_bdytide_2ddta=.true.
+            !
+            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
@@ -268,36 +256,59 @@
       END DO ! loop on ib_bdy
 
-      IF( nn_timing == 1 ) CALL timing_stop('tide_init')
-
-   END SUBROUTINE tide_init
-
-
-   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE tide_update  ***
+      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
+
+   END SUBROUTINE bdytide_init
+
+   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. 
       !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter
-!!gm doctor jit ==> kit
-      TYPE(OBC_INDEX), INTENT( in )  ::   idx     ! OBC indices
-      TYPE(OBC_DATA),  INTENT(inout) ::   dta     ! OBC external data
-      TYPE(TIDES_DATA),INTENT( in )  ::   td      ! tidal harmonics data
-      INTEGER,INTENT(in),OPTIONAL    ::   jit     ! Barotropic timestep counter (for timesplitting option)
-      INTEGER,INTENT( in ), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
-                                                       ! is present then units = subcycle timesteps.
-                                                       ! time_offset = 0 => get data at "now" time level
-                                                       ! time_offset = -1 => get data at "before" time level
-                                                       ! time_offset = +1 => get data at "after" time level
-                                                       ! etc.
+      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
+      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices
+      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data
+      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data
+      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
+                                                        ! is present then units = subcycle timesteps.
+                                                        ! time_offset = 0  => get data at "now"    time level
+                                                        ! time_offset = -1 => get data at "before" time level
+                                                        ! time_offset = +1 => get data at "after"  time level
+                                                        ! etc.
       !!
-      INTEGER                          :: itide, igrd, ib      ! dummy loop indices
-      INTEGER                          :: time_add             ! time offset in units of timesteps
-      REAL(wp)                         :: z_arg, z_sarg      
-      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
-      !!----------------------------------------------------------------------
-
-      IF( nn_timing == 1 ) CALL timing_start('tide_update')
+      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, zramp      
+      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
+      !!----------------------------------------------------------------------
+
+      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
+      IF ( PRESENT(jit) ) THEN
+        IF ( jit /= 1 ) zflag=0
+      ENDIF
+
+      IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN
+        !
+        kt_tide = kt
+        !
+        IF(lwp) THEN
+           WRITE(numout,*)
+           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
+           WRITE(numout,*) '~~~~~~~~~~~~~~ '
+        ENDIF
+        !
+        CALL tide_init_elevation ( idx, td )
+        CALL tide_init_velocities( idx, td )
+        !
+      ENDIF 
 
       time_add = 0
@@ -306,429 +317,144 @@
       ENDIF
          
-      ! Note tide phase speeds are in deg/hour, so we need to convert the
-      ! elapsed time in seconds to hours by dividing by 3600.0
       IF( PRESENT(jit) ) THEN  
-         z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
+         z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) )
       ELSE                              
-         z_arg = (kt+time_add) * rdt * rad / 3600.0
+         z_arg = ((kt-kt_tide)+time_add) * rdt
       ENDIF
 
-      DO itide = 1, td%ncpt
-         z_sarg = z_arg * td%speed(itide)
+      ! 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)
          z_cost(itide) = COS( z_sarg )
          z_sist(itide) = SIN( z_sarg )
       END DO
 
-      DO itide = 1, td%ncpt
-         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)
-            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2)
+      DO itide = 1, nb_harmo
+         igrd=1                              ! SSH on tracer grid
+         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)
-            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2)
+         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)
-            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2)
+         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
-
-!!gm  doctor naming of dummy argument variables!!!   and all local variables
-
-   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
-      !!----------------------------------------------------------------------
-      !!                   ***  SUBROUTINE uvset  ***
-      !!                     
-      !! ** Purpose : - adjust tidal forcing for date factors
-      !!                
-      !!----------------------------------------------------------------------
-      implicit none
-      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
-      INTEGER, INTENT( in ) ::   iday    ! start time days
-      INTEGER, INTENT( in ) ::   imnth   ! start time month
-      INTEGER, INTENT( in ) ::   iyr     ! start time year
-      !!
-!!gm  nc is jptides_max ????
-      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents
-      CHARACTER(len=8), DIMENSION(nc) :: cname
-      INTEGER                 ::   year, vd, ivdy, ndc, i, k
-      REAL(wp)                ::   ss, h, p, en, p1, rtd
-      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction
-      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction
-      REAL(wp), DIMENSION(nc) ::   u, v, zig
-      !!
-      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   &
-         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   &
-         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      /
-      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   &
-         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   &
-         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   /
-      !!----------------------------------------------------------------------
-!
-!     ihs             -  start time gmt on ...
-!     iday/imnth/iyr  -  date   e.g.   12/10/87
-!
-      CALL vday(iday,imnth,iyr,ivdy)
-      vd = ivdy
-!
-!rp   note change of year number for d. blackman shpen
-!rp   if(iyr.ge.1000) year=iyr-1900
-!rp   if(iyr.lt.1000) year=iyr
-      year = iyr
-!
-!.....year  =  year of required data
-!.....vd    =  day of required data..set up for 0000gmt day year
-      ndc = nc
-!.....ndc   =  number of constituents allowed
-!!gm use rpi ?
-      rtd = 360.0 / 6.2831852
-      DO i = 1, ndc
-         zig(i) = zig(i)*rtd
-         ! sigo(i)= zig(i)
+   END SUBROUTINE bdytide_update
+
+   SUBROUTINE tide_init_elevation( idx, td )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init_elevation  ***
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
+      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.
+   
+      ilen0(1) =  SIZE(td%ssh0(:,1,1))
+
+      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
+
+      DO itide = 1, nb_harmo
+         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 , ilen0(igrd)
+            mod_tide(ib)=mod_tide(ib)*ftide(itide)
+            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+         ENDDO
+         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))
+         ENDDO
       END DO
 
-!!gm try to avoid RETURN  in F90
-      IF( year == 0 )   RETURN
-      
-         CALL shpen( year, vd, ss, h , p , en, p1 )
-         CALL ufset( p   , en, u , f )
-         CALL vset ( ss  , h , p , en, p1, v )
-         !
-         DO k = 1, nc
-            z_vplu(k) = v(k) + u(k)
-            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
-            DO WHILE( z_vplu(k) < 0    )
-               z_vplu(k) = z_vplu(k) + 360.0
-            END DO
-            DO WHILE( z_vplu(k) > 360. )
-               z_vplu(k) = z_vplu(k) - 360.0
-            END DO
-         END DO
-         !
-      END SUBROUTINE uvset
-
-
-      SUBROUTINE vday( iday, imnth, iy, ivdy )
-      !!----------------------------------------------------------------------
-      !!                   *** SUBROUTINE vday  ***
-      !!                  
-      !! ** Purpose : - adjust tidal forcing for date factors
-      !!                
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ????
-      INTEGER, INTENT(  out) ::   ivdy              ! ???
-      !! 
-      INTEGER ::   iyr
-      !!------------------------------------------------------------------------------
-
-!!gm   nday_year in day mode is the variable compiuted here, no?
-!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year
-
-      !calculate day number in year from day/month/year
-       if(imnth.eq.1) ivdy=iday
-       if(imnth.eq.2) ivdy=iday+31
-       if(imnth.eq.3) ivdy=iday+59
-       if(imnth.eq.4) ivdy=iday+90
-       if(imnth.eq.5) ivdy=iday+120
-       if(imnth.eq.6) ivdy=iday+151
-       if(imnth.eq.7) ivdy=iday+181
-       if(imnth.eq.8) ivdy=iday+212
-       if(imnth.eq.9) ivdy=iday+243
-       if(imnth.eq.10) ivdy=iday+273
-       if(imnth.eq.11) ivdy=iday+304
-       if(imnth.eq.12) ivdy=iday+334
-       iyr=iy
-       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
-       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
-       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
-      !
-   END SUBROUTINE vday
-
-!!doctor norme for dummy arguments
-
-   SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE shpen  ***
-      !!                   
-      !! ** Purpose : - calculate astronomical arguments for tides
-      !!                this version from d. blackman 30 nove 1990
-      !!
-      !!----------------------------------------------------------------------
-!!gm add INTENT in, out or inout....    DOCTOR name....
-!!gm please do not use variable name with a single letter:  impossible to search in a code
-      INTEGER  ::   year,vd
-      REAL(wp) ::   s,h,p,en,p1
-      !!
-      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd
-      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
-      !!
-      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           &
-         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    &
-         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    &
-         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    &
-         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    &
-         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    &
-         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    &
-         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    &
-         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    &
-         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    &
-         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    &
-         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    &
-         &       52.57 /
-      !!----------------------------------------------------------------------
-
-      cycle = 360.0
-      ilc = 0
-      icent = year / 100
-      yr = year - icent * 100
-      t = icent - 20
-!
-!     for the following equations
-!     time origin is fixed at 00 hr of jan 1st,2000.
-!     see notes by cartwright
-!
-!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90)
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      it = icent - 20
-      if (it) 1,2,2
-    1 iday = it/4 -it
-      go to 3
-    2 iday = (it+3)/4 - it
-!
-!     t is in julian century
-!     correction in gegorian calander where only century year divisible
-!     by 4 is leap year.
-!
-    3 continue
-!
-      td = 0.0
-!
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      if (yr) 4,5,4
-!
-    4 iyd = 365*yr
-      ild = (yr-1)/4
-      if((icent - (icent/4)*4) .eq. 0) ilc = 1
-      td = iyd + ild + ilc
-!
-    5 td = td + iday + vd -1.0 - 0.5
-      t = t + (td/36525.0)
-!
-      ipos=year-1899
-      if (ipos .lt. 0) go to 7
-      if (ipos .gt. 83) go to 6
-!
-      delta = (delt(ipos+1)+delt(ipos))/2.0
-      go to 7
-!
-    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
-!
-    7 deltat = delta * 1.0e-6
-!
-!!gm   precision of the computation   :  example for s it should be replace by:
-!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results
-      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
-      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat
-      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat
-      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat
-      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t
-      !
-      nn = s / cycle
-      s = s - nn * cycle
-      IF( s < 0.e0 )   s = s + cycle
-      !
-      nn = h / cycle
-      h = h - cycle * nn
-      IF( h < 0.e0 )   h = h + cycle
-      !
-      nn = p / cycle
-      p = p - cycle * nn
-      IF( p < 0.e0)   p = p + cycle
-      !
-      nn = en / cycle
-      en = en - cycle * nn
-      IF( en < 0.e0 )   en = en + cycle
-      en = cycle - en
-      !
-      nn = p1 / cycle
-      p1 = p1 - nn * cycle
-      !
-   END SUBROUTINE shpen
-
-
-   SUBROUTINE ufset( p, cn, b, a )
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE ufset  ***
-      !!                    
-      !! ** Purpose : - calculate nodal parameters for the tides
-      !!                
-      !!----------------------------------------------------------------------
-!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!gm  nc is jptides_max ????
-      integer nc
-      parameter (nc=15)
-      REAL(wp) p,cn
-      !!
-      
-!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z"
-      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
-      REAL(wp) ::   a(nc), b(nc)
-      INTEGER  ::   ndc, k
-      !!----------------------------------------------------------------------
-
-      ndc = nc
-
-!    a=f       ,  b =u
-!    t is  zero as compared to tifa.
-!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module
-      rad = 6.2831852d0/360.0
-      pw = p  * rad
-      nw = cn * rad
-      w1 = cos(   nw )
-      w2 = cos( 2*nw )
-      w3 = cos( 3*nw )
-      w4 = sin(   nw )
-      w5 = sin( 2*nw )
-      w6 = sin( 3*nw )
-      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   &
-         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw )
-      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   &
-         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw )
-      !
-      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
-      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
-      !   q1
-      a(2) = a(1)
-      b(2) = b(1)
-      !   o1
-      a(3) = 1.0
-      b(3) = 0.0
-      !   p1
-      a(4) = 1.0
-      b(4) = 0.0
-      !   s1
-      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
-      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
-      !   k1
-      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
-      b(6) = -0.0374*w4
-      !  2n2
-      a(7) = a(6)
-      b(7) = b(6)
-      !  mu2
-      a(8) = a(6)
-      b(8) = b(6)
-      !   n2
-      a(9) = a(6)
-      b(9) = b(6)
-      !  nu2
-      a(10) = a(6)
-      b(10) = b(6)
-      !   m2
-      a(11) = SQRT( w7 * w7 + w8 * w8 )
-      b(11) = ATAN( w8 / w7 )
-!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ????
-      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992
-      !   l2
-      a(12) = 1.0
-      b(12) = 0.0
-      !   t2
-      a(13)= a(12)
-      b(13)= b(12)
-      !   s2
-      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
-      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
-      !   k2
-      a(15) = a(6)*a(6)
-      b(15) = 2*b(6)
-      !   m4
-!!gm  old coding,  remove GOTO and label of lines
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      DO 40 k = 1,ndc
-         b(k) = b(k)/rad
-32       if (b(k)) 34,35,35
-34       b(k) = b(k) + 360.0
-         go to 32
-35       if (b(k)-360.0) 40,37,37
-37       b(k) = b(k)-360.0
-         go to 35
-40    continue
-      !
-   END SUBROUTINE ufset
-   
-
-   SUBROUTINE vset( s,h,p,en,p1,v)
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE vset  ***
-      !!                   
-      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
-      !!                
-      !!----------------------------------------------------------------------
-!!gm  doctor naming of dummy argument variables!!!   and all local variables
-!!gm  nc is jptides_max ????
-!!gm  en argument is not used: suppress it ?
-      integer nc
-      parameter (nc=15)
-      real(wp) s,h,p,en,p1
-      real(wp)   v(nc)
-      !!
-      integer ndc, k
-      !!----------------------------------------------------------------------
-
-      ndc = nc
-      !   v s  are computed here.
-      v(1) =-3*s +h +p +270      ! Q1
-      v(2) =-2*s +h +270.0       ! O1
-      v(3) =-h +270              ! P1
-      v(4) =180                  ! S1
-      v(5) =h +90.0              ! K1
-      v(6) =-4*s +2*h +2*p       ! 2N2
-      v(7) =-4*(s-h)             ! MU2
-      v(8) =-3*s +2*h +p         ! N2
-      v(9) =-3*s +4*h -p         ! MU2
-      v(10) =-2*s +2*h           ! M2
-      v(11) =-s +2*h -p +180     ! L2 
-      v(12) =-h +p1              ! T2
-      v(13) =0                   ! S2
-      v(14) =h+h                 ! K2
-      v(15) =2*v(10)             ! M4
-!
-!!gm  old coding,  remove GOTO and label of lines
-!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
-      do 72 k = 1, ndc
-69       if( v(k) )   70,71,71
-70       v(k) = v(k)+360.0
-         go to 69
-71       if( v(k) - 360.0 )   72,73,73
-73       v(k) = v(k)-360.0
-         go to 71
-72    continue
-      !
-   END SUBROUTINE vset
-
+      DEALLOCATE(mod_tide,phi_tide)
+
+   END SUBROUTINE tide_init_elevation
+
+   SUBROUTINE tide_init_velocities( idx, td )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE tide_init_elevation  ***
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
+      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(ilen0(igrd)),phi_tide(ilen0(igrd)))
+
+      DO itide = 1, nb_harmo
+         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, ilen0(igrd)
+            mod_tide(ib)=mod_tide(ib)*ftide(itide)
+            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+         ENDDO
+         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))
+         ENDDO
+      END DO
+
+      DEALLOCATE(mod_tide,phi_tide)
+
+      igrd=3                                 ! V grid.
+
+      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
+
+      DO itide = 1, nb_harmo
+         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, ilen0(igrd)
+            mod_tide(ib)=mod_tide(ib)*ftide(itide)
+            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+         ENDDO
+         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))
+         ENDDO
+      END DO
+
+      DEALLOCATE(mod_tide,phi_tide)
+
+  END SUBROUTINE tide_init_velocities
 #else
    !!----------------------------------------------------------------------
    !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
    !!----------------------------------------------------------------------
-!!gm  are you sure we need to define filtide and tide_cpt ?
-   CHARACTER(len=80), PUBLIC               ::   filtide                !: Filename root for tidal input files
-   CHARACTER(len=4 ), PUBLIC, DIMENSION(1) ::   tide_cpt               !: Names of tidal components used.
-
 CONTAINS
-   SUBROUTINE tide_init                ! Empty routine
-   END SUBROUTINE tide_init
-   SUBROUTINE tide_data                ! Empty routine
-   END SUBROUTINE tide_data
-   SUBROUTINE tide_update( kt, kit )   ! Empty routine
-      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit
-   END SUBROUTINE tide_update
+   SUBROUTINE bdytide_init             ! Empty routine
+      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?'
+   END SUBROUTINE bdytide_init
+   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_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3583)
@@ -23,8 +23,10 @@
    USE in_out_manager  ! I/O manager
 
+
    IMPLICIT NONE
    PRIVATE
 
    PUBLIC bdy_tra      ! routine called in tranxt.F90 
+   PUBLIC bdy_tra_dmp  ! routine called in step.F90 
 
    !!----------------------------------------------------------------------
@@ -53,8 +55,18 @@
          CASE(jp_frs)
             CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(2)
+            CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(3)
+            CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE(4)
+            CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
          CASE DEFAULT
             CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
          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
@@ -90,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 )
       !
@@ -97,5 +107,170 @@
       !
    END SUBROUTINE bdy_tra_frs
-   
+  
+   SUBROUTINE bdy_tra_spe( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_frs  ***
+      !!                    
+      !! ** Purpose : Apply a specified value for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_spe')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            tsa(ii,ij,ik,jp_tem) = dta%tem(ib,ik) * tmask(ii,ij,ik)
+            tsa(ii,ij,ik,jp_sal) = dta%sal(ib,ik) * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_spe')
+      !
+   END SUBROUTINE bdy_tra_spe
+
+   SUBROUTINE bdy_tra_nmn( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_nmn  ***
+      !!                    
+      !! ** Purpose : Duplicate the value for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_nmn')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ! search the sense of the gradient
+            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
+            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
+            IF ( zcoef1+zcoef2 == 0) THEN
+               ! corner
+               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
+               tsa(ii,ij,ik,jp_tem) = tsa(ii-1,ij  ,ik,jp_tem) * tmask(ii-1,ij  ,ik) + &
+                 &                    tsa(ii+1,ij  ,ik,jp_tem) * tmask(ii+1,ij  ,ik) + &
+                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + &
+                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik)
+               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + &
+                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + &
+                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + &
+                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik)
+               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+            ELSE
+               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik)
+               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik)
+            ENDIF
+         END DO
+      END DO
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_nmn')
+      !
+   END SUBROUTINE bdy_tra_nmn
+
+   SUBROUTINE bdy_tra_rnf( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_rnf  ***
+      !!                    
+      !! ** Purpose : Apply the runoff values for tracers at open boundaries:
+      !!                  - specified to 0.1 PSU for the salinity
+      !!                  - duplicate the value for the temperature
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_rnf')
+      !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+            jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+            tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii,ij,ik)
+            tsa(ii,ij,ik,jp_sal) =                        0.1 * tmask(ii,ij,ik)
+         END DO
+      END DO
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_rnf')
+      !
+   END SUBROUTINE bdy_tra_rnf
+
+   SUBROUTINE bdy_tra_dmp( kt )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra_dmp  ***
+      !!                    
+      !! ** Purpose : Apply damping for tracers at open boundaries.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,         INTENT(in) ::   kt
+      !! 
+      REAL(wp) ::   zwgt           ! boundary weight
+      REAL(wp) ::   zta, zsa, ztime
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! 2D addresses
+      INTEGER  ::   ib_bdy         ! Loop index
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_tra_dmp')
+      !
+      DO ib_bdy=1, nb_bdy
+         IF ( ln_tra_dmp(ib_bdy) ) THEN
+            igrd = 1                       ! Everything is at T-points here
+            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)
+               zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
+               DO ik = 1, jpkm1
+                  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
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp')
+      !
+   END SUBROUTINE bdy_tra_dmp
+ 
 #else
    !!----------------------------------------------------------------------
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 3583)
@@ -403,5 +403,5 @@
          IF( lk_obc )   CALL obc_dta_bt ( kt, jn   )
          IF( lk_bdy )   CALL bdy_dta ( kt, jit=jn, time_offset=+1 )
-         IF ( ln_tide_pot ) CALL upd_tide( kt, jn )
+         IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn )
 
          !                                                !* after ssh_e
@@ -453,5 +453,5 @@
                   ENDIF
                   ! add tidal astronomical forcing
-                  IF ( ln_tide_pot ) THEN 
+                  IF ( ln_tide_pot .AND. lk_tide ) THEN 
                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
@@ -503,5 +503,5 @@
                   ENDIF
                   ! add tidal astronomical forcing
-                  IF ( ln_tide_pot ) THEN
+                  IF ( ln_tide_pot .AND. lk_tide ) THEN
                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
@@ -550,5 +550,5 @@
                   ENDIF
                   ! add tidal astronomical forcing
-                  IF ( ln_tide_pot ) THEN
+                  IF ( ln_tide_pot .AND. lk_tide ) THEN
                   zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)
                   zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 3583)
@@ -667,5 +667,5 @@
       !!----------------------------------------------------------------------
 #if defined key_bdy
-      USE bdy_oce, ONLY:  dta_global         ! workspace to read in global data arrays
+      USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays
 #endif 
 
@@ -681,13 +681,9 @@
       INTEGER                                 ::   ilendta  ! length of data in file
       INTEGER                                 ::   idvar    ! variable ID
-      INTEGER                                 ::   ib, ik   ! loop counters
+      INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
       INTEGER                                 ::   ierr
-      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data
       !!---------------------------------------------------------------------
             
-#if defined key_bdy
-      dta_read => dta_global
-#endif
-
       ipi = SIZE( dta, 1 )
       ipj = 1
@@ -696,6 +692,19 @@
       idvar   = iom_varid( num, clvar )
       ilendta = iom_file(num)%dimsz(1,idvar)
+
+#if defined key_bdy
+      ipj = iom_file(num)%dimsz(2,idvar)
+      IF (ipj == 1) THEN ! we assume that this is a structured open boundary file
+         dta_read => dta_global
+      ELSE
+         dta_read => dta_global2
+      ENDIF
+#endif
+
+      IF ((lwp).AND.(ipi>ilendta*ipj)) CALL ctl_stop('fld_map problem with file dimension')
+
       IF(lwp) WRITE(numout,*) 'Dim size for ',TRIM(clvar),' is ', ilendta
       IF(lwp) WRITE(numout,*) 'Number of levels for ',TRIM(clvar),' is ', ipk
+
 
       SELECT CASE( ipk )
@@ -706,9 +715,19 @@
       END SELECT
       !
-      DO ib = 1, ipi
-         DO ik = 1, ipk
-            dta(ib,1,ik) =  dta_read(map(ib),1,ik)
+      IF (ipj==1) THEN
+         DO ib = 1, ipi
+            DO ik = 1, ipk
+               dta(ib,1,ik) =  dta_read(map(ib),1,ik)
+            END DO
          END DO
-      END DO
+      ELSE ! we assume that this is a structured open boundary file
+         DO ib = 1, ipi
+            jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
+            ji=map(ib)-(jj-1)*ilendta
+            DO ik = 1, ipk
+               dta(ib,1,ik) =  dta_read(ji,jj,ik)
+            END DO
+         END DO
+      ENDIF
 
    END SUBROUTINE fld_map
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 3583)
@@ -66,10 +66,9 @@
       !!
       INTEGER            ::   ierror  ! local integer 
-      REAL(wp)           ::   zpref   ! local scalar
       !!
       CHARACTER(len=100) ::  cn_dir   ! Root directory for location of ssr files
       TYPE(FLD_N)        ::  sn_apr   ! informations about the fields to be read
       !!
-      NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr
+      NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr, rpref, ln_apr_obc
       !!----------------------------------------------------------------------
       !
@@ -113,8 +112,7 @@
          !
          !                                            !* control check
-         IF( ln_apr_obc  )   &
-            CALL ctl_stop( 'sbc_apr: inverse barometer added to OBC ssh data not yet implemented ' )
-         IF( ln_apr_obc .AND. .NOT. lk_obc )   &
-            CALL ctl_stop( 'sbc_apr: add inverse barometer to OBC requires to use key_obc' )
+         IF ( ln_apr_obc  ) THEN
+            IF(lwp) WRITE(numout,*) '         Inverse barometer added to OBC ssh data'
+         ENDIF
          IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts )   &
             CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' )
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90	(revision 3583)
@@ -14,5 +14,5 @@
   USE daymod
   USE dynspg_oce
-  USE tide_mod
+  USE tideini
   USE iom
 
@@ -21,22 +21,9 @@
 
   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro
-  LOGICAL, PUBLIC :: ln_tide_pot = .false.
+
 #if defined key_tide
 
   LOGICAL, PUBLIC, PARAMETER ::   lk_tide  = .TRUE.
-
-  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide 
-
-  REAL(wp), ALLOCATABLE, DIMENSION(:) ::  &
-       v0tide,      &
-       utide,       &
-       ftide
-
   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot,phi_pot
-
-  INTEGER, PUBLIC :: nb_harmo
-  INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide
-  INTEGER, PUBLIC :: nn_tide, kt_tide
-
   !!---------------------------------------------------------------------------------
   !!   OPA 9.0 , LODYC-IPSL  (2003)
@@ -51,65 +38,34 @@
     !! * Arguments
     INTEGER, INTENT( in ) ::   kt     ! ocean time-step
-    !! * Local declarations
-    INTEGER  :: jk,ji
-    CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
     !!----------------------------------------------------------------------
 
-    NAMELIST/nam_tide/ln_tide_pot,nb_harmo,clname,nn_tide
+    IF ( kt == nit000 .AND. .NOT. lk_dynspg_ts )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' )
 
-    IF ( kt == nit000 ) THEN
+    IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN
+      !
+      kt_tide = kt
 
-       IF( .NOT. lk_dynspg_ts  )  CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' )
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'sbc_tide : (re)Initialization of the tidal potential at kt=',kt
+         WRITE(numout,*) '~~~~~~~ '
+      ENDIF
 
-    ! Read Namelist nam_tide
+      IF(lwp) THEN
+         IF ( kt == nit000 ) WRITE(numout,*) 'Apply astronomical potential : ln_tide_pot =', ln_tide_pot
+         CALL flush(numout)
+      ENDIF
 
-    nn_tide=INT(rday/rdt)
+      IF ( kt == nit000 ) ALLOCATE(amp_pot(jpi,jpj,nb_harmo))
+      IF ( kt == nit000 ) ALLOCATE(phi_pot(jpi,jpj,nb_harmo))
+      IF ( kt == nit000 ) ALLOCATE(pot_astro(jpi,jpj))
 
-    CALL tide_init_Wave
+      amp_pot(:,:,:) = 0.e0
+      phi_pot(:,:,:) = 0.e0
+      pot_astro(:,:) = 0.e0
 
-    REWIND ( numnam )
-    READ   ( numnam, nam_tide )
-
-    IF(lwp) THEN
-       WRITE(numout,*)
-       WRITE(numout,*) 'sbc_tide : Initialization of the tidal components'
-       WRITE(numout,*) '~~~~~~~ '
+      IF ( ln_tide_pot ) CALL tide_init_potential
+      !
     ENDIF
-
-    IF(lwp) THEN
-       WRITE(numout,*) '        Namelist nam_tide'
-       WRITE(numout,*) '        Apply astronomical potential : ln_tide_pot =', ln_tide_pot
-       WRITE(numout,*) '        nb_harmo    = ', nb_harmo
-       CALL flush(numout)
-    ENDIF
-
-    ALLOCATE(ntide     (nb_harmo))
-    DO jk=1,nb_harmo
-       DO ji=1,jpmax_harmo
-          IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN
-             ntide(jk) = ji
-             EXIT
-          END IF
-       END DO
-    END DO
-    ALLOCATE(omega_tide(nb_harmo))
-    ALLOCATE(v0tide    (nb_harmo))
-    ALLOCATE(utide     (nb_harmo))
-    ALLOCATE(ftide     (nb_harmo))
-    ALLOCATE(amp_pot(jpi,jpj,nb_harmo))
-    ALLOCATE(phi_pot(jpi,jpj,nb_harmo))
-    ALLOCATE(pot_astro(jpi,jpj))
-    ENDIF
-
-    IF ( MOD( kt - 1, nn_tide ) == 0 ) THEN
-      kt_tide = kt
-      CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo)
-    ENDIF
-
-    amp_pot(:,:,:) = 0.e0
-    phi_pot(:,:,:) = 0.e0
-    pot_astro(:,:) = 0.e0
-
-    IF (ln_tide_pot          ) CALL tide_init_potential
 
   END SUBROUTINE sbc_tide
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 3583)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/tideini.F90	(revision 3583)
@@ -0,0 +1,126 @@
+MODULE tideini
+  !!=================================================================================
+  !!                       ***  MODULE  tideini  ***
+  !! Initialization of tidal forcing
+  !! History :  9.0  !  07  (O. Le Galloudec)  Original code
+  !!=================================================================================
+  !! * Modules used
+  USE oce             ! ocean dynamics and tracers variables
+  USE dom_oce         ! ocean space and time domain
+  USE in_out_manager  ! I/O units
+  USE ioipsl          ! NetCDF IPSL library
+  USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+  USE phycst
+  USE daymod
+  USE dynspg_oce
+  USE tide_mod
+  USE iom
+
+  IMPLICIT NONE
+  PUBLIC
+
+  REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::  &
+       omega_tide,  &
+       v0tide,      &
+       utide,       &
+       ftide
+
+  LOGICAL, PUBLIC :: ln_tide_pot = .false., ln_tide_ramp = .false.
+  REAL(wp), PUBLIC :: rdttideramp 
+  INTEGER, PUBLIC :: nb_harmo
+  INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide
+  INTEGER, PUBLIC :: kt_tide
+
+  !!---------------------------------------------------------------------------------
+  !!   OPA 9.0 , LODYC-IPSL  (2003)
+  !!---------------------------------------------------------------------------------
+
+CONTAINS
+   
+  SUBROUTINE tide_init ( kt )
+    !!----------------------------------------------------------------------
+    !!                 ***  ROUTINE tide_init  ***
+    !!----------------------------------------------------------------------      
+    !! * Local declarations
+    INTEGER  :: ji, jk
+    INTEGER, INTENT( in ) ::   kt     ! ocean time-step
+    CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
+    !
+    NAMELIST/nam_tide/ln_tide_pot, ln_tide_ramp, rdttideramp, clname
+    !!----------------------------------------------------------------------
+
+    IF ( kt == nit000 ) THEN
+       !
+       IF(lwp) THEN
+          WRITE(numout,*)
+          WRITE(numout,*) 'tide_init : Initialization of the tidal components'
+          WRITE(numout,*) '~~~~~~~~~ '
+       ENDIF
+       !
+       CALL tide_init_Wave
+       !
+       clname(:)=''
+       !
+       ! Read Namelist nam_tide
+       REWIND ( numnam )
+       READ   ( numnam, nam_tide )
+       !
+       nb_harmo=0
+       DO jk=1,jpmax_harmo
+          DO ji=1,jpmax_harmo
+             IF(TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN
+                nb_harmo=nb_harmo+1
+             ENDIF
+          END DO
+       ENDDO
+       !
+       IF(lwp) THEN
+          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
+       !
+       ALLOCATE(ntide(nb_harmo))
+       DO jk=1,nb_harmo
+          DO ji=1,jpmax_harmo
+             IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN
+                ntide(jk) = ji
+                EXIT
+             END IF
+          END DO
+       END DO
+       !
+       ALLOCATE(omega_tide(nb_harmo))
+       ALLOCATE(v0tide    (nb_harmo))
+       ALLOCATE(utide     (nb_harmo))
+       ALLOCATE(ftide     (nb_harmo))
+       kt_tide = kt
+       !
+    ENDIF
+
+    IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN
+       !
+       IF(lwp) THEN
+          WRITE(numout,*)
+          WRITE(numout,*) 'tide_ini : Update of the tidal components at kt=',kt
+          WRITE(numout,*) '~~~~~~~~ '
+       ENDIF
+       CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo)
+       DO jk =1,nb_harmo
+         IF(lwp) WRITE(numout,*) Wave(ntide(jk))%cname_tide,utide(jk),ftide(jk),v0tide(jk),omega_tide(jk)
+         call flush(numout)
+       END DO
+       !
+       kt_tide = kt
+       !
+    ENDIF
+
+  END SUBROUTINE tide_init
+   
+END MODULE tideini
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90	(revision 3583)
@@ -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_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 3583)
@@ -46,8 +46,9 @@
    USE mppini          ! shared/distributed memory setting (mpp_init routine)
    USE domain          ! domain initialization             (dom_init routine)
+   USE tideini         ! tidal components initialization   (tide_ini routine)
    USE obcini          ! open boundary cond. initialization (obc_ini routine)
    USE bdyini          ! open boundary cond. initialization (bdy_init routine)
    USE bdydta          ! open boundary cond. initialization (bdy_dta_init routine)
-   USE bdytides        ! open boundary cond. initialization (tide_init routine)
+   USE bdytides        ! open boundary cond. initialization (bdytide_init routine)
    USE istate          ! initial state setting          (istate_init routine)
    USE ldfdyn          ! lateral viscosity setting      (ldfdyn_init routine)
@@ -75,4 +76,6 @@
    USE mod_ioclient
 #endif
+   USE sbctide, ONLY: lk_tide
+
 
    IMPLICIT NONE
@@ -314,13 +317,16 @@
 
       IF( lk_obc        )   CALL     obc_init   ! Open boundaries 
+
+                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
+
+      IF( lk_tide       )   CALL tide_init( nit000 )    ! Initialisation of the tidal harmonics
+
       IF( lk_bdy        )   CALL     bdy_init       ! Open boundaries initialisation
       IF( lk_bdy        )   CALL     bdy_dta_init   ! Open boundaries initialisation of external data arrays
-      IF( lk_bdy        )   CALL     tide_init      ! Open boundaries initialisation of tidal harmonic forcing
+      IF( lk_bdy        )   CALL     bdytide_init   ! Open boundaries initialisation of tidal harmonic forcing
 
                             CALL flush(numout)
                             CALL dyn_nept_init  ! simplified form of Neptune effect
                             CALL flush(numout)
-
-                            CALL  istate_init   ! ocean initial state (Dynamics and tracers)
 
       !                                     ! Ocean physics
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 3583)
@@ -95,4 +95,5 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
+      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
@@ -187,4 +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         )   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
@@ -219,4 +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           )   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
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 3583)
@@ -55,5 +55,8 @@
 
    USE bdy_par          ! for lk_bdy
+   USE bdy_oce          ! for dmp logical
    USE bdydta           ! open boundary condition data     (bdy_dta routine)
+   USE bdytra           ! bdy cond. for tracers            (bdy_tra routine)
+   USE bdydyn3d         ! bdy cond. for baroclinic vel.  (bdy_dyn3d routine)
 
    USE sshwzv           ! vertical velocity and ssh        (ssh_wzv routine)
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_MERCATOR_CLUSTER
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_MERCATOR_CLUSTER	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/BATCH_TEMPLATE/batch-ifort_MERCATOR_CLUSTER	(revision 3583)
@@ -4,8 +4,8 @@
 #PBS -e sette.$PBS_JOBID.err
 #PBS -o sette.$PBS_JOBID.out
-#PBS -l nodes=2:ppn=8
-#PBS -q multi
+#PBS -l nodes=NODES:ppn=NPROCSNODE
+#PBS -q QUEUE
 #PBS -l walltime=03:00:00
-#PBS -l mem=24gb
+#PBS -l mem=MEMgb
 
 #
@@ -22,5 +22,6 @@
 # Local settings for machine IBM Power6 (VARGAS at IDRIS France)
 #
-export MPIRUN="mpiexec -n $OCEANCORES"
+#cbr export MPIRUN="mpiexec -n $OCEANCORES"
+export MPIRUN="mpirun -np $OCEANCORES"
 
 #
@@ -50,5 +51,4 @@
 #
 cd ${EXE_DIR}
-
   echo Running on host `hostname`
   echo Time is `date`
@@ -60,4 +60,5 @@
   echo "Running time ${MPIRUN} ./opa"
      time ${MPIRUN} ./opa
+#cbr     mpirun -np $OCEANCORES ./opa
   else
   echo "Running time ./opa"
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/param.cfg
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/param.cfg	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/param.cfg	(revision 3583)
@@ -1,4 +1,4 @@
 #- forcing files storing 
-FORCING_DIR=~/FORCING
+FORCING_DIR=/data/rd_exchange/cbricaud/NEMO_3.4_dev2012/FORCING
 #- input files storing 
 INPUT_DIR=${CONFIG_DIR}/${NEW_CONF}/EXP00
Index: /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/prepare_job.sh
===================================================================
--- /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/prepare_job.sh	(revision 3582)
+++ /branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/SETTE/prepare_job.sh	(revision 3583)
@@ -169,5 +169,5 @@
 					# number of processes required is an integer multiple of 4
 					#
-					NB_NODES=$( echo $NB_PROC | awk '{print $1 / 4}')
+					NB_NODES=$( echo $NB_PROC | awk '{print $1  $1 / 4}')
 				else
 					#
@@ -178,4 +178,27 @@
 	       			fi
 				;;
+                        ifort_MERCATOR_CLUSTER)
+                                echo NB_PROCS ${NB_PROC}
+                                echo NB_NODES ${NB_NODES}
+                                echo  NB_PROC ${NB_PROC} 
+                                if [ ${NB_PROC} -eq 1 ] ; then
+                                   NB_NODES=1
+                                   QUEUE=monoproc
+                                   NB_PROC_NODE=${NB_PROC} 
+                                else
+                                   if [ ${NB_PROC} -le 8 ] ; then 
+                                      NB_NODES=1
+                                      QUEUE=mono
+                                      NB_PROC_NODE=${NB_PROC} 
+                                   else
+                                      NB_NODES=$( echo $NB_PROC | awk '{print $1 - $1 % 8}'  | awk '{print $1 / 8 }')
+                                      QUEUE=multi
+                                      NB_PROC_NODE=8
+                                  fi    
+                                fi    
+                                echo NB_PROCS     ${NB_PROC}
+                                echo NB_NODES     ${NB_NODES}
+                                echo NB_PROC_NODE ${NB_PROC_NODE} 
+				;;
 			*)
 				NB_NODES=${NB_PROC}
@@ -186,5 +209,9 @@
 # Pass settings into job file by using sed to edit predefined strings
 #
-        cat ${SETTE_DIR}/job_batch_template | sed -e"s/NODES/${NB_NODES}/" -e"s/NPROCS/${NB_PROC}/" \
+         ((mem=5*NB_PROC))
+         echo NB_NODES NB_PROC QUEUE ${NB_NODES} ${NB_PROC} ${QUEUE}
+        cat ${SETTE_DIR}/job_batch_template | sed -e"s/NODES/${NB_NODES}/" -e"s/NPROCSNODE/${NB_PROC_NODE}/" \
+             -e"s/NPROCS/${NB_PROC}/" \
+             -e"s/QUEUE/${QUEUE}/" -e"s/MEM/${mem}/" \
              -e"s:DEF_SETTE_DIR:${SETTE_DIR}:" -e"s:DEF_INPUT_DIR:${INPUT_DIR}:" \
              -e"s:DEF_EXE_DIR:${EXE_DIR}:" \
