Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90	(revision 2800)
@@ -150,7 +150,7 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_vvl = .FALSE.   !: fixed grid flag
 #endif
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur  , hvr    !: inverse of u and v-points ocean depth (1/m)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu   , hv     !: depth at u- and v-points (meters)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_0 , hv_0   !: refernce depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur  , hvr    !: inverse of u and v-points ocean depth (1/m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu   , hv     !: depth at u- and v-points (meters)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_0 , hv_0   !: refernce depth at u- and v-points (meters)
 
    INTEGER, PUBLIC ::   nla10              !: deepest    W level Above  ~10m (nlb10 - 1)
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 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 2800)
@@ -30,6 +30,5 @@
    USE domvvl          ! variable volume
    USE obc_oce         ! ocean open boundary conditions
-   USE obcdyn3d        ! open boundary condition for baroclinic velocities
-   USE obcdyn2d        ! open boundary condition for barotropic variables
+   USE obcdyn          ! open boundary condition for baroclinic velocities
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
    USE in_out_manager  ! I/O manager
@@ -155,16 +154,13 @@
       IF( .NOT. lk_dynspg_flt ) THEN
 
-         CALL obc_dyn3d( kt )
-         !
-    !!!! ENDA'S FIX: NEED TO THINK ABOUT THIS !!!!
-         CALL obc_dta( kt+1, jit=0 )
-         CALL obc_dyn2d( kt, sshn_b )
-         !
-!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
-         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
-         !
-         IF( ln_vol_cst )   CALL obc_vol( kt )
-         !
-         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
+         CALL obc_dyn( kt )
+
+!!$!!gm ERROR - potential BUG: sshn should not be modified at this stage !!   ssh_nxt not alrady called
+!!$         CALL lbc_lnk( sshn, 'T', 1. )         ! Boundary conditions on sshn
+!!$         !
+!!$         IF( ln_vol_cst )   CALL obc_vol( kt )
+!!$         !
+!!$         IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh      : ', mask1=tmask )
+
       ENDIF
       !
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 2800)
@@ -33,5 +33,5 @@
    USE solpcg          ! preconditionned conjugate gradient solver
    USE solsor          ! Successive Over-relaxation solver
-   USE obcdyn3d        ! ocean open boundary condition (obc_dyn3d routines)
+   USE obcdyn          ! ocean open boundary condition on dynamics
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
    USE cla             ! cross land advection
@@ -180,6 +180,6 @@
 
 #if defined key_obc
-      CALL obc_dyn3d( kt )    ! Update velocities on each open boundary
-      CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
+      CALL obc_dyn( kt )      ! Update velocities on each open boundary
+      CALL obc_vol( kt )      ! Correction of the barotropic component velocity to control the volume of the system
 #endif
 #if defined key_agrif
Index: /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 2800)
@@ -34,9 +34,9 @@
 
   !                                                                         !!! Time splitting scheme (key_dynspg_ts) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_e, ssha_e   ! sea surface heigth (now, after, average)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ua_e  , va_e     ! barotropic velocities (after)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hu_e  , hv_e     ! now ocean depth ( = Ho+sshn_e )
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hur_e , hvr_e    ! inverse of hu_e and hv_e
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshn_b       ! before field without time-filter
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_e, ssha_e   ! sea surface heigth (now, after, average)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   ua_e  , va_e     ! barotropic velocities (after)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   hu_e  , hv_e     ! now ocean depth ( = Ho+sshn_e )
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   hur_e , hvr_e    ! inverse of hu_e and hv_e
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:), TARGET ::   sshn_b           ! before field without time-filter
 
    !!----------------------------------------------------------------------
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 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 2800)
@@ -466,5 +466,14 @@
          !                                                !  -----------------------
                                                                ! OBC open boundaries
