Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 3116)
@@ -7,4 +7,5 @@
    !!            3.0  !  2008-04  (NEMO team)  add in the reference version     
    !!            3.3  !  2010-09  (D. Storkey) add ice boundary conditions
+   !!            3.4  !  2011     (D. Storkey, J. Chanut) OBC-BDY merge
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -19,30 +20,69 @@
    PUBLIC
 
+   TYPE, PUBLIC ::   OBC_INDEX    !: Indices and weights which define the open boundary
+      INTEGER,          DIMENSION(jpbgrd) ::  nblen
+      INTEGER,          DIMENSION(jpbgrd) ::  nblenrim
+      INTEGER, POINTER, DIMENSION(:,:)   ::  nbi
+      INTEGER, POINTER, DIMENSION(:,:)   ::  nbj
+      INTEGER, POINTER, DIMENSION(:,:)   ::  nbr
+      INTEGER, POINTER, DIMENSION(:,:)   ::  nbmap
+      REAL   , POINTER, DIMENSION(:,:)   ::  nbw
+      REAL   , POINTER, DIMENSION(:)     ::  flagu
+      REAL   , POINTER, DIMENSION(:)     ::  flagv
+   END TYPE OBC_INDEX
+
+   TYPE, PUBLIC ::   OBC_DATA     !: Storage for external data
+      REAL, POINTER, DIMENSION(:)     ::  ssh
+      REAL, POINTER, DIMENSION(:)     ::  u2d
+      REAL, POINTER, DIMENSION(:)     ::  v2d
+      REAL, POINTER, DIMENSION(:,:)   ::  u3d
+      REAL, POINTER, DIMENSION(:,:)   ::  v3d
+      REAL, POINTER, DIMENSION(:,:)   ::  tem
+      REAL, POINTER, DIMENSION(:,:)   ::  sal
+#if defined key_lim2
+      REAL, POINTER, DIMENSION(:)     ::  frld
+      REAL, POINTER, DIMENSION(:)     ::  hicif
+      REAL, POINTER, DIMENSION(:)     ::  hsnif
+#endif
+   END TYPE OBC_DATA
+
    !!----------------------------------------------------------------------
    !! Namelist variables
    !!----------------------------------------------------------------------
-   CHARACTER(len=80) ::   cn_mask        !: Name of unstruct. bdy mask file
-   CHARACTER(len=80) ::   cn_dta_frs_T   !: Name of unstruct. bdy data file at T points for FRS conditions
-   CHARACTER(len=80) ::   cn_dta_frs_U   !: Name of unstruct. bdy data file at U points for FRS conditions
-   CHARACTER(len=80) ::   cn_dta_frs_V   !: Name of unstruct. bdy data file at V points for FRS conditions
-   CHARACTER(len=80) ::   cn_dta_fla_T   !: Name of unstruct. bdy data file at T points for Flather scheme
-   CHARACTER(len=80) ::   cn_dta_fla_U   !: Name of unstruct. bdy data file at U points for Flather scheme
-   CHARACTER(len=80) ::   cn_dta_fla_V   !: Name of unstruct. bdy data file at V points for Flather scheme
+   CHARACTER(len=80), DIMENSION(jp_bdy) ::   cn_coords_file !: Name of bdy coordinates file
+   CHARACTER(len=80)                    ::   cn_mask_file   !: Name of bdy mask file
    !
-   LOGICAL ::   ln_tides = .false.    !: =T apply tidal harmonic forcing along open boundaries
-   LOGICAL ::   ln_vol  = .false.     !: =T volume correction             
-   LOGICAL ::   ln_mask = .false.     !: =T read bdymask from file
-   LOGICAL ::   ln_clim = .false.     !: =T bdy data files contain  1 time dump  (-->bdy forcing will be constant) 
-   !                                  !                         or 12 months     (-->bdy forcing will be cyclic) 
-   LOGICAL ::   ln_dyn_fla  = .false. !: =T Flather boundary conditions on barotropic velocities
-   LOGICAL ::   ln_dyn_frs  = .false. !: =T FRS boundary conditions on velocities
-   LOGICAL ::   ln_tra_frs  = .false. !: =T FRS boundary conditions on tracers (T and S)
-   LOGICAL ::   ln_ice_frs  = .false. !: =T FRS boundary conditions on seaice (leads fraction, ice depth, snow depth)
+   LOGICAL, DIMENSION(jp_bdy) ::   ln_coords_file           !: =T read bdy coordinates from file; 
+   !                                                        !: =F read bdy coordinates from namelist
+   LOGICAL                    ::   ln_mask_file             !: =T read bdymask from file
+   LOGICAL                    ::   ln_vol                   !: =T volume correction             
    !
-   INTEGER ::   nn_rimwidth = 7       !: boundary rim width
-   INTEGER ::   nn_dtactl   = 1	     !: = 0 use the initial state as bdy dta ; = 1 read it in a NetCDF file
-   INTEGER ::   nn_volctl   = 1       !: = 0 the total volume will have the variability of the surface Flux E-P 
-   !                                  !  = 1 the volume will be constant during all the integration.
+   INTEGER                    ::   nb_bdy                   !: number of open boundary sets
+   INTEGER, DIMENSION(jp_bdy) ::   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.
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d                 ! Choice of boundary condition for barotropic variables (U,V,SSH)
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn2d_dta      	    !: = 0 use the initial state as bdy 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_bdy) ::   nn_dyn3d                 ! Choice of boundary condition for baroclinic velocities 
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dyn3d_dta      	    !: = 0 use the initial state as bdy dta ; 
+                                                            !: = 1 read it in a NetCDF file
+   INTEGER, DIMENSION(jp_bdy) ::   nn_tra                   ! Choice of boundary condition for active tracers (T and S)
+   INTEGER, DIMENSION(jp_bdy) ::   nn_tra_dta      	    !: = 0 use the initial state as bdy dta ; 
+                                                            !: = 1 read it in a NetCDF file
+#if defined key_lim2
+   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2              ! Choice of boundary condition for sea ice variables 
+   INTEGER, DIMENSION(jp_bdy) ::   nn_ice_lim2_dta          !: = 0 use the initial state as bdy dta ; 
+                                                            !: = 1 read it in a NetCDF file
+#endif
+   !
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp2d_in              ! Damping timescale (days) for 2D solution for inward radiation or FRS 
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp2d_out             ! Damping timescale (days) for 2D solution for outward radiation 
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp3d_in              ! Damping timescale (days) for 3D solution for inward radiation or FRS 
+   INTEGER, DIMENSION(jp_bdy) ::   nn_dmp3d_out             ! Damping timescale (days) for 3D solution for outward radiation
 
+   
    !!----------------------------------------------------------------------
    !! Global variables
@@ -52,31 +92,21 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bdyvmask   !: Mask defining computational domain at V-points
 
+   REAL(wp)                                    ::   bdysurftot !: Lateral surface of unstructured open boundary
+
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pssh       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   phur       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   phvr       !: Pointers for barotropic fields 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pu2d       !: 
+   REAL(wp), POINTER, DIMENSION(:,:)           ::   pv2d       !: 
+
    !!----------------------------------------------------------------------
-   !! Unstructured open boundary data variables
+   !! open boundary data variables
    !!----------------------------------------------------------------------
-   INTEGER, DIMENSION(jpbgrd) ::   nblen    = 0           !: Size of bdy data on a proc for each grid type
-   INTEGER, DIMENSION(jpbgrd) ::   nblenrim = 0           !: Size of bdy data on a proc for first rim ind
-   INTEGER, DIMENSION(jpbgrd) ::   nblendta = 0           !: Size of bdy data in file
 
-   INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbi, nbj        !: i and j indices of bdy dta
-   INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbr             !: Discrete distance from rim points
-   INTEGER, DIMENSION(jpbdim,jpbgrd) ::   nbmap           !: Indices of data in file for data in memory 
-    
-   REAL(wp) ::   bdysurftot	                            !: Lateral surface of unstructured open boundary
-
-   REAL(wp), DIMENSION(jpbdim)        ::   flagu, flagv   !: Flag for normal velocity compnt for velocity components
-   REAL(wp), DIMENSION(jpbdim,jpbgrd) ::   nbw            !: Rim weights of bdy data
-
-   REAL(wp), DIMENSION(jpbdim)     ::   sshbdy            !: Now clim of bdy sea surface height (Flather)
-   REAL(wp), DIMENSION(jpbdim)     ::   ubtbdy, vbtbdy    !: Now clim of bdy barotropic velocity components
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tbdy  , sbdy      !: Now clim of bdy temperature and salinity  
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ubdy  , vbdy	  !: Now clim of bdy velocity components
-   REAL(wp), DIMENSION(jpbdim) ::   sshtide               !: Tidal boundary array : SSH
-   REAL(wp), DIMENSION(jpbdim) ::   utide, vtide          !: Tidal boundary array : U and V
-#if defined key_lim2
-   REAL(wp), DIMENSION(jpbdim) ::   frld_bdy    !: now ice leads fraction climatology   
-   REAL(wp), DIMENSION(jpbdim) ::   hicif_bdy   !: Now ice  thickness climatology
-   REAL(wp), DIMENSION(jpbdim) ::   hsnif_bdy   !: now snow thickness
-#endif
+   INTEGER,  DIMENSION(jp_bdy)                     ::   nn_dta            !: =0 => *all* data is set to initial conditions
+                                                                          !: =1 => some data to be read in from data files
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays
+   TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process)
+   TYPE(OBC_DATA) , DIMENSION(jp_bdy)              ::   dta_bdy           !: bdy external data (local process)
 
    !!----------------------------------------------------------------------
@@ -94,7 +124,6 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( bdytmask(jpi,jpj) , tbdy(jpbdim,jpk) , sbdy(jpbdim,jpk) ,     &
-         &      bdyumask(jpi,jpj) , ubdy(jpbdim,jpk) ,                        &
-         &      bdyvmask(jpi,jpj) , vbdy(jpbdim,jpk) ,                    STAT=bdy_oce_alloc )
+      ALLOCATE( bdytmask(jpi,jpj) , bdyumask(jpi,jpj), bdyvmask(jpi,jpj),                    &  
+         &      STAT=bdy_oce_alloc )
          !
       IF( lk_mpp             )   CALL mpp_sum ( bdy_oce_alloc )
@@ -112,2 +141,3 @@
    !!======================================================================
 END MODULE bdy_oce
+
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 3116)
@@ -17,8 +17,12 @@
 
    LOGICAL, PUBLIC, PARAMETER ::   lk_bdy  = .TRUE.   !: Unstructured Ocean Boundary Condition flag
-   INTEGER, PUBLIC, PARAMETER ::   jpbdta  = 20000    !: Max length of bdy field in file
-   INTEGER, PUBLIC, PARAMETER ::   jpbdim  = 20000    !: Max length of bdy field on a processor
+   INTEGER, PUBLIC, PARAMETER ::   jp_bdy  = 10       !: Maximum number of bdy sets
    INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000     !: Max number of time dumps per file
-   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 6	      !: Number of horizontal grid types used  (T, u, v, f)
+   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 3	      !: Number of horizontal grid types used  (T, U, V)
+
+   !! Flags for choice of schemes
+   INTEGER, PUBLIC, PARAMETER ::   jp_none         = 0        !: Flag for no open boundary condition
+   INTEGER, PUBLIC, PARAMETER ::   jp_frs          = 1        !: Flag for Flow Relaxation Scheme
+   INTEGER, PUBLIC, PARAMETER ::   jp_flather      = 2        !: Flag for Flather
 #else
    !!----------------------------------------------------------------------
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3116)
@@ -10,19 +10,20 @@
    !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations 
    !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
+   !!            3.4  ????????????????
    !!----------------------------------------------------------------------
 #if defined key_bdy
    !!----------------------------------------------------------------------
-   !!   'key_bdy'                     Unstructured Open Boundary Conditions
-   !!----------------------------------------------------------------------
-   !!   bdy_dta_frs    : read u, v, t, s data along open boundaries
-   !!   bdy_dta_fla : read depth-mean velocities and elevation along open boundaries        
+   !!   'key_bdy'                     Open Boundary Conditions
+   !!----------------------------------------------------------------------
+   !!    bdy_dta        : read external data along open boundaries from file
+   !!    bdy_dta_init   : initialise arrays etc for reading of external data
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
    USE phycst          ! physical constants
-   USE bdy_oce         ! ocean open boundary conditions
+   USE bdy_oce         ! ocean open boundary conditions  
    USE bdytides        ! tidal forcing at boundaries
-   USE iom
-   USE ioipsl
+   USE fldread         ! read input fields
+   USE iom             ! IOM library
    USE in_out_manager  ! I/O logical units
 #if defined key_lim2
@@ -33,30 +34,18 @@
    PRIVATE
 
-   PUBLIC   bdy_dta_frs      ! routines called by step.F90
-   PUBLIC   bdy_dta_fla 
-   PUBLIC   bdy_dta_alloc    ! routine called by bdy_init.F90
-
-   INTEGER ::   numbdyt, numbdyu, numbdyv                      ! logical units for T-, U-, & V-points data file, resp.
-   INTEGER ::   ntimes_bdy                                     ! exact number of time dumps in data files
-   INTEGER ::   nbdy_b, nbdy_a                                 ! record of bdy data file for before and after time step
-   INTEGER ::   numbdyt_bt, numbdyu_bt, numbdyv_bt             ! logical unit for T-, U- & V-points data file, resp.
-   INTEGER ::   ntimes_bdy_bt                                  ! exact number of time dumps in data files
-   INTEGER ::   nbdy_b_bt, nbdy_a_bt                           ! record of bdy data file for before and after time step
-
-   INTEGER, DIMENSION (jpbtime) ::   istep, istep_bt           ! time array in seconds in each data file
-
-   REAL(wp) ::  zoffset                                        ! time offset between time origin in file & start time of model run
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tbdydta, sbdydta   ! time interpolated values of T and S bdy data   
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ubdydta, vbdydta   ! time interpolated values of U and V bdy data 
-   REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta ! Arrays used for time interpolation of bdy data   
-   REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            ! bdy data of ssh
-
-#if defined key_lim2
-   REAL(wp), DIMENSION(jpbdim,2)     ::   frld_bdydta          ! }
-   REAL(wp), DIMENSION(jpbdim,2)     ::   hicif_bdydta         ! } Arrays used for time interp. of ice bdy data 
-   REAL(wp), DIMENSION(jpbdim,2)     ::   hsnif_bdydta         ! }
-#endif
-
+   PUBLIC   bdy_dta          ! routine called by step.F90 and dynspg_ts.F90
+   PUBLIC   bdy_dta_init     ! routine called by nemogcm.F90
+
+   INTEGER, ALLOCATABLE, DIMENSION(:)   ::   nb_bdy_fld        ! Number of fields to update for each boundary set.
+   INTEGER                              ::   nb_bdy_fld_sum    ! Total number of fields to update for all boundary sets.
+
+   LOGICAL,           DIMENSION(jp_bdy) ::   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)
@@ -66,1165 +55,636 @@
 CONTAINS
 
-  FUNCTION bdy_dta_alloc()
-     !!----------------------------------------------------------------------
-     USE lib_mpp, ONLY: ctl_warn, mpp_sum
-     !
-     INTEGER :: bdy_dta_alloc
-     !!----------------------------------------------------------------------
-     !
-     ALLOCATE(tbdydta(jpbdim,jpk,2), sbdydta(jpbdim,jpk,2), &
-              ubdydta(jpbdim,jpk,2), vbdydta(jpbdim,jpk,2), Stat=bdy_dta_alloc)
-
-     IF( lk_mpp           ) CALL mpp_sum ( bdy_dta_alloc )
-     IF(bdy_dta_alloc /= 0) CALL ctl_warn('bdy_dta_alloc: failed to allocate arrays')
-
-   END FUNCTION bdy_dta_alloc
-
-
-   SUBROUTINE bdy_dta_frs( kt )
+      SUBROUTINE bdy_dta( kt, jit, time_offset )
       !!----------------------------------------------------------------------
-      !!                   ***  SUBROUTINE bdy_dta_frs  ***
+      !!                   ***  SUBROUTINE bdy_dta  ***
       !!                    
-      !! ** Purpose :   Read unstructured boundary data for FRS condition.
+      !! ** Purpose :   Update external data for open boundary conditions
       !!
-      !! ** Method  :   At the first timestep, read in boundary data for two
-      !!                times from the file and time-interpolate. At other 
-      !!                timesteps, check to see if we need another time from 
-      !!                the file. If so read it in. Time interpolate.
+      !! ** Method  :   Use fldread.F90
+      !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index (for timesplitting option, otherwise zero)
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: wrk_2d_22, wrk_2d_23   ! 2D workspace
       !!
-      CHARACTER(LEN=80), DIMENSION(3) ::   clfile               ! names of input files
-      CHARACTER(LEN=70 )              ::   clunits              ! units attribute of time coordinate
-      LOGICAL ::   lect                                         ! flag for reading
-      INTEGER ::   it, ib, ik, igrd                             ! dummy loop indices
-      INTEGER ::   igrd_start, igrd_end                         ! start and end of loops on igrd
-      INTEGER ::   idvar                                        ! netcdf var ID
-      INTEGER ::   iman, i15, imois                             ! Time variables for monthly clim forcing
-      INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv
-      INTEGER ::   itimer, totime
-      INTEGER ::   ii, ij                                       ! array addresses
-      INTEGER ::   ipi, ipj, ipk, inum                          ! local integers (NetCDF read)
-      INTEGER ::   iyear0, imonth0, iday0
-      INTEGER ::   ihours0, iminutes0, isec0
-      INTEGER ::   iyear, imonth, iday, isecs
-      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files
-      REAL(wp) ::   dayfrac, zxy, zoffsett
-      REAL(wp) ::   zoffsetu, zoffsetv
-      REAL(wp) ::   dayjul0, zdayjulini
-      REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files
-      REAL(wp), DIMENSION(jpbdta,1,jpk) ::   zdta               ! temporary array for data fields
+      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_bdy, jfld, jstart, jend, ib, ii, ij, ik, igrd  ! local indices
+      INTEGER,          DIMENSION(jpbgrd) ::   ilen1 
+      INTEGER, POINTER, DIMENSION(:)      ::   nblen, nblenrim  ! short cuts
+      !!
       !!---------------------------------------------------------------------------
 
