Index: trunk/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- trunk/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 1058)
+++ trunk/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 1125)
@@ -1,14 +1,20 @@
 MODULE bdydta
-   !!==============================================================================
-   !!                            ***  MODULE bdydta  ***
+   !!======================================================================
+   !!                       ***  MODULE bdydta  ***
    !! Open boundary data : read the data for the unstructured open boundaries.
-   !!==============================================================================
-#if defined key_bdy || defined key_bdy_tides
-   !!------------------------------------------------------------------------------
-   !!   'key_bdy'         :                   Unstructured Open Boundary Conditions
-   !!------------------------------------------------------------------------------
-   !!   bdy_dta           : read u, v, t, s data along each open boundary
-   !!------------------------------------------------------------------------------
-   !! * Modules used
+   !!======================================================================
+   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
+   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
+   !!             -   !  2007-07  (D. Storkey) add bdy_dta_bt
+   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
+   !!----------------------------------------------------------------------
+#if defined key_bdy
+   !!----------------------------------------------------------------------
+   !!   'key_bdy'                     Unstructured Open Boundary Conditions
+   !!----------------------------------------------------------------------
+   !!   bdy_dta    : read u, v, t, s data along open boundaries
+   !!   bdy_dta_bt : read depth-mean velocities and elevation along open
+   !!                boundaries        
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
@@ -24,82 +30,621 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC bdy_dta                         ! routines called by step.F90
-   PUBLIC bdy_dta_bt 
+   PUBLIC   bdy_dta      ! routines called by step.F90
+   PUBLIC   bdy_dta_bt 
+
+   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 model 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 model 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), DIMENSION(jpbdim,jpk,2) ::   tbdydta, sbdydta     !: time interpolated values of T and S bdy data   
+   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   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
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! $Id: $ 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bdy_dta( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  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.
+      !!----------------------------------------------------------------------
+      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                          ! 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), DIMENSION(jpbtime)      ::   zstepr             ! REAL time array from data files
+      REAL(wp), DIMENSION(jpbdta,1,jpk) ::   zdta               ! temporary array for data fields
+      !!---------------------------------------------------------------------------
+
+      IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_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 variables
+ 
+      iman = INT( raamo )      ! Number of months in a year
+
+      i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 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 : Initialize unstructured boundary data'
+         IF(lwp) WRITE(numout,*)    '~~~~~~~' 
+
+         IF     ( nbdy_dta == 0 ) THEN
+            !
+            IF(lwp) WRITE(numout,*) '          Bdy data are taken from initial conditions'
+            !
+         ELSEIF (nbdy_dta == 1) THEN
+            !
+            IF(lwp) WRITE(numout,*) '          Bdy data are read in netcdf files'
+            !
+            dayfrac = adatrj  - FLOAT( itimer ) / 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) = filbdy_data_T
+            clfile(2) = filbdy_data_U
+            clfile(3) = filbdy_data_V
+            !                                                  
+            ! how many files are we to read in?
+            igrd_start = 1
+            igrd_end   = 3
+            IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN
+               ! No T-grid file.
+               igrd_start = 2
+            ELSEIF ( .NOT. ln_bdy_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), 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) , 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
+
+               ! 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.NE.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
+               !
+               IF    ( igrd == 1 ) THEN
+                 ntimes_bdyt = ntimes_bdy
+                 zoffsett = zoffset
+                 istept(:) = INT( zstepr(:) + zoffset )
+               ELSEIF(igrd == 2 ) THEN
+                 ntimes_bdyu = ntimes_bdy
+                 zoffsetu = zoffset
+                 istepu(:) = INT( zstepr(:) + zoffset )
+               ELSEIF(igrd == 3 ) THEN
+                 ntimes_bdyv = ntimes_bdy
+                 zoffsetv = zoffset
+                 istepv(:) = INT( zstepr(:) + zoffset )
+               ENDIF
+               !
+            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_bdy_clim ) THEN
+              CALL ctl_stop( 'There is only one time dump in data files',   &
+                 &           'Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
+            ENDIF
+
+            IF( ln_bdy_clim ) THEN
+              IF( ntimes_bdy /= 1 .AND. ntimes_bdy /= 12 ) THEN
+                 CALL ctl_stop( 'For climatological boundary forcing (ln_bdy_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
+            !
+            WRITE(numout,*) 'Time offset is ',zoffset
+            WRITE(numout,*) 'First record to read is ',nbdy_b
+
+         ENDIF ! endif (nbdy_dta == 1)
+
+
+         ! 1.2  Read first record in file if necessary (ie if nbdy_dta == 1)
+         ! *****************************************************************
+
+         IF( nbdy_dta == 0) THEN      ! boundary data arrays are filled with initial conditions
+            !
+            IF (ln_bdy_tra_frs) THEN
+              igrd = 1            ! T-points data 
+              DO ib = 1, nblen(igrd)
+                ii = nbi(ib,igrd)
+                ij = nbj(ib,igrd)
+                DO ik = 1, jpkm1
+                  tbdy(ib,ik) = tn(ii, ij, ik)
+                  sbdy(ib,ik) = sn(ii, ij, ik)
+                ENDDO
+              END DO
+            ENDIF
+
+            IF(ln_bdy_dyn_frs) THEN
+              igrd = 2            ! U-points data 
+              DO ib = 1, nblen(igrd)
+                ii = nbi(ib,igrd)
+                ij = nbj(ib,igrd)
+                DO ik = 1, jpkm1
+                  ubdy(ib,ik) = un(ii, ij, ik)
+                ENDDO
+              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)
+                ENDDO
+              END DO
+            ENDIF
+            !
+         ELSEIF( nbdy_dta == 1 ) THEN    ! Set first record in the climatological case:   
+            !
+            IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN
+               nbdy_a = 1
+            ELSEIF( ln_bdy_clim .AND. ntimes_bdy == iman ) THEN
+               nbdy_b = 0
+               nbdy_a = imois
+            ELSE
+               nbdy_a = nbdy_b
+            ENDIF
    
-   !! * local modules variables
-   INTEGER ::   &
-      bdynumt,         &  ! logical unit for T-points data file
-      bdynumu,         &  ! logical unit for U-points data file
-      bdynumv,         &  ! logical unit for V-points data file
-      kbdy,            &  ! Exact number of time dumps in data files
-      nbdy1,           &  ! Time record in bdy data file BEFORE the model current time
-      nbdy2               ! Time record in bdy data file AFTER the model current time
-
-   !! * local modules variables for barotropic variables
-   INTEGER ::   &
-      bdynumt_bt,         &  ! logical unit for T-points data file
-      bdynumu_bt,         &  ! logical unit for U-points data file
-      bdynumv_bt,         &  ! logical unit for V-points data file
-      kbdy_bt,            &  ! Exact number of time dumps in data files
-      nbdy1_bt,           &  ! Time record in bdy data file BEFORE the model current time
-      nbdy2_bt               ! Time record in bdy data file AFTER the model current time
-
-   INTEGER, DIMENSION (jpbtime) ::   & 
-      istep, istep_bt        ! time array in seconds in each data file
-
-   REAL(wp) ::  &
-      offset              ! time offset between time origin in file and start time of model run
-
-   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   &
-      tbdydta, sbdydta,	& !: Arrays used for time interpolation of bdy data   
-      ubdydta, vbdydta
-
-   REAL(wp), DIMENSION(jpbdim,2) ::   &
-      ubtbdydta, vbtbdydta,	& !: Arrays used for time interpolation of bdy data   
-      sshbdydta                   !: bdy data of ssh
-
-   !!---------------------------------------------------------------------------------
-   !!   OPA 9.0 , LODYC-IPSL  (2003)
-   !!---------------------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE bdy_dta( kt )
+            ! Read first record:
+            ipj  = 1
+            ipk  = jpk
+            igrd = 1
+            ipi  = nblendta(igrd)
+
+            IF(ln_bdy_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
+               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
+               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_bdy_tra_frs
+ 
+            IF(ln_bdy_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
+               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
+               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_bdy_dyn_frs
+
+
+            IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN
+               ! First data time is after start of run
+               ! Put first value in both time levels
+               nbdy_b = nbdy_a
+               IF(ln_bdy_tra_frs) THEN
+                 tbdydta(:,:,1) = tbdydta(:,:,2)
+                 sbdydta(:,:,1) = sbdydta(:,:,2)
+               ENDIF
+               IF(ln_bdy_dyn_frs) THEN
+                 ubdydta(:,:,1) = ubdydta(:,:,2)
+                 vbdydta(:,:,1) = vbdydta(:,:,2)
+               ENDIF
+            END IF
+
+         END IF ! nbdy_dta == 0/1
+ 
+         ! In the case of constant boundary forcing fill bdy arrays once for all
+         IF ((ln_bdy_clim).AND.(ntimes_bdy==1)) THEN
+            IF(ln_bdy_tra_frs) THEN
+              tbdy  (:,:) = tbdydta  (:,:,2)
+              sbdy  (:,:) = sbdydta  (:,:,2)
+            ENDIF
+            IF(ln_bdy_dyn_frs) THEN
+              ubdy  (:,:) = ubdydta  (:,:,2)
+              vbdy  (:,:) = vbdydta  (:,:,2)
+            ENDIF
+
+            IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )
+            IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu )
+            IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv )
+         END IF
+
+      ENDIF                                            ! End if nit000
+
+
+      !                                                !---------------------!
+      !                                                !  at each time step  !
+      !                                                !---------------------!
+
+      IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN 
+         !
+         ! Read one more record if necessary
+         !**********************************
+
+        IF( ln_bdy_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_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN
+
+           IF ( nbdy_a < ntimes_bdy ) THEN
+              nbdy_b = nbdy_a
+              nbdy_a = nbdy_a + 1
+              lect  =.true.
+           ELSE
+              ! We have reached the end of the file
+              ! put the last data time into both time levels
+              nbdy_b = nbdy_a
+              IF(ln_bdy_tra_frs) THEN
+                tbdydta(:,:,1) =  tbdydta(:,:,2)
+                sbdydta(:,:,1) =  sbdydta(:,:,2)
+              ENDIF
+              IF(ln_bdy_dyn_frs) THEN
+                ubdydta(:,:,1) =  ubdydta(:,:,2)
+                vbdydta(:,:,1) =  vbdydta(:,:,2)
+              ENDIF
+            END IF ! nbdy_a < ntimes_bdy
+
+        END IF
+         
+        IF( lect ) THEN
+           ! Swap arrays
+           IF(ln_bdy_tra_frs) THEN
+             tbdydta(:,:,1) =  tbdydta(:,:,2)
+             sbdydta(:,:,1) =  sbdydta(:,:,2)
+           ENDIF
+           IF(ln_bdy_dyn_frs) THEN
+             ubdydta(:,:,1) =  ubdydta(:,:,2)
+             vbdydta(:,:,1) =  vbdydta(:,:,2)
+           ENDIF
+ 
+           ! read another set
+           ipj  = 1
+           ipk  = jpk
+
+           IF(ln_bdy_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_bdy_tra_frs
+
+           IF(ln_bdy_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_bdy_dyn_frs
+
+           !
+           IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy_b ',nbdy_b
+           IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy_a ',nbdy_a
+           IF (.NOT.ln_bdy_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_bdy_clim ) THEN   ;   zxy = FLOAT( nday                  ) / FLOAT( nobis(nbdy_b) ) + 0.5 - i15
+       ELSE                     ;   zxy = FLOAT( istep(nbdy_b) - itimer ) / FLOAT( istep(nbdy_b) - istep(nbdy_a) )
+       END IF
+
+          IF(ln_bdy_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_bdy_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
+
+      END IF                       !end if ((nbdy_dta==1).AND.(ntimes_bdy>1))
+    
+
+      !                                                !---------------------!
+      !                                                !     last call       !
+      !                                                !---------------------!
+      IF( kt == nitend ) THEN
+          IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )              ! Closing of the 3 files
+          IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu )
+          IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv )
+      ENDIF
+      !
+      ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs
+
+   END SUBROUTINE bdy_dta
+
+
+   SUBROUTINE bdy_dta_bt( kt, jit )
       !!---------------------------------------------------------------------------
-      !!                      ***  SUBROUTINE bdy_dta  ***
+      !!                      ***  SUBROUTINE bdy_dta_bt  ***
       !!                    
-      !! ** Purpose :   Read unstructured boundary data
+      !! ** Purpose :   Read unstructured boundary data for Flather condition
       !!
-      !! ** Method  :   
+      !! ** 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
+      !                                      ! (for timesplitting option, otherwise zero)
       !!
-      !! History :  OPA 9.0 ! 05-01 (J. Chanut, A. Sellar) Original
-      !!                 "  ! 07-01 (D. Storkey) Update to use IOM module.
-      !!---------------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
-                                             ! (for timesplitting option, otherwise zero)
-
-      !! * Local declarations
       LOGICAL ::   lect                      ! flag for reading
-      INTEGER ::   jt, jb, jk, jgrd          ! dummy loop indices
+      INTEGER ::   it, ib, igrd              ! dummy loop indices
       INTEGER ::   idvar                     ! netcdf var ID
       INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing
-      INTEGER ::   kbdyt, kbdyu, kbdyv
+      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 ::   kyear, kmonth, kday, ksecs
+      INTEGER ::   iyear, imonth, iday, isecs
       INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files
-      REAL(wp) ::   dayfrac, zxy, offsett
-      REAL(wp) ::   offsetu, offsetv
-      REAL(wp) ::   dayjul0, zdayjulini
-      REAL(wp), DIMENSION(jpbtime)      ::   istepr             ! REAL time array from data files
-      REAL(wp), DIMENSION(jpbdta,1,jpk) ::   pdta3              ! temporary array for data fields
-      CHARACTER(LEN=80), DIMENSION(3)   ::   bdyfile
+      REAL(wp) ::   dayfrac, zxy, zoffsett
+      REAL(wp) ::   zoffsetu, zoffsetv
+      REAL(wp) ::   dayjul0, zdayjulini, zntotime
+      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(3)   ::   clfile
       CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate
-
       !!---------------------------------------------------------------------------
+
+!!gm   add here the same style as in bdy_dta
+!!gm      clearly bdy_dta_bt and bdy_dta  can be combined...   
+!!gm      too many things duplicated in the read of data...   simplification can be done
 
       ! -------------------- !
@@ -111,4 +656,5 @@
       ! Some time variables for monthly climatological forcing:
       ! *******************************************************
+ !!gm  here  use directely daymod variables
  
       iman  = INT( raamo ) ! Number of months in a year
@@ -123,25 +669,38 @@
       ! *********************************************
 
-      itimer = (kt-nit000+1)*rdt      ! current time in seconds for interpolation 
-
-      ! -------------------- !
-      ! 1.   First call only !
-      ! -------------------- !
-
-      IF ( kt == nit000 ) THEN
-
-        istep(:) = 0
-
-        nbdy1=0
-        nbdy2=0
-
-      ! 1.1  Get time information from bdy data file
-      ! ********************************************
+      itimer = ((kt-1)-nit000+1)*rdt      ! current time in seconds for interpolation 
+      itimer = itimer + jit*rdtbt         ! in non-climatological case
+
+      IF ( ln_bdy_tides ) THEN
+
+         ! -------------------------------------!
+         ! Update BDY fields with tidal forcing !
+         ! -------------------------------------!  
+
+         CALL tide_update( kt, jit ) 
+  
+      ENDIF
+
+      IF ( ln_bdy_dyn_fla ) THEN
+
+         ! -------------------------------------!
+         ! Update BDY fields with model data    !
+         ! -------------------------------------!  
+
+      !                                                !-------------------!
+      IF( kt == nit000 ) 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 :Initialize unstructured boundary data'
+        IF(lwp) WRITE(numout,*)    'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.'
         IF(lwp) WRITE(numout,*)    '~~~~~~~' 
 
-        IF (nbdy_dta == 0) THEN
+        IF( nbdy_dta == 0 ) THEN
           IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions'
 
@@ -154,13 +713,21 @@
                                            ! necessary time dumps in file are included
 
-          bdyfile(1) = filbdy_data_T
-          bdyfile(2) = filbdy_data_U
-          bdyfile(3) = filbdy_data_V
-
-          DO jgrd = 1,3
-
-            CALL iom_open( bdyfile(jgrd), inum )
-            CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits ) 
-            CALL iom_close( inum )
+          clfile(1) = filbdy_data_bt_T
+          clfile(2) = filbdy_data_bt_U
+          clfile(3) = filbdy_data_bt_V
+
+          DO igrd = 1,3
+
+            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 
@@ -170,15 +737,14 @@
             CALL ymds2ju(iyear0, imonth0, iday0, real(isec0), dayjul0)
             ! Compute model initialization time 