-         IF( lk_obc .OR. ln_tides )   CALL obc_dyn2d( kt, sshn_e ) 
+#if defined key_obc
+         pssh => sshn_e
+         phur => hur_e
+         phvr => hvr_e
+         pu2d => ua_e
+         pv2d => va_e
+
+         IF( lk_obc )   CALL obc_dyn2d( kt ) 
+#endif
+
          !
          CALL lbc_lnk( ua_e  , 'U', -1. )                      ! local domain boundaries 
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 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obc_oce.F90	(revision 2800)
@@ -88,4 +88,10 @@
    REAL(wp)                                    ::   obcsurftot !: Lateral surface of unstructured open boundary
 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !: 
+
    !!----------------------------------------------------------------------
    !! open boundary data variables
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 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 2800)
@@ -22,5 +22,5 @@
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
-   USE obc_oce         ! ocean open boundary conditions
+   USE obc_oce         ! ocean open boundary conditions  
    USE obctides        ! tidal forcing at boundaries
    USE fldread         ! read input fields
@@ -40,8 +40,12 @@
    INTEGER                              ::   nb_obc_fld_sum    ! Total number of fields to update for all boundary sets.
 
+   LOGICAL,           DIMENSION(jp_obc) ::   ln_full_vel_array ! =T => full velocities in 3D boundary conditions
+                                                               ! =F => baroclinic velocities in 3D boundary conditions
+
    TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:), TARGET ::   bf        ! structure of input fields (file informations, fields read)
 
    TYPE(MAP_POINTER), ALLOCATABLE, DIMENSION(:) :: nbmap_ptr   ! array of pointers to nbmap
 
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
@@ -60,21 +64,130 @@
       !!                
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: wrk_2d_9, wrk_2d_10      ! 2D workspace
+      !!
       INTEGER, INTENT( in )           ::   kt    ! ocean time-step index 
       INTEGER, INTENT( in ), OPTIONAL ::   jit   ! subcycle time-step index (for timesplitting option)
       !!
-      INTEGER     ::  ib_obc, jfld, jstart, jend            ! local indices
-      INTEGER, POINTER, DIMENSION(:)  ::   nblen, nblenrim  ! short cuts
+      INTEGER     ::  ib_obc, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices
+      INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
+      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
       !!
       !!---------------------------------------------------------------------------