-
-      IF( ln_dyn_frs .OR. ln_tra_frs    &
-         &               .OR. ln_ice_frs ) THEN  ! If these are both false then this routine does nothing
-
-      ! -------------------- !
-      !    Initialization    !
-      ! -------------------- !
-
-      lect   = .false.           ! If true, read a time record
-
-      ! Some time variables for monthly climatological forcing:
-      ! *******************************************************
-
-!!gm  here  use directely daymod calendar variables
- 
-      iman = INT( raamo )      ! Number of months in a year
-
-      i15 = INT( 2*REAL( nday, wp ) / ( REAL( nmonth_len(nmonth), wp ) + 0.5 ) )
-      ! i15=0 if the current day is in the first half of the month, else i15=1
-
-      imois = nmonth + i15 - 1            ! imois is the first month record
-      IF( imois == 0 )   imois = iman
-
-      ! Time variable for non-climatological forcing:
-      ! *********************************************
-      itimer = (kt-nit000+1)*rdt      ! current time in seconds for interpolation 
-
-
-      !                                                !-------------------!
-      IF( kt == nit000 ) THEN                          !  First call only  !
-         !                                             !-------------------!
-         istep(:) = 0
-         nbdy_b   = 0
-         nbdy_a   = 0
-
-         ! Get time information from bdy data file
-         ! ***************************************
-
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*)    'bdy_dta_frs : Initialize unstructured boundary data'
-         IF(lwp) WRITE(numout,*)    '~~~~~~~' 
-
-         IF     ( nn_dtactl == 0 ) THEN
-            !
-            IF(lwp) WRITE(numout,*) '          Bdy data are taken from initial conditions'
-            !
-         ELSEIF (nn_dtactl == 1) THEN
-            !
-            IF(lwp) WRITE(numout,*) '          Bdy data are read in netcdf files'
-            !
-            dayfrac = adatrj  - REAL( itimer, wp ) / 86400.   ! day fraction at time step kt-1
-            dayfrac = dayfrac - INT ( dayfrac )               !
-            totime  = ( nitend - nit000 + 1 ) * rdt           ! Total time of the run to verify that all the
-            !                                                 ! necessary time dumps in file are included
-            !
-            clfile(1) = cn_dta_frs_T
-            clfile(2) = cn_dta_frs_U
-            clfile(3) = cn_dta_frs_V
-            !                                                  
-            ! how many files are we to read in?
-            igrd_start = 1
-            igrd_end   = 3
-            IF(.NOT. ln_tra_frs .AND. .NOT. ln_ice_frs) THEN       ! No T-grid file.
-               igrd_start = 2
-            ELSEIF ( .NOT. ln_dyn_frs ) THEN                           ! No U-grid or V-grid file.
-               igrd_end   = 1         
-            ENDIF
-
-            DO igrd = igrd_start, igrd_end                     !  loop over T, U & V grid  !
-               !                                               !---------------------------!
-               CALL iom_open( clfile(igrd), inum )
-               CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy, cdunits=clunits ) 
-
-               SELECT CASE( igrd )
-                  CASE (1)   ;   numbdyt = inum
-                  CASE (2)   ;   numbdyu = inum
-                  CASE (3)   ;   numbdyv = inum
-               END SELECT
-
-               ! Calculate time offset 
-               READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0
-               ! Convert time origin in file to julian days 
-               isec0 = isec0 + ihours0*60.*60. + iminutes0*60.
-               CALL ymds2ju(iyear0, imonth0, iday0, REAL(isec0, wp), dayjul0)
-               ! Compute model initialization time 
-               iyear  = ndastp / 10000
-               imonth = ( ndastp - iyear * 10000 ) / 100
-               iday   = ndastp - iyear * 10000 - imonth * 100
-               isecs  = dayfrac * 86400
-               CALL ymds2ju(iyear, imonth, iday, REAL(isecs, wp) , zdayjulini)
-               ! offset from initialization date:
-               zoffset = (dayjul0-zdayjulini)*86400
-               !
-7000           FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
-
-               !! TO BE DONE... Check consistency between calendar from file 
-               !! (available optionally from iom_gettime) and calendar in model 
-               !! when calendar in model available outside of IOIPSL.
-
-               IF(lwp) WRITE(numout,*) 'number of times: ',ntimes_bdy
-               IF(lwp) WRITE(numout,*) 'offset: ',zoffset
-               IF(lwp) WRITE(numout,*) 'totime: ',totime
-               IF(lwp) WRITE(numout,*) 'zstepr: ',zstepr(1:ntimes_bdy)
-
-               ! Check that there are not too many times in the file. 
-               IF( ntimes_bdy > jpbtime ) THEN
-                  WRITE(ctmp1,*) 'Check file: ', clfile(igrd), 'jpbtime= ', jpbtime, ' ntimes_bdy= ', ntimes_bdy
-                  CALL ctl_stop( 'Number of time dumps in files exceed jpbtime parameter', ctmp1 )
-               ENDIF
-
-               ! Check that time array increases:
-               it = 1
-               DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 ) 
-                  it = it + 1
-               END DO
-               !
-               IF( it /= ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN
-                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
-                     CALL ctl_stop( 'Time array in unstructured boundary data files',   &
-                        &           'does not continuously increase.'               , ctmp1 )
-               ENDIF
-               !
-               ! Check that times in file span model run time:
-               IF( zstepr(1) + zoffset > 0 ) THEN
-                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
-                     CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 )
-               END IF
-               IF( zstepr(ntimes_bdy) + zoffset < totime ) THEN
-                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
-                     CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 )
-               END IF
-               !
-               SELECT CASE( igrd )
-                  CASE (1)
-                    ntimes_bdyt = ntimes_bdy
-                    zoffsett = zoffset
-                    istept(:) = INT( zstepr(:) + zoffset )
-                    numbdyt = inum
-                  CASE (2)
-                    ntimes_bdyu = ntimes_bdy
-                    zoffsetu = zoffset
-                    istepu(:) = INT( zstepr(:) + zoffset )
-                    numbdyu = inum
-                  CASE (3)
-                    ntimes_bdyv = ntimes_bdy
-                    zoffsetv = zoffset
-                    istepv(:) = INT( zstepr(:) + zoffset )
-                    numbdyv = inum
-               END SELECT
-               !
-            END DO                                         ! end loop over T, U & V grid 
-
-            IF (igrd_start == 1 .and. igrd_end == 3) THEN
-               ! Only test differences if we are reading in 3 files
-               ! Verify time consistency between files  
-               IF( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN
-                  CALL ctl_stop( 'Bdy data files must have the same number of time dumps',   &
-                  &           'Multiple time frequencies not implemented yet'  )
-               ENDIF
-               ntimes_bdy = ntimes_bdyt
-               !
-               IF( zoffsetu /= zoffsett .OR. zoffsetv /= zoffsett ) THEN
-                  CALL ctl_stop( 'Bdy data files must have the same time origin',   &
-                  &           'Multiple time frequencies not implemented yet' )
-               ENDIF
-               zoffset = zoffsett
-            ENDIF
-
-            IF( igrd_start == 1 ) THEN   ;   istep(:) = istept(:)
-            ELSE                         ;   istep(:) = istepu(:)
-            ENDIF
-
-            ! Check number of time dumps:              
-            IF( ntimes_bdy == 1 .AND. .NOT. ln_clim ) THEN
-              CALL ctl_stop( 'There is only one time dump in data files',   &
-                 &           'Choose ln_clim=.true. in namelist for constant bdy forcing.' )
-            ENDIF
-
-            IF( ln_clim ) THEN
-              IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN
-                 CALL ctl_stop( 'For climatological boundary forcing (ln_clim=.true.),',   &
-                    &           'bdy data files must contain 1 or 12 time dumps.' )
-              ELSEIF( ntimes_bdy ==  1 ) THEN
-                IF(lwp) WRITE(numout,*)
-                IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
-              ELSEIF( ntimes_bdy == 12 ) THEN
-                IF(lwp) WRITE(numout,*)
-                IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
-              ENDIF
-            ENDIF
-
-            ! Find index of first record to read (before first model time). 
-            it = 1
-            DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 )
-               it = it + 1
-            END DO
-            nbdy_b = it
-            !
-            IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset
-            IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b
-
-         ENDIF ! endif (nn_dtactl == 1)
-
-
-         ! 1.2  Read first record in file if necessary (ie if nn_dtactl == 1)
-         ! *****************************************************************
-
-         IF( nn_dtactl == 0 ) THEN      ! boundary data arrays are filled with initial conditions
-            !
-            IF (ln_tra_frs) THEN
-               igrd = 1            ! T-points data 
-               DO ib = 1, nblen(igrd)
-                  ii = nbi(ib,igrd)
-                  ij = nbj(ib,igrd)
+      IF(wrk_in_use(2, 22,23) ) THEN
+         CALL ctl_stop('bdy_dta: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
+      END IF
+
+      ! Initialise data arrays once for all from initial conditions where required
+      !---------------------------------------------------------------------------
+      IF( kt .eq. nit000 .and. .not. PRESENT(jit) ) THEN
+
+         ! Calculate depth-mean currents
+         !-----------------------------
+         pu2d => wrk_2d_22
+         pu2d => wrk_2d_23
+
+         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_bdy = 1, nb_bdy
+
+            nblen => idx_bdy(ib_bdy)%nblen
+            nblenrim => idx_bdy(ib_bdy)%nblenrim
+
+            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 
+               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(:) = nblen(:)
+               ELSE
+                  ilen1(:) = nblenrim(:)
+               ENDIF
+               igrd = 1
+               DO ib = 1, ilen1(igrd)
+                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  dta_bdy(ib_bdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+               END DO 
+               igrd = 2
+               DO ib = 1, ilen1(igrd)
+                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  dta_bdy(ib_bdy)%u2d(ib) = pu2d(ii,ij) * umask(ii,ij,1)         
+               END DO 
+               igrd = 3
+               DO ib = 1, ilen1(igrd)
+                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  dta_bdy(ib_bdy)%v2d(ib) = pv2d(ii,ij) * vmask(ii,ij,1)         
+               END DO 
+            ENDIF
+
+            IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 
+               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(:) = nblen(:)
+               ELSE
+                  ilen1(:) = nblenrim(:)
+               ENDIF
+               igrd = 2 
+               DO ib = 1, ilen1(igrd)
                   DO ik = 1, jpkm1
-                     tbdy(ib,ik) = tsn(ii,ij,ik,jp_tem)
-                     sbdy(ib,ik) = tsn(ii,ij,ik,jp_sal)
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - pu2d(ii,ij) ) * umask(ii,ij,ik)         
+                  END DO
+               END DO 
+               igrd = 3 
+               DO ib = 1, ilen1(igrd)
+                  DO ik = 1, jpkm1
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - pv2d(ii,ij) ) * vmask(ii,ij,ik)         
+                     END DO
+               END DO 
+            ENDIF
+
+            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN 
+               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(:) = nblen(:)
+               ELSE
+                  ilen1(:) = nblenrim(:)
+               ENDIF
+               igrd = 1                       ! Everything is at T-points here
+               DO ib = 1, ilen1(igrd)
+                  DO ik = 1, jpkm1
+                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                     dta_bdy(ib_bdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
+                     dta_bdy(ib_bdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
+                  END DO
+               END DO 
+            ENDIF
+
+#if defined key_lim2
+            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 
+               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(:) = nblen(:)
+               ELSE
+                  ilen1(:) = nblenrim(:)
+               ENDIF
+               igrd = 1                       ! Everything is at T-points here
+               DO ib = 1, ilen1(igrd)
+                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  dta_bdy(ib_bdy)%frld(ib) = frld(ii,ij) * tmask(ii,ij,1)         
+                  dta_bdy(ib_bdy)%hicif(ib) = hicif(ii,ij) * tmask(ii,ij,1)         
+                  dta_bdy(ib_bdy)%hsnif(ib) = hsnif(ii,ij) * tmask(ii,ij,1)         
+               END DO 
+            ENDIF
+#endif
+
+         ENDDO ! ib_bdy
+
+      ENDIF ! kt .eq. nit000
+
+      ! update external data from files
+      !--------------------------------
+     
+      jstart = 1
+      DO ib_bdy = 1, nb_bdy   
+         IF( nn_dta(ib_bdy) .eq. 1 ) THEN ! skip this bit if no external data required
+      
+            IF( PRESENT(jit) ) THEN
+               ! Update barotropic boundary conditions only
+               ! jit is optional argument for fld_read and tide_update
+               IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
+                  IF( nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     dta_bdy(ib_bdy)%ssh(:) = 0.0
+                     dta_bdy(ib_bdy)%u2d(:) = 0.0
+                     dta_bdy(ib_bdy)%v2d(:) = 0.0
+                  ENDIF
+                  IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data
+                     jend = jstart + 2
+                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset )
+                  ENDIF
+                  IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                     CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), jit=jit, time_offset=time_offset )
+                  ENDIF
+               ENDIF
+            ELSE
+               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                  dta_bdy(ib_bdy)%ssh(:) = 0.0
+                  dta_bdy(ib_bdy)%u2d(:) = 0.0
+                  dta_bdy(ib_bdy)%v2d(:) = 0.0
+               ENDIF
+               IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data
+                  jend = jstart + nb_bdy_fld(ib_bdy) - 1
+                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )
+               ENDIF
+               IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                  CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )
+               ENDIF
+            ENDIF
+            jstart = jend+1
+
+            ! If full velocities in boundary data then split into barotropic and baroclinic data
+            ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same
+            ! time as the dynspg_ts option). 
+
+            IF( ln_full_vel_array(ib_bdy) .and.                                             & 
+           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 
+
+               igrd = 2                      ! zonal velocity
+               dta_bdy(ib_bdy)%u2d(:) = 0.0
+               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
+                  DO ik = 1, jpkm1
+                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) &
+              &                                + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik)
+                  END DO
+                  dta_bdy(ib_bdy)%u2d(ib) =  dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij)
+                  DO ik = 1, jpkm1
+                     dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 
                   END DO
                END DO
-            ENDIF
-
-            IF(ln_dyn_frs) THEN
-               igrd = 2            ! U-points data 
-               DO ib = 1, nblen(igrd)
-                  ii = nbi(ib,igrd)
-                  ij = nbj(ib,igrd)
+
+               igrd = 3                      ! meridional velocity
+               dta_bdy(ib_bdy)%v2d(:) = 0.0
+               DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+                  ii   = idx_bdy(ib_bdy)%nbi(ib,igrd)
+                  ij   = idx_bdy(ib_bdy)%nbj(ib,igrd)
                   DO ik = 1, jpkm1
-                     ubdy(ib,ik) = un(ii, ij, ik)
+                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) &
+              &                                + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik)
+                  END DO
+                  dta_bdy(ib_bdy)%v2d(ib) =  dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij)
+                  DO ik = 1, jpkm1
+                     dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 
                   END DO
                END DO
