Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2865)
@@ -31,5 +31,5 @@
    USE obc_oce         ! ocean open boundary conditions
    USE obcdta          ! ocean open boundary conditions
-   USE obcdyn3d        ! ocean open boundary conditions
+   USE obcdyn          ! ocean open boundary conditions
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
    USE in_out_manager  ! I/O manager
@@ -153,7 +153,6 @@
 # if defined key_obc
       !                                !* OBC open boundaries
-      IF( .NOT. lk_dynspg_flt ) THEN
-
-         CALL obc_dyn3d( kt )
+      IF( lk_dynspg_exp ) CALL obc_dyn( kt )
+      IF( lk_dynspg_ts )  CALL obc_dyn( kt, dyn3d_only=.true. )
 
 !!$!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
@@ -163,6 +162,4 @@
 !!$         !
 !!$         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
-
-      ENDIF
       !
 # endif
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2865)
@@ -355,5 +355,5 @@
          !                                                !* Update the forcing (OBC and tides)
          !                                                !  ------------------
-         IF( lk_obc )   CALL obc_dta ( kt, jit=jn   )
+         IF( lk_obc )   CALL obc_dta ( kt, jit=jn, time_offset=+1 )
 
          !                                                !* after ssh_e
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 2865)
@@ -57,13 +57,7 @@
    LOGICAL                    ::   ln_mask_file             !: =T read obcmask from file
    LOGICAL                    ::   ln_vol                   !: =T volume correction             
-   LOGICAL, DIMENSION(jp_obc) ::   ln_clim                  !: =T obc data files contain climatological data (time-cyclic)
    !
    INTEGER                    ::   nb_obc                   !: number of open boundary sets
-   INTEGER, DIMENSION(jp_obc) ::   nn_rimwidth              !: boundary rim width
-   INTEGER, DIMENSION(jp_obc) ::   nn_dtactl      	    !: = 0 use the initial state as obc dta ; 
-                                                            !: = 1 read it in a NetCDF file
-   INTEGER, DIMENSION(jp_obc) ::   nn_tides                 !: = 0 no tidal harmonic forcing
-                                                            !: = 1 apply ONLY tidal harmonic forcing for barotropic solution
-                                                            !: = 2 ADD tidal harmonic forcing to other barotropic boundary data
+   INTEGER, DIMENSION(jp_obc) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
    !                                                        !  = 1 the volume will be constant during all the integration.
@@ -71,4 +65,6 @@
    INTEGER, DIMENSION(jp_obc) ::   nn_dyn2d_dta      	    !: = 0 use the initial state as obc dta ; 
                                                             !: = 1 read it in a NetCDF file
+                                                            !: = 2 read tidal harmonic forcing from a NetCDF file
+                                                            !: = 3 read external data AND tidal harmonic forcing from NetCDF files
    INTEGER, DIMENSION(jp_obc) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities 
    INTEGER, DIMENSION(jp_obc) ::   nn_dyn3d_dta      	    !: = 0 use the initial state as obc dta ; 
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 2865)
@@ -55,5 +55,5 @@
 CONTAINS
 
-      SUBROUTINE obc_dta( kt, jit )
+      SUBROUTINE obc_dta( kt, jit, time_offset )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE obc_dta  ***
@@ -69,4 +69,10 @@
       INTEGER, INTENT( in )           ::   kt    ! ocean time-step index 
       INTEGER, INTENT( in ), OPTIONAL ::   jit   ! subcycle time-step index (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     ::  ib_obc, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices
@@ -202,31 +208,31 @@
             IF( PRESENT(jit) ) THEN
                ! Update barotropic boundary conditions only
-               ! jit is optional argument for fld_read
-               IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) THEN
-                  IF( nn_tides(ib_obc) .eq. 1 ) THEN
+               ! jit is optional argument for fld_read and tide_update
+               IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN
+                  IF( nn_dyn2d_dta(ib_obc) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
                      dta_obc(ib_obc)%ssh(:) = 0.0
                      dta_obc(ib_obc)%u2d(:) = 0.0
                      dta_obc(ib_obc)%v2d(:) = 0.0
-                  ELSE
+                  ENDIF
+                  IF( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .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 )
+                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset )
                   ENDIF
-               ENDIF
-               IF( nn_tides(ib_obc) .gt. 0 ) THEN
-                  CALL tide_update( kt=kt, jit=jit, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) )
+                  IF( nn_dyn2d_dta(ib_obc) .ge. 2 ) THEN ! update tidal harmonic forcing
+                     CALL tide_update( kt=kt, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc), jit=jit, time_offset=time_offset )
+                  ENDIF
                ENDIF
             ELSE
-               IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN 
-                  jend = jstart + nb_obc_fld(ib_obc) - 1
-                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), timeshift=1 )
-               ENDIF
-               IF( nn_tides(ib_obc) .eq. 1 ) THEN
+               IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
                   dta_obc(ib_obc)%ssh(:) = 0.0
                   dta_obc(ib_obc)%u2d(:) = 0.0
                   dta_obc(ib_obc)%v2d(:) = 0.0
                ENDIF