+
+      IF(wrk_in_use(2, 9,10) ) THEN
+         CALL ctl_stop('obc_dta: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
+      END IF
 
       ! for nn_dtactl = 0, initialise data arrays once for all
       ! from initial conditions
       !-------------------------------------------------------
-      IF( kt .eq. 1 .and. .not. PRESENT(jit) ) THEN
-
+      IF( kt .eq. nit000 .and. .not. PRESENT(jit) ) THEN
+
+         ! Calculate depth-mean currents
+         !-----------------------------
+         pu2d => wrk_2d_9
+         pu2d => wrk_2d_10
+
+         pu2d(:,:) = 0.e0
+         pv2d(:,:) = 0.e0
+
+         DO ik = 1, jpkm1   !! Vertically integrated momentum trends
+             pu2d(:,:) = pu2d(:,:) + fse3u(:,:,ik) * umask(:,:,ik) * un(:,:,ik)
+             pv2d(:,:) = pv2d(:,:) + fse3v(:,:,ik) * vmask(:,:,ik) * vn(:,:,ik)
+         END DO
+         pu2d(:,:) = pu2d(:,:) * hur(:,:)
+         pv2d(:,:) = pv2d(:,:) * hvr(:,:)
+         
          DO ib_obc = 1, nb_obc
             IF( nn_dtactl(ib_obc) .eq. 0 ) THEN
 
-               !!! TO BE DONE !!!
+               nblen => idx_obc(ib_obc)%nblen
+               nblenrim => idx_obc(ib_obc)%nblenrim
+
+               IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN 
+                  IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
+                     ilen1(:) = nblen(:)
+                  ELSE
+                     ilen1(:) = nblenrim(:)
+                  ENDIF
+                  igrd = 1
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                     ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                     dta_obc(ib_obc)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                  END DO 
+                  igrd = 2
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                     ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                     dta_obc(ib_obc)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1)         
+                  END DO 
+                  igrd = 3
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                     ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                     dta_obc(ib_obc)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1)         
+                  END DO 
+               ENDIF
+
+               IF( nn_dyn3d(ib_obc) .gt. 0 ) THEN 
+                  IF( nn_dyn3d(ib_obc) .eq. jp_frs ) THEN
+                     ilen1(:) = nblen(:)
+                  ELSE
+                     ilen1(:) = nblenrim(:)
+                  ENDIF
+                  igrd = 2 
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                        ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                        dta_obc(ib_obc)%u3d(ib,ik) =  ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)         
+                     END DO
+                  END DO 
+                  igrd = 3 
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                        ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                        dta_obc(ib_obc)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik)         
+                     END DO
+                  END DO 
+               ENDIF
+
+               IF( nn_tra(ib_obc) .gt. 0 ) THEN 
+                  IF( nn_tra(ib_obc) .eq. jp_frs ) THEN
+                     ilen1(:) = nblen(:)
+                  ELSE
+                     ilen1(:) = nblenrim(:)
+                  ENDIF
+                  igrd = 1                       ! Everything is at T-points here
+                  DO ib = 1, ilen1(igrd)
+                     DO ik = 1, jpkm1
+                        ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                        ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                        dta_obc(ib_obc)%tem(ib,ik) = tn(ii,ij,ik) * tmask(ii,ij,ik)         
+                        dta_obc(ib_obc)%sal(ib,ik) = sn(ii,ij,ik) * tmask(ii,ij,ik)         
+                     END DO
+                  END DO 
+               ENDIF
+
+#if defined key_lim2
+               IF( nn_ice_lim2(ib_obc) .gt. 0 ) THEN 
+                  IF( nn_ice_lim2(ib_obc) .eq. jp_frs ) THEN
+                     ilen1(:) = nblen(:)
+                  ELSE
+                     ilen1(:) = nblenrim(:)
+                  ENDIF
+                  igrd = 1                       ! Everything is at T-points here
+                  DO ib = 1, ilen1(igrd)
+                     ii = idx_obc(ib_obc)%nbi(ib,igrd)
+                     ij = idx_obc(ib_obc)%nbj(ib,igrd)
+                     dta_obc(ib_obc)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)         
+                     dta_obc(ib_obc)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)         
+                     dta_obc(ib_obc)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)         
+                  END DO 
+               ENDIF
+#endif
 
             ENDIF
@@ -103,6 +216,46 @@
             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_obc) ) THEN 
+
+               igrd = 2                      ! zonal velocity
+               dta_obc(ib_obc)%u2d(:) = 0.0
+               DO ib = 1, idx_obc(ib_obc)%nblen(igrd)
+                  ii   = idx_obc(ib_obc)%nbi(ib,igrd)
+                  ij   = idx_obc(ib_obc)%nbj(ib,igrd)
+                  DO ik = 1, jpkm1
+                     dta_obc(ib_obc)%u2d(ib) = dta_obc(ib_obc)%u2d(ib) &
+              &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_obc(ib_obc)%u3d(ib,ik)
+                  END DO
+                  dta_obc(ib_obc)%u2d(ib) =  dta_obc(ib_obc)%u2d(ib) * hur(ii,ij)
+                  DO ik = 1, jpkm1
+                     dta_obc(ib_obc)%u3d(ib,ik) = dta_obc(ib_obc)%u3d(ib,ik) - dta_obc(ib_obc)%u2d(ib) 
+                  END DO
+               END DO
+
+               igrd = 3                      ! meridional velocity
+               dta_obc(ib_obc)%v2d(:) = 0.0
+               DO ib = 1, idx_obc(ib_obc)%nblen(igrd)
+                  ii   = idx_obc(ib_obc)%nbi(ib,igrd)
+                  ij   = idx_obc(ib_obc)%nbj(ib,igrd)
+                  DO ik = 1, jpkm1
+                     dta_obc(ib_obc)%v2d(ib) = dta_obc(ib_obc)%v2d(ib) &
+              &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_obc(ib_obc)%v3d(ib,ik)
+                  END DO
+                  dta_obc(ib_obc)%v2d(ib) =  dta_obc(ib_obc)%v2d(ib) * hvr(ii,ij)
+                  DO ik = 1, jpkm1
+                     dta_obc(ib_obc)%v3d(ib,ik) = dta_obc(ib_obc)%v3d(ib,ik) - dta_obc(ib_obc)%v2d(ib) 
+                  END DO
+               END DO
+    
+            ENDIF
+
          END IF ! nn_dtactl(ib_obc) = 1
       END DO  ! ib_obc