-               !
-               igrd = 3            ! V-points data 
-               DO ib = 1, nblen(igrd)            
-                  ii = nbi(ib,igrd)
-                  ij = nbj(ib,igrd)
-                  DO ik = 1, jpkm1
-                     vbdy(ib,ik) = vn(ii, ij, ik)
-                  END DO
-               END DO
-            ENDIF
-            !
-#if defined key_lim2
-            IF( ln_ice_frs ) THEN
-               igrd = 1            ! T-points data
-               DO ib = 1, nblen(igrd)
-                  frld_bdy (ib) =  frld(nbi(ib,igrd), nbj(ib,igrd))
-                  hicif_bdy(ib) = hicif(nbi(ib,igrd), nbj(ib,igrd))
-                  hsnif_bdy(ib) = hsnif(nbi(ib,igrd), nbj(ib,igrd))
-               END DO
-            ENDIF
-#endif
-         ELSEIF( nn_dtactl == 1 ) THEN    ! Set first record in the climatological case:   
-            !
-            IF( ln_clim .AND. ntimes_bdy == 1 ) THEN
-               nbdy_a = 1
-            ELSEIF( ln_clim .AND. ntimes_bdy == iman ) THEN
-               nbdy_b = 0
-               nbdy_a = imois
+    
+            ENDIF
+
+         END IF ! nn_dta(ib_bdy) = 1
+      END DO  ! ib_bdy
+
+      IF(wrk_not_released(2, 22,23) )    CALL ctl_stop('bdy_dta: ERROR: failed to release workspace arrays.')
+
+      END SUBROUTINE bdy_dta
+
+
+      SUBROUTINE bdy_dta_init
+      !!----------------------------------------------------------------------
+      !!                   ***  SUBROUTINE bdy_dta_init  ***
+      !!                    
+      !! ** Purpose :   Initialise arrays for reading of external data 
+      !!                for open boundary conditions
+      !!
+      !! ** Method  :   Use fldread.F90
+      !!                
+      !!----------------------------------------------------------------------
+      USE dynspg_oce, ONLY: lk_dynspg_ts
+      !!
+      INTEGER     ::  ib_bdy, jfld, jstart, jend, ierror  ! local indices
+      !!
+      CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
+      CHARACTER(len=100), DIMENSION(nb_bdy)  ::   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 bdy 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(:)     ::   ibdy           ! bdy set for a particular jfld
+      INTEGER, ALLOCATABLE, DIMENSION(:)     ::   igrid         ! index for grid type (1,2,3 = T,U,V)
+      INTEGER, POINTER, DIMENSION(:)         ::   nblen, nblenrim  ! short cuts
+      TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   blf_i         !  array of namelist information structures
+      TYPE(FLD_N) ::   bn_tem, bn_sal, bn_u3d, bn_v3d   ! 
+      TYPE(FLD_N) ::   bn_ssh, bn_u2d, bn_v2d           ! informations about the fields to be read
+#if defined key_lim2
+      TYPE(FLD_N) ::   bn_frld, bn_hicif, bn_hsnif      !
+#endif
+      NAMELIST/nambdy_dta/ cn_dir, bn_tem, bn_sal, bn_u3d, bn_v3d, bn_ssh, bn_u2d, bn_v2d 
+#if defined key_lim2
+      NAMELIST/nambdy_dta/ bn_frld, bn_hicif, bn_hsnif
+#endif
+      NAMELIST/nambdy_dta/ ln_full_vel
+      !!---------------------------------------------------------------------------
+
+      ! Set nn_dta
+      DO ib_bdy = 1, nb_bdy
+         nn_dta(ib_bdy) = MAX(  nn_dyn2d_dta(ib_bdy)       &
+                               ,nn_dyn3d_dta(ib_bdy)       &
+                               ,nn_tra_dta(ib_bdy)         &
+#if defined key_ice_lim2
+                               ,nn_ice_lim2_dta(ib_bdy)    &
+#endif
+                              )
+         IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1
+      END DO
+
+      ! Work out upper bound of how many fields there are to read in and allocate arrays
+      ! ---------------------------------------------------------------------------
+      ALLOCATE( nb_bdy_fld(nb_bdy) )
+      nb_bdy_fld(:) = 0
+      DO ib_bdy = 1, nb_bdy         
+         IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN
+            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
+         ENDIF
+         IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) THEN
+            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
+         ENDIF
+         IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1  ) THEN
+            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 2
+         ENDIF
+#if defined key_lim2
+         IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1  ) THEN
+            nb_bdy_fld(ib_bdy) = nb_bdy_fld(ib_bdy) + 3
+         ENDIF
+#endif               
+      ENDDO            
+
+      nb_bdy_fld_sum = SUM( nb_bdy_fld )
+
+      ALLOCATE( bf(nb_bdy_fld_sum), STAT=ierror )
+      IF( ierror > 0 ) THEN   
+         CALL ctl_stop( 'bdy_dta: unable to allocate bf structure' )   ;   RETURN  
+      ENDIF
+      ALLOCATE( blf_i(nb_bdy_fld_sum), STAT=ierror )
+      IF( ierror > 0 ) THEN   
+         CALL ctl_stop( 'bdy_dta: unable to allocate blf_i structure' )   ;   RETURN  
+      ENDIF
+      ALLOCATE( nbmap_ptr(nb_bdy_fld_sum), STAT=ierror )
+      IF( ierror > 0 ) THEN   
+         CALL ctl_stop( 'bdy_dta: unable to allocate nbmap_ptr structure' )   ;   RETURN  
+      ENDIF
+      ALLOCATE( ilen1(nb_bdy_fld_sum), ilen3(nb_bdy_fld_sum) ) 
+      ALLOCATE( ibdy(nb_bdy_fld_sum) ) 
+      ALLOCATE( igrid(nb_bdy_fld_sum) ) 
+
+      ! Read namelists
+      ! --------------
+      REWIND(numnam)
+      jfld = 0 
+      DO ib_bdy = 1, nb_bdy         
+         IF( nn_dta(ib_bdy) .eq. 1 ) THEN
+            ! 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  !
+            !                    !  name       !  (hours)  !   name           !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs     !
+            bn_ssh     = FLD_N(  'bdy_ssh'     ,    24     ,  'sossheig'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_u2d     = FLD_N(  'bdy_vel2d_u' ,    24     ,  'vobtcrtx'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_v2d     = FLD_N(  'bdy_vel2d_v' ,    24     ,  'vobtcrty'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_u3d     = FLD_N(  'bdy_vel3d_u' ,    24     ,  'vozocrtx'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_v3d     = FLD_N(  'bdy_vel3d_v' ,    24     ,  'vomecrty'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_tem     = FLD_N(  'bdy_tem'     ,    24     ,  'votemper'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_sal     = FLD_N(  'bdy_sal'     ,    24     ,  'vosaline'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+#if defined key_lim2
+            bn_frld    = FLD_N(  'bdy_frld'    ,    24     ,  'ildsconc'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_hicif   = FLD_N(  'bdy_hicif'   ,    24     ,  'iicethic'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+            bn_hsnif   = FLD_N(  'bdy_hsnif'   ,    24     ,  'isnothic'    ,  .false.   , .false. ,   'yearly'  , ''       , ''        )
+#endif
+
+            ! Important NOT to rewind here.
+            READ( numnam, nambdy_dta )
+
+            cn_dir_array(ib_bdy) = cn_dir
+            ln_full_vel_array(ib_bdy) = ln_full_vel
+
+            IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts )  THEN
+               CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',&
+            &                  'with dynspg_ts option' )   ;   RETURN  
+            ENDIF             
+
+            nblen => idx_bdy(ib_bdy)%nblen
+            nblenrim => idx_bdy(ib_bdy)%nblenrim
+
+            ! Only read in necessary fields for this set.
+            ! Important that barotropic variables come first.
+            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
+
+               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_ssh
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 1
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+                  ilen3(jfld) = 1
+               ENDIF
+
+               IF( .not. ln_full_vel_array(ib_bdy) ) THEN
+
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_u2d
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 2
+                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+                     ilen1(jfld) = nblen(igrid(jfld))
+                  ELSE
+                     ilen1(jfld) = nblenrim(igrid(jfld))
+                  ENDIF
+                  ilen3(jfld) = 1
+
+                  jfld = jfld + 1
+                  blf_i(jfld) = bn_v2d
+                  ibdy(jfld) = ib_bdy
+                  igrid(jfld) = 3
+                  IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+                     ilen1(jfld) = nblen(igrid(jfld))
+                  ELSE
+                     ilen1(jfld) = nblenrim(igrid(jfld))
+                  ENDIF
+                  ilen3(jfld) = 1
+
+               ENDIF
+
+            ENDIF
+
+            ! baroclinic velocities
+            IF( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ) .or. &
+           &      ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.  &
+           &        ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_u3d
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 2
+               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = jpk
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_v3d
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 3
+               IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = jpk
+
+            ENDIF
+
+            ! temperature and salinity
+            IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 1 ) THEN
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_tem
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 1
+               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = jpk
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_sal
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 1
+               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = jpk
+
+            ENDIF
+
+#if defined key_lim2
+            ! sea ice
+            IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 1 ) THEN
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_frld
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 1
+               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = 1
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_hicif
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 1
+               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = 1
+
+               jfld = jfld + 1
+               blf_i(jfld) = bn_hsnif
+               ibdy(jfld) = ib_bdy
+               igrid(jfld) = 1
+               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
+                  ilen1(jfld) = nblen(igrid(jfld))
+               ELSE
+                  ilen1(jfld) = nblenrim(igrid(jfld))
+               ENDIF
+               ilen3(jfld) = 1
+
+            ENDIF
+#endif
+            ! Recalculate field counts
+            !-------------------------
+            nb_bdy_fld_sum = 0
+            IF( ib_bdy .eq. 1 ) THEN 
+               nb_bdy_fld(ib_bdy) = jfld
+               nb_bdy_fld_sum     = jfld              
             ELSE
-               nbdy_a = nbdy_b
-            ENDIF
-   
-            ! Read first record:
-            ipj  = 1
-            ipk  = jpk
-            igrd = 1
-            ipi  = nblendta(igrd)
-
-            IF(ln_tra_frs) THEN
-               !
-               igrd = 1                                           ! Temperature
-               IF( nblendta(igrd) <=  0 ) THEN 
-                  idvar = iom_varid( numbdyt, 'votemper' )
-                  nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
-               ENDIF
-               IF(lwp) WRITE(numout,*) 'Dim size for votemper is ', nblendta(igrd)
-               ipi = nblendta(igrd)
-               CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
-               !
-               DO ib = 1, nblen(igrd)
-                  DO ik = 1, jpkm1
-                     tbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
-                  END DO
-               END DO
-               !
-               igrd = 1                                           ! salinity
-               IF( nblendta(igrd) .le. 0 ) THEN 
-                  idvar = iom_varid( numbdyt, 'vosaline' )
-                  nblendta(igrd) = iom_file(numbdyt)%dimsz(1,idvar)
-               ENDIF
-               IF(lwp) WRITE(numout,*) 'Dim size for vosaline is ', nblendta(igrd)
-               ipi = nblendta(igrd)
-               CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
-               !
-               DO ib = 1, nblen(igrd)
-                  DO ik = 1, jpkm1
-                     sbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
-                  END DO
-               END DO
-            ENDIF  ! ln_tra_frs
- 
-            IF( ln_dyn_frs ) THEN
-               !
-               igrd = 2                                           ! u-velocity
-               IF ( nblendta(igrd) .le. 0 ) THEN 
-                 idvar = iom_varid( numbdyu,'vozocrtx' )
-                 nblendta(igrd) = iom_file(numbdyu)%dimsz(1,idvar)
-               ENDIF
-               IF(lwp) WRITE(numout,*) 'Dim size for vozocrtx is ', nblendta(igrd)
-               ipi = nblendta(igrd)
-               CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
-               DO ib = 1, nblen(igrd)
-                  DO ik = 1, jpkm1
-                     ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
-                  END DO
-               END DO
-               !
-               igrd = 3                                           ! v-velocity
-               IF ( nblendta(igrd) .le. 0 ) THEN 
-                 idvar = iom_varid( numbdyv,'vomecrty' )
-                 nblendta(igrd) = iom_file(numbdyv)%dimsz(1,idvar)
-               ENDIF
-               IF(lwp) WRITE(numout,*) 'Dim size for vomecrty is ', nblendta(igrd)
-               ipi = nblendta(igrd)
-               CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
-               DO ib = 1, nblen(igrd)
-                  DO ik = 1, jpkm1
-                     vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
-                  END DO
-               END DO
-            ENDIF ! ln_dyn_frs
-
-#if defined key_lim2
-            IF( ln_ice_frs ) THEN
-              !
-              igrd=1                                              ! leads fraction
-              IF(lwp) WRITE(numout,*) 'Dim size for ildsconc is ',nblendta(igrd)
-              ipi=nblendta(igrd)
-              CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a )
-              DO ib=1, nblen(igrd)
-                frld_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
-              END DO
-              !
-              igrd=1                                              ! ice thickness
-              IF(lwp) WRITE(numout,*) 'Dim size for iicethic is ',nblendta(igrd)
-              ipi=nblendta(igrd)
-              CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a )
-              DO ib=1, nblen(igrd)
-                hicif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
-              END DO
-              !
-              igrd=1                                              ! snow thickness
-              IF(lwp) WRITE(numout,*) 'Dim size for isnowthi is ',nblendta(igrd)
-              ipi=nblendta(igrd)
-              CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a )
-              DO ib=1, nblen(igrd)
-                hsnif_bdydta(ib,2) =  zdta(nbmap(ib,igrd),1,1)
-              END DO
-            ENDIF ! just if ln_ice_frs is set
-#endif
-
-            IF( .NOT.ln_clim .AND. istep(1) > 0 ) THEN     ! First data time is after start of run
-               nbdy_b = nbdy_a                                 ! Put first value in both time levels
-               IF( ln_tra_frs ) THEN
-                 tbdydta(:,:,1) = tbdydta(:,:,2)
-                 sbdydta(:,:,1) = sbdydta(:,:,2)
-               ENDIF
-               IF( ln_dyn_frs ) THEN
-                 ubdydta(:,:,1) = ubdydta(:,:,2)
-                 vbdydta(:,:,1) = vbdydta(:,:,2)
-               ENDIF
-#if defined key_lim2
-               IF( ln_ice_frs ) THEN
-                  frld_bdydta (:,1) =  frld_bdydta(:,2)
-                  hicif_bdydta(:,1) = hicif_bdydta(:,2)
-                  hsnif_bdydta(:,1) = hsnif_bdydta(:,2)
-               ENDIF
-#endif
-            END IF
-            !
-         END IF   ! nn_dtactl == 0/1
- 
-         ! In the case of constant boundary forcing fill bdy arrays once for all
-         IF( ln_clim .AND. ntimes_bdy == 1 ) THEN
-            IF( ln_tra_frs ) THEN
-               tbdy  (:,:) = tbdydta  (:,:,2)
-               sbdy  (:,:) = sbdydta  (:,:,2)
-            ENDIF
-            IF( ln_dyn_frs) THEN
-               ubdy  (:,:) = ubdydta  (:,:,2)
-               vbdy  (:,:) = vbdydta  (:,:,2)
-            ENDIF
-#if defined key_lim2
-            IF( ln_ice_frs ) THEN
-               frld_bdy (:) = frld_bdydta (:,2)
-               hicif_bdy(:) = hicif_bdydta(:,2)
-               hsnif_bdy(:) = hsnif_bdydta(:,2)
-            ENDIF
-#endif
-
-            IF( ln_tra_frs .OR. ln_ice_frs) CALL iom_close( numbdyt )
-            IF( ln_dyn_frs                    ) CALL iom_close( numbdyu )
-            IF( ln_dyn_frs                    ) CALL iom_close( numbdyv )
-         END IF
-         !
-      ENDIF                                            ! End if nit000
-
-
-      !                                                !---------------------!
-      IF( nn_dtactl == 1 .AND. ntimes_bdy > 1 ) THEN    !  at each time step  !
-         !                                             !---------------------!
-         ! Read one more record if necessary
-         !**********************************
-
-         IF( ln_clim .AND. imois /= nbdy_b ) THEN      ! remember that nbdy_b=0 for kt=nit000
-            nbdy_b = imois
-            nbdy_a = imois + 1
-            nbdy_b = MOD( nbdy_b, iman )   ;   IF( nbdy_b == 0 ) nbdy_b = iman
-            nbdy_a = MOD( nbdy_a, iman )   ;   IF( nbdy_a == 0 ) nbdy_a = iman
-            lect=.true.
-         ELSEIF( .NOT.ln_clim .AND. itimer >= istep(nbdy_a) ) THEN
-
-            IF( nbdy_a < ntimes_bdy ) THEN
-               nbdy_b = nbdy_a
-               nbdy_a = nbdy_a + 1
-               lect  =.true.
+               nb_bdy_fld(ib_bdy) = jfld - nb_bdy_fld_sum
+               nb_bdy_fld_sum = nb_bdy_fld_sum + nb_bdy_fld(ib_bdy)
+            ENDIF
+
+         ENDIF ! nn_dta .eq. 1
+      ENDDO ! ib_bdy
+
+
+      DO jfld = 1, nb_bdy_fld_sum
+         ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) )
+         IF( blf_i(jfld)%ln_tint ) ALLOCATE( bf(jfld)%fdta(ilen1(jfld),1,ilen3(jfld),2) )
+         nbmap_ptr(jfld)%ptr => idx_bdy(ibdy(jfld))%nbmap(:,igrid(jfld))
+      ENDDO
+
+      ! fill bf with blf_i and control print
+      !-------------------------------------
+      jstart = 1
+      DO ib_bdy = 1, nb_bdy
+         jend = jstart + nb_bdy_fld(ib_bdy) - 1
+         CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta', 'open boundary conditions', 'nambdy_dta' )
+         jstart = jend + 1
+      ENDDO
+
+      ! Initialise local boundary data arrays
+      ! nn_xxx_dta=0 : allocate space - will be filled from initial conditions later
+      ! nn_xxx_dta=1 : point to "fnow" arrays
+      !-------------------------------------
+
+      jfld = 0
+      DO ib_bdy=1, nb_bdy
+
+         nblen => idx_bdy(ib_bdy)%nblen
+         nblenrim => idx_bdy(ib_bdy)%nblenrim
+
+         IF (nn_dyn2d(ib_bdy) .gt. 0) THEN
+            IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN
+               IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  ilen0(1:3) = nblenrim(1:3)
+               ENDIF
+               ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) )
+               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) )
+               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) )
             ELSE