-            kyear  = ndastp / 10000
-            kmonth = ( ndastp - kyear * 10000 ) / 100
-            kday   = ndastp - kyear * 10000 - kmonth * 100
-            ksecs  = dayfrac * 86400
-            CALL ymds2ju(kyear, kmonth, kday, real(ksecs) , zdayjulini)
-            ! offset from initialization date:
-            offset = (dayjul0-zdayjulini)*86400
+            iyear  = ndastp / 10000
+            imonth = ( ndastp - iyear * 10000 ) / 100
+            iday   = ndastp - iyear * 10000 - imonth * 100
+            isecs  = dayfrac * 86400
+            CALL ymds2ju(iyear, imonth, iday, real(isecs) , 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 
@@ -186,124 +752,120 @@
             !! when calendar in model available outside of IOIPSL.
 
-            IF(lwp) WRITE(numout,*) 'kdby (number of times): ',kbdy
-            IF(lwp) WRITE(numout,*) 'offset: ',offset
-            IF(lwp) WRITE(numout,*) 'totime: ',totime
-            IF(lwp) WRITE(numout,*) 'istepr: ',istepr
-
             ! Check that there are not too many times in the file. 
-            IF (kbdy > jpbtime) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter'
-              IF(lwp) WRITE(numout,*) ' ==========  jpbtime= ',jpbtime,' kbdy= ', kbdy             
-              IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd)
-              nstop = nstop + 1
+            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
+               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_bdy_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 == 1 ) CALL ctl_stop( &
+                    'There is only one time dump in data files', &
+                    'Set ln_bdy_clim=.true. in namelist for constant bdy forcing.' )
+
+               zinterval_s = zstepr(2) - zstepr(1)
+               zinterval_e = zstepr(ntimes_bdy) - zstepr(ntimes_bdy-1)
+
+               IF ( zstepr(1) - zinterval_s / 2.0 > 0 ) THEN	          
+                  IF(lwp) WRITE(numout,*) 'First bdy time relative to nit000:', zstepr(1)
+                  IF(lwp) WRITE(numout,*) 'Interval between first two times: ', zinterval_s
+                  CALL ctl_stop( 'First data time is after start of run', & 
+                       'by more than half a meaning period', &
+                       'Check file: ' // TRIM(clfile(igrd)) )
+               END IF
+
+               IF ( zstepr(ntimes_bdy) + zinterval_e / 2.0 < zntotime ) THEN
+                  IF(lwp) WRITE(numout,*) 'Last bdy time relative to nit000:', zstepr(ntimes_bdy)
+                  IF(lwp) WRITE(numout,*) 'Interval between last two times: ', zinterval_e
+                  CALL ctl_stop( 'Last data time is before end of run', & 
+                       'by more than half a meaning period', &
+                       'Check file: ' // TRIM(clfile(igrd))  )
+               END IF
+
+            END IF ! .NOT. ln_bdy_clim
+
+            IF ( igrd .EQ. 1) THEN
+              ntimes_bdyt = ntimes_bdy_bt
+              zoffsett = zoffset
+              istept(:) = INT( zstepr(:) + zoffset )
+            ELSE IF (igrd .EQ. 2) THEN
+              ntimes_bdyu = ntimes_bdy_bt
+              zoffsetu = zoffset
+              istepu(:) = INT( zstepr(:) + zoffset )
+            ELSE IF (igrd .EQ. 3) THEN
+              ntimes_bdyv = ntimes_bdy_bt
+              zoffsetv = zoffset
+              istepv(:) = INT( zstepr(:) + zoffset )
             ENDIF
 
-            ! Check that time array increases:
-        
-            jt=1
-            DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy-1)))
-              jt=jt+1
-            END DO
-
-            IF ((jt.NE.kbdy).AND.(kbdy.GT.1)) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files'
-              IF(lwp) WRITE(numout,*) ' ==========  does not continuously increase. Check file: '         
-              IF(lwp) WRITE(numout,*) bdyfile(jgrd)
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            ENDIF
-
-            ! Check that times in file span model run time:
- 
-            IF (istepr(1)+offset > 0) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time'
-              IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd)
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            END IF
-
-            IF (istepr(kbdy)+offset < totime) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time'
-              IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd) 
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            END IF
-
-            IF ( jgrd .EQ. 1) THEN
-              kbdyt = kbdy
-              offsett = offset
-              istept(:) = INT( istepr(:) + offset )
-            ELSE IF (jgrd .EQ. 2) THEN
-              kbdyu = kbdy
-              offsetu = offset
-              istepu(:) = INT( istepr(:) + offset )
-            ELSE IF (jgrd .EQ. 3) THEN
-              kbdyv = kbdy
-              offsetv = offset
-              istepv(:) = INT( istepr(:) + offset )
-            ENDIF
-
           ENDDO
 
       ! Verify time consistency between files  
 