+
+      IF(wrk_not_released(2, 9,10) )    CALL ctl_stop('obc_dta: ERROR: failed to release workspace arrays.')
 
       END SUBROUTINE obc_dta
@@ -119,9 +272,14 @@
       !!                
       !!----------------------------------------------------------------------
+      USE dynspg_oce, ONLY: lk_dynspg_ts
+      !!
       INTEGER     ::  ib_obc, jfld, jstart, jend, ierror  ! local indices
       !!
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
       CHARACTER(len=100), DIMENSION(nb_obc)  ::   cn_dir_array  ! Root directory for location of data files
+      LOGICAL                                ::   ln_full_vel   ! =T => full velocities in 3D boundary data
+                                                                ! =F => baroclinic velocities in 3D boundary data
       INTEGER                                ::   ilen_global   ! Max length required for global obc dta arrays
+      INTEGER,              DIMENSION(jpbgrd) ::  ilen0         ! size of local arrays
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ilen1, ilen3  ! size of 1st and 3rd dimensions of local arrays
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   iobc           ! obc set for a particular jfld
@@ -138,8 +296,9 @@
       NAMELIST/namobc_dta/ bn_frld, bn_hicif, bn_hsnif
 #endif
+      NAMELIST/namobc_dta/ ln_full_vel
       !!---------------------------------------------------------------------------
 
-      ! Work out how many fields there are to read in and allocate arrays
-      ! -----------------------------------------------------------------
+      ! Work out upper bound of how many fields there are to read in and allocate arrays
+      ! ---------------------------------------------------------------------------
       ALLOCATE( nb_obc_fld(nb_obc) )
       nb_obc_fld(:) = 0
@@ -189,4 +348,5 @@
             ! set file information
             cn_dir = './'        ! directory in which the model is executed
+            ln_full_vel = .false.
             ! ... default values (NB: frequency positive => hours, negative => months)
             !                    !  file       ! frequency !  variable        ! time intep !  clim   ! 'yearly' or ! weights  ! rotation  !
@@ -209,4 +369,10 @@
 
             cn_dir_array(ib_obc) = cn_dir
+            ln_full_vel_array(ib_obc) = ln_full_vel
+
+            IF( ln_full_vel_array(ib_obc) .and. lk_dynspg_ts )  THEN
+               CALL ctl_stop( 'obc_dta_init: ERROR, cannot specify full velocities in boundary data',&
+            &                  'with dynspg_ts option' )   ;   RETURN  
+            ENDIF             
 
             nblen => idx_obc(ib_obc)%nblen
@@ -217,41 +383,44 @@
             IF( nn_dyn2d(ib_obc) .gt. 0 ) THEN 
 