-               ! We have reached the end of the file
-               ! put the last data time into both time levels
-               nbdy_b = nbdy_a
-               IF(ln_tra_frs) THEN
-                  tbdydta(:,:,1) =  tbdydta(:,:,2)
-                  sbdydta(:,:,1) =  sbdydta(:,:,2)
-               ENDIF
-               IF(ln_dyn_frs) THEN
-                  ubdydta(:,:,1) =  ubdydta(:,:,2)
-                  vbdydta(:,:,1) =  vbdydta(:,:,2)
-               ENDIF
-#if defined key_lim2
-               IF(ln_ice_frs) THEN
-                  frld_bdydta (:,1) =  frld_bdydta (:,2)
-                  hicif_bdydta(:,1) =  hicif_bdydta(:,2)
-                  hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
-               ENDIF
-#endif
-            END IF ! nbdy_a < ntimes_bdy
-            !
-        END IF
-         
-        IF( lect ) THEN           ! Swap arrays
-           IF( ln_tra_frs ) THEN
-             tbdydta(:,:,1) =  tbdydta(:,:,2)
-             sbdydta(:,:,1) =  sbdydta(:,:,2)
-           ENDIF
-           IF( ln_dyn_frs ) THEN
-             ubdydta(:,:,1) =  ubdydta(:,:,2)
-             vbdydta(:,:,1) =  vbdydta(:,:,2)
-           ENDIF
-#if defined key_lim2
-           IF( ln_ice_frs ) THEN
-             frld_bdydta (:,1) =  frld_bdydta (:,2)
-             hicif_bdydta(:,1) =  hicif_bdydta(:,2)
-             hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
-           ENDIF
-#endif 
-           ! read another set
-           ipj  = 1
-           ipk  = jpk
-
-           IF( ln_tra_frs ) THEN
-              ! 
-              igrd = 1                                   ! temperature
-              ipi  = nblendta(igrd)
-              CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
-              DO ib = 1, nblen(igrd)
-                 DO ik = 1, jpkm1
-                    tbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
-                 END DO
-              END DO
-              !
-              igrd = 1                                   ! salinity
-              ipi  = nblendta(igrd)
-              CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
-              DO ib = 1, nblen(igrd)
-                 DO ik = 1, jpkm1
-                    sbdydta(ib,ik,2) = zdta(nbmap(ib,igrd),1,ik)
-                 END DO
-              END DO
-           ENDIF ! ln_tra_frs
-
-           IF(ln_dyn_frs) THEN
-              !
-              igrd = 2                                   ! u-velocity
-              ipi  = nblendta(igrd)
-              CALL iom_get ( numbdyu, jpdom_unknown,'vozocrtx',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
-              DO ib = 1, nblen(igrd)
-                DO ik = 1, jpkm1
-                  ubdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
-                END DO
-              END DO
-              !
-              igrd = 3                                   ! v-velocity
-              ipi  = nblendta(igrd)
-              CALL iom_get ( numbdyv, jpdom_unknown,'vomecrty',zdta(1:ipi,1:ipj,1:ipk),nbdy_a )
-              DO ib = 1, nblen(igrd)
-                 DO ik = 1, jpkm1
-                    vbdydta(ib,ik,2) =  zdta(nbmap(ib,igrd),1,ik)
-                 END DO
-              END DO
-           ENDIF ! ln_dyn_frs
-           !
-#if defined key_lim2
-           IF(ln_ice_frs) THEN
-             !
-             igrd = 1                                    ! ice concentration
-             ipi=nblendta(igrd)
-             CALL iom_get ( numbdyt, jpdom_unknown,'ildsconc',zdta(1:ipi,:,1),nbdy_a )
-             DO ib=1, nblen(igrd)
-               frld_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
-             END DO
-             !
-             igrd=1                                      ! ice thickness
-             ipi=nblendta(igrd)
-             CALL iom_get ( numbdyt, jpdom_unknown,'iicethic',zdta(1:ipi,:,1),nbdy_a )
-             DO ib=1, nblen(igrd)
-               hicif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
-             END DO
-             !
-             igrd=1                                      ! snow thickness
-             ipi=nblendta(igrd)
-             CALL iom_get ( numbdyt, jpdom_unknown,'isnowthi',zdta(1:ipi,:,1),nbdy_a )
-             DO ib=1, nblen(igrd)
-               hsnif_bdydta(ib,2) =  zdta( nbmap(ib,igrd), 1, 1 )
-             END DO
-           ENDIF ! ln_ice_frs
-#endif
-           !
-           IF(lwp) WRITE(numout,*) 'bdy_dta_frs : first record file used nbdy_b ',nbdy_b
-           IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a ',nbdy_a
-           IF (.NOT.ln_clim) THEN
-              IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep(nbdy_b)
-              IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer
-              IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy_a)
-           ENDIF
-           !
-       ENDIF ! end lect=.true.
-
-
-       ! Interpolate linearly
-       ! ********************
-       ! 
-       IF( ln_clim ) THEN   ;   zxy = REAL( nday                   ) / REAL( nmonth_len(nbdy_b) ) + 0.5 - i15
-       ELSEIF( istep(nbdy_b) == istep(nbdy_a) ) THEN 
-                                    zxy = 0.0_wp
-       ELSE                     ;   zxy = REAL( istep(nbdy_b) - itimer ) / REAL( istep(nbdy_b) - istep(nbdy_a) )
-       END IF
-
-          IF(ln_tra_frs) THEN
-             igrd = 1                                   ! temperature & salinity
-             DO ib = 1, nblen(igrd)
-               DO ik = 1, jpkm1
-                 tbdy(ib,ik) = zxy * tbdydta(ib,ik,2) + (1.-zxy) * tbdydta(ib,ik,1)
-                 sbdy(ib,ik) = zxy * sbdydta(ib,ik,2) + (1.-zxy) * sbdydta(ib,ik,1)
-               END DO
-             END DO
-          ENDIF
-
-          IF(ln_dyn_frs) THEN
-             igrd = 2                                   ! u-velocity
-             DO ib = 1, nblen(igrd)
-               DO ik = 1, jpkm1
-                 ubdy(ib,ik) = zxy * ubdydta(ib,ik,2) + (1.-zxy) * ubdydta(ib,ik,1)   
-               END DO
-             END DO
-             !
-             igrd = 3                                   ! v-velocity
-             DO ib = 1, nblen(igrd)
-               DO ik = 1, jpkm1
-                 vbdy(ib,ik) = zxy * vbdydta(ib,ik,2) + (1.-zxy) * vbdydta(ib,ik,1)   
-               END DO
-             END DO
-          ENDIF
-
-#if defined key_lim2
-          IF(ln_ice_frs) THEN
-            igrd=1
-            DO ib=1, nblen(igrd)
-               frld_bdy(ib) = zxy *  frld_bdydta(ib,2) + (1.-zxy) *  frld_bdydta(ib,1)
-              hicif_bdy(ib) = zxy * hicif_bdydta(ib,2) + (1.-zxy) * hicif_bdydta(ib,1)
-              hsnif_bdy(ib) = zxy * hsnif_bdydta(ib,2) + (1.-zxy) * hsnif_bdydta(ib,1)
-            END DO
-          ENDIF ! just if ln_ice_frs is true
-#endif
-
-      END IF                       !end if ((nn_dtactl==1).AND.(ntimes_bdy>1))
-    
-
-      !                                                !---------------------!
-      !                                                !     last call       !
-      !                                                !---------------------!
-      IF( kt == nitend ) THEN
-          IF(ln_tra_frs .or. ln_ice_frs) CALL iom_close( numbdyt )              ! Closing of the 3 files
-          IF(ln_dyn_frs) CALL iom_close( numbdyu )
-          IF(ln_dyn_frs) CALL iom_close( numbdyv )
-      ENDIF
-      !
-      ENDIF ! ln_dyn_frs .OR. ln_tra_frs
-      !
-   END SUBROUTINE bdy_dta_frs
-
-
-   SUBROUTINE bdy_dta_fla( kt, jit, icycl )
-      !!---------------------------------------------------------------------------
-      !!                      ***  SUBROUTINE bdy_dta_fla  ***
-      !!                    
-      !! ** Purpose :   Read unstructured boundary data for Flather condition
-      !!
-      !! ** Method  :  At the first timestep, read in boundary data for two
-      !!               times from the file and time-interpolate. At other 
-      !!               timesteps, check to see if we need another time from 
-      !!               the file. If so read it in. Time interpolate.
-      !!---------------------------------------------------------------------------
-!!gm DOCTOR names :   argument integer :  start with "k"
-      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
-      INTEGER, INTENT( in ) ::   jit         ! barotropic time step index
-      INTEGER, INTENT( in ) ::   icycl       ! number of cycles need for final file close
-      !                                      ! (for timesplitting option, otherwise zero)
-      !!
-      LOGICAL ::   lect                      ! flag for reading
-      INTEGER ::   it, ib, igrd              ! dummy loop indices
-      INTEGER ::   idvar                     ! netcdf var ID
-      INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing
-      INTEGER ::   ntimes_bdyt, ntimes_bdyu, ntimes_bdyv
-      INTEGER ::   itimer, totime
-      INTEGER ::   ipi, ipj, ipk, inum       ! temporary integers (NetCDF read)
-      INTEGER ::   iyear0, imonth0, iday0
-      INTEGER ::   ihours0, iminutes0, isec0
-      INTEGER ::   iyear, imonth, iday, isecs
-      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files
-      REAL(wp) ::   dayfrac, zxy, zoffsett
-      REAL(wp) ::   zoffsetu, zoffsetv
-      REAL(wp) ::   dayjul0, zdayjulini
-      REAL(wp) ::   zinterval_s, zinterval_e                    ! First and last interval in time axis
-      REAL(wp), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files
-      REAL(wp), DIMENSION(jpbdta,1)     ::   zdta               ! temporary array for data fields
-      CHARACTER(LEN=80), DIMENSION(6)   ::   clfile
-      CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate
-      !!---------------------------------------------------------------------------
-
-!!gm   add here the same style as in bdy_dta_frs
-!!gm      clearly bdy_dta_fla and bdy_dta_frs  can be combined...   
-!!gm      too many things duplicated in the read of data...   simplification can be done
-
-      ! -------------------- !
-      !    Initialization    !
-      ! -------------------- !
-
-      lect   = .false.           ! If true, read a time record
-
-      ! Some time variables for monthly climatological forcing:
-      ! *******************************************************
- !!gm  here  use directely daymod variables
- 
-      iman  = INT( raamo ) ! Number of months in a year
-
-      i15 = INT( 2*REAL( nday, wp ) / ( REAL( nmonth_len(nmonth), wp ) + 0.5 ) )
-      ! i15=0 if the current day is in the first half of the month, else i15=1
-
-      imois = nmonth + i15 - 1            ! imois is the first month record
-      IF( imois == 0 ) imois = iman
-
-      ! Time variable for non-climatological forcing:
-      ! *********************************************
-
-      itimer = ((kt-1)-nit000+1)*rdt                      ! current time in seconds for interpolation 
-      itimer = itimer + jit*rdt/REAL(nn_baro,wp)      ! in non-climatological case
-
-      IF ( ln_tides ) THEN
-
-         ! -------------------------------------!
-         ! Update BDY fields with tidal forcing !
-         ! -------------------------------------!  
-
-         CALL tide_update( kt, jit ) 
-  
-      ENDIF
-
-      IF ( ln_dyn_fla ) THEN
-
-         ! -------------------------------------!
-         ! Update BDY fields with model data    !
-         ! -------------------------------------!  
-
-      !                                                !-------------------!
-      IF( kt == nit000 .and. jit ==2 ) THEN            !  First call only  !
-         !                                             !-------------------!
-         istep_bt(:) = 0
-         nbdy_b_bt    = 0
-         nbdy_a_bt    = 0
-
-         ! Get time information from bdy data file
-         ! ***************************************
-
-        IF(lwp) WRITE(numout,*)
-        IF(lwp) WRITE(numout,*)    'bdy_dta_fla :Initialize unstructured boundary data for barotropic variables.'
-        IF(lwp) WRITE(numout,*)    '~~~~~~~' 
-
-        IF( nn_dtactl == 0 ) THEN
-          IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions'
-
-        ELSEIF (nn_dtactl == 1) THEN
-          IF(lwp) WRITE(numout,*)  'Bdy data are read in netcdf files'
-
-          dayfrac = adatrj  - REAL(itimer,wp)/86400. ! day fraction at time step kt-1
-          dayfrac = dayfrac - INT (dayfrac)          !
-          totime = (nitend-nit000+1)*rdt             ! Total time of the run to verify that all the
-                                                     ! necessary time dumps in file are included
-
-          clfile(4) = cn_dta_fla_T
-          clfile(5) = cn_dta_fla_U
-          clfile(6) = cn_dta_fla_V
-
-          DO igrd = 4,6
-
-            CALL iom_open( clfile(igrd), inum )
-            CALL iom_gettime( inum, zstepr, kntime=ntimes_bdy_bt, cdunits=clunits ) 
-
-            SELECT CASE( igrd )
-               CASE (4) 
-                  numbdyt_bt = inum
-               CASE (5) 
-                  numbdyu_bt = inum
-               CASE (6) 
-                  numbdyv_bt = inum
-            END SELECT
-
-            ! Calculate time offset 
-            READ(clunits,7000) iyear0, imonth0, iday0, ihours0, iminutes0, isec0
-            ! Convert time origin in file to julian days 
-            isec0 = isec0 + ihours0*60.*60. + iminutes0*60.
-            CALL ymds2ju(iyear0, imonth0, iday0, REAL(isec0, wp), dayjul0)
-            ! Compute model initialization time 
-            iyear  = ndastp / 10000
-            imonth = ( ndastp - iyear * 10000 ) / 100
-            iday   = ndastp - iyear * 10000 - imonth * 100
-            isecs  = dayfrac * 86400
-            CALL ymds2ju(iyear, imonth, iday, REAL(isecs, wp) , zdayjulini)
-            ! zoffset from initialization date:
-            zoffset = (dayjul0-zdayjulini)*86400
-            !
-
-7000 FORMAT('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
-
-            !! TO BE DONE... Check consistency between calendar from file 
-            !! (available optionally from iom_gettime) and calendar in model 
-            !! when calendar in model available outside of IOIPSL.
-
-            ! Check that there are not too many times in the file. 
-            IF (ntimes_bdy_bt > jpbtime) CALL ctl_stop( &
-                 'Number of time dumps in bdy file exceed jpbtime parameter', &
-                 'Check file:' // TRIM(clfile(igrd))  )
-
-            ! Check that time array increases (or interp will fail):
-            DO it = 2, ntimes_bdy_bt
-               IF ( zstepr(it-1) >= zstepr(it) ) THEN
-                  CALL ctl_stop('Time array in unstructured boundary data file', &
-                       'does not continuously increase.',               &
-                       'Check file:' // TRIM(clfile(igrd))  )
-                  EXIT
-               END IF
-            END DO
-
-            IF ( .NOT. ln_clim ) THEN
-               ! Check that times in file span model run time:
-
-               ! Note: the fields may be time means, so we allow nit000 to be before
-               ! first time in the file, provided that it falls inside the meaning
-               ! period of the first field.  Until we can get the meaning period
-               ! from the file, use the interval between fields as a proxy.
-               ! If nit000 is before the first time, use the value at first time
-               ! instead of extrapolating.  This is done by putting time 1 into
-               ! both time levels.
-               ! The same applies to the last time level: see setting of lect below.
-
-               IF ( ntimes_bdy_bt == 1 ) CALL ctl_stop( &
-                    'There is only one time dump in data files', &
-                    'Set ln_clim=.true. in namelist for constant bdy forcing.' )
-
-               zinterval_s = zstepr(2) - zstepr(1)
-               zinterval_e = zstepr(ntimes_bdy_bt) - zstepr(ntimes_bdy_bt-1)
-
-               IF( zstepr(1) + zoffset > 0 ) THEN
-                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
-                     CALL ctl_stop( 'First time dump in bdy file is after model initial time', ctmp1 )
-               END IF
-               IF( zstepr(ntimes_bdy_bt) + zoffset < totime ) THEN
-                     WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
-                     CALL ctl_stop( 'Last time dump in bdy file is before model final time', ctmp1 )
-               END IF
-            END IF ! .NOT. ln_clim
-
-            IF ( igrd .EQ. 4) THEN
-              ntimes_bdyt = ntimes_bdy_bt
-              zoffsett = zoffset
-              istept(:) = INT( zstepr(:) + zoffset )
-            ELSE IF (igrd .EQ. 5) THEN
-              ntimes_bdyu = ntimes_bdy_bt
-              zoffsetu = zoffset
-              istepu(:) = INT( zstepr(:) + zoffset )
-            ELSE IF (igrd .EQ. 6) THEN
-              ntimes_bdyv = ntimes_bdy_bt
-              zoffsetv = zoffset
-              istepv(:) = INT( zstepr(:) + zoffset )
-            ENDIF
-
-          ENDDO
-
-      ! Verify time consistency between files  
-
-          IF ( ntimes_bdyu /= ntimes_bdyt .OR. ntimes_bdyv /= ntimes_bdyt ) THEN
-             CALL ctl_stop( &
-             'Time axis lengths differ between bdy data files', &
-             'Multiple time frequencies not implemented yet' )
-          ELSE
-            ntimes_bdy_bt = ntimes_bdyt
-          ENDIF
-
-          IF (zoffsetu.NE.zoffsett .OR. zoffsetv.NE.zoffsett) THEN
-            CALL ctl_stop( & 
-            'Bdy data files must have the same time origin', &
-            'Multiple time frequencies not implemented yet'  )
-          ENDIF
-          zoffset = zoffsett
-
-      !! Check that times are the same in the three files... HERE.
-          istep_bt(:) = istept(:)
-
-      ! Check number of time dumps:              
-          IF (ln_clim) THEN
-            SELECT CASE ( ntimes_bdy_bt )
-            CASE( 1 )
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
-              IF(lwp) WRITE(numout,*)             
-            CASE( 12 )
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
-              IF(lwp) WRITE(numout,*) 
-            CASE DEFAULT
-              CALL ctl_stop( &
-                'For climatological boundary forcing (ln_clim=.true.),',&
-                'bdy data files must contain 1 or 12 time dumps.' )
-            END SELECT
-          ENDIF
-
-      ! Find index of first record to read (before first model time). 
-
-          it=1
-          DO WHILE ( ((istep_bt(it+1)) <= 0 ).AND.(it.LE.(ntimes_bdy_bt-1)))
-            it=it+1
-          END DO
-          nbdy_b_bt = it
-
-          IF(lwp) WRITE(numout,*) 'Time offset is ',zoffset
-          IF(lwp) WRITE(numout,*) 'First record to read is ',nbdy_b_bt
-
-        ENDIF ! endif (nn_dtactl == 1)
-
-      ! 1.2  Read first record in file if necessary (ie if nn_dtactl == 1)
-      ! *****************************************************************
-
-        IF ( nn_dtactl == 0) THEN
-          ! boundary data arrays are filled with initial conditions
-          igrd = 5            ! U-points data 
-          DO ib = 1, nblen(igrd)              
-            ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1)
-          END DO
-
-          igrd = 6            ! V-points data 
-          DO ib = 1, nblen(igrd)              
-            vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1)
-          END DO
-
-          igrd = 4            ! T-points data 
-          DO ib = 1, nblen(igrd)              
-            sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd))
-          END DO
-
-        ELSEIF (nn_dtactl == 1) THEN
- 
-        ! Set first record in the climatological case:   
-          IF ((ln_clim).AND.(ntimes_bdy_bt==1)) THEN
-            nbdy_a_bt = 1
-          ELSEIF ((ln_clim).AND.(ntimes_bdy_bt==iman)) THEN
-            nbdy_b_bt = 0
-            nbdy_a_bt = imois
-          ELSE
-            nbdy_a_bt = nbdy_b_bt
-          END IF
- 
-         ! Open Netcdf files:
-
-          CALL iom_open ( cn_dta_fla_T, numbdyt_bt )
-          CALL iom_open ( cn_dta_fla_U, numbdyu_bt )
-          CALL iom_open ( cn_dta_fla_V, numbdyv_bt )
-
-         ! Read first record:
-          ipj=1
-          igrd=4
-          ipi=nblendta(igrd)
-
-          ! ssh
-          igrd=4
-          IF ( nblendta(igrd) .le. 0 ) THEN 
-            idvar = iom_varid( numbdyt_bt,'sossheig' )
-            nblendta(igrd) = iom_file(numbdyt_bt)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for sossheig is ',nblendta(igrd)
-          ipi=nblendta(igrd)
-
-          CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt )
-
-          DO ib=1, nblen(igrd)
-            sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
-          END DO
- 
-          ! u-velocity
-          igrd=5
-          IF ( nblendta(igrd) .le. 0 ) THEN 
-            idvar = iom_varid( numbdyu_bt,'vobtcrtx' )
-            nblendta(igrd) = iom_file(numbdyu_bt)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(igrd)
-          ipi=nblendta(igrd)
-
-          CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt )
-
-          DO ib=1, nblen(igrd)
-            ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
-          END DO
-
-          ! v-velocity
-          igrd=6
-          IF ( nblendta(igrd) .le. 0 ) THEN 
-            idvar = iom_varid( numbdyv_bt,'vobtcrty' )
-            nblendta(igrd) = iom_file(numbdyv_bt)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(igrd)
-          ipi=nblendta(igrd)
-
-          CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt )
-
-          DO ib=1, nblen(igrd)
-            vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
-          END DO
-
-        END IF
- 
-        ! In the case of constant boundary forcing fill bdy arrays once for all
-        IF ((ln_clim).AND.(ntimes_bdy_bt==1)) THEN
-
-          ubtbdy  (:) = ubtbdydta  (:,2)
-          vbtbdy  (:) = vbtbdydta  (:,2)
-          sshbdy  (:) = sshbdydta  (:,2)
-
-          CALL iom_close( numbdyt_bt )
-          CALL iom_close( numbdyu_bt )
-          CALL iom_close( numbdyv_bt )
-
-        END IF
-
-      ENDIF ! End if nit000
-
-      ! -------------------- !
-      ! 2. At each time step !
-      ! -------------------- !
-
-      IF ((nn_dtactl==1).AND.(ntimes_bdy_bt>1)) THEN 
-
-      ! 2.1 Read one more record if necessary
-      !**************************************
-
-        IF ( (ln_clim).AND.(imois/=nbdy_b_bt) ) THEN ! remember that nbdy_b_bt=0 for kt=nit000
-         nbdy_b_bt = imois
-         nbdy_a_bt = imois+1
-         nbdy_b_bt = MOD( nbdy_b_bt, iman )
-         IF( nbdy_b_bt == 0 ) nbdy_b_bt = iman
-         nbdy_a_bt = MOD( nbdy_a_bt, iman )
-         IF( nbdy_a_bt == 0 ) nbdy_a_bt = iman
-         lect=.true.
-
-        ELSEIF ((.NOT.ln_clim).AND.(itimer >= istep_bt(nbdy_a_bt))) THEN
-          nbdy_b_bt=nbdy_a_bt
-          nbdy_a_bt=nbdy_a_bt+1
-          lect=.true.
-        END IF
-         
-        IF (lect) THEN
-
-        ! Swap arrays
-          sshbdydta(:,1) =  sshbdydta(:,2)
-          ubtbdydta(:,1) =  ubtbdydta(:,2)
-          vbtbdydta(:,1) =  vbtbdydta(:,2)
- 
-        ! read another set
-
-          ipj=1
-          ipk=jpk
-          igrd=4
-          ipi=nblendta(igrd)
-
-          
-          ! ssh
-          igrd=4
-          ipi=nblendta(igrd)
-
-          CALL iom_get ( numbdyt_bt, jpdom_unknown,'sossheig',zdta(1:ipi,1:ipj),nbdy_a_bt )
-
-          DO ib=1, nblen(igrd)
-            sshbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
-          END DO
-
-          ! u-velocity
-          igrd=5
-          ipi=nblendta(igrd)
-
-          CALL iom_get ( numbdyu_bt, jpdom_unknown,'vobtcrtx',zdta(1:ipi,1:ipj),nbdy_a_bt )
-
-          DO ib=1, nblen(igrd)
-            ubtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
-          END DO
-
-          ! v-velocity
-          igrd=6
-          ipi=nblendta(igrd)
-
-          CALL iom_get ( numbdyv_bt, jpdom_unknown,'vobtcrty',zdta(1:ipi,1:ipj),nbdy_a_bt )
-
-          DO ib=1, nblen(igrd)
-            vbtbdydta(ib,2) =  zdta(nbmap(ib,igrd),1)
-          END DO
-
-
-         IF(lwp) WRITE(numout,*) 'bdy_dta_fla : first record file used nbdy_b_bt ',nbdy_b_bt
-         IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a_bt ',nbdy_a_bt
-         IF (.NOT.ln_clim) THEN
-           IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy_b_bt)
-           IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer
-           IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt)
-         ENDIF
-        END IF ! end lect=.true.
-
-
-      ! 2.2   Interpolate linearly:
-      ! ***************************
-    
-        IF (ln_clim) THEN
-          zxy = REAL( nday, wp ) / REAL( nmonth_len(nbdy_b_bt), wp ) + 0.5 - i15
-        ELSE          
-          zxy = REAL(istep_bt(nbdy_b_bt)-itimer, wp) / REAL(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt), wp)
-        END IF
-
-          igrd=4
-          DO ib=1, nblen(igrd)
-            sshbdy(ib) = zxy      * sshbdydta(ib,2) + &
-                       (1.-zxy) * sshbdydta(ib,1)   
-          END DO
-
-          igrd=5
-          DO ib=1, nblen(igrd)
-            ubtbdy(ib) = zxy      * ubtbdydta(ib,2) + &
-                         (1.-zxy) * ubtbdydta(ib,1)   
-          END DO
-
-          igrd=6
-          DO ib=1, nblen(igrd)
-            vbtbdy(ib) = zxy      * vbtbdydta(ib,2) + &
-                         (1.-zxy) * vbtbdydta(ib,1)   
-          END DO
-
-
-      END IF !end if ((nn_dtactl==1).AND.(ntimes_bdy_bt>1))
-    
-      ! ------------------- !
-      ! Last call kt=nitend !
-      ! ------------------- !
-
-      ! Closing of the 3 files
-      IF( kt == nitend   .and. jit == icycl ) THEN
-          CALL iom_close( numbdyt_bt )
-          CALL iom_close( numbdyu_bt )
-          CALL iom_close( numbdyv_bt )
-      ENDIF
-
-      ENDIF ! ln_dyn_frs
-
-      END SUBROUTINE bdy_dta_fla
-
+               IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN
+                  jfld = jfld + 1
+                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)
+               ENDIF
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%u2d => bf(jfld)%fnow(:,1,1)
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%v2d => bf(jfld)%fnow(:,1,1)
+            ENDIF
+         ENDIF
+
+         IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN
+            IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN
+               ilen0(1:3) = nblen(1:3)
+            ELSE
+               ilen0(1:3) = nblenrim(1:3)
+            ENDIF
+            ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) )
+            ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) )
+         ENDIF
+         IF ( ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 1 ).or. &
+           &  ( ln_full_vel_array(ib_bdy) .and. nn_dyn2d(ib_bdy) .gt. 0 .and.   &
+           &    ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) ) THEN
+            jfld = jfld + 1
+            dta_bdy(ib_bdy)%u3d => bf(jfld)%fnow(:,1,:)
+            jfld = jfld + 1
+            dta_bdy(ib_bdy)%v3d => bf(jfld)%fnow(:,1,:)
+         ENDIF
+
+         IF (nn_tra(ib_bdy) .gt. 0) THEN
+            IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN
+               IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  ilen0(1:3) = nblenrim(1:3)
+               ENDIF
+               ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) )
+               ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) )
+            ELSE
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%tem => bf(jfld)%fnow(:,1,:)
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%sal => bf(jfld)%fnow(:,1,:)
+            ENDIF
+         ENDIF
+
+#if defined key_lim2
+         IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN
+            IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN
+               IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN
+                  ilen0(1:3) = nblen(1:3)
+               ELSE
+                  ilen0(1:3) = nblenrim(1:3)
+               ENDIF
+               ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) )
+               ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) )
+               ALLOCATE( dta_bdy(ib_bdy)%hsnif(ilen0(1)) )
+            ELSE
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%frld  => bf(jfld)%fnow(:,1,1)
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%hicif => bf(jfld)%fnow(:,1,1)
+               jfld = jfld + 1
+               dta_bdy(ib_bdy)%hsnif => bf(jfld)%fnow(:,1,1)
+            ENDIF
+         ENDIF
+#endif
+
+      ENDDO ! ib_bdy 
+
+      END SUBROUTINE bdy_dta_init
 
 #else
    !!----------------------------------------------------------------------