-          IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN
-            IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps'
-            IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet'
-            nstop = nstop + 1
+          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
-          kbdy = kbdyt
-
-          IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN
-            IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin'
-            IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet'
-            nstop = nstop + 1
+
+          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
-          offset = offsett
+          zoffset = zoffsett
 
       !! Check that times are the same in the three files... HERE.
-          istep(:) = istept(:)
+          istep_bt(:) = istept(:)
 
       ! Check number of time dumps:              
-          IF ((kbdy.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files'
-            IF(lwp) WRITE(numout,*) ' ==========  Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.'            
-            IF(lwp) WRITE(numout,*)
-            nstop = nstop + 1
-          ENDIF
-
           IF (ln_bdy_clim) THEN
-            IF ((kbdy.NE.1).AND.(kbdy.NE.12)) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),'
-              IF(lwp) WRITE(numout,*) ' ==========  bdy data files must contain 1 or 12 time dumps.'         
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            ELSEIF (kbdy.EQ.1 ) 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,*)             
-            ELSEIF (kbdy.EQ.12) THEN
+            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,*) 
-            ENDIF
+            CASE DEFAULT
+              CALL ctl_stop( &
+                'For climatological boundary forcing (ln_bdy_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). 
 
-          jt=1
-          DO WHILE ( ((istep(jt+1)) <= 0 ).AND.(jt.LE.(kbdy-1)))
-            jt=jt+1
-          END DO
-          nbdy1 = jt
-
-          WRITE(numout,*) 'Time offset is ',offset
-          WRITE(numout,*) 'First record to read is ',nbdy1
+          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
+
+          WRITE(numout,*) 'Time offset is ',zoffset
+          WRITE(numout,*) 'First record to read is ',nbdy_b_bt
 
         ENDIF ! endif (nbdy_dta == 1)
@@ -314,120 +876,85 @@
         IF ( nbdy_dta == 0) THEN
           ! boundary data arrays are filled with initial conditions
-          DO jk = 1, jpkm1 
-
-#if ! defined key_barotropic
-            jgrd = 1            ! T-points data 
-            DO jb = 1, nblen(jgrd)              
-              tbdy(jb,jk) = tn(nbi(jb,jgrd), nbj(jb,jgrd), jk)
-              sbdy(jb,jk) = sn(nbi(jb,jgrd), nbj(jb,jgrd), jk)
-            END DO
-#endif
-
-            jgrd = 2            ! U-points data 
-            DO jb = 1, nblen(jgrd)              
-              ubdy(jb,jk) = un(nbi(jb,jgrd), nbj(jb,jgrd), jk)
-            END DO
-
-            jgrd = 3            ! V-points data 
-            DO jb = 1, nblen(jgrd)              
-              vbdy(jb,jk) = vn(nbi(jb,jgrd), nbj(jb,jgrd), jk)
-            END DO
-          END DO 
+          igrd = 2            ! U-points data 
+          DO ib = 1, nblen(igrd)              
+            ubtbdy(ib) = un(nbi(ib,igrd), nbj(ib,igrd), 1)
+          END DO
+
+          igrd = 3            ! V-points data 
+          DO ib = 1, nblen(igrd)              
+            vbtbdy(ib) = vn(nbi(ib,igrd), nbj(ib,igrd), 1)
+          END DO
+
+          igrd = 1            ! T-points data 
+          DO ib = 1, nblen(igrd)              
+            sshbdy(ib) = sshn(nbi(ib,igrd), nbj(ib,igrd))
+          END DO
 
         ELSEIF (nbdy_dta == 1) THEN
  
         ! Set first record in the climatological case:   
-          IF ((ln_bdy_clim).AND.(kbdy==1)) THEN
-            nbdy2 = 1
-          ELSEIF ((ln_bdy_clim).AND.(kbdy==iman)) THEN
-            nbdy1 = 0
-            nbdy2 = imois
+          IF ((ln_bdy_clim).AND.(ntimes_bdy_bt==1)) THEN
+            nbdy_a_bt = 1
+          ELSEIF ((ln_bdy_clim).AND.(ntimes_bdy_bt==iman)) THEN
+            nbdy_b_bt = 0
+            nbdy_a_bt = imois
           ELSE
-            nbdy2 = nbdy1
+            nbdy_a_bt = nbdy_b_bt
           END IF
  
          ! Open Netcdf files:
 
-          CALL iom_open ( filbdy_data_T, bdynumt )
-          CALL iom_open ( filbdy_data_U, bdynumu )
-          CALL iom_open ( filbdy_data_V, bdynumv )
+          CALL iom_open ( filbdy_data_bt_T, numbdyt_bt )
+          CALL iom_open ( filbdy_data_bt_U, numbdyu_bt )
+          CALL iom_open ( filbdy_data_bt_V, numbdyv_bt )
 
          ! Read first record:
           ipj=1
-          ipk=jpk
-          jgrd=1
-          ipi=nblendta(jgrd)
-
-#if ! defined key_barotropic
-          !temperature
-
-          jgrd=1
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumt,'votemper' )
-            nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar)
+          igrd=1
+          ipi=nblendta(igrd)
+
+          ! ssh
+          igrd=1
+          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 votemper is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              tbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-
-          ! salinity
-
-          jgrd=1
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumt,'vosaline' )
-            nblendta(jgrd) = iom_file(bdynumt)%dimsz(1,idvar)
+          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=2
+          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 vosaline is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              sbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-#endif
-          
-          ! u-velocity
-
-          jgrd=2
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumu,'vozocrtx' )
-            nblendta(jgrd) = iom_file(bdynumu)%dimsz(1,idvar)
+          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=3
+          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 vozocrtx is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              ubdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-
-          ! v-velocity
-          jgrd=3
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumv,'vomecrty' )
-            nblendta(jgrd) = iom_file(bdynumv)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for vomecrty is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              vbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
+          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
 