-               jfld = jfld + 1
-               blf_i(jfld) = bn_ssh
-               iobc(jfld) = ib_obc
-               igrid(jfld) = 1
-               IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
-               ilen3(jfld) = 1
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_u2d
-               iobc(jfld) = ib_obc
-               igrid(jfld) = 2
-               IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
-               ilen3(jfld) = 1
-
-               jfld = jfld + 1
-               blf_i(jfld) = bn_v2d
-               iobc(jfld) = ib_obc
-               igrid(jfld) = 3
-               IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(jfld) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(jfld) = nblenrim(igrid(jfld))
-               ENDIF
-               ilen3(jfld) = 1
+               IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_ssh
+                  iobc(jfld) = ib_obc
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+                  ilen3(jfld) = 1
+               ENDIF
+
+               IF( .not. ln_full_vel_array(ib_obc) ) THEN
+
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_u2d
+                  iobc(jfld) = ib_obc
+                  igrid(jfld) = 2
+                  IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
+                     ilen1(jfld) = nblen(igrid(jfld))
+                  ELSE
+                     ilen1(jfld) = nblenrim(igrid(jfld))
+                  ENDIF
+                  ilen3(jfld) = 1
+
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_v2d
+                  iobc(jfld) = ib_obc
+                  igrid(jfld) = 3
+                  IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
+                     ilen1(jfld) = nblen(igrid(jfld))
+                  ELSE
+                     ilen1(jfld) = nblenrim(igrid(jfld))
+                  ENDIF
+                  ilen3(jfld) = 1
+
+               ENDIF
 
             ENDIF
 
             ! baroclinic velocities
-            IF( nn_dyn3d(ib_obc) .gt. 0 ) THEN
+            IF( nn_dyn3d(ib_obc) .gt. 0 .or. &
+                  ( ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) ) THEN
 
                jfld = jfld + 1
@@ -345,10 +514,18 @@
             ENDIF
 #endif
+            ! Recalculate field counts
+            !-------------------------
+            nb_obc_fld_sum = 0
+            IF( ib_obc .eq. 1 ) THEN 
+               nb_obc_fld(ib_obc) = jfld
+               nb_obc_fld_sum     = jfld              
+            ELSE
+               nb_obc_fld(ib_obc) = jfld - nb_obc_fld_sum
+               nb_obc_fld_sum = nb_obc_fld_sum + nb_obc_fld(ib_obc)
+            ENDIF
+
          ENDIF ! nn_dtactl .eq. 1
       ENDDO ! ib_obc
 
-      IF( jfld .ne. nb_obc_fld_sum ) THEN
-         CALL ctl_stop( 'obc_dta: error in initialisation: jpfld .ne. nb_obc_fld_sum' )   ;   RETURN  
-      ENDIF
 
       DO jfld = 1, nb_obc_fld_sum
@@ -385,46 +562,40 @@
             IF (nn_dyn2d(ib_obc) .gt. 0) THEN
                IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(1) = nblen(1)