-   !!   Dummy module                   NO Unstruct Open Boundary Conditions
+   !!   Dummy module                   NO Open Boundary Conditions
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE bdy_dta_frs( kt )              ! Empty routine
-      WRITE(*,*) 'bdy_dta_frs: You should not have seen this print! error?', kt
-   END SUBROUTINE bdy_dta_frs
-   SUBROUTINE bdy_dta_fla( kt, kit, icycle )      ! Empty routine
-      WRITE(*,*) 'bdy_dta_frs: You should not have seen this print! error?', kt, kit
-   END SUBROUTINE bdy_dta_fla
+   SUBROUTINE bdy_dta( kt, jit, time_offset ) ! Empty routine
+      INTEGER, INTENT( in )           ::   kt    
+      INTEGER, INTENT( in ), OPTIONAL ::   jit   
+      INTEGER, INTENT( in ), OPTIONAL ::   time_offset
+      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_dta
+   SUBROUTINE bdy_dta_init()                  ! Empty routine
+      WRITE(*,*) 'bdy_dta_init: You should not have seen this print! error?'
+   END SUBROUTINE bdy_dta_init
 #endif
 
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 3116)
@@ -15,14 +15,14 @@
    !!   'key_bdy' :                    Unstructured Open Boundary Condition
    !!----------------------------------------------------------------------
-   !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary
-   !!   bdy_dyn_fla    : Flather condition for barotropic solution
+   !!   bdy_dyn3d        : apply open boundary conditions to baroclinic velocities
+   !!   bdy_dyn3d_frs    : apply Flow Relaxation Scheme
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers 
    USE dom_oce         ! ocean space and time domain
+   USE dynspg_oce      
    USE bdy_oce         ! ocean open boundary conditions
-   USE dynspg_oce      ! for barotropic variables
-   USE phycst          ! physical constants
+   USE bdydyn2d        ! open boundary conditions for barotropic solution
+   USE bdydyn3d        ! open boundary conditions for baroclinic velocities
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE bdytides        ! for tidal harmonic forcing at boundary
    USE in_out_manager  !
 
@@ -30,9 +30,8 @@
    PRIVATE
 
-   PUBLIC   bdy_dyn_frs   ! routine called in dynspg_flt (free surface case ONLY)
-# if defined key_dynspg_exp || defined key_dynspg_ts
-   PUBLIC   bdy_dyn_fla   ! routine called in dynspg_flt (free surface case ONLY)
-# endif
+   PUBLIC   bdy_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)
@@ -42,141 +41,80 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn_frs( kt )
+   SUBROUTINE bdy_dyn( kt, dyn3d_only )
       !!----------------------------------------------------------------------
-      !!                  ***  SUBROUTINE bdy_dyn_frs  ***
+      !!                  ***  SUBROUTINE bdy_dyn  ***
       !!
-      !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the  
-      !!                case of unstructured open boundaries.
+      !! ** Purpose : - Wrapper routine for bdy_dyn2d and bdy_dyn3d.
       !!
-      !! 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   ! Main time step counter
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: wrk_2d_7, wrk_2d_8      ! 2D workspace
       !!
