Index: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 2715)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3294)
@@ -10,19 +10,22 @@
    !!            3.3  !  2010-09  (E.O'Dea) modifications for Shelf configurations 
    !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
+   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
    !!----------------------------------------------------------------------
 #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 wrk_nemo        ! Memory Allocation
+   USE timing          ! Timing
    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 +36,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 +57,640 @@
 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.
-      !!
-      !! ** 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.
+      !! ** Purpose :   Update external data for open boundary conditions
+      !!
+      !! ** Method  :   Use fldread.F90
+      !!                
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index (for timesplitting option, otherwise zero)
-      !!
-      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( nn_timing == 1 ) CALL timing_start('bdy_dta')
+
+      ! Initialise data arrays once for all from initial conditions where required
+      !---------------------------------------------------------------------------
+      IF( kt .eq. nit000 .and. .not. PRESENT(jit) ) THEN
+
+         ! Calculate depth-mean currents
+         !-----------------------------
+         CALL wrk_alloc(jpi,jpj,pu2d,pv2d) 
+
+         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) = tn(ii,ij,ik)
-                     sbdy(ib,ik) = sn(ii,ij,ik)
+                     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
+
+         CALL wrk_dealloc(jpi,jpj,pu2d,pv2d) 
+
+      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( nn_timing == 1 ) CALL timing_stop('bdy_dta')
+
+      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
+      !!---------------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init')
+
+      ! 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_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 
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_init')
+
+      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
 