@@ -435,15 +962,14 @@
  
         ! In the case of constant boundary forcing fill bdy arrays once for all
-        IF ((ln_bdy_clim).AND.(kbdy==1)) THEN
-#if ! defined key_barotropic
-          tbdy  (:,:) = tbdydta  (:,:,2)
-          sbdy  (:,:) = sbdydta  (:,:,2)
-#endif
-          ubdy  (:,:) = ubdydta  (:,:,2)
-          vbdy  (:,:) = vbdydta  (:,:,2)
-
-          CALL iom_close( bdynumt )
-          CALL iom_close( bdynumu )
-          CALL iom_close( bdynumv )
+        IF ((ln_bdy_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
 
@@ -454,534 +980,21 @@
       ! -------------------- !
 
-      IF ((nbdy_dta==1).AND.(kbdy>1)) THEN 
+      IF ((nbdy_dta==1).AND.(ntimes_bdy_bt>1)) THEN 
 
       ! 2.1 Read one more record if necessary
       !**************************************
 
-        IF ( (ln_bdy_clim).AND.(imois/=nbdy1) ) THEN ! remember that nbdy1=0 for kt=nit000
-         nbdy1 = imois
-         nbdy2 = imois+1
-         nbdy1 = MOD( nbdy1, iman )
-         IF( nbdy1 == 0 ) nbdy1 = iman
-         nbdy2 = MOD( nbdy2, iman )
-         IF( nbdy2 == 0 ) nbdy2 = iman
+        IF ( (ln_bdy_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_bdy_clim).AND.(itimer >= istep(nbdy2))) THEN
-          nbdy1=nbdy2
-          nbdy2=nbdy2+1
-          lect=.true.
-        END IF
-         
-        IF (lect) THEN
-
-        ! Swap arrays
-#if ! defined key_barotropic
-          tbdydta(:,:,1) =  tbdydta(:,:,2)
-          sbdydta(:,:,1) =  sbdydta(:,:,2)
-#endif
-          ubdydta(:,:,1) =  ubdydta(:,:,2)
-          vbdydta(:,:,1) =  vbdydta(:,:,2)
- 
-        ! read another set
-
-          ipj=1
-          ipk=jpk
-          jgrd=1
-          ipi=nblendta(jgrd)
-
-#if ! defined key_barotropic
-          !temperature
- 
-          jgrd=1
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumt, jpdom_unknown,'votemper',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              tbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-
-          ! salinity
-
-          jgrd=1
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumt, jpdom_unknown,'vosaline',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              sbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-#endif
-          
-          ! u-velocity
-
-          jgrd=2
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumu, jpdom_unknown,'vozocrtx',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              ubdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-
-          ! v-velocity
-          jgrd=3
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumv, jpdom_unknown,'vomecrty',pdta3(1:ipi,1:ipj,1:ipk),nbdy2 )
-
-          DO jb=1, nblen(jgrd)
-            DO jk=1,jpkm1
-              vbdydta(jb,jk,2) =  pdta3(nbmap(jb,jgrd),1,jk)
-            END DO
-          END DO
-
-         IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy1 ',nbdy1
-         IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy2 ',nbdy2
-         IF (.NOT.ln_bdy_clim) THEN
-           IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep(nbdy1)
-           IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer
-           IF(lwp) WRITE(numout,*) 'second record time (s): ', istep(nbdy2)
-         ENDIF
-        END IF ! end lect=.true.
-
-
-      ! 2.2   Interpolate linearly:
-      ! ***************************
-    
-        IF (ln_bdy_clim) THEN
-          zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1) ) + 0.5 - i15
-        ELSE          
-          zxy = FLOAT(istep(nbdy1)-itimer) / FLOAT(istep(nbdy1)-istep(nbdy2))
-        END IF
-
-#if ! defined key_barotropic
-          jgrd=1
-          DO jb=1, nblen(jgrd)
-            DO jk=1, jpkm1
-              tbdy(jb,jk) = zxy      * tbdydta(jb,jk,2) + &
-                            (1.-zxy) * tbdydta(jb,jk,1)
-              sbdy(jb,jk) = zxy      * sbdydta(jb,jk,2) + &
-                            (1.-zxy) * sbdydta(jb,jk,1)
-            END DO
-          END DO
-#endif
-
-          jgrd=2
-          DO jb=1, nblen(jgrd)
-            DO jk=1, jpkm1
-              ubdy(jb,jk) = zxy      * ubdydta(jb,jk,2) + &
-                            (1.-zxy) * ubdydta(jb,jk,1)   
-            END DO
-          END DO
-
-          jgrd=3
-          DO jb=1, nblen(jgrd)
-            DO jk=1, jpkm1
-              vbdy(jb,jk) = zxy      * vbdydta(jb,jk,2) + &
-                            (1.-zxy) * vbdydta(jb,jk,1)   
-            END DO
-          END DO
-
-      END IF !end if ((nbdy_dta==1).AND.(kbdy>1))
-    
-      ! ------------------- !
-      ! Last call kt=nitend !
-      ! ------------------- !
-
-      ! Closing of the 3 files
-      IF( kt == nitend ) THEN
-          CALL iom_close( bdynumt )
-          CALL iom_close( bdynumu )
-          CALL iom_close( bdynumv )
-      ENDIF
-
-   END SUBROUTINE bdy_dta
-
-   SUBROUTINE bdy_dta_bt( kt, jit )
-      !!---------------------------------------------------------------------------
-      !!                      ***  SUBROUTINE bdy_dta_bt  ***
-      !!                    
-      !! ** Purpose :   Read unstructured boundary data for barotropic variables
-      !!
-      !! ** Method  :   
-      !!
-      !! History : OPA 9.0  ! 07-07 (D. Storkey) Original
-      !!---------------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
-      INTEGER, INTENT( in ) ::   jit         ! barotropic time step index
-                                             ! (for timesplitting option, otherwise zero)
-
-      !! * Local declarations
-      LOGICAL ::   lect                      ! flag for reading
-      INTEGER ::   jt, jb, jgrd              ! dummy loop indices
-      INTEGER ::   idvar                     ! netcdf var ID
-      INTEGER ::   iman, i15, imois          ! Time variables for monthly clim forcing
-      INTEGER ::   kbdyt, kbdyu, kbdyv
-      INTEGER ::   itimer, totime
-      INTEGER ::   ipi, ipj, ipk, inum       ! temporary integers (NetCDF read)
-      INTEGER ::   iyear0, imonth0, iday0
-      INTEGER ::   ihours0, iminutes0, isec0
-      INTEGER ::   kyear, kmonth, kday, ksecs
-      INTEGER, DIMENSION(jpbtime) ::   istept, istepu, istepv   ! time arrays from data files
-      REAL(wp) ::   dayfrac, zxy, offsett
-      REAL(wp) ::   offsetu, offsetv
-      REAL(wp) ::   dayjul0, zdayjulini
-      REAL(wp), DIMENSION(jpbtime)      ::   istepr             ! REAL time array from data files
-      REAL(wp), DIMENSION(jpbdta,1)     ::   pdta2              ! temporary array for data fields
-      REAL(wp), DIMENSION(jpbdta,1,jpk) ::   pdta3              ! temporary array for data fields
-      CHARACTER(LEN=80), DIMENSION(3)   ::   bdyfile
-      CHARACTER(LEN=70 )                ::   clunits            ! units attribute of time coordinate
-
-      !!---------------------------------------------------------------------------
-
-      ! -------------------- !
-      !    Initialization    !
-      ! -------------------- !
-
-      lect   = .false.           ! If true, read a time record
-
-      ! Some time variables for monthly climatological forcing:
-      ! *******************************************************
- 
-      iman  = INT( raamo ) ! Number of months in a year
-
-      i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 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*rdtbt         ! in non-climatological case
-
-      ! -------------------------------------!
-      ! Update BDY fields with tidal forcing !
-      ! -------------------------------------!  
-
-      IF ( lk_bdy_tides ) CALL tide_update( kt, jit ) 
-  
-      IF ( lk_bdy ) THEN 
-
-      ! -------------------- !
-      ! 1.   First call only !
-      ! -------------------- !
-
-      IF ( kt == nit000 ) THEN
-
-        istep_bt(:) = 0
-
-        nbdy1_bt=0
-        nbdy2_bt=0
-
-      ! 1.1  Get time information from bdy data file
-      ! ********************************************
-
-        IF(lwp) WRITE(numout,*)
-        IF(lwp) WRITE(numout,*)    'bdy_dta_bt :Initialize unstructured boundary data for barotropic variables.'
-        IF(lwp) WRITE(numout,*)    '~~~~~~~' 
-
-        IF (nbdy_dta == 0) THEN
-          IF(lwp) WRITE(numout,*)  'Bdy data are taken from initial conditions'
-
-        ELSEIF (nbdy_dta == 1) THEN
-          IF(lwp) WRITE(numout,*)  'Bdy data are read in netcdf files'
-
-          dayfrac = adatrj-float(itimer)/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
-
-          bdyfile(1) = filbdy_data_bt_T
-          bdyfile(2) = filbdy_data_bt_U
-          bdyfile(3) = filbdy_data_bt_V
-
-          DO jgrd = 1,3
-
-            CALL iom_open( bdyfile(jgrd), inum )
-            CALL iom_gettime( inum, istepr, kntime=kbdy, cdunits=clunits ) 
-            CALL iom_close( inum )
-
-            ! 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), dayjul0)
-            ! Compute model initialization time 
-            kyear  = ndastp / 10000
-            kmonth = ( ndastp - kyear * 10000 ) / 100
-            kday   = ndastp - kyear * 10000 - kmonth * 100
-            ksecs  = dayfrac * 86400
-            CALL ymds2ju(kyear, kmonth, kday, real(ksecs) , zdayjulini)
-            ! offset from initialization date:
-            offset = (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,*) 'kdby (number of times): ',kbdy_bt
-            IF(lwp) WRITE(numout,*) 'offset: ',offset
-            IF(lwp) WRITE(numout,*) 'totime: ',totime
-            IF(lwp) WRITE(numout,*) 'istepr: ',istepr
-
-            ! Check that there are not too many times in the file. 
-            IF (kbdy_bt > jpbtime) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : Number of time dumps in files exceed jpbtime parameter'
-              IF(lwp) WRITE(numout,*) ' ==========  jpbtime= ',jpbtime,' kbdy_bt= ', kbdy_bt             
-              IF(lwp) WRITE(numout,*) ' Check file: ', bdyfile(jgrd)
-              nstop = nstop + 1
-            ENDIF
-
-            ! Check that time array increases:
-        
-            jt=1
-            DO WHILE ((istepr(jt+1).GT.istepr(jt)).AND.(jt.LE.(kbdy_bt-1)))
-              jt=jt+1
-            END DO
-
-            IF ((jt.NE.kbdy_bt).AND.(kbdy_bt.GT.1)) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : Time array in unstructured boundary data files'
-              IF(lwp) WRITE(numout,*) ' ==========  does not continuously increase. Check file: '         
-              IF(lwp) WRITE(numout,*) bdyfile(jgrd)
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            ENDIF
-
-            ! Check that times in file span model run time:
- 
-            IF (istepr(1)+offset > 0) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : First time dump in bdy file is after model initial time'
-              IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd)
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            END IF
-
-            IF (istepr(kbdy_bt)+offset < totime) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : Last time dump in bdy file is before model final time'
-              IF(lwp) WRITE(numout,*) ' ==========  Check file: ',bdyfile(jgrd) 
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            END IF
-
-            IF ( jgrd .EQ. 1) THEN
-              kbdyt = kbdy_bt
-              offsett = offset
-              istept(:) = INT( istepr(:) + offset )
-            ELSE IF (jgrd .EQ. 2) THEN
-              kbdyu = kbdy_bt
-              offsetu = offset
-              istepu(:) = INT( istepr(:) + offset )
-            ELSE IF (jgrd .EQ. 3) THEN
-              kbdyv = kbdy_bt
-              offsetv = offset
-              istepv(:) = INT( istepr(:) + offset )
-            ENDIF
-
-          ENDDO
-
-      ! Verify time consistency between files  
-
-          IF (kbdyu.NE.kbdyt .OR. kbdyv.NE.kbdyt) THEN
-            IF(lwp) WRITE(numout,*) 'Bdy data files must have the same number of time dumps'
-            IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet'
-            nstop = nstop + 1
-          ENDIF
-          kbdy_bt = kbdyt
-
-          IF (offsetu.NE.offsett .OR. offsetv.NE.offsett) THEN
-            IF(lwp) WRITE(numout,*) 'Bdy data files must have the same time origin'
-            IF(lwp) WRITE(numout,*) 'Multiple time frequencies not implemented yet'
-            nstop = nstop + 1
-          ENDIF
-          offset = offsett
-
-      !! Check that times are the same in the three files... HERE.
-          istep_bt(:) = istept(:)
-
-      ! Check number of time dumps:              
-          IF ((kbdy_bt.EQ.1).AND.(.NOT.ln_bdy_clim)) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) ' E R R O R : There is only one time dump in data files'
-            IF(lwp) WRITE(numout,*) ' ==========  Choose ln_bdy_clim=.true. in namelist for constant bdy forcing.'            
-            IF(lwp) WRITE(numout,*)
-            nstop = nstop + 1
-          ENDIF
-
-          IF (ln_bdy_clim) THEN
-            IF ((kbdy_bt.NE.1).AND.(kbdy_bt.NE.12)) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) ' E R R O R : For climatological boundary forcing (ln_bdy_clim=.true.),'
-              IF(lwp) WRITE(numout,*) ' ==========  bdy data files must contain 1 or 12 time dumps.'         
-              IF(lwp) WRITE(numout,*)
-              nstop = nstop + 1
-            ELSEIF (kbdy_bt.EQ.1 ) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) 'We assume constant boundary forcing from bdy data files'
-              IF(lwp) WRITE(numout,*)             
-            ELSEIF (kbdy_bt.EQ.12) THEN
-              IF(lwp) WRITE(numout,*)
-              IF(lwp) WRITE(numout,*) 'We assume monthly (and cyclic) boundary forcing from bdy data files'
-              IF(lwp) WRITE(numout,*) 
-            ENDIF
-          ENDIF
-
-      ! Find index of first record to read (before first model time). 
-
-          jt=1
-          DO WHILE ( ((istep_bt(jt+1)) <= 0 ).AND.(jt.LE.(kbdy_bt-1)))
-            jt=jt+1
-          END DO
-          nbdy1_bt = jt
-
-          WRITE(numout,*) 'Time offset is ',offset
-          WRITE(numout,*) 'First record to read is ',nbdy1_bt
-
-        ENDIF ! endif (nbdy_dta == 1)
-
-      ! 1.2  Read first record in file if necessary (ie if nbdy_dta == 1)
-      ! *****************************************************************
-
-        IF ( nbdy_dta == 0) THEN
-          ! boundary data arrays are filled with initial conditions
-          jgrd = 2            ! U-points data 
-          DO jb = 1, nblen(jgrd)              
-            ubtbdy(jb) = un(nbi(jb,jgrd), nbj(jb,jgrd), 1)
-          END DO
-
-          jgrd = 3            ! V-points data 
-          DO jb = 1, nblen(jgrd)              
-            vbtbdy(jb) = vn(nbi(jb,jgrd), nbj(jb,jgrd), 1)
-          END DO
-
-          jgrd = 1            ! T-points data 
-          DO jb = 1, nblen(jgrd)              
-            sshbdy(jb) = sshn(nbi(jb,jgrd), nbj(jb,jgrd))
-          END DO
-
-        ELSEIF (nbdy_dta == 1) THEN
- 
-        ! Set first record in the climatological case:   
-          IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN
-            nbdy2_bt = 1
-          ELSEIF ((ln_bdy_clim).AND.(kbdy_bt==iman)) THEN
-            nbdy1_bt = 0
-            nbdy2_bt = imois
-          ELSE
-            nbdy2_bt = nbdy1_bt
-          END IF
- 
-         ! Open Netcdf files:
-
-          CALL iom_open ( filbdy_data_bt_T, bdynumt_bt )
-          CALL iom_open ( filbdy_data_bt_U, bdynumu_bt )
-          CALL iom_open ( filbdy_data_bt_V, bdynumv_bt )
-
-         ! Read first record:
-          ipj=1
-          jgrd=1
-          ipi=nblendta(jgrd)
-
-          ! ssh
-          jgrd=1
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumt_bt,'sossheig' )
-            nblendta(jgrd) = iom_file(bdynumt_bt)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for sossheig is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt )
-
-          DO jb=1, nblen(jgrd)
-            sshbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1)
-          END DO
- 
-          ! u-velocity
-          jgrd=2
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumu_bt,'vobtcrtx' )
-            nblendta(jgrd) = iom_file(bdynumu_bt)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for vobtcrtx is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt )
-
-          DO jb=1, nblen(jgrd)
-            ubtbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1)
-          END DO
-
-          ! v-velocity
-          jgrd=3
-          IF ( nblendta(jgrd) .le. 0 ) THEN 
-            idvar = iom_varid( bdynumv_bt,'vobtcrty' )
-            nblendta(jgrd) = iom_file(bdynumv_bt)%dimsz(1,idvar)
-          ENDIF
-          WRITE(numout,*) 'Dim size for vobtcrty is ',nblendta(jgrd)
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt )
-
-          DO jb=1, nblen(jgrd)
-            vbtbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1)
-          END DO
-
-        END IF
- 
-        ! In the case of constant boundary forcing fill bdy arrays once for all
-        IF ((ln_bdy_clim).AND.(kbdy_bt==1)) THEN
-
-          ubtbdy  (:) = ubtbdydta  (:,2)
-          vbtbdy  (:) = vbtbdydta  (:,2)
-          sshbdy  (:) = sshbdydta  (:,2)
-
-          CALL iom_close( bdynumt_bt )
-          CALL iom_close( bdynumu_bt )
-          CALL iom_close( bdynumv_bt )
-
-        END IF
-
-      ENDIF ! End if nit000
-
-      ! -------------------- !
-      ! 2. At each time step !
-      ! -------------------- !
-
-      IF ((nbdy_dta==1).AND.(kbdy_bt>1)) THEN 
-
-      ! 2.1 Read one more record if necessary
-      !**************************************
-
-        IF ( (ln_bdy_clim).AND.(imois/=nbdy1_bt) ) THEN ! remember that nbdy1_bt=0 for kt=nit000
-         nbdy1_bt = imois
-         nbdy2_bt = imois+1
-         nbdy1_bt = MOD( nbdy1_bt, iman )
-         IF( nbdy1_bt == 0 ) nbdy1_bt = iman
-         nbdy2_bt = MOD( nbdy2_bt, iman )
-         IF( nbdy2_bt == 0 ) nbdy2_bt = iman
-         lect=.true.
-
-        ELSEIF ((.NOT.ln_bdy_clim).AND.(itimer >= istep_bt(nbdy2_bt))) THEN
-          nbdy1_bt=nbdy2_bt
-          nbdy2_bt=nbdy2_bt+1
+        ELSEIF ((.NOT.ln_bdy_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
@@ -998,45 +1011,45 @@
           ipj=1
           ipk=jpk
-          jgrd=1
-          ipi=nblendta(jgrd)
+          igrd=1
+          ipi=nblendta(igrd)
 
           
           ! ssh
-          jgrd=1
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumt_bt, jpdom_unknown,'sossheig',pdta2(1:ipi,1:ipj),nbdy2_bt )
-
-          DO jb=1, nblen(jgrd)
-            sshbdydta(jb,2) =  pdta2(nbmap(jb,jgrd),1)
+          igrd=1
+          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
-          jgrd=2
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumu_bt, jpdom_unknown,'vobtcrtx',pdta2(1:ipi,1:ipj),nbdy2_bt )
-
-          DO jb=1, nblen(jgrd)
-            ubtbdydta(jb,2) =  pdta3(nbmap(jb,jgrd),1)
+          igrd=2
+          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
-          jgrd=3
-          ipi=nblendta(jgrd)
-
-          CALL iom_get ( bdynumv_bt, jpdom_unknown,'vobtcrty',pdta2(1:ipi,1:ipj),nbdy2_bt )
-
-          DO jb=1, nblen(jgrd)
-            vbtbdydta(jb,2) =  pdta3(nbmap(jb,jgrd),1)
-          END DO
-
-
-         IF(lwp) WRITE(numout,*) 'bdy_dta : first record file used nbdy1_bt ',nbdy1_bt
-         IF(lwp) WRITE(numout,*) '~~~~~~~~  last  record file used nbdy2_bt ',nbdy2_bt
+          igrd=3
+          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 : 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_bdy_clim) THEN
-           IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy1_bt)
+           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(nbdy2_bt)
+           IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy_a_bt)
          ENDIF
         END IF ! end lect=.true.