-      INTEGER  ::   jb, jk         ! dummy loop indices
-      INTEGER  ::   ii, ij, igrd   ! local integers
-      REAL(wp) ::   zwgt           ! boundary weight
-      !!----------------------------------------------------------------------
-      !
-      IF(ln_dyn_frs) THEN       ! If this is false, then this routine does nothing. 
-         !
-         IF( kt == nit000 ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'bdy_dyn_frs : Flow Relaxation Scheme on momentum'
-            IF(lwp) WRITE(numout,*) '~~~~~~~'
-         ENDIF
-         !
-         igrd = 2                      ! Relaxation of zonal velocity
-         DO jb = 1, nblen(igrd)
-            DO jk = 1, jpkm1
-               ii   = nbi(jb,igrd)
-               ij   = nbj(jb,igrd)
-               zwgt = nbw(jb,igrd)
-               ua(ii,ij,jk) = ( ua(ii,ij,jk) * ( 1.- zwgt ) + ubdy(jb,jk) * zwgt ) * umask(ii,ij,jk)
-            END DO
-         END DO
-         !
-         igrd = 3                      ! Relaxation of meridional velocity
-         DO jb = 1, nblen(igrd)
-            DO jk = 1, jpkm1
-               ii   = nbi(jb,igrd)
-               ij   = nbj(jb,igrd)
-               zwgt = nbw(jb,igrd)
-               va(ii,ij,jk) = ( va(ii,ij,jk) * ( 1.- zwgt ) + vbdy(jb,jk) * zwgt ) * vmask(ii,ij,jk)
-            END DO
-         END DO 
-         CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
-         !
-      ENDIF ! ln_dyn_frs
-      !
-   END SUBROUTINE bdy_dyn_frs
+      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  
+      !!
 
+      IF(wrk_in_use(2, 7,8) ) THEN
+         CALL ctl_stop('bdy_dyn: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
+      END IF
 
-# if defined   key_dynspg_exp   ||   defined key_dynspg_ts
-   !!----------------------------------------------------------------------
-   !!   'key_dynspg_exp'        OR              explicit sea surface height
-   !!   'key_dynspg_ts '                  split-explicit sea surface height
-   !!----------------------------------------------------------------------
-   
-!! Option to use Flather with dynspg_flt not coded yet...
+      ll_dyn2d = .true.
+      ll_dyn3d = .true.
 
-   SUBROUTINE bdy_dyn_fla( pssh )
-      !!----------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdy_dyn_fla  ***
-      !!             
-      !!              - Apply Flather boundary conditions on normal barotropic velocities 
-      !!                (ln_dyn_fla=.true. or ln_tides=.true.)
-      !!
-      !! ** WARNINGS about FLATHER implementation:
-      !!1. According to Palma and Matano, 1998 "after ssh" is used. 
-      !!   In ROMS and POM implementations, it is "now ssh". In the current 
-      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable. 
-      !!   So I use "before ssh" in the following.
-      !!
-      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of 
-      !!   fact, the model ssh just inside the dynamical boundary is used (the outside  
-      !!   ssh in the code is not updated).
-      !!
-      !! References:  Flather, R. A., 1976: A tidal model of the northwest European
-      !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh
+      IF( PRESENT(dyn3d_only) ) THEN
+         IF( dyn3d_only ) ll_dyn2d = .false.
+      ENDIF
 
-      INTEGER  ::   jb, igrd                         ! dummy loop indices
-      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
-      REAL(wp) ::   zcorr                            ! Flather correction
-      REAL(wp) ::   zforc                            ! temporary scalar
-      !!----------------------------------------------------------------------
+      !-------------------------------------------------------
+      ! Set pointers
+      !-------------------------------------------------------
 
-      ! ---------------------------------!
-      ! Flather boundary conditions     :!
-      ! ---------------------------------! 
-     
-      IF(ln_dyn_fla .OR. ln_tides) THEN ! If these are both false, then this routine does nothing. 
+      pssh => sshn
+      phur => hur
+      phvr => hvr
+      pu2d => wrk_2d_7
+      pv2d => wrk_2d_8
 
-         ! Fill temporary array with ssh data (here spgu):
-         igrd = 4
-         spgu(:,:) = 0.0
-         DO jb = 1, nblenrim(igrd)
-            ii = nbi(jb,igrd)
-            ij = nbj(jb,igrd)
-            IF( ln_dyn_fla ) spgu(ii, ij) = sshbdy(jb)
-            IF( ln_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb)
-         END DO
-         !
-         igrd = 5      ! Flather bc on u-velocity; 
-         !             ! remember that flagu=-1 if normal velocity direction is outward
-         !             ! I think we should rather use after ssh ?
-         DO jb = 1, nblenrim(igrd)
-            ii  = nbi(jb,igrd)
-            ij  = nbj(jb,igrd) 
-            iim1 = ii + MAX( 0, INT( flagu(jb) ) )   ! T pts i-indice inside the boundary
-            iip1 = ii - MIN( 0, INT( flagu(jb) ) )   ! T pts i-indice outside the boundary 
-            !
-            zcorr = - flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
-            zforc = ubtbdy(jb) + utide(jb)
-            ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
-         END DO
-         !
-         igrd = 6      ! Flather bc on v-velocity
-         !             ! remember that flagv=-1 if normal velocity direction is outward
-         DO jb = 1, nblenrim(igrd)
-            ii  = nbi(jb,igrd)
-            ij  = nbj(jb,igrd) 
-            ijm1 = ij + MAX( 0, INT( flagv(jb) ) )   ! T pts j-indice inside the boundary
-            ijp1 = ij - MIN( 0, INT( flagv(jb) ) )   ! T pts j-indice outside the boundary 
-            !
-            zcorr = - flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
-            zforc = vbtbdy(jb) + vtide(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. )   !
-         !
-      ENDIF ! ln_dyn_fla .or. ln_tides
-      !
-   END SUBROUTINE bdy_dyn_fla
-#endif
+      !-------------------------------------------------------
+      ! 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
+      !-------------------------------------------------------
+
+      IF( ll_dyn2d ) CALL bdy_dyn2d( kt )
+
+      IF( ll_dyn3d ) CALL bdy_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('bdy_dyn: ERROR: failed to release workspace arrays.')
+
+   END SUBROUTINE bdy_dyn
 
 #else
@@ -185,11 +123,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE bdy_dyn_frs( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
-   END SUBROUTINE bdy_dyn_frs
-   SUBROUTINE bdy_dyn_fla( pssh )    ! Empty routine
-      REAL :: pssh(:,:)
-      WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?', pssh(1,1)
-   END SUBROUTINE bdy_dyn_fla
+   SUBROUTINE bdy_dyn( kt )      ! Empty routine
+      WRITE(*,*) 'bdy_dyn: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_dyn
 #endif
 
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3116)
@@ -10,4 +10,6 @@
    !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
    !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
+   !!            3.4  !  2011     (D. Storkey, J. Chanut) OBC-BDY merge
+   !!                 !  --- Renamed bdyini.F90 -> bdyini.F90 ---
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -19,8 +21,5 @@
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain
-   USE obc_par         ! ocean open boundary conditions
    USE bdy_oce         ! unstructured open boundary conditions
-   USE bdydta, ONLY: bdy_dta_alloc ! open boundary data
-   USE bdytides        ! tides at open boundaries initialization (tide_init routine)
    USE in_out_manager  ! I/O units
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -52,104 +51,581 @@
       !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
       !!----------------------------------------------------------------------      
-      INTEGER  ::   ii, ij, ik, igrd, ib, ir   ! dummy loop indices
-      INTEGER  ::   icount, icountr, ib_len, ibr_max   ! local integers
-      INTEGER  ::   iw, ie, is, in, inum, id_dummy     !   -       -
-      INTEGER  ::   igrd_start, igrd_end               !   -       -
-      REAL(wp) ::   zefl, zwfl, znfl, zsfl              ! local scalars
-      INTEGER, DIMENSION (2)             ::   kdimsz
-      INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
-      INTEGER, DIMENSION(jpbdta, jpbgrd) ::   nbrdta           ! Discrete distance from rim points
-      REAL(wp), DIMENSION(jpidta,jpjdta) ::   zmask            ! global domain mask
-      REAL(wp), DIMENSION(jpbdta,1)      ::   zdta             ! temporary array 
-      CHARACTER(LEN=80),DIMENSION(6)     ::   clfile
+      ! namelist variables
+      !-------------------
+      INTEGER, PARAMETER          :: jp_nseg = 100
+      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
+      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft
+      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft
+      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft
+      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft
+
+      ! local variables
+      !-------------------
+      INTEGER  ::   ib_bdy, 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         !   -       -
+      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
+      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
+      REAL   , POINTER  ::  flagu, flagv                   !    -   -
+      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
+      INTEGER, DIMENSION (2)                ::   kdimsz
+      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays 
+      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
+      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
+      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask
+      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
+      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid
       !!
-      NAMELIST/nambdy/cn_mask, cn_dta_frs_T, cn_dta_frs_U, cn_dta_frs_V,   &
-         &            cn_dta_fla_T, cn_dta_fla_U, cn_dta_fla_V,            &
-         &            ln_tides, ln_clim, ln_vol, ln_mask,                  &
-         &            ln_dyn_fla, ln_dyn_frs, ln_tra_frs,ln_ice_frs,       &
-         &            nn_dtactl, nn_rimwidth, nn_volctl
+      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
+         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
+         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         &  
+#if defined key_lim2
+         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
+#endif
+         &             ln_vol, nn_volctl,                                  &
+         &             nn_rimwidth, nn_dmp2d_in, nn_dmp2d_out,             &
+         &             nn_dmp3d_in, nn_dmp3d_out
+      !!
+      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             &
+                             nbdysegw, jpiwob, jpjwdt, jpjwft,             &
+                             nbdysegn, jpjnob, jpindt, jpinft,             &
+                             nbdysegs, jpjsob, jpisdt, jpisft
+
       !!----------------------------------------------------------------------
 
+      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
+
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries'
+      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
       IF(lwp) WRITE(numout,*) '~~~~~~~~'
       !
-      !                                      ! allocate bdy_oce arrays
-      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
-      IF( bdy_dta_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate dta arrays' )
 
       IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
-         &                               ' and unstructured open boundary condition are not compatible' )
-
-      IF( lk_obc      )   CALL ctl_stop( 'Straight open boundaries,',   &
-         &                               ' and unstructured open boundaries are not compatible' )
-
-      ! ---------------------------
-      REWIND( numnam )                    ! Read namelist parameters
+         &                               ' and general open boundary condition are not compatible' )
+
+      cgrid= (/'t','u','v'/)
+
+      ! -----------------------------------------
+      ! Initialise and read namelist parameters
+      ! -----------------------------------------
+
+      nb_bdy            = 0
+      ln_coords_file(:) = .false.
+      cn_coords_file(:) = ''
+      ln_mask_file      = .false.
+      cn_mask_file(:)   = ''
+      nn_dyn2d(:)       = 0
+      nn_dyn2d_dta(:)   = -1  ! uninitialised flag
+      nn_dyn3d(:)       = 0
+      nn_dyn3d_dta(:)   = -1  ! uninitialised flag
+      nn_tra(:)         = 0
+      nn_tra_dta(:)     = -1  ! uninitialised flag
+#if defined key_lim2
+      nn_ice_lim2(:)    = 0
+      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
+#endif
+      ln_vol            = .false.
+      nn_volctl         = -1  ! uninitialised flag
+      nn_rimwidth(:)    = -1  ! uninitialised flag
+      nn_dmp2d_in(:)    = -1  ! uninitialised flag
+      nn_dmp2d_out(:)   = -1  ! uninitialised flag
+      nn_dmp3d_in(:)    = -1  ! uninitialised flag
+      nn_dmp3d_out(:)   = -1  ! uninitialised flag
+
+      REWIND( numnam )                    
       READ  ( numnam, nambdy )
+
+      ! -----------------------------------------
+      ! Check and write out namelist parameters
+      ! -----------------------------------------
 
       !                                   ! control prints
       IF(lwp) WRITE(numout,*) '         nambdy'
 
-      !                                         ! check type of data used (nn_dtactl value)
-      IF(lwp) WRITE(numout,*) 'nn_dtactl =', nn_dtactl      
-      IF(lwp) WRITE(numout,*)
-      SELECT CASE( nn_dtactl )                   ! 
-      CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
-      CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
-      CASE DEFAULT   ;   CALL ctl_stop( 'nn_dtactl must be 0 or 1' )
-      END SELECT
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nn_rimwidth = ', nn_rimwidth
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) '      nn_volctl = ', nn_volctl
-
-      IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
-         SELECT CASE ( nn_volctl )
+      IF( nb_bdy .eq. 0 ) THEN 
+        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
+      ELSE
+        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
+      ENDIF
+
+      DO ib_bdy = 1,nb_bdy
+        IF(lwp) WRITE(numout,*) ' ' 
+        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
+
+        IF( ln_coords_file(ib_bdy) ) THEN
+           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
+        ELSE
+           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+
+        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
+        SELECT CASE( nn_dyn2d(ib_bdy) )                  
+          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
+        END SELECT
+        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
+           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   ! 
+              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
+              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
+              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
+        IF(lwp) WRITE(numout,*)
+
+        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
+        SELECT CASE( nn_dyn3d(ib_bdy) )                  
+          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
+        END SELECT
+        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
+           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   ! 
+              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
+              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
+              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
+           END SELECT
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+
+        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
+        SELECT CASE( nn_tra(ib_bdy) )                  
+          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
+        END SELECT
+        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
+           SELECT CASE( nn_tra_dta(ib_bdy) )                   ! 
+              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
+              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
+              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
+           END SELECT
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+
+#if defined key_lim2
+        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
+        SELECT CASE( nn_ice_lim2(ib_bdy) )                  
+          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
+          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
+        END SELECT
+        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN 
+           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   ! 
+              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
+              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
+              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
+           END SELECT
+        ENDIF
+        IF(lwp) WRITE(numout,*)
+#endif
+
+        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)
+        IF(lwp) WRITE(numout,*)
+
+      ENDDO
+
+     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
+       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
+       IF(lwp) WRITE(numout,*)
+       SELECT CASE ( nn_volctl )
          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
          CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
-         END SELECT
-         IF(lwp) WRITE(numout,*)
-      ELSE
-         IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries'
-         IF(lwp) WRITE(numout,*)
-      ENDIF
-
-      IF( ln_tides ) THEN
-        IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*)
-      ENDIF
-
-      IF( ln_dyn_fla ) THEN
-        IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*)
-      ENDIF
-
-      IF( ln_dyn_frs ) THEN
-        IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*)
-      ENDIF
-
-      IF( ln_tra_frs ) THEN
-        IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*)
-      ENDIF
-
-      IF( ln_ice_frs ) THEN
-        IF(lwp) WRITE(numout,*) 'FRS condition on ice fields at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*)
-      ENDIF
-
-      IF( ln_tides )   CALL tide_init      ! Read tides namelist 
-
-
-      ! Read arrays defining unstructured open boundaries
+       END SELECT
+       IF(lwp) WRITE(numout,*)
+     ELSE
+       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
+       IF(lwp) WRITE(numout,*)
+     ENDIF
+
       ! -------------------------------------------------
+      ! Initialise indices arrays for open boundaries
+      ! -------------------------------------------------
+
+      ! Work out global dimensions of boundary data
+      ! ---------------------------------------------
+      REWIND( numnam )                    
+      DO ib_bdy = 1, nb_bdy
+
+         jpbdta = 1
+         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
+ 
+            ! No REWIND here because may need to read more than one nambdy_index namelist.
+            READ  ( numnam, nambdy_index )
+
+            ! Automatic boundary definition: if nbdysegX = -1
+            ! set boundary to whole side of model domain.
+            IF( nbdysege == -1 ) THEN
+               nbdysege = 1
+               jpieob(1) = jpiglo - 1
+               jpjedt(1) = 2
+               jpjeft(1) = jpjglo - 1
+            ENDIF
+            IF( nbdysegw == -1 ) THEN
+               nbdysegw = 1
+               jpiwob(1) = 2
+               jpjwdt(1) = 2
+               jpjwft(1) = jpjglo - 1
+            ENDIF
+            IF( nbdysegn == -1 ) THEN
+               nbdysegn = 1
+               jpjnob(1) = jpjglo - 1
+               jpindt(1) = 2
+               jpinft(1) = jpiglo - 1
+            ENDIF
+            IF( nbdysegs == -1 ) THEN
+               nbdysegs = 1
+               jpjsob(1) = 2
+               jpisdt(1) = 2
+               jpisft(1) = jpiglo - 1
+            ENDIF
+
+            nblendta(:,ib_bdy) = 0
+            DO iseg = 1, nbdysege
+               igrd = 1
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
+               igrd = 3
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)               
+            ENDDO
+            DO iseg = 1, nbdysegw
+               igrd = 1
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
+               igrd = 3
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)               
+            ENDDO
+            DO iseg = 1, nbdysegn
+               igrd = 1
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)               
+               igrd = 3
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
+            ENDDO
+            DO iseg = 1, nbdysegs
+               igrd = 1
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
+               igrd = 2
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)
+               igrd = 3
+               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
+            ENDDO
+
+            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
+            jpbdta = MAXVAL(nblendta(:,ib_bdy))               
+
+
+         ELSE            ! Read size of arrays in boundary coordinates file.
+
+
+            CALL iom_open( cn_coords_file(ib_bdy), inum )
+            jpbdta = 1
+            DO igrd = 1, jpbgrd
+               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz )  
+               nblendta(igrd,ib_bdy) = kdimsz(1)
+               jpbdta = MAX(jpbdta, kdimsz(1))
+            ENDDO
+
+         ENDIF 
+
+      ENDDO ! ib_bdy
+
+      ! Allocate arrays
+      !---------------
+      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
+         &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
+
+      ALLOCATE( dta_global(jpbdta, 1, jpk) )
+
+      ! Calculate global boundary index arrays or read in from file
+      !------------------------------------------------------------
+      REWIND( numnam )                    
+      DO ib_bdy = 1, nb_bdy
+
+         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters
+
+            ! No REWIND here because may need to read more than one nambdy_index namelist.
+            READ  ( numnam, nambdy_index )
+
+            ! Automatic boundary definition: if nbdysegX = -1
+            ! set boundary to whole side of model domain.
+            IF( nbdysege == -1 ) THEN
+               nbdysege = 1
+               jpieob(1) = jpiglo - 1
+               jpjedt(1) = 2
+               jpjeft(1) = jpjglo - 1
+            ENDIF
+            IF( nbdysegw == -1 ) THEN
+               nbdysegw = 1
+               jpiwob(1) = 2
+               jpjwdt(1) = 2
+               jpjwft(1) = jpjglo - 1
+            ENDIF
+            IF( nbdysegn == -1 ) THEN
+               nbdysegn = 1
+               jpjnob(1) = jpjglo - 1
+               jpindt(1) = 2
+               jpinft(1) = jpiglo - 1
+            ENDIF
+            IF( nbdysegs == -1 ) THEN
+               nbdysegs = 1
+               jpjsob(1) = 2
+               jpisdt(1) = 2
+               jpisft(1) = jpiglo - 1
+            ENDIF
+
+            ! ------------ T points -------------
+            igrd = 1  
+            icount = 0
+            DO ir = 1, nn_rimwidth(ib_bdy)
+               ! east
+               DO iseg = 1, nbdysege
+                  DO ij = jpjedt(iseg), jpjeft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
+                     nbjdta(icount, igrd, ib_bdy) = ij
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! west
+               DO iseg = 1, nbdysegw
+                  DO ij = jpjwdt(iseg), jpjwft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+                     nbjdta(icount, igrd, ib_bdy) = ij
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! north
+               DO iseg = 1, nbdysegn
+                  DO ii = jpindt(iseg), jpinft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = ii
+                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! south
+               DO iseg = 1, nbdysegs
+                  DO ii = jpisdt(iseg), jpisft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = ii
+                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            ! ------------ U points -------------
+            igrd = 2  
+            icount = 0
+            DO ir = 1, nn_rimwidth(ib_bdy)
+               ! east
+               DO iseg = 1, nbdysege
+                  DO ij = jpjedt(iseg), jpjeft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
+                     nbjdta(icount, igrd, ib_bdy) = ij
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! west
+               DO iseg = 1, nbdysegw
+                  DO ij = jpjwdt(iseg), jpjwft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+                     nbjdta(icount, igrd, ib_bdy) = ij
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! north
+               DO iseg = 1, nbdysegn
+                  DO ii = jpindt(iseg), jpinft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = ii
+                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! south
+               DO iseg = 1, nbdysegs
+                  DO ii = jpisdt(iseg), jpisft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = ii
+                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+            ENDDO
+
+            ! ------------ V points -------------
+            igrd = 3  
+            icount = 0
+            DO ir = 1, nn_rimwidth(ib_bdy)
+               ! east
+               DO iseg = 1, nbdysege
+                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
+                     nbjdta(icount, igrd, ib_bdy) = ij
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! west
+               DO iseg = 1, nbdysegw
+                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
+                     nbjdta(icount, igrd, ib_bdy) = ij
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! north
+               DO iseg = 1, nbdysegn
+                  DO ii = jpindt(iseg), jpinft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = ii
+                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+               ! south
+               DO iseg = 1, nbdysegs
+                  DO ii = jpisdt(iseg), jpisft(iseg)
+                     icount = icount + 1
+                     nbidta(icount, igrd, ib_bdy) = ii
+                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir + 1
+                     nbrdta(icount, igrd, ib_bdy) = ir
+                  ENDDO
+               ENDDO
+            ENDDO
+
+         ELSE            ! Read global index arrays from boundary coordinates file.
+
+            DO igrd = 1, jpbgrd
+               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+               DO ii = 1,nblendta(igrd,ib_bdy)
+                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+               END DO
+               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+               DO ii = 1,nblendta(igrd,ib_bdy)
+                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+               END DO
+               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
+               DO ii = 1,nblendta(igrd,ib_bdy)
+                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
+               END DO
+
+               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
+               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
+               IF (ibr_max < nn_rimwidth(ib_bdy))   &
+                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
+
+            END DO
+            CALL iom_close( inum )
+
+         ENDIF 
+
+      ENDDO 
+
+      ! Work out dimensions of boundary data on each processor
+      ! ------------------------------------------------------
+     
+      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
+      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
+      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
+      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
+
+      DO ib_bdy = 1, nb_bdy
+         DO igrd = 1, jpbgrd
+            icount  = 0
+            icountr = 0
+            idx_bdy(ib_bdy)%nblen(igrd)    = 0
+            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
+            DO ib = 1, nblendta(igrd,ib_bdy)
+               ! check that data is in correct order in file
+               ibm1 = MAX(1,ib-1)
+               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
+                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
+                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
+                                   'A utility for re-ordering boundary coordinates and data files exists in CDFTOOLS')
+                  ENDIF    
+               ENDIF
+               ! check if point is in local domain
+               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
+                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
+                  !
+                  icount = icount  + 1
+                  !
+                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
+               ENDIF
+            ENDDO
+            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
+            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc        
+         ENDDO  ! igrd
+
+         ! Allocate index arrays for this boundary set
+         !--------------------------------------------
+         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
+         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
+         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
+         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
+
+         ! Dispatch mapping indices and discrete distances on each processor
+         ! -----------------------------------------------------------------
+
+         DO igrd = 1, jpbgrd
+            icount  = 0
+            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
+            DO ir=1, nn_rimwidth(ib_bdy)
+               DO ib = 1, nblendta(igrd,ib_bdy)
+                  ! check if point is in local domain and equals ir
+                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
+                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
+                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
+                     !
+                     icount = icount  + 1
+                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
+                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
+                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
+                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDDO 
+
+         ! Compute rim weights for FRS scheme
+         ! ----------------------------------
+         DO igrd = 1, jpbgrd
+            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
+               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
+               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
+!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic
+!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
+            END DO
+         END DO 
+
+      ENDDO
+
+      ! ------------------------------------------------------
+      ! Initialise masks and find normal/tangential directions
+      ! ------------------------------------------------------
 
       ! Read global 2D mask at T-points: bdytmask