-               IF( nn_tides(ib_obc) .gt. 0 ) THEN
-                  !!! THINK ABOUT kt, jit VALUES !!!
-                  CALL tide_update( kt=kt, jit=0, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc) )
+               IF( nb_obc_fld(ib_obc) .gt. 0 ) THEN ! update external data
+                  jend = jstart + nb_obc_fld(ib_obc) - 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_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .ge. 2 ) THEN ! update tidal harmonic forcing
+                  CALL tide_update( kt=kt, idx=idx_obc(ib_obc), dta=dta_obc(ib_obc), td=tides(ib_obc), time_offset=time_offset )
                ENDIF
             ENDIF
@@ -238,5 +244,5 @@
 
             IF( ln_full_vel_array(ib_obc) .and.                                             & 
-           &    ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn3d_dta(ib_obc) .eq. 1 ) ) THEN 
+           &    ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 .or. nn_dyn3d_dta(ib_obc) .eq. 1 ) ) THEN 
 
                igrd = 2                      ! zonal velocity
@@ -326,4 +332,5 @@
 #endif
                               )
+         IF(nn_dta(ib_obc) .gt. 1) nn_dta(ib_obc) = 1
       END DO
 
@@ -333,5 +340,5 @@
       nb_obc_fld(:) = 0
       DO ib_obc = 1, nb_obc         
-         IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) THEN
+         IF( nn_dyn2d(ib_obc) .gt. 0 .and. ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) THEN
             nb_obc_fld(ib_obc) = nb_obc_fld(ib_obc) + 3
          ENDIF
@@ -408,7 +415,7 @@
             ! Only read in necessary fields for this set.
             ! Important that barotropic variables come first.
-            IF( nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) THEN 
-
-               IF( nn_dyn2d(ib_obc) .ne. jp_frs .and. nn_tides(ib_obc) .ne. 1) THEN
+            IF( nn_dyn2d(ib_obc) .gt. 0 .and. ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) THEN 
+
+               IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
@@ -419,5 +426,5 @@
                ENDIF
 
-               IF( .not. ln_full_vel_array(ib_obc) .and. nn_tides(ib_obc) .ne. 1 ) THEN
+               IF( .not. ln_full_vel_array(ib_obc) ) THEN
 
                   jfld = jfld + 1
@@ -449,5 +456,6 @@
             ! baroclinic velocities
             IF( ( nn_dyn3d(ib_obc) .gt. 0 .and. nn_dyn3d_dta(ib_obc) .eq. 1 ) .or. &
-                  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 .and. nn_dyn2d_dta(ib_obc) .eq. 1 ) ) THEN
+           &      ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 .and.  &
+           &        ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) ) THEN
 
                jfld = jfld + 1
@@ -583,5 +591,5 @@
 
          IF (nn_dyn2d(ib_obc) .gt. 0) THEN
-            IF( nn_dyn2d_dta(ib_obc) .eq. 0 .or. ln_full_vel_array(ib_obc) .or. nn_tides(ib_obc) .eq. 1 ) THEN
+            IF( nn_dyn2d_dta(ib_obc) .eq. 0 .or. nn_dyn2d_dta(ib_obc) .eq. 2 .or. ln_full_vel_array(ib_obc) ) THEN
                IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
                   ilen0(1:3) = nblen(1:3)
@@ -614,5 +622,6 @@
          ENDIF
          IF ( ( nn_dyn3d(ib_obc) .gt. 0 .and. nn_dyn3d_dta(ib_obc) .eq. 1 ).or. &
-           &  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) ) THEN
+           &  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 .and.   &
+           &    ( nn_dyn2d_dta(ib_obc) .eq. 1 .or. nn_dyn2d_dta(ib_obc) .eq. 3 ) ) ) THEN
             jfld = jfld + 1
             dta_obc(ib_obc)%u3d => bf(jfld)%fnow(:,1,:)
@@ -646,7 +655,7 @@
                   ilen0(1:3) = nblenrim(1:3)
                ENDIF
-               ALLOCATE( dta_obc(ib_obc)%ssh(ilen0(1)) )
-               ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(1)) )
-               ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(1)) )
+               ALLOCATE( dta_obc(ib_obc)%frld(ilen0(1)) )
+               ALLOCATE( dta_obc(ib_obc)%hicif(ilen0(1)) )
+               ALLOCATE( dta_obc(ib_obc)%hsnif(ilen0(1)) )
             ELSE
                jfld = jfld + 1
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90	(revision 2865)
@@ -41,5 +41,5 @@
 CONTAINS
 
-   SUBROUTINE obc_dyn( kt )
+   SUBROUTINE obc_dyn( kt, dyn3d_only )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE obc_dyn  ***
@@ -51,7 +51,9 @@
       USE wrk_nemo, ONLY: wrk_2d_7, wrk_2d_8      ! 2D workspace
       !!