@@ -1047,29 +1060,29 @@
     
         IF (ln_bdy_clim) THEN
-          zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1_bt) ) + 0.5 - i15
+          zxy = FLOAT( nday ) / FLOAT( nobis(nbdy_b_bt) ) + 0.5 - i15
         ELSE          
-          zxy = FLOAT(istep_bt(nbdy1_bt)-itimer) / FLOAT(istep_bt(nbdy1_bt)-istep_bt(nbdy2_bt))
+          zxy = FLOAT(istep_bt(nbdy_b_bt)-itimer) / FLOAT(istep_bt(nbdy_b_bt)-istep_bt(nbdy_a_bt))
         END IF
 
-          jgrd=1
-          DO jb=1, nblen(jgrd)
-            sshbdy(jb) = zxy      * sshbdydta(jb,2) + &
-                       (1.-zxy) * sshbdydta(jb,1)   
-          END DO
-
-          jgrd=2
-          DO jb=1, nblen(jgrd)
-            ubtbdy(jb) = zxy      * ubtbdydta(jb,2) + &
-                         (1.-zxy) * ubtbdydta(jb,1)   
-          END DO
-
-          jgrd=3
-          DO jb=1, nblen(jgrd)
-            vbtbdy(jb) = zxy      * vbtbdydta(jb,2) + &
-                         (1.-zxy) * vbtbdydta(jb,1)   
-          END DO
-
-
-      END IF !end if ((nbdy_dta==1).AND.(kbdy_bt>1))
+          igrd=1
+          DO ib=1, nblen(igrd)
+            sshbdy(ib) = zxy      * sshbdydta(ib,2) + &
+                       (1.-zxy) * sshbdydta(ib,1)   
+          END DO
+
+          igrd=2
+          DO ib=1, nblen(igrd)
+            ubtbdy(ib) = zxy      * ubtbdydta(ib,2) + &
+                         (1.-zxy) * ubtbdydta(ib,1)   
+          END DO
+
+          igrd=3
+          DO ib=1, nblen(igrd)
+            vbtbdy(ib) = zxy      * vbtbdydta(ib,2) + &
+                         (1.-zxy) * vbtbdydta(ib,1)   
+          END DO
+
+
+      END IF !end if ((nbdy_dta==1).AND.(ntimes_bdy_bt>1))
     
       ! ------------------- !