-      ! *****************************************
+      ! -----------------------------------------
       ! bdytmask = 1  on the computational domain AND on open boundaries
       !          = 0  elsewhere   
@@ -158,6 +634,6 @@
          zmask(         :                ,:) = 0.e0
          zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0          
-      ELSE IF( ln_mask ) THEN
-         CALL iom_open( cn_mask, inum )
+      ELSE IF( ln_mask_file ) THEN
+         CALL iom_open( cn_mask_file, inum )
          CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
          CALL iom_close( inum )
@@ -184,154 +660,4 @@
 
 
-      ! Read discrete distance and mapping indices
-      ! ******************************************
-      nbidta(:,:) = 0.e0
-      nbjdta(:,:) = 0.e0
-      nbrdta(:,:) = 0.e0
-
-      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
-         icount = 0
-         DO ir = 1, nn_rimwidth                  ! Define west boundary (from ii=2 to ii=1+nn_rimwidth):
-            DO ij = 3, jpjglo-2
-               icount = icount + 1
-               nbidta(icount,:) = ir + 1 + (jpizoom-1)
-               nbjdta(icount,:) = ij     + (jpjzoom-1) 
-               nbrdta(icount,:) = ir
-            END DO
-         END DO
-         !
-         DO ir = 1, nn_rimwidth                  ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nn_rimwidth):
-            DO ij=3,jpjglo-2
-               icount = icount + 1
-               nbidta(icount,:) = jpiglo-ir + (jpizoom-1)
-               nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points
-               nbjdta(icount,:) = ij + (jpjzoom-1)
-               nbrdta(icount,:) = ir
-            END DO
-         END DO
-         !       
-      ELSE            ! Read indices and distances in unstructured boundary data files 
-         !
-         IF( ln_tides ) THEN             ! Read tides input files for preference in case there are no bdydata files
-            clfile(4) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc'
-            clfile(5) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc'
-            clfile(6) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc'
-         ENDIF
-         IF( ln_dyn_fla .AND. .NOT. ln_tides ) THEN 
-            clfile(4) = cn_dta_fla_T
-            clfile(5) = cn_dta_fla_U
-            clfile(6) = cn_dta_fla_V
-         ENDIF
-
-         IF( ln_tra_frs ) THEN 
-            clfile(1) = cn_dta_frs_T
-            IF( .NOT. ln_dyn_frs ) THEN 
-               clfile(2) = cn_dta_frs_T     ! Dummy read re read T file for sake of 6 files
-               clfile(3) = cn_dta_frs_T     !
-            ENDIF
-         ENDIF          
-         IF( ln_dyn_frs ) THEN 
-            IF( .NOT. ln_tra_frs )   clfile(1) = cn_dta_frs_U      ! Dummy Read 
-            clfile(2) = cn_dta_frs_U
-            clfile(3) = cn_dta_frs_V 
-         ENDIF
-
-         !                                   ! how many files are we to read in?
-         IF(ln_tides .OR. ln_dyn_fla)   igrd_start = 4
-         !
-         IF(ln_tra_frs    ) THEN   ;   igrd_start = 1
-         ELSEIF(ln_dyn_frs) THEN   ;   igrd_start = 2
-         ENDIF
-         !
-         IF( ln_tra_frs   )   igrd_end = 1
-         !
-         IF(ln_dyn_fla .OR. ln_tides) THEN   ;   igrd_end = 6
-         ELSEIF( ln_dyn_frs             ) THEN   ;   igrd_end = 3
-         ENDIF
-
-         DO igrd = igrd_start, igrd_end
-            CALL iom_open( clfile(igrd), inum )
-            id_dummy = iom_varid( inum, 'nbidta', kdimsz=kdimsz )  
-            IF(lwp) WRITE(numout,*) 'kdimsz : ',kdimsz
-            ib_len = kdimsz(1)
-            IF( ib_len > jpbdta)   CALL ctl_stop(  'Boundary data array in file too long.',                  &
-                &                                  'File :', TRIM(clfile(igrd)),'increase parameter jpbdta.' )
-
-            CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) )
-            DO ii = 1,ib_len
-               nbidta(ii,igrd) = INT( zdta(ii,1) )
-            END DO
-            CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) )
-            DO ii = 1,ib_len
-               nbjdta(ii,igrd) = INT( zdta(ii,1) )
-            END DO
-            CALL iom_get( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) )
-            DO ii = 1,ib_len
-               nbrdta(ii,igrd) = INT( zdta(ii,1) )
-            END DO
-            CALL iom_close( inum )
-
-            IF( igrd < 4) THEN            ! Check that rimwidth in file is big enough for Frs case(barotropic is one):
-               ibr_max = MAXVAL( nbrdta(:,igrd) )
-               IF(lwp) WRITE(numout,*)
-               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
-               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth
-               IF (ibr_max < nn_rimwidth)   CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file' )
-            ENDIF !Check igrd < 4
-            !
-         END DO
-         !
-      ENDIF 
-
-      ! Dispatch mapping indices and discrete distances on each processor
-      ! *****************************************************************
-     
-      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
-      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
-      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
-      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
-
-      DO igrd = igrd_start, igrd_end
-         icount  = 0
-         icountr = 0
-         nblen   (igrd) = 0
-         nblenrim(igrd) = 0
-         nblendta(igrd) = 0
-         DO ir=1, nn_rimwidth
-            DO ib = 1, jpbdta
-               ! check if point is in local domain and equals ir
-               IF(  nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND.   &
-                  & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND.   &
-                  & nbrdta(ib,igrd) == ir  ) THEN
-                  !
-                  icount = icount  + 1
-                  !
-                  IF( ir == 1 )   icountr = icountr+1
-                  IF (icount > jpbdim) THEN
-                     IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'
-                     nstop = nstop + 1
-                  ELSE
-                     nbi(icount, igrd)  = nbidta(ib,igrd)- mig(1)+1
-                     nbj(icount, igrd)  = nbjdta(ib,igrd)- mjg(1)+1
-                     nbr(icount, igrd)  = nbrdta(ib,igrd)
-                     nbmap(icount,igrd) = ib
-                  ENDIF            
-               ENDIF
-            END DO
-         END DO
-         nblenrim(igrd) = icountr !: length of rim boundary data on each proc
-         nblen   (igrd) = icount  !: length of boundary data on each proc        
-      END DO 
-
-      ! Compute rim weights
-      ! -------------------
-      DO igrd = igrd_start, igrd_end
-         DO ib = 1, nblen(igrd)
-            nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 )                     ! tanh formulation
-!           nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr(ib,igrd))/FLOAT(nn_rimwidth))**2      ! quadratic
-!           nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr(ib,igrd))/FLOAT(nn_rimwidth)          ! linear
-         END DO
-      END DO 
-   
       ! Mask corrections
       ! ----------------
@@ -361,22 +687,30 @@
       ! bdy masks and bmask are now set to zero on boundary points:
       igrd = 1       ! In the free surface case, bmask is at T-points
-      DO ib = 1, nblenrim(igrd)     
-        bmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
-      END DO
+      DO ib_bdy = 1, nb_bdy
+        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
+          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
+        ENDDO
+      ENDDO
       !
       igrd = 1
-      DO ib = 1, nblenrim(igrd)      
-        bdytmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
-      END DO
+      DO ib_bdy = 1, nb_bdy
+        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)      
+          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
+        ENDDO
+      ENDDO
       !
       igrd = 2
-      DO ib = 1, nblenrim(igrd)
-        bdyumask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
-      END DO
+      DO ib_bdy = 1, nb_bdy
+        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
+        ENDDO
+      ENDDO
       !
       igrd = 3
-      DO ib = 1, nblenrim(igrd)
-        bdyvmask(nbi(ib,igrd), nbj(ib,igrd)) = 0.e0
-      END DO
+      DO ib_bdy = 1, nb_bdy
+        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
+        ENDDO
+      ENDDO
 
       ! Lateral boundary conditions
@@ -384,20 +718,24 @@
       CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
 
-      IF( ln_vol .OR. ln_dyn_fla ) THEN      ! Indices and directions of rim velocity components
-         !
+      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
+
+         idx_bdy(ib_bdy)%flagu(:) = 0.e0
+         idx_bdy(ib_bdy)%flagv(:) = 0.e0
+         icount = 0 
+
          !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
          !flagu =  0 : u is tangential
          !flagu =  1 : u is normal to the boundary and is direction is inward
-         icount = 0 
-         flagu(:) = 0.e0
- 
+  
          igrd = 2      ! u-component 
-         DO ib = 1, nblenrim(igrd)  
-            zefl=bdytmask(nbi(ib,igrd)  , nbj(ib,igrd))
-            zwfl=bdytmask(nbi(ib,igrd)+1, nbj(ib,igrd))
-            IF( zefl + zwfl ==2 ) THEN
-               icount = icount +1
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
+            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+            zefl = bdytmask(nbi  ,nbj)
+            zwfl = bdytmask(nbi+1,nbj)
+            IF( zefl + zwfl == 2 ) THEN
+               icount = icount + 1
             ELSE
-               flagu(ib)=-zefl+zwfl
+               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
             ENDIF
          END DO
@@ -406,14 +744,15 @@
          !flagv =  0 : u is tangential
          !flagv =  1 : u is normal to the boundary and is direction is inward
-         flagv(:) = 0.e0
 
          igrd = 3      ! v-component
-         DO ib = 1, nblenrim(igrd)  
-            znfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd))
-            zsfl = bdytmask(nbi(ib,igrd), nbj(ib,igrd)+1)
-            IF( znfl + zsfl ==2 ) THEN
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
+            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
+            znfl = bdytmask(nbi,nbj  )
+            zsfl = bdytmask(nbi,nbj+1)
+            IF( znfl + zsfl == 2 ) THEN
                icount = icount + 1
             ELSE
-               flagv(ib) = -znfl + zsfl
+               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
             END IF
          END DO
@@ -422,5 +761,5 @@
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
-               ' are not boundary points. Check nbi, nbj, indices.'
+               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
             IF(lwp) WRITE(numout,*) ' ========== '
             IF(lwp) WRITE(numout,*)
@@ -428,5 +767,5 @@
          ENDIF 
     
-      ENDIF
+      ENDDO
 
       ! Compute total lateral surface for volume correction:
@@ -435,46 +774,41 @@
       IF( ln_vol ) THEN  
          igrd = 2      ! Lateral surface at U-points
-         DO ib = 1, nblenrim(igrd)
-            bdysurftot = bdysurftot + hu     (nbi(ib,igrd)  ,nbj(ib,igrd))                      &
-               &                    * e2u    (nbi(ib,igrd)  ,nbj(ib,igrd)) * ABS( flagu(ib) )   &
-               &                    * tmask_i(nbi(ib,igrd)  ,nbj(ib,igrd))                      &
-               &                    * tmask_i(nbi(ib,igrd)+1,nbj(ib,igrd))                   
-         END DO
+         DO ib_bdy = 1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               flagu => idx_bdy(ib_bdy)%flagu(ib)
+               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
+                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
+                  &                    * tmask_i(nbi  , nbj)                           &
+                  &                    * tmask_i(nbi+1, nbj)                   
+            ENDDO
+         ENDDO
 
          igrd=3 ! Add lateral surface at V-points
-         DO ib = 1, nblenrim(igrd)
-            bdysurftot = bdysurftot + hv     (nbi(ib,igrd),nbj(ib,igrd)  )                      &
-               &                    * e1v    (nbi(ib,igrd),nbj(ib,igrd)  ) * ABS( flagv(ib) )   &
-               &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)  )                      &
-               &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1)
-         END DO
+         DO ib_bdy = 1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               nbj => idx_bdy(ib_bdy)%nbi(ib,igrd)
+               flagv => idx_bdy(ib_bdy)%flagv(ib)
+               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
+                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
+                  &                    * tmask_i(nbi, nbj  )                           &
+                  &                    * tmask_i(nbi, nbj+1)
+            ENDDO
+         ENDDO
          !
          IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
       END IF   
-
-      ! Initialise bdy data arrays
-      ! --------------------------
-      tbdy(:,:) = 0.e0
-      sbdy(:,:) = 0.e0
-      ubdy(:,:) = 0.e0
-      vbdy(:,:) = 0.e0
-      sshbdy(:) = 0.e0
-      ubtbdy(:) = 0.e0
-      vbtbdy(:) = 0.e0
-#if defined key_lim2
-      frld_bdy(:) = 0.e0
-      hicif_bdy(:) = 0.e0
-      hsnif_bdy(:) = 0.e0
-#endif
-
-      ! Read in tidal constituents and adjust for model start time
-      ! ----------------------------------------------------------
-      IF( ln_tides )   CALL tide_data
       !
+      ! Tidy up
+      !--------
+      DEALLOCATE(nbidta, nbjdta, nbrdta)
+
    END SUBROUTINE bdy_init
 
 #else
    !!---------------------------------------------------------------------------------
-   !!   Dummy module                                   NO unstructured open boundaries
+   !!   Dummy module                                   NO open boundaries
    !!---------------------------------------------------------------------------------
 CONTAINS
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3116)
@@ -11,9 +11,8 @@
 #if defined key_bdy
    !!----------------------------------------------------------------------
-   !!   'key_bdy'     Unstructured Open Boundary Condition
+   !!   'key_bdy'     Open Boundary Condition
    !!----------------------------------------------------------------------
    !!   PUBLIC
-   !!      tide_init     : read of namelist
-   !!      tide_data     : read in and initialisation of tidal constituents at boundary
+   !!      tide_init     : read of namelist and initialisation of tidal harmonics data
    !!      tide_update   : calculation of tidal forcing at each timestep
    !!   PRIVATE
@@ -33,25 +32,23 @@
    USE bdy_oce         ! ocean open boundary conditions
    USE daymod          ! calendar
+   USE fldread, ONLY: fld_map
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tide_init     ! routine called in bdyini
-   PUBLIC   tide_data     ! routine called in bdyini
+   PUBLIC   tide_init     ! routine called in nemo_init
    PUBLIC   tide_update   ! routine called in bdydyn
 
-   LOGICAL, PUBLIC            ::   ln_tide_date          !: =T correct tide phases and amplitude for model start date
-   INTEGER, PUBLIC, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents
-   INTEGER, PUBLIC            ::   ntide                 !: Actual number of tidal constituents
-
-   CHARACTER(len=80), PUBLIC                         ::   filtide    !: Filename root for tidal input files
-   CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used.
-
-   INTEGER , PUBLIC, DIMENSION(jptides_max) ::   nindx        !: ???
-   REAL(wp), PUBLIC, DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
-   
-   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   ! Tidal constituents : SSH
-   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     ! Tidal constituents : U
-   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     ! Tidal constituents : V
+   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
+      INTEGER                                ::   ncpt       !: Actual number of tidal components
+      REAL(wp), POINTER, DIMENSION(:)        ::   speed      !: Phase speed of tidal constituent (deg/hr)
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V
+   END TYPE TIDES_DATA
+
+   INTEGER, PUBLIC, PARAMETER                  ::   jptides_max = 15      !: Max number of tidal contituents
+
+   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET ::   tides                 !: External tidal harmonics data
    
    !!----------------------------------------------------------------------
@@ -66,8 +63,25 @@
       !!                    ***  SUBROUTINE tide_init  ***
       !!                     
-      !! ** Purpose : - Read in namelist for tides
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER ::   itide                  ! dummy loop index 
+      !! ** Purpose : - Read in namelist for tides and initialise external
+      !!                tidal harmonics data
+      !!
+      !!----------------------------------------------------------------------
+      !! namelist variables
+      !!-------------------
+      LOGICAL                                   ::   ln_tide_date !: =T correct tide phases and amplitude for model start date
+      CHARACTER(len=80)                         ::   filtide      !: Filename root for tidal input files
+      CHARACTER(len= 4), DIMENSION(jptides_max) ::   tide_cpt     !: Names of tidal components used.
+      REAL(wp),          DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
+      !!
+      INTEGER, DIMENSION(jptides_max)           ::   nindx              !: index of pre-set tidal components
+      INTEGER                                   ::   ib_bdy, itide, ib  !: dummy loop indices
+      INTEGER                                   ::   inum, igrd
+      INTEGER, DIMENSION(3)                     ::   ilen0              !: length of boundary data (from OBC arrays)
+      CHARACTER(len=80)                         ::   clfile             !: full file name for tidal input file 
+      REAL(wp)                                  ::   z_arg, z_atde, z_btde, z1t, z2t           
+      REAL(wp),DIMENSION(jptides_max)           ::   z_vplu, z_ftc
+      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read           !: work space to read in tidal harmonics data
+      !!
+      TYPE(TIDES_DATA),  POINTER                ::   td                 !: local short cut   
       !!
       NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed
@@ -78,262 +92,239 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~'
 