-                  ilen1(2) = nblen(2)
-                  ilen1(3) = nblen(3)
-               ELSE
-                  ilen1(1) = nblenrim(1)
-                  ilen1(2) = nblenrim(2)
-                  ilen1(3) = nblenrim(3)
-               ENDIF
-               ALLOCATE( dta_obc(ib_obc)%ssh(ilen1(1)) )
-               ALLOCATE( dta_obc(ib_obc)%u2d(ilen1(2)) )
-               ALLOCATE( dta_obc(ib_obc)%v2d(ilen1(3)) )
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  ilen0(1:3) = nblenrim(1:3)
+               ENDIF
+               ALLOCATE( dta_obc(ib_obc)%ssh(ilen0(1)) )
+               ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(2)) )
+               ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(3)) )
             ENDIF
             IF (nn_dyn3d(ib_obc) .gt. 0) THEN
                IF( nn_dyn3d(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(2) = nblen(2)
-                  ilen1(3) = nblen(3)
-               ELSE
-                  ilen1(2) = nblenrim(2)
-                  ilen1(3) = nblenrim(3)
-               ENDIF
-               ALLOCATE( dta_obc(ib_obc)%u3d(ilen1(2),jpk) )
-               ALLOCATE( dta_obc(ib_obc)%v3d(ilen1(3),jpk) )
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  ilen0(1:3) = nblenrim(1:3)
+               ENDIF
+               ALLOCATE( dta_obc(ib_obc)%u3d(ilen0(2),jpk) )
+               ALLOCATE( dta_obc(ib_obc)%v3d(ilen0(3),jpk) )
             ENDIF
             IF (nn_tra(ib_obc) .gt. 0) THEN
                IF( nn_tra(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(1) = nblen(1)
-               ELSE
-                  ilen1(1) = nblenrim(1)
-               ENDIF
-               ALLOCATE( dta_obc(ib_obc)%tem(ilen1(1),jpk) )
-               ALLOCATE( dta_obc(ib_obc)%sal(ilen1(1),jpk) )
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  ilen0(1:3) = nblenrim(1:3)
+               ENDIF
+               ALLOCATE( dta_obc(ib_obc)%tem(ilen0(1),jpk) )
+               ALLOCATE( dta_obc(ib_obc)%sal(ilen0(1),jpk) )
             ENDIF
 #if defined key_lim2
             IF (nn_ice_lim2(ib_obc) .gt. 0) THEN
                IF( nn_ice_lim2(ib_obc) .eq. jp_frs ) THEN
-                  ilen1(1) = nblen(igrid(jfld))
-               ELSE
-                  ilen1(1) = nblenrim(igrid(jfld))
-               ENDIF
-               ALLOCATE( dta_obc(ib_obc)%ssh(ilen1(1)) )
-               ALLOCATE( dta_obc(ib_obc)%u2d(ilen1(1)) )
-               ALLOCATE( dta_obc(ib_obc)%v2d(ilen1(1)) )
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  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)) )
             ENDIF
 #endif
@@ -436,12 +607,29 @@
             !-----------------------------------------------------------
             IF (nn_dyn2d(ib_obc) .gt. 0) THEN
-               jfld = jfld + 1
-               dta_obc(ib_obc)%ssh => bf(jfld)%fnow(:,1,1)
-               jfld = jfld + 1
-               dta_obc(ib_obc)%u2d => bf(jfld)%fnow(:,1,1)
-               jfld = jfld + 1
-               dta_obc(ib_obc)%v2d => bf(jfld)%fnow(:,1,1)
-            ENDIF
-            IF (nn_dyn3d(ib_obc) .gt. 0) THEN
+               IF( nn_dyn2d(ib_obc) .ne. jp_frs ) THEN
+                  jfld = jfld + 1
+                  dta_obc(ib_obc)%ssh => bf(jfld)%fnow(:,1,1)
+               ENDIF
+               IF( ln_full_vel_array(ib_obc) ) THEN
+                  ! In this case we need space but we aren't reading it 
+                  ! directly from the external file. 
+                  IF( nn_dyn2d(ib_obc) .eq. jp_frs ) THEN
+                     ilen0(2) = nblen(2)
+                     ilen0(3) = nblen(3)
+                  ELSE
+                     ilen0(2) = nblenrim(2)
+                     ilen0(3) = nblenrim(3)
+                  ENDIF
+                  ALLOCATE( dta_obc(ib_obc)%u2d(ilen0(2)) )
+                  ALLOCATE( dta_obc(ib_obc)%v2d(ilen0(3)) )
+               ELSE
+                  jfld = jfld + 1
+                  dta_obc(ib_obc)%u2d => bf(jfld)%fnow(:,1,1)
+                  jfld = jfld + 1
+                  dta_obc(ib_obc)%v2d => bf(jfld)%fnow(:,1,1)
+               ENDIF
+            ENDIF
+            IF (nn_dyn3d(ib_obc) .gt. 0 .or. &
+              &   ln_full_vel_array(ib_obc) .and. nn_dyn2d(ib_obc) .gt. 0 ) THEN
                jfld = jfld + 1
                dta_obc(ib_obc)%u3d => bf(jfld)%fnow(:,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 2800)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn.F90	(revision 2800)
@@ -0,0 +1,122 @@
+MODULE obcdyn
+   !!======================================================================
+   !!                       ***  MODULE  obcdyn  ***
+   !! Unstructured Open Boundary Cond. :   Flow relaxation scheme on velocities
+   !!======================================================================
+   !! History :  1.0  !  2005-02  (J. Chanut, A. Sellar)  Original code
+   !!             -   !  2007-07  (D. Storkey) Move Flather implementation to separate routine.
+   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
+   !!            3.2  !  2008-04  (R. Benshila) consider velocity instead of transport 
+   !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations 
+   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
+   !!----------------------------------------------------------------------
+#if defined key_obc 
+   !!----------------------------------------------------------------------
+   !!   'key_obc' :                    Unstructured Open Boundary Condition
+   !!----------------------------------------------------------------------
+   !!   obc_dyn3d        : apply open boundary conditions to baroclinic velocities
+   !!   obc_dyn3d_frs    : apply Flow Relaxation Scheme
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE dynspg_oce      
+   USE obc_oce         ! ocean open boundary conditions
+   USE obcdyn2d        ! open boundary conditions for barotropic solution
+   USE obcdyn3d        ! open boundary conditions for baroclinic velocities
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  !
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   obc_dyn     ! routine called in dynspg_flt (if lk_dynspg_flt) or 
+                        ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp)
+
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: obcdyn.F90 2528 2010-12-27 17:33:53Z rblod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE obc_dyn( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE obc_dyn  ***
+      !!
+      !! ** Purpose : - Wrapper routine for obc_dyn2d and obc_dyn3d.
+      !!
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: wrk_2d_7, wrk_2d_8      ! 2D workspace
+      !!
+      INTEGER, INTENT( in ) :: kt     ! Main time step counter
+      !!
+      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter
+      !!
+
+      IF(wrk_in_use(2, 7,8) ) THEN
+         CALL ctl_stop('obc_dyn: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
+      END IF
+
+      !-------------------------------------------------------
+      ! Set pointers
+      !-------------------------------------------------------
+
+      pssh => sshn_b
+      phur => hur
+      phvr => hvr
+      pu2d => wrk_2d_7
+      pv2d => wrk_2d_8
+
+      !-------------------------------------------------------
+      ! Split velocities into barotropic and baroclinic parts
+      !-------------------------------------------------------
+
+      pu2d(:,:) = 0.e0
+      pv2d(:,:) = 0.e0
+      DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+      END DO
+      pu2d(:,:) = pu2d(:,:) * phur(:,:)
+      pv2d(:,:) = pv2d(:,:) * phvr(:,:)
+      DO jk = 1 , jpkm1
+         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:)
+         va(:,:,jk) = va(:,:,jk) - pv2d(:,:)
+      END DO
+
+      !-------------------------------------------------------
+      ! Apply boundary conditions to barotropic and baroclinic
+      ! parts separately
+      !-------------------------------------------------------
+
+      CALL obc_dyn2d( kt )
+
+      CALL obc_dyn3d( kt )
+
+      !-------------------------------------------------------
+      ! Recombine velocities
+      !-------------------------------------------------------
+
+      DO jk = 1 , jpkm1
+         ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk)
+         va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk)
+      END DO
+
+      IF(wrk_not_released(2, 7,8) )    CALL ctl_stop('obc_dyn: ERROR: failed to release workspace arrays.')
+
+   END SUBROUTINE obc_dyn
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                   NO Unstruct Open Boundary Conditions
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE obc_dyn( kt )      ! Empty routine
+      WRITE(*,*) 'obc_dyn: You should not have seen this print! error?', kt
+   END SUBROUTINE obc_dyn
+#endif
+
+   !!======================================================================
+END MODULE obcdyn
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 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcdyn2d.F90	(revision 2800)
@@ -39,5 +39,5 @@
 CONTAINS
 