-      INTEGER, INTENT( in ) :: kt     ! Main time step counter
+      INTEGER, INTENT( in )           :: kt               ! Main time step counter
+      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities
       !!
       INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter
+      LOGICAL               :: ll_dyn2d, ll_dyn3d  
       !!
 
@@ -60,9 +62,16 @@
       END IF
 
+      ll_dyn2d = .true.
+      ll_dyn3d = .true.
+
+      IF( PRESENT(dyn3d_only) ) THEN
+         IF( dyn3d_only ) ll_dyn2d = .false.
+      ENDIF
+
       !-------------------------------------------------------
       ! Set pointers
       !-------------------------------------------------------
 
-      pssh => sshn_b
+      pssh => sshn
       phur => hur
       phvr => hvr
@@ -92,7 +101,7 @@
       !-------------------------------------------------------
 
-      CALL obc_dyn2d( kt )
+      IF( ll_dyn2d ) CALL obc_dyn2d( kt )
 
-      CALL obc_dyn3d( kt )
+      IF( ll_dyn3d ) CALL obc_dyn3d( kt )
 
       !-------------------------------------------------------
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90	(revision 2865)
@@ -55,5 +55,5 @@
             CYCLE
          CASE(jp_frs)
-            CALL obc_dyn2d_frs( idx_obc(ib_obc), dta_obc(ib_obc), kt )
+            CALL obc_dyn2d_frs( idx_obc(ib_obc), dta_obc(ib_obc) )
          CASE(jp_flather)
             CALL obc_dyn2d_fla( idx_obc(ib_obc), dta_obc(ib_obc) )
@@ -65,5 +65,5 @@
    END SUBROUTINE obc_dyn2d
 
-   SUBROUTINE obc_dyn2d_frs( idx, dta, kt )
+   SUBROUTINE obc_dyn2d_frs( idx, dta )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE obc_dyn2d_frs  ***
@@ -76,5 +76,4 @@
       !!               topography. Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER,         INTENT(in) ::   kt
       TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
       TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcini.F90	(revision 2865)
@@ -51,5 +51,16 @@
       !! ** Input   :  obc_init.nc, input file for unstructured open boundaries
       !!----------------------------------------------------------------------      
-      INTEGER  ::   ib_obc, ii, ij, ik, igrd, ib, ir       ! dummy loop indices
+      ! namelist variables
+      !-------------------
+      INTEGER, PARAMETER          :: jp_nseg = 100
+      INTEGER                     :: nobcsege, nobcsegw, nobcsegn, nobcsegs 
+      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
+
+      ! local variables
+      !-------------------
+      INTEGER  ::   ib_obc, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
       INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
       INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
@@ -72,7 +83,13 @@
          &             nn_ice_lim2, nn_ice_lim2_dta,                       &
 #endif
-         &             nn_tides, ln_vol, ln_clim, nn_volctl,               &
+         &             ln_vol, nn_volctl,                                  &
          &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             &
          &             nn_dmp3d_in, nn_dmp3d_out
+      !!
+      NAMELIST/namobc_index/ nobcsege, jpieob, jpjedt, jpjeft,             &
+                             nobcsegw, jpiwob, jpjwdt, jpjwft,             &
+                             nobcsegn, jpjnob, jpindt, jpinft,             &
+                             nobcsegs, jpjsob, jpisdt, jpisft
+
       !!----------------------------------------------------------------------
 
@@ -109,6 +126,4 @@
 #endif
       ln_vol            = .false.
-      ln_clim(:)        = .false.
-      nn_tides(:)       =  0  ! default to no tidal forcing
       nn_volctl         = -1  ! uninitialised flag
       nn_rimwidth(:)    = -1  ! uninitialised flag
@@ -137,4 +152,11 @@
         IF(lwp) WRITE(numout,*) ' ' 
         IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_obc,'------' 
+
+        IF( ln_coords_file(ib_obc) ) THEN
+           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_obc))
+        ELSE
+           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
+        ENDIF
+        IF(lwp) WRITE(numout,*)
 
         IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
@@ -149,5 +171,7 @@
               CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for obc data'        
               CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be 0 or 1' )
+              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
+              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
+              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
            END SELECT
         ENDIF
@@ -201,20 +225,6 @@
 #endif
 
-        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth
+        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_obc)
         IF(lwp) WRITE(numout,*)
-
-        SELECT CASE( nn_tides(ib_obc) ) 
-        CASE(0)
-          IF(lwp) WRITE(numout,*) 'No tidal harmonic forcing'
-          IF(lwp) WRITE(numout,*)
-        CASE(1)
-          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ONLY for barotropic solution'
-          IF(lwp) WRITE(numout,*)
-        CASE(2)
-          IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing ADDED to other barotropic boundary conditions'
-          IF(lwp) WRITE(numout,*)
-        CASE DEFAULT
-          CALL ctl_stop( 'obc_ini: ERROR: incorrect value for nn_tides ' )
-        END SELECT
 
       ENDDO
@@ -240,4 +250,5 @@
       ! Work out global dimensions of boundary data
       ! ---------------------------------------------