-      ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
-      ln_tide_date = .false.
-      filtide(:) = ''
-      tide_speed(:) = 0.0
-      tide_cpt(:) = ''
-      REWIND( numnam )                                ! Read namelist parameters
-      READ  ( numnam, nambdy_tide )
-      !                                               ! Count number of components specified
-      ntide = jptides_max
-      DO itide = 1, jptides_max
-        IF( tide_cpt(itide) == '' ) THEN
-           ntide = itide-1
-           exit
-        ENDIF
-      END DO
-
-      !                                               ! find constituents in standard list
-      DO itide = 1, ntide
-         nindx(itide) = 0
-         IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
-         IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
-         IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
-         IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
-         IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
-         IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
-         IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
-         IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
-         IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
-         IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
-         IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
-         IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
-         IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
-         IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
-         IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
-         IF( nindx(itide) == 0  .AND. lwp ) THEN
-            WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
-            CALL ctl_warn( ctmp1 )
-         ENDIF
-      END DO
-      !                                               ! Parameter control and print
-      IF( ntide < 1 ) THEN 
-         CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' )
-      ELSE
-         IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
-         IF(lwp) WRITE(numout,*) '             tidal components specified ', ntide
-         IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:ntide)
-         IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
-         IF(lwp) WRITE(numout,*) '                ', tide_speed(1:ntide)
-      ENDIF
-
-      ! Initialisation of tidal harmonics arrays
-      sshtide(:) = 0.e0
-      utide  (:) = 0.e0
-      vtide  (:) = 0.e0
-      !
+      REWIND(numnam)
+      DO ib_bdy = 1, nb_bdy
+         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
+
+            td => tides(ib_bdy)
+
+            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
+            ln_tide_date = .false.
+            filtide(:) = ''
+            tide_speed(:) = 0.0
+            tide_cpt(:) = ''
+
+            ! Don't REWIND here - may need to read more than one of these namelists.
+            READ  ( numnam, nambdy_tide )
+            !                                               ! Count number of components specified
+            td%ncpt = 0
+            DO itide = 1, jptides_max
+              IF( tide_cpt(itide) /= '' ) THEN
+                 td%ncpt = td%ncpt + 1
+              ENDIF
+            END DO
+
+            ! Fill in phase speeds from namelist
+            ALLOCATE( td%speed(td%ncpt) )
+            td%speed = tide_speed(1:td%ncpt)
+
+            ! Find constituents in standard list
+            DO itide = 1, td%ncpt
+               nindx(itide) = 0
+               IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
+               IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
+               IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
+               IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
+               IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
+               IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
+               IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
+               IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
+               IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
+               IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
+               IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
+               IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
+               IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
+               IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
+               IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
+               IF( nindx(itide) == 0  .AND. lwp ) THEN
+                  WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
+                  CALL ctl_warn( ctmp1 )
+               ENDIF
+            END DO
+            !                                               ! Parameter control and print
+            IF( td%ncpt < 1 ) THEN 
+               CALL ctl_stop( '          Did not find any tidal components in namelist nambdy_tide' )
+            ELSE
+               IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
+               IF(lwp) WRITE(numout,*) '             tidal components specified ', td%ncpt
+               IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:td%ncpt)
+               IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
+               IF(lwp) WRITE(numout,*) '                ', tide_speed(1:td%ncpt)
+            ENDIF
+
+
+            ! Allocate space for tidal harmonics data - 
+            ! get size from OBC data arrays
+            ! ---------------------------------------
+
+            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 
+            ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) )
+
+            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 
+            ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) )
+
+            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 
+            ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) )
+
+            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) )
+
+
+            ! Open files and read in tidal forcing data
+            ! -----------------------------------------
+
+            DO itide = 1, td%ncpt
+               !                                                              ! SSH fields
+               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
+               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               CALL iom_open( clfile, inum )
+               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+               td%ssh(:,itide,1) = dta_read(1:ilen0(1),1,1)
+               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
+               td%ssh(:,itide,2) = dta_read(1:ilen0(1),1,1)
+               CALL iom_close( inum )
+               !                                                              ! U fields
+               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
+               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               CALL iom_open( clfile, inum )
+               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+               td%u(:,itide,1) = dta_read(1:ilen0(2),1,1)
+               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
+               td%u(:,itide,2) = dta_read(1:ilen0(2),1,1)
+               CALL iom_close( inum )
+               !                                                              ! V fields
+               clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
+               IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
+               CALL iom_open( clfile, inum )
+               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+               td%v(:,itide,1) = dta_read(1:ilen0(3),1,1)
+               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
+               td%v(:,itide,2) = dta_read(1:ilen0(3),1,1)
+               CALL iom_close( inum )
+               !
+            END DO ! end loop on tidal components
+
+            IF( ln_tide_date ) THEN      ! correct for date factors
+
+!! used nmonth, nyear and nday from daymod....
+               ! Calculate date corrects for 15 standard consituents
+               ! This is the initialisation step, so nday, nmonth, nyear are the 
+               ! initial date/time of the integration.
+                 print *, nday,nmonth,nyear
+                 nyear  = int(ndate0 / 10000  )                          ! initial year
+                 nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month
+                 nday   = int(ndate0 - nyear * 10000 - nmonth * 100)
+
+               CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
+
+               IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
+
+               DO itide = 1, td%ncpt       ! loop on tidal components
+                  !
+                  IF( nindx(itide) /= 0 ) THEN
+!!gm use rpi  and rad global variable  
+                     z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
+                     z_atde=z_ftc(nindx(itide))*cos(z_arg)
+                     z_btde=z_ftc(nindx(itide))*sin(z_arg)
+                     IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide),   &
+                        &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
+                  ELSE
+                     z_atde = 1.0_wp
+                     z_btde = 0.0_wp
+                  ENDIF
+                  !                                         !  elevation         
+                  igrd = 1
+                  DO ib = 1, ilen0(igrd)                
+                     z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2)
+                     z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1)
+                     td%ssh(ib,itide,1) = z1t
+                     td%ssh(ib,itide,2) = z2t
+                  END DO
+                  !                                         !  u       
+                  igrd = 2
+                  DO ib = 1, ilen0(igrd)                
+                     z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2)
+                     z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1)
+                     td%u(ib,itide,1) = z1t
+                     td%u(ib,itide,2) = z2t
+                  END DO
+                  !                                         !  v       
+                  igrd = 3
+                  DO ib = 1, ilen0(igrd)                
+                     z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2)
+                     z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1)
+                     td%v(ib,itide,1) = z1t
+                     td%v(ib,itide,2) = z2t
+                  END DO
+                  !
+               END DO   ! end loop on tidal components
+               !
+            ENDIF ! date correction
+            !
+         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
+         !
+      END DO ! loop on ib_bdy
+
    END SUBROUTINE tide_init
 
 
-   SUBROUTINE tide_data
-      !!----------------------------------------------------------------------
-      !!                    ***  SUBROUTINE tide_data  ***
-      !!                    
-      !! ** Purpose : - Read in tidal harmonics data and adjust for the start 
-      !!                time of the model run. 
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER ::   itide, igrd, ib        ! dummy loop indices
-      CHARACTER(len=80) :: clfile         ! full file name for tidal input file 
-      INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read)
-      INTEGER, DIMENSION(6) :: lendta=0   ! length of data in the file (note may be different from nblendta!)
-      REAL(wp) ::  z_arg, z_atde, z_btde, z1t, z2t           
-      REAL(wp), DIMENSION(jpbdta,1) ::   zdta   ! temporary array for data fields
-      REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc
-      !!------------------------------------------------------------------------------
-
-      ! Open files and read in tidal forcing data
-      ! -----------------------------------------
-
-      ipj = 1
-
-      DO itide = 1, ntide
-         !                                                              ! SSH fields
-         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
-         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-         CALL iom_open( clfile, inum )
-         igrd = 4
-         IF( nblendta(igrd) <= 0 ) THEN 
-            idvar = iom_varid( inum,'z1' )
-            IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar)
-            nblendta(igrd) = iom_file(inum)%dimsz(1,idvar)
-            WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd)
-         ENDIF
-         ipi = nblendta(igrd)
-         CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) )
-         DO ib = 1, nblenrim(igrd)
-            ssh1(ib,itide) = zdta(nbmap(ib,igrd),1)
-         END DO
-         CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) )
-         DO ib = 1, nblenrim(igrd)
-            ssh2(ib,itide) = zdta(nbmap(ib,igrd),1)
-         END DO
-         CALL iom_close( inum )
-         !
-         !                                                              ! U fields
-         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
-         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
-         CALL iom_open( clfile, inum )
-         igrd = 5
-         IF( lendta(igrd) <= 0 ) THEN 
-            idvar = iom_varid( inum,'u1' )
-            lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
-            WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd)
-         ENDIF
-         ipi = lendta(igrd)
-         CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) )
-         DO ib = 1, nblenrim(igrd)
-            u1(ib,itide) = zdta(nbmap(ib,igrd),1)
-         END DO
-         CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) )
-         DO ib = 1, nblenrim(igrd)
-            u2(ib,itide) = zdta(nbmap(ib,igrd),1)
-         END DO
-         CALL iom_close( inum )
-         !
-         !                                                              ! V fields
-         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
-         if(lwp) write(numout,*) 'Reading data from file ', clfile
-         CALL iom_open( clfile, inum )
-         igrd = 6
-         IF( lendta(igrd) <= 0 ) THEN 
-            idvar = iom_varid( inum,'v1' )
-            lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
-            WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd)
-         ENDIF
-         ipi = lendta(igrd)
-         CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) )
-         DO ib = 1, nblenrim(igrd)
-            v1(ib,itide) = zdta(nbmap(ib,igrd),1)
-         END DO
-         CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) )
-         DO ib=1, nblenrim(igrd)
-            v2(ib,itide) = zdta(nbmap(ib,igrd),1)
-         END DO
-         CALL iom_close( inum )
-         !
-      END DO ! end loop on tidal components
-
-      IF( ln_tide_date ) THEN      ! correct for date factors
-
-!! used nmonth, nyear and nday from daymod....
-         ! Calculate date corrects for 15 standard consituents
-         ! This is the initialisation step, so nday, nmonth, nyear are the 
-         ! initial date/time of the integration.
-           print *, nday,nmonth,nyear
-           nyear  = int(ndate0 / 10000  )                           ! initial year
-           nmonth = int((ndate0 - nyear * 10000 ) / 100 )          ! initial month
-           nday   = int(ndate0 - nyear * 10000 - nmonth * 100)
-
-         CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
-
-         IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
-
-         DO itide = 1, ntide       ! loop on tidal components
-            !
-            IF( nindx(itide) /= 0 ) THEN
-!!gm use rpi  and rad global variable  
-               z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
-               z_atde=z_ftc(nindx(itide))*cos(z_arg)
-               z_btde=z_ftc(nindx(itide))*sin(z_arg)
-               IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide),   &
-                  &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
-            ELSE
-               z_atde = 1.0_wp
-               z_btde = 0.0_wp
-            ENDIF
-            !                                         !  elevation         
-            igrd = 4
-            DO ib = 1, nblenrim(igrd)                
-               z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide)
-               z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide)
-               ssh1(ib,itide) = z1t
-               ssh2(ib,itide) = z2t
-            END DO
-            !                                         !  u       
-            igrd = 5
-            DO ib = 1, nblenrim(igrd)                
-               z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide)
-               z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide)
-               u1(ib,itide) = z1t
-               u2(ib,itide) = z2t
-            END DO
-            !                                         !  v       
-            igrd = 6
-            DO ib = 1, nblenrim(igrd)                
-               z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide)
-               z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide)
-               v1(ib,itide) = z1t
-               v2(ib,itide) = z2t
-            END DO
-            !
-         END DO   ! end loop on tidal components
-         !
-      ENDIF ! date correction
-      !
-   END SUBROUTINE tide_data
-
-
-   SUBROUTINE tide_update ( kt, jit )
+   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE tide_update  ***
       !!                
-      !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays. 
+      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
       !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter
+      INTEGER, INTENT( in )          ::   kt      ! Main timestep counter
 !!gm doctor jit ==> kit
-      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option)
-      !!
-      INTEGER  ::   itide, igrd, ib      ! dummy loop indices
-      REAL(wp) ::   z_arg, z_sarg            !            
+      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
 
-      DO itide = 1, ntide
-         z_sarg = z_arg * tide_speed(itide)
+      DO itide = 1, td%ncpt
+         z_sarg = z_arg * td%speed(itide)
          z_cost(itide) = COS( z_sarg )
          z_sist(itide) = SIN( z_sarg )
       END DO
 
-      ! summing of tidal constituents into BDY arrays
-      sshtide(:) = 0.0
-      utide (:) = 0.0
-      vtide (:) = 0.0
-      !
-      DO itide = 1, ntide
-         igrd=4                              ! SSH on tracer grid.
-         DO ib = 1, nblenrim(igrd)
-            sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide)
-            !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide)
+      DO itide = 1, td%ncpt
+         igrd=1                              ! SSH on tracer grid.
+         DO ib = 1, idx%nblenrim(igrd)
+            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)
+            !    if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2)
          END DO
-         igrd=5                              ! U grid
-         DO ib=1, nblenrim(igrd)
-            utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide)
-            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide)
+         igrd=2                              ! U grid
+         DO ib=1, idx%nblenrim(igrd)
+            dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide)
+            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2)
          END DO
-         igrd=6                              ! V grid
-         DO ib=1, nblenrim(igrd)
-            vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide)
-            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide)
+         igrd=3                              ! V grid
+         DO ib=1, idx%nblenrim(igrd)
+            dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide)
+            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2)
          END DO
       END DO
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 3116)
@@ -11,9 +11,11 @@
    !!   'key_bdy'                     Unstructured Open Boundary Conditions
    !!----------------------------------------------------------------------
-   !!   bdy_tra_frs        : Relaxation of tracers on unstructured open boundaries
+   !!   bdy_tra            : Apply open boundary conditions to T and S
+   !!   bdy_tra_frs        : Apply Flow Relaxation Scheme
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain variables 
    USE bdy_oce         ! ocean open boundary conditions
+   USE bdydta, ONLY:   bf
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! I/O manager
@@ -22,5 +24,5 @@
    PRIVATE
 
-   PUBLIC bdy_tra_frs     ! routine called in tranxt.F90 
+   PUBLIC bdy_tra      ! routine called in tranxt.F90 
 
    !!----------------------------------------------------------------------
@@ -31,14 +33,40 @@
 CONTAINS
 
-   SUBROUTINE bdy_tra_frs( kt )
+   SUBROUTINE bdy_tra( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d  ***
+      !!
+      !! ** Purpose : - Apply open boundary conditions for baroclinic velocities
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) :: kt     ! Main time step counter
+      !!
+      INTEGER               :: ib_bdy ! Loop index
+
+      DO ib_bdy=1, nb_bdy
+
+         SELECT CASE( nn_tra(ib_bdy) )
+         CASE(jp_none)
+            CYCLE
+         CASE(jp_frs)
+            CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt )
+         CASE DEFAULT
+            CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
+         END SELECT
+      ENDDO
+
+   END SUBROUTINE bdy_tra
+
+   SUBROUTINE bdy_tra_frs( idx, dta, kt )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_tra_frs  ***
       !!                    
-      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers in the  
-      !!              case of unstructured open boundaries.
+      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries.
       !! 
       !! Reference : Engedahl H., 1995, Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt
+      INTEGER,         INTENT(in) ::   kt
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data
       !! 
       REAL(wp) ::   zwgt           ! boundary weight
@@ -47,28 +75,20 @@
       !!----------------------------------------------------------------------
       !
-      IF(ln_tra_frs) THEN       ! If this is false, then this routine does nothing. 
-         !
-         IF( kt == nit000 ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'bdy_tra_frs : Flow Relaxation Scheme for tracers'
-            IF(lwp) WRITE(numout,*) '~~~~~~~'
-         ENDIF
-         !
-         igrd = 1                       ! Everything is at T-points here
-         DO ib = 1, nblen(igrd)
-            DO ik = 1, jpkm1
-               ii = nbi(ib,igrd)
-               ij = nbj(ib,igrd)
-               zwgt = nbw(ib,igrd)
-               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)         
-               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)
-            END DO
-         END DO 
-         !                                              ! Boundary points should be updated
-         CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )     
-         CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    
-         !
-      ENDIF ! ln_tra_frs
       !
+      igrd = 1                       ! Everything is at T-points here
+      DO ib = 1, idx%nblen(igrd)
+         DO ik = 1, jpkm1
+            ii = idx%nbi(ib,igrd)
+            ij = idx%nbj(ib,igrd)
+            zwgt = idx%nbw(ib,igrd)
+            tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta%tem(ib,ik) - tsa(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik)         
+            tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta%sal(ib,ik) - tsa(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik)
+         END DO
+      END DO 
+      !
+      CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. )   ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. )    ! Boundary points should be updated
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+   !
    END SUBROUTINE bdy_tra_frs
    
@@ -78,7 +98,7 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE bdy_tra_frs(kt)      ! Empty routine
-      WRITE(*,*) 'bdy_tra_frs: You should not have seen this print! error?', kt
-   END SUBROUTINE bdy_tra_frs
+   SUBROUTINE bdy_tra(kt)      ! Empty routine
+      WRITE(*,*) 'bdy_tra: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_tra
 #endif
 
Index: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 3115)
+++ branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 3116)
@@ -71,6 +71,7 @@
       !!
       INTEGER  ::   ji, jj, jk, jb, jgrd
-      INTEGER  ::   ii, ij
+      INTEGER  ::   ib_bdy, ii, ij
       REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
+      TYPE(OBC_INDEX), POINTER :: idx
       !!-----------------------------------------------------------------------------
 
@@ -91,19 +92,24 @@
       ! ------------------------------------------------
       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)
+      DO ib_bdy = 1, nb_bdy
+         idx => idx_bdy(ib_bdy)
+
+         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
-      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) 
+         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
@@ -118,21 +124,26 @@
       ! -------------------------------------------------------------
       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)
+      DO ib_bdy = 1, nb_bdy
+         idx => idx_bdy(ib_bdy)
+
+         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
-      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)
+         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