-   SUBROUTINE obc_dyn2d( kt, pssh )
+   SUBROUTINE obc_dyn2d( kt )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE obc_dyn2d  ***
@@ -47,5 +47,4 @@
       !!----------------------------------------------------------------------
       INTEGER,                      INTENT(in) ::   kt   ! Main time step counter
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh
       !!
       INTEGER                                  ::   ib_obc ! Loop counter
@@ -56,6 +55,8 @@
          CASE(jp_none)
             CYCLE
+         CASE(jp_frs)
+            CALL obc_dyn2d_frs( idx_obc(ib_obc), dta_obc(ib_obc), kt )
          CASE(jp_flather)
-            CALL obc_dyn2d_fla( idx_obc(ib_obc), dta_obc(ib_obc), pssh )
+            CALL obc_dyn2d_fla( idx_obc(ib_obc), dta_obc(ib_obc) )
          CASE DEFAULT
             CALL ctl_stop( 'obc_dyn3d : unrecognised option for open boundaries for barotropic variables' )
@@ -65,7 +66,50 @@
    END SUBROUTINE obc_dyn2d
 
-   SUBROUTINE obc_dyn2d_fla( idx, dta, pssh )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE obc_dyn_fla  ***
+   SUBROUTINE obc_dyn2d_frs( idx, dta, kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE obc_dyn2d_frs  ***
+      !!
+      !! ** Purpose : - Apply the Flow Relaxation Scheme for barotropic velocities
+      !!                at open boundaries.
+      !!
+      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in 
+      !!               a three-dimensional baroclinic ocean model with realistic
+      !!               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
+      !!
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      !!----------------------------------------------------------------------
+      !
+      !
+      igrd = 2                      ! Relaxation of zonal velocity
+      DO jb = 1, idx%nblen(igrd)
+         ii   = idx%nbi(jb,igrd)
+         ij   = idx%nbj(jb,igrd)
+         zwgt = idx%nbw(jb,igrd)
+         pu2d(ii,ij) = ( pu2d(ii,ij) + zwgt * ( dta%u2d(jb) - pu2d(ii,ij) ) ) * umask(ii,ij,1)
+      END DO
+      !
+      igrd = 3                      ! Relaxation of meridional velocity
+      DO jb = 1, idx%nblen(igrd)
+         ii   = idx%nbi(jb,igrd)
+         ij   = idx%nbj(jb,igrd)
+         zwgt = idx%nbw(jb,igrd)
+         pv2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1)
+      END DO 
+      CALL lbc_lnk( pu2d, 'U', -1. ) 
+      CALL lbc_lnk( pv2d, 'V', -1. )   ! Boundary points should be updated
+      !
+
+   END SUBROUTINE obc_dyn2d_frs
+
+
+   SUBROUTINE obc_dyn2d_fla( idx, dta )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE obc_dyn2d_fla  ***
       !!             
       !!              - Apply Flather boundary conditions on normal barotropic velocities 