@@ -1079,10 +1092,10 @@
       ! Closing of the 3 files
       IF( kt == nitend ) THEN
-          CALL iom_close( bdynumt_bt )
-          CALL iom_close( bdynumu_bt )
-          CALL iom_close( bdynumv_bt )
+          CALL iom_close( numbdyt_bt )
+          CALL iom_close( numbdyu_bt )
+          CALL iom_close( numbdyv_bt )
       ENDIF
 
-      ENDIF ! lk_bdy
+      ENDIF ! ln_bdy_dyn_frs
 
       END SUBROUTINE bdy_dta_bt
@@ -1090,16 +1103,13 @@
 
 #else
-   !!------------------------------------------------------------------------------
-   !!   default option:           Dummy module NO Unstruct Open Boundary Conditions
-   !!------------------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   Dummy module                   NO Unstruct Open Boundary Conditions
+   !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE bdy_dta( kt )             ! Dummy routine
-      INTEGER, INTENT (in) :: kt
+   SUBROUTINE bdy_dta( kt )              ! Empty routine
       WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dta
-
-   SUBROUTINE bdy_dta_bt( kt, jit )             ! Dummy routine
-      INTEGER, INTENT (in) :: kt, jit
-      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, jit
+   SUBROUTINE bdy_dta_bt( kt, kit )      ! Empty routine
+      WRITE(*,*) 'bdy_dta: You should not have seen this print! error?', kt, kit
    END SUBROUTINE bdy_dta_bt
 #endif