+      REWIND( numnam )                    
       DO ib_obc = 1, nb_obc
 
@@ -245,7 +256,70 @@
          IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Work out size of global arrays from namelist parameters
  
-
-           !! 1. Read parameters from namelist
-           !! 2. Work out global size of boundary data arrays nblendta and jpbdta
+            ! No REWIND here because may need to read more than one namobc_index namelist.
+            READ  ( numnam, namobc_index )
+
+            ! Automatic boundary definition: if nobcsegX = -1
+            ! set boundary to whole side of model domain.
+            IF( nobcsege == -1 ) THEN
+               nobcsege = 1
+               jpieob(1) = jpiglo - 1
+               jpjedt(1) = 2
+               jpjeft(1) = jpjglo - 1
+            ENDIF
+            IF( nobcsegw == -1 ) THEN
+               nobcsegw = 1
+               jpiwob(1) = 2
+               jpjwdt(1) = 2
+               jpjwft(1) = jpjglo - 1
+            ENDIF
+            IF( nobcsegn == -1 ) THEN
+               nobcsegn = 1
+               jpjnob(1) = jpjglo - 1
+               jpindt(1) = 2
+               jpinft(1) = jpiglo - 1
+            ENDIF
+            IF( nobcsegs == -1 ) THEN
+               nobcsegs = 1
+               jpjsob(1) = 2
+               jpisdt(1) = 2
+               jpisft(1) = jpiglo - 1
+            ENDIF
+
+            nblendta(:,ib_obc) = 0
+            DO iseg = 1, nobcsege
+               igrd = 1
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjeft(iseg) - jpjedt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjeft(iseg) - jpjedt(iseg) + 1               
+               igrd = 3
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjeft(iseg) - jpjedt(iseg)               
+            ENDDO
+            DO iseg = 1, nobcsegw
+               igrd = 1
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjwft(iseg) - jpjwdt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjwft(iseg) - jpjwdt(iseg) + 1               
+               igrd = 3
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpjwft(iseg) - jpjwdt(iseg)               
+            ENDDO
+            DO iseg = 1, nobcsegn
+               igrd = 1
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpinft(iseg) - jpindt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpinft(iseg) - jpindt(iseg)               
+               igrd = 3
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpinft(iseg) - jpindt(iseg) + 1
+            ENDDO
+            DO iseg = 1, nobcsegs
+               igrd = 1
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpisft(iseg) - jpisdt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpisft(iseg) - jpisdt(iseg)
+               igrd = 3
+               nblendta(igrd,ib_obc) = nblendta(igrd,ib_obc) + jpisft(iseg) - jpisdt(iseg) + 1               
+            ENDDO
+
+            nblendta(:,ib_obc) = nblendta(:,ib_obc) * nn_rimwidth(ib_obc)
+            jpbdta = MAXVAL(nblendta(:,ib_obc))               
 
 
@@ -263,5 +337,5 @@
          ENDIF 
 
-      ENDDO
+      ENDDO ! ib_obc
 
       ! Allocate arrays
@@ -274,9 +348,164 @@
       ! Calculate global boundary index arrays or read in from file
       !------------------------------------------------------------
+      REWIND( numnam )                    
       DO ib_obc = 1, nb_obc
 
          IF( .NOT. ln_coords_file(ib_obc) ) THEN ! Calculate global index arrays from namelist parameters
 