@@ -100,4 +144,6 @@
      
 !!!!!!!!!!!! SOME WORK TO DO ON THE TIDES HERE !!!!!!!!!!!!!
+
+!!! REPLACE spgu with nemo_wrk work space
 
       ! Fill temporary array with ssh data (here spgu):
@@ -120,8 +166,8 @@
          iip1 = ii - MIN( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice outside the boundary 
          !
-         zcorr = - idx%flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
+         zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
 !!$         zforc = dta%u2d(jb) + utide(jb)
          zforc = dta%u2d(jb)
-         ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
+         pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
       END DO
       !
@@ -134,11 +180,11 @@
          ijp1 = ij - MIN( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice outside the boundary 
          !
-         zcorr = - idx%flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
+         zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
 !!$         zforc = dta%v2d(jb) + vtide(jb)
          zforc = dta%v2d(jb)
-         va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
-      END DO
-      CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated
-      CALL lbc_lnk( va_e, 'V', -1. )   !
+         pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
+      END DO
+      CALL lbc_lnk( pu2d, 'U', -1. )   ! Boundary points should be updated
+      CALL lbc_lnk( pv2d, 'V', -1. )   !
       !
       !
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 2799)
+++ /branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2800)
@@ -134,5 +134,8 @@
       !!---------------------------------------------------------------------
       ! Note that shifting time to be centrered in the middle of sbc time step impacts only nsec_* variables of the calendar 
-      IF( present(timeshift) ) THEN
+      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
       ELSE