-            !! Calculate global index arrays from namelist parameters
+            ! No REWIND here because may need to read more than one namobc_index namelist.
+            READ  ( numnam, namobc_index )
+
+            ! Automatic boundary definition: if nobcsegX = -1
+            ! set boundary to whole side of model domain.
+            IF( nobcsege == -1 ) THEN
+               nobcsege = 1
+               jpieob(1) = jpiglo - 1
+               jpjedt(1) = 2
+               jpjeft(1) = jpjglo - 1
+            ENDIF
+            IF( nobcsegw == -1 ) THEN
+               nobcsegw = 1
+               jpiwob(1) = 2
+               jpjwdt(1) = 2
+               jpjwft(1) = jpjglo - 1
+            ENDIF
+            IF( nobcsegn == -1 ) THEN
+               nobcsegn = 1
+               jpjnob(1) = jpjglo - 1
+               jpindt(1) = 2
+               jpinft(1) = jpiglo - 1
+            ENDIF
+            IF( nobcsegs == -1 ) THEN
+               nobcsegs = 1
+               jpjsob(1) = 2
+               jpisdt(1) = 2
+               jpisft(1) = jpiglo - 1
+            ENDIF
+
+            ! ------------ T points -------------
+            igrd = 1  
+            icount = 0
+            DO ir = 1, nn_rimwidth(ib_obc)
+               ! east
+               DO iseg = 1, nobcsege
+                  DO ij = jpjedt(iseg), jpjeft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = jpieob(iseg) - ir + 1
+                     nbjdta(icount, igrd, ib_obc) = ij
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! west
+               DO iseg = 1, nobcsegw
+                  DO ij = jpjwdt(iseg), jpjwft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = jpiwob(iseg) + ir - 1
+                     nbjdta(icount, igrd, ib_obc) = ij
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! north
+               DO iseg = 1, nobcsegn
+                  DO ii = jpindt(iseg), jpinft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = ii
+                     nbjdta(icount, igrd, ib_obc) = jpjnob(iseg) - ir + 1
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! south
+               DO iseg = 1, nobcsegs
+                  DO ii = jpisdt(iseg), jpisft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = ii
+                     nbjdta(icount, igrd, ib_obc) = jpjsob(iseg) + ir + 1
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            ! ------------ U points -------------
+            igrd = 2  
+            icount = 0
+            DO ir = 1, nn_rimwidth(ib_obc)
+               ! east
+               DO iseg = 1, nobcsege
+                  DO ij = jpjedt(iseg), jpjeft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = jpieob(iseg) - ir
+                     nbjdta(icount, igrd, ib_obc) = ij
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! west
+               DO iseg = 1, nobcsegw
+                  DO ij = jpjwdt(iseg), jpjwft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = jpiwob(iseg) + ir - 1
+                     nbjdta(icount, igrd, ib_obc) = ij
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! north
+               DO iseg = 1, nobcsegn
+                  DO ii = jpindt(iseg), jpinft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = ii
+                     nbjdta(icount, igrd, ib_obc) = jpjnob(iseg) - ir + 1
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! south
+               DO iseg = 1, nobcsegs
+                  DO ii = jpisdt(iseg), jpisft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = ii
+                     nbjdta(icount, igrd, ib_obc) = jpjsob(iseg) + ir + 1
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            ! ------------ V points -------------
+            igrd = 3  
+            icount = 0
+            DO ir = 1, nn_rimwidth(ib_obc)
+               ! east
+               DO iseg = 1, nobcsege
+                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = jpieob(iseg) - ir + 1
+                     nbjdta(icount, igrd, ib_obc) = ij
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! west
+               DO iseg = 1, nobcsegw
+                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = jpiwob(iseg) + ir - 1
+                     nbjdta(icount, igrd, ib_obc) = ij
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! north
+               DO iseg = 1, nobcsegn
+                  DO ii = jpindt(iseg), jpinft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = ii
+                     nbjdta(icount, igrd, ib_obc) = jpjnob(iseg) - ir
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+               ! south
+               DO iseg = 1, nobcsegs
+                  DO ii = jpisdt(iseg), jpisft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_obc) = ii
+                     nbjdta(icount, igrd, ib_obc) = jpjsob(iseg) + ir + 1
+                     nbrdta(icount, igrd, ib_obc) = ir
+                  ENDDO
+               ENDDO
+            ENDDO
 
          ELSE            ! Read global index arrays from boundary coordinates file.
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctides.F90	(revision 2865)
@@ -94,5 +94,5 @@
       REWIND(numnam)
       DO ib_obc = 1, nb_obc
-         IF( nn_tides(ib_obc) .gt. 0 ) THEN
+         IF( nn_dyn2d_dta(ib_obc) .ge. 2 ) THEN
 
             td => tides(ib_obc)
@@ -260,5 +260,5 @@
             ENDIF ! date correction
             !
-         ENDIF ! nn_tides .gt. 0
+         ENDIF ! nn_dyn2d_dta(ib_obc) .ge. 2
          !
       END DO ! loop on ib_obc
@@ -267,5 +267,5 @@
 
 
-   SUBROUTINE tide_update ( kt, jit, idx, dta, td )
+   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE tide_update  ***
@@ -276,21 +276,32 @@
       INTEGER, INTENT( in )          ::   kt      ! Main timestep counter
 !!gm doctor jit ==> kit
-      INTEGER, INTENT( in )          ::   jit     ! Barotropic timestep counter (for timesplitting option)
       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                          :: 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
       !!----------------------------------------------------------------------
 
+      time_add = 0
+      IF( PRESENT(time_offset) ) THEN
+         time_add = time_offset
+      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( jit == 0 ) THEN  
-         z_arg = kt * rdt * rad / 3600.0
-      ELSE                              ! we are in a barotropic subcycle (for timesplitting option)
-!         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0
-         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
+      IF( PRESENT(jit) ) THEN  
+         z_arg = ( (kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
+      ELSE                              
+         z_arg = (kt+time_add) * rdt * rad / 3600.0
       ENDIF
 
@@ -301,5 +312,4 @@
       END DO
 
-      !
       DO itide = 1, td%ncpt
          igrd=1                              ! SSH on tracer grid.
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctra.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctra.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obctra.F90	(revision 2865)
@@ -52,5 +52,5 @@
             CALL obc_tra_frs( idx_obc(ib_obc), dta_obc(ib_obc), kt )
          CASE DEFAULT
-            CALL ctl_stop( 'obc_tra : unrecognised option for open boundaries for T an S' )
+            CALL ctl_stop( 'obc_tra : unrecognised option for open boundaries for T and S' )
          END SELECT
       ENDDO
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90	(revision 2865)
@@ -10,149 +10,160 @@
    !!----------------------------------------------------------------------
 #if   defined key_obc   &&   defined key_dynspg_flt
-!!$   !!----------------------------------------------------------------------
-!!$   !!   'key_obc'            AND      unstructured open boundary conditions
-!!$   !!   'key_dynspg_flt'                              filtered free surface
-!!$   !!----------------------------------------------------------------------
-!!$   USE oce             ! ocean dynamics and tracers 
-!!$   USE dom_oce         ! ocean space and time domain 
-!!$   USE phycst          ! physical constants
-!!$   USE obc_oce         ! ocean open boundary conditions
-!!$   USE lib_mpp         ! for mppsum
-!!$   USE in_out_manager  ! I/O manager
-!!$   USE sbc_oce         ! ocean surface boundary conditions
-!!$
-!!$   IMPLICIT NONE
-!!$   PRIVATE
-!!$
-!!$   PUBLIC obc_vol        ! routine called by dynspg_flt.h90
-!!$
-!!$   !! * Substitutions
-!!$#  include "domzgr_substitute.h90"
-!!$   !!----------------------------------------------------------------------
-!!$   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-!!$   !! $Id$ 
-!!$   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-!!$   !!----------------------------------------------------------------------
-!!$CONTAINS
-!!$
-!!$   SUBROUTINE obc_vol( kt )
-!!$      !!----------------------------------------------------------------------
-!!$      !!                      ***  ROUTINE obcvol  ***
-!!$      !!
-!!$      !! ** Purpose :   This routine is called in dynspg_flt to control 
-!!$      !!      the volume of the system. A correction velocity is calculated
-!!$      !!      to correct the total transport through the unstructured OBC. 
-!!$      !!      The total depth used is constant (H0) to be consistent with the 
-!!$      !!      linear free surface coded in OPA 8.2
-!!$      !!
-!!$      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
-!!$      !!      the total transport through all open boundaries (trans_obc) minus
-!!$      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral 
-!!$      !!      surface (obcsurftot) of the unstructured boundary. 
-!!$      !!         zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot)
-!!$      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
-!!$      !!                       over all the domain in m3/s at each time step.
-!!$      !!      z_cflxemp < 0 when precipitation dominate
-!!$      !!      z_cflxemp > 0 when evaporation dominate
-!!$      !!
-!!$      !!      There are 2 options (user's desiderata): 
-!!$      !!         1/ The volume changes according to E-P, this is the default
-!!$      !!            option. In this case the cumulate E-P flux are setting to
-!!$      !!            zero (z_cflxemp=0) to calculate the correction velocity. So
-!!$      !!            it will only balance the flux through open boundaries.
-!!$      !!            (set nn_volctl to 0 in tne namelist for this option)
-!!$      !!         2/ The volume is constant even with E-P flux. In this case
-!!$      !!            the correction velocity must balance both the flux 
-!!$      !!            through open boundaries and the ones through the free
-!!$      !!            surface. 
-!!$      !!            (set nn_volctl to 1 in tne namelist for this option)
-!!$      !!----------------------------------------------------------------------
-!!$      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-!!$      !!
-!!$      INTEGER  ::   ji, jj, jk, jb, jgrd
-!!$      INTEGER  ::   ii, ij
-!!$      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
-!!$      !!-----------------------------------------------------------------------------
-!!$
-!!$      IF( ln_vol ) THEN
-!!$
-!!$      IF( kt == nit000 ) THEN 
-!!$         IF(lwp) WRITE(numout,*)
-!!$         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC'
-!!$         IF(lwp) WRITE(numout,*)'~~~~~~~'
-!!$      END IF 
-!!$
-!!$      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
-!!$      ! -----------------------------------------------------------------------
-!!$      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
-!!$      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
-!!$
-!!$      ! Transport through the unstructured open boundary
-!!$      ! ------------------------------------------------
-!!$      zubtpecor = 0.e0
-!!$      jgrd = 2                               ! cumulate u component contribution first 
-!!$      DO jb = 1, nblenrim(jgrd)
-!!$         DO jk = 1, jpkm1
-!!$            ii = nbi(jb,jgrd)
-!!$            ij = nbj(jb,jgrd)
-!!$            zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
-!!$         END DO
-!!$      END DO
-!!$      jgrd = 3                               ! then add v component contribution
-!!$      DO jb = 1, nblenrim(jgrd)
-!!$         DO jk = 1, jpkm1
-!!$            ii = nbi(jb,jgrd)
-!!$            ij = nbj(jb,jgrd)
-!!$            zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
-!!$         END DO
-!!$      END DO
-!!$      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
-!!$
-!!$      ! The normal velocity correction
-!!$      ! ------------------------------
-!!$      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot 
-!!$      ELSE                   ;   zubtpecor =   zubtpecor             / obcsurftot
-!!$      END IF
-!!$
-!!$      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
-!!$      ! -------------------------------------------------------------
-!!$      ztranst = 0.e0
-!!$      jgrd = 2                               ! correct u component
-!!$      DO jb = 1, nblenrim(jgrd)
-!!$         DO jk = 1, jpkm1
-!!$            ii = nbi(jb,jgrd)
-!!$            ij = nbj(jb,jgrd)
-!!$            ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk)
-!!$            ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
-!!$         END DO
-!!$      END DO
-!!$      jgrd = 3                              ! correct v component
-!!$      DO jb = 1, nblenrim(jgrd)
-!!$         DO jk = 1, jpkm1
-!!$            ii = nbi(jb,jgrd)
-!!$            ij = nbj(jb,jgrd)
-!!$            va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk)
-!!$            ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
-!!$         END DO
-!!$      END DO
-!!$      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
-!!$ 
-!!$      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
-!!$      ! ------------------------------------------------------
-!!$      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-!!$         IF(lwp) WRITE(numout,*)
-!!$         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
-!!$         IF(lwp) WRITE(numout,*)'~~~~~~~ '
-!!$         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
-!!$         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', obcsurftot, '(m2)'
-!!$         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
-!!$         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
-!!$      END IF 
-!!$      !
-!!$      END IF ! ln_vol
-!!$
-!!$   END SUBROUTINE obc_vol
-!!$
-!!$#else
+   !!----------------------------------------------------------------------
+   !!   'key_obc'            AND      unstructured open boundary conditions
+   !!   'key_dynspg_flt'                              filtered free surface
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain 
+   USE phycst          ! physical constants
+   USE obc_oce         ! ocean open boundary conditions
+   USE lib_mpp         ! for mppsum
+   USE in_out_manager  ! I/O manager
+   USE sbc_oce         ! ocean surface boundary conditions
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC obc_vol        ! routine called by dynspg_flt.h90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE obc_vol( kt )
+      !!----------------------------------------------------------------------
+      !!                      ***  ROUTINE obcvol  ***
+      !!
+      !! ** Purpose :   This routine is called in dynspg_flt to control 
+      !!      the volume of the system. A correction velocity is calculated
+      !!      to correct the total transport through the unstructured OBC. 
+      !!      The total depth used is constant (H0) to be consistent with the 
+      !!      linear free surface coded in OPA 8.2
+      !!
+      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
+      !!      the total transport through all open boundaries (trans_obc) minus
+      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral 
+      !!      surface (obcsurftot) of the unstructured boundary. 
+      !!         zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot)
+      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
+      !!                       over all the domain in m3/s at each time step.
+      !!      z_cflxemp < 0 when precipitation dominate
+      !!      z_cflxemp > 0 when evaporation dominate
+      !!
+      !!      There are 2 options (user's desiderata): 
+      !!         1/ The volume changes according to E-P, this is the default
+      !!            option. In this case the cumulate E-P flux are setting to
+      !!            zero (z_cflxemp=0) to calculate the correction velocity. So
+      !!            it will only balance the flux through open boundaries.
+      !!            (set nn_volctl to 0 in tne namelist for this option)
+      !!         2/ The volume is constant even with E-P flux. In this case
+      !!            the correction velocity must balance both the flux 
+      !!            through open boundaries and the ones through the free
+      !!            surface. 
+      !!            (set nn_volctl to 1 in tne namelist for this option)
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !!
+      INTEGER  ::   ji, jj, jk, jb, jgrd
+      INTEGER  ::   ib_obc, ii, ij
+      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
+      TYPE(OBC_INDEX), POINTER :: idx
+      !!-----------------------------------------------------------------------------
+
+      IF( ln_vol ) THEN
+
+      IF( kt == nit000 ) THEN 
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC'
+         IF(lwp) WRITE(numout,*)'~~~~~~~'
+      END IF 
+
+      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
+      ! -----------------------------------------------------------------------
+      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
+      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
+
+      ! Transport through the unstructured open boundary
+      ! ------------------------------------------------
+      zubtpecor = 0.e0
+      DO ib_obc = 1, nb_obc
+         idx => idx_obc(ib_obc)
+
+         jgrd = 2                               ! cumulate u component contribution first 
+         DO jb = 1, idx%nblenrim(jgrd)
+            DO jk = 1, jpkm1
+               ii = idx%nbi(jb,jgrd)
+               ij = idx%nbj(jb,jgrd)
+               zubtpecor = zubtpecor + idx%flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+            END DO
+         END DO
+         jgrd = 3                               ! then add v component contribution
+         DO jb = 1, idx%nblenrim(jgrd)
+            DO jk = 1, jpkm1
+               ii = idx%nbi(jb,jgrd)
+               ij = idx%nbj(jb,jgrd)
+               zubtpecor = zubtpecor + idx%flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
+            END DO
+         END DO
+
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
+
+      ! The normal velocity correction
+      ! ------------------------------
+      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot 
+      ELSE                   ;   zubtpecor =   zubtpecor             / obcsurftot
+      END IF
+
+      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
+      ! -------------------------------------------------------------
+      ztranst = 0.e0
+      DO ib_obc = 1, nb_obc
+         idx => idx_obc(ib_obc)
+
+         jgrd = 2                               ! correct u component
+         DO jb = 1, idx%nblenrim(jgrd)
+            DO jk = 1, jpkm1
+               ii = idx%nbi(jb,jgrd)
+               ij = idx%nbj(jb,jgrd)
+               ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb) * zubtpecor * umask(ii,ij,jk)
+               ztranst = ztranst + idx%flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+            END DO
+         END DO
+         jgrd = 3                              ! correct v component
+         DO jb = 1, idx%nblenrim(jgrd)
+            DO jk = 1, jpkm1
+               ii = idx%nbi(jb,jgrd)
+               ij = idx%nbj(jb,jgrd)
+               va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb) * zubtpecor * vmask(ii,ij,jk)
+               ztranst = ztranst + idx%flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
+            END DO
+         END DO
+
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
+ 
+      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
+      ! ------------------------------------------------------
+      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
+         IF(lwp) WRITE(numout,*)'~~~~~~~ '
+         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
+         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', obcsurftot, '(m2)'
+         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
+         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
+      END IF 
+      !
+      END IF ! ln_vol
+
+   END SUBROUTINE obc_vol
+
+#else
    !!----------------------------------------------------------------------
    !!   Dummy module                   NO Unstruct Open Boundary Conditions
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2865)
@@ -104,5 +104,5 @@
 CONTAINS
 
-   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, jit, timeshift )
+   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, jit, time_offset )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_read  ***
@@ -121,5 +121,8 @@
       TYPE(MAP_POINTER),INTENT(in), OPTIONAL, DIMENSION(:) ::   map   ! global-to-local mapping index
       INTEGER  , INTENT(in   ), OPTIONAL     ::   jit       ! subcycle timestep for timesplitting option
-      INTEGER  , INTENT(in   ), OPTIONAL     ::   timeshift ! provide fields at time other than "now"
+      INTEGER  , INTENT(in   ), OPTIONAL     ::   time_offset ! provide fields at time other than "now"
+                                                              ! time_offset = -1 => fields at "before" time level
+                                                              ! time_offset = +1 => fields at "after" time levels
+                                                              ! etc.
       !!
       INTEGER  ::   imf        ! size of the structure sd
@@ -128,4 +131,5 @@
       INTEGER  ::   isecend    ! number of second since Jan. 1st 00h of nit000 year at nitend
       INTEGER  ::   isecsbc    ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step
+      INTEGER  ::   time_add   ! local time_offset variable
       LOGICAL  ::   llnxtyr    ! open next year  file?
       LOGICAL  ::   llnxtmth   ! open next month file?
@@ -143,12 +147,15 @@
       ENDIF
 
+      time_add = 0
+      IF( PRESENT(time_offset) ) THEN
+         time_add = time_offset
+      ENDIF
+         
       ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 
       IF( present(jit) ) THEN 
          ! ignore kn_fsbc in this case
-         isecsbc = nsec_year + nsec1jan000 + jit*rdt/REAL(nn_baro,wp)
-      ELSE IF( present(timeshift) ) THEN
-         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + timeshift * rdttra(1)  ! middle of sbc time step
+         isecsbc = nsec_year + nsec1jan000 + (jit+time_add)*rdt/REAL(nn_baro,wp) 
       ELSE
-         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1))  ! middle of sbc time step
+         isecsbc = nsec_year + nsec1jan000 + NINT(0.5 * REAL(kn_fsbc - 1,wp) * rdttra(1)) + time_add * rdttra(1)  ! middle of sbc time step
       ENDIF
       imf = SIZE( sd )
@@ -263,6 +270,7 @@
                   clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f7.2,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', records b/a: ', i4.4, '/', i4.4, ' (days ', f7.2,'/', f7.2, ')')"
-                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &
+                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &            
                      & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
+                  WRITE(numout, *) 'time_add is : ',time_add
                ENDIF
                ! temporal interpolation weights
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 2865)
@@ -35,8 +35,8 @@
    !! free surface                                      !  before  ! now    ! after  !
    !! ------------                                      !  fields  ! fields ! trends !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::            sshf_n          !: sea surface height at f-point [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshb   , sshn   , ssha   !: sea surface height at t-point [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   sshu_b , sshu_n , sshu_a !: sea surface height at u-point [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   sshv_b , sshv_n , sshv_a !: sea surface height at u-point [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::            sshf_n          !: sea surface height at f-point [m]
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   spgu, spgv               !: horizontal surface pressure gradient
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2864)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2865)
@@ -97,5 +97,5 @@
       IF( lk_dtasal  )   CALL dta_sal( kstp )         ! update 3D salinity data
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
-      IF( lk_obc     )   CALL obc_dta( kstp )         ! update dynamic and tracer data at open boundaries
+      IF( lk_obc     )   CALL obc_dta( kstp, time_offset=+1 ) ! update dynamic and tracer data at open boundaries
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
