Index: /trunk/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 911)
@@ -0,0 +1,112 @@
+MODULE bdy_oce
+   !!======================================================================
+   !!                       ***  MODULE bdy_oce   ***
+   !! Unstructured Open Boundary Cond. :   define related variables
+   !!======================================================================
+#if defined key_bdy || key_bdy_tides
+   !!----------------------------------------------------------------------
+   !!   'key_bdy' :                    Unstructured Open Boundary Condition
+   !!----------------------------------------------------------------------
+   !! History :
+   !!  9.0   01/05   (J. Chanut, A. Sellar)  Original code
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_oce         ! ocean parameters
+   USE bdy_par         ! Unstructured boundary parameters
+
+   IMPLICIT NONE
+   PUBLIC
+
+   !!----------------------------------------------------------------------
+   !! Namelist variables
+   !!----------------------------------------------------------------------
+
+   CHARACTER(len=80) :: &
+      filbdy_mask, &         !: Name of unstruct. bdy mask file
+      filbdy_data_T, &       !: Name of unstruct. bdy data file at T points
+      filbdy_data_U, &       !: Name of unstruct. bdy data file at U points
+      filbdy_data_V, &       !: Name of unstruct. bdy data file at V points
+      filbdy_data_bt_T, &    !: Name of unstruct. bdy data file at T points for barotropic variables
+      filbdy_data_bt_U, &    !: Name of unstruct. bdy data file at U points for barotropic variables
+      filbdy_data_bt_V       !: Name of unstruct. bdy data file at V points for barotropic variables
+
+   LOGICAL ::  &
+      ln_bdy_clim =.false.!: if true, we assume that bdy data files contain 
+                          !  1 time dump (-->bdy forcing will be constant) 
+                          !  or 12 months  (-->bdy forcing will be cyclic) 
+   LOGICAL ::  &
+      ln_bdy_fla =.false. !: if true, flather boundary conditions
+
+   LOGICAL ::  &
+      ln_bdy_vol =.false. !: if true, volume correction             
+
+   LOGICAL ::  &
+      ln_bdy_mask = .false. !: if true, read bdymask from file
+
+   INTEGER ::  &
+      nb_rimwidth = 7,	& !: boundary rim width
+      nbdy_dta = 1	  !: = 0 use the initial state as bdy dta
+                          !: = 1 read bdy data in netcdf file
+   INTEGER ::  & 
+      volbdy = 1          !: = 0 the total volume will have the variability 
+                          !       of the surface Flux E-P else (volbdy = 1) 
+                          !       the volume will be constant
+                          !  = 1 the volume will be constant during all the 
+                          !       integration.
+
+   !!----------------------------------------------------------------------
+   !! Global variables
+   !!----------------------------------------------------------------------
+
+   REAL(wp), DIMENSION(jpi, jpj) ::  & !:
+      bdytmask, &  !: Mask defining computational domain at T-points
+      bdyumask, &  !: Mask defining computational domain at U-points
+      bdyvmask     !: Mask defining computational domain at V-points
+
+   !!----------------------------------------------------------------------
+   !! Unstructured open boundary data variables
+   !!----------------------------------------------------------------------
+
+   INTEGER, DIMENSION(jpbgrd) ::  &
+      nblen  = 0,	& !: Size of bdy data on a proc for each grid type
+      nblenrim = 0,	& !: Size of bdy data on a proc for first rim ind
+      nblendta = 0        !: Size of bdy data in file
+
+   INTEGER, DIMENSION(jpbdim, jpbgrd) ::  & 
+      nbi, nbj,		& !: i and j indices of bdy dta
+      nbr,		& !: Discrete distance from rim points
+      nbmap		  !: Indices of data in file for data in memory 
+    
+   REAL(wp) :: &
+      bdysurftot	  !: Lateral surface of unstructured open boundary
+
+   REAL(wp), DIMENSION(jpbdim) ::  &
+      flagu,		& !: Flag for normal velocity compnt for u-velocity
+      flagv		  !: Flag for normal velocity compnt for v-velocity
+
+   REAL(wp), DIMENSION(jpbdim, jpbgrd) ::  &
+      nbw		  !: Rim weights of bdy data
+
+   REAL(wp), DIMENSION(jpbdim) ::  &
+      sshbdy,		& !: Now clim of bdy sea surface height (Flather)
+      ubtbdy, vbtbdy	  !: Now clim of bdy barotropic velocity components
+
+   REAL(wp), DIMENSION(jpbdim,jpk) ::  &
+      tbdy, sbdy,	& !: Now clim of bdy temperature and salinity  
+      ubdy, vbdy	  !: Now clim of bdy velocity components
+
+#if defined key_bdy_tides
+   REAL(wp), DIMENSION(jpbdim) ::  &
+      sshtide,          & !: Tidal boundary array : SSH
+      utide,            & !: Tidal boundary array : U
+      vtide	          !: Tidal boundary array : V
+#endif
+		    
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                                       Empty module
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE bdy_oce
Index: /trunk/NEMO/OPA_SRC/BDY/bdy_par.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 911)
@@ -0,0 +1,47 @@
+MODULE bdy_par
+   !!======================================================================
+   !!                  ***  MODULE bdy_par   ***
+   !! Unstructured Open Boundary Cond. :   define related parameters
+   !!======================================================================
+#if defined key_bdy || key_bdy_tides
+   !!----------------------------------------------------------------------
+   !!   'key_bdy' :                    Unstructured Open Boundary Condition
+   !!----------------------------------------------------------------------
+   !! history :
+   !!  9.0   01/05   (J. Chanut, A. Sellar)  Original code
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE par_oce         ! ocean parameters
+
+   IMPLICIT NONE
+   PUBLIC
+
+#if defined key_bdy
+   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .TRUE. !: Unstructured Ocean 
+                                                   !Boundary Condition flag
+#else
+   ! tides only at boundaries
+   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .FALSE. !: Unstructured Ocean 
+                                                    !Boundary Condition flag
+#endif
+
+   !!----------------------------------------------------------------------
+   !! Unstructured open boundary parameters
+   !!----------------------------------------------------------------------
+
+   INTEGER, PARAMETER ::   & 
+      jpbdta     = 5000,  & !: Max length of bdy field in file
+      jpbdim     = 5000,  & !: Max length of bdy field on a processor
+      jpbtime    = 1000,  & !: Max number of time dumps per file
+      jpbgrd     = 3	    !: Number of horizontal grid types used 
+                            !  (T, u, v, f)
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :                         NO open boundary condition
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy = .FALSE.!: Unstructured Ocean 
+                                                   !Boundary Condition flag
+#endif
+
+   !!======================================================================
+END MODULE bdy_par
Index: /trunk/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 911)
@@ -0,0 +1,1108 @@
+MODULE bdydta
+   !!==============================================================================
+   !!                            ***  MODULE bdydta  ***
+   !! Open boundary data : read the data for the unstructured open boundaries.
+   !!==============================================================================
+#if defined key_bdy || key_bdy_tides
+   !!------------------------------------------------------------------------------
+   !!   'key_bdy'         :                   Unstructured Open Boundary Conditions
+   !!------------------------------------------------------------------------------
+   !!   bdy_dta           : read u, v, t, s data along each open boundary
+   !!------------------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers
+   USE dom_oce         ! ocean space and time domain
+   USE daymod          ! calendar
+   USE phycst          ! physical constants
+   USE bdy_oce         ! ocean open boundary conditions
+   USE bdytides        ! tidal forcing at boundaries
+   USE iom
+   USE ioipsl
+   USE in_out_manager  ! I/O logical units
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC bdy_dta                         ! routines called by step.F90
+   PUBLIC bdy_dta_bt 
+   
+   !! * 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 )
+      !!---------------------------------------------------------------------------
+      !!                      ***  SUBROUTINE bdy_dta  ***
+      !!                    
+      !! ** Purpose :   Read unstructured boundary data
+      !!
+      !! ** Method  :   
+      !!
+      !! 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 ::   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,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-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
+      ! ********************************************
+
+        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
+
+          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 )
+
+            ! 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
+            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
+            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
+          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
+          ENDIF
+          offset = offsett
+
+      !! Check that times are the same in the three files... HERE.
+          istep(:) = 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
+              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
+              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(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
+
+        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
+          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 
+
+        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
+          ELSE
+            nbdy2 = nbdy1
+          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 )
+
+         ! 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)
+          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)
+          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)
+          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
+          END DO
+
+        END IF
+ 
+        ! 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 )
+        END IF
+
+      ENDIF ! End if nit000
+
+      ! -------------------- !
+      ! 2. At each time step !
+      ! -------------------- !
+
+      IF ((nbdy_dta==1).AND.(kbdy>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
+         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
+          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
+          jgrd=1
+          ipi=nblendta(jgrd)
+
+          
+          ! 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)
+          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)
+          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
+         IF (.NOT.ln_bdy_clim) THEN
+           IF(lwp) WRITE(numout,*) 'first  record time (s): ', istep_bt(nbdy1_bt)
+           IF(lwp) WRITE(numout,*) 'model time (s)        : ', itimer
+           IF(lwp) WRITE(numout,*) 'second record time (s): ', istep_bt(nbdy2_bt)
+         ENDIF
+        END IF ! end lect=.true.
+
+
+      ! 2.2   Interpolate linearly:
+      ! ***************************
+    
+        IF (ln_bdy_clim) THEN
+          zxy = FLOAT( nday ) / FLOAT( nobis(nbdy1_bt) ) + 0.5 - i15
+        ELSE          
+          zxy = FLOAT(istep_bt(nbdy1_bt)-itimer) / FLOAT(istep_bt(nbdy1_bt)-istep_bt(nbdy2_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))
+    
+      ! ------------------- !
+      ! Last call kt=nitend !
+      ! ------------------- !
+
+      ! Closing of the 3 files
+      IF( kt == nitend ) THEN
+          CALL iom_close( bdynumt_bt )
+          CALL iom_close( bdynumu_bt )
+          CALL iom_close( bdynumv_bt )
+      ENDIF
+
+      ENDIF ! lk_bdy
+
+      END SUBROUTINE bdy_dta_bt
+
+
+#else
+   !!------------------------------------------------------------------------------
+   !!   default option:           Dummy module NO Unstruct Open Boundary Conditions
+   !!------------------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_dta( kt )             ! Dummy routine
+      INTEGER, INTENT (in) :: kt
+      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
+   END SUBROUTINE bdy_dta_bt
+#endif
+
+   !!==============================================================================
+END MODULE bdydta
Index: /trunk/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 911)
@@ -0,0 +1,207 @@
+MODULE bdydyn
+   !!=================================================================================
+   !!                       ***  MODULE  bdydyn  ***
+   !! Ocean dynamics:   Flow relaxation scheme of velocities on unstruc. open boundary
+   !!=================================================================================
+
+   !!---------------------------------------------------------------------------------
+   !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary
+   !!   bdy_dyn_fla    : Flather condition for barotropic solution
+   !!---------------------------------------------------------------------------------
+#if defined key_bdy || defined key_bdy_tides
+   !!---------------------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE bdy_oce         ! ocean open boundary conditions
+   USE dynspg_oce      ! for barotropic variables
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE bdytides        ! for tidal harmonic forcing at boundary
+   USE in_out_manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC bdy_dyn_frs  ! routine called in dynspg_flt (free surface case ONLY)
+#if defined key_dynspg_exp || defined key_dynspg_ts
+   PUBLIC bdy_dyn_fla  ! routine called in dynspg_flt (free surface case ONLY)
+#endif
+
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bdy_dyn_frs ( kt )
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE bdy_dyn_frs
+      !!                     ***********************
+      !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the  
+      !!                case of unstructured open boundaries.
+      !!
+      !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in 
+      !!               a three-dimensional baroclinic ocean model with realistic
+      !!               topography. Tellus, 365-382.
+      !!
+      !! History :
+      !!   9.0  !  05-02 (J. Chanut, A. Sellar) Original
+      !!        !  07-07 (D. Storkey) Move Flather implementation to separate routine.
+      !!------------------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt    ! Main time step counter
+
+      !! * Local declarations
+      INTEGER ::   jb, jk, jgrd        ! dummy loop indices
+      INTEGER ::   ii, ij              ! 2D addresses
+      REAL(wp) :: zwgt                 ! boundary weight
+
+      !!------------------------------------------------------------------------------
+      !!  OPA 9.0, LODYC-IPSL (2003)
+      !!------------------------------------------------------------------------------
+
+      ! ---------------------------!
+      ! FRS on the total velocity :!
+      ! ---------------------------!  
+  
+      jgrd=2 !: Relaxation of zonal velocity
+      DO jb = 1, nblen(jgrd)
+        DO jk = 1, jpkm1
+          ii = nbi(jb,jgrd)
+          ij = nbj(jb,jgrd)
+          zwgt = nbw(jb,jgrd)
+
+          ua(ii,ij,jk) = ( ua(ii,ij,jk)*(1.-zwgt) + ubdy(jb,jk)*zwgt ) &
+                                                        * umask(ii,ij,jk)
+        END DO
+      END DO
+
+      jgrd=3 !: Relaxation of meridional velocity
+      DO jb = 1, nblen(jgrd)
+        DO jk = 1, jpkm1
+          ii = nbi(jb,jgrd)
+          ij = nbj(jb,jgrd)
+          zwgt = nbw(jb,jgrd)
+
+          va(ii,ij,jk) = ( va(ii,ij,jk)*(1.-zwgt) + vbdy(jb,jk)*zwgt ) &
+                                                        * vmask(ii,ij,jk)
+        END DO
+      END DO 
+
+      CALL lbc_lnk( ua, 'U', 1. ) ! Boundary points should be updated
+      CALL lbc_lnk( va, 'V', 1. ) !
+
+   END SUBROUTINE bdy_dyn_frs
+
+#if defined key_dynspg_exp || defined key_dynspg_ts
+!! Option to use Flather with dynspg_flt not coded yet...
+   SUBROUTINE bdy_dyn_fla
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE bdy_dyn_fla
+      !!                     ***********************
+      !!              - Apply Flather boundary conditions on normal barotropic velocities 
+      !!                (ln_bdy_fla=.true.)
+      !!
+      !! ** WARNINGS about FLATHER implementation:
+      !!1. According to Palma and Matano, 1998 "after ssh" is used. 
+      !!   In ROMS and POM implementations, it is "now ssh". In the current 
+      !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable. 
+      !!   So I use "before ssh" in the following.
+      !!
+      !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of 
+      !!   fact, the model ssh just inside the dynamical boundary is used (the outside  
+      !!   ssh in the code is not updated).
+      !!
+      !!             - Flather, R. A., 1976: A tidal model of the northwest European
+      !!               continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.     
+      !! History :
+      !!   9.0  !  05-02 (J. Chanut, A. Sellar) Original
+      !!        !  07-07 (D. Storkey) Flather algorithm in separate routine.
+      !!------------------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   jb, jk, jgrd, ji, jj, jim1, jip1, jjm1, jjp1   ! dummy loop indices
+      INTEGER ::   ii, ij                    ! 2D addresses
+      REAL(wp) ::  corr ! Flather correction
+      REAL(wp) :: zwgt                       ! boundary weight
+      REAL(wp) :: elapsed
+
+      !!------------------------------------------------------------------------------
+      !!  OPA 9.0, LODYC-IPSL (2003)
+      !!------------------------------------------------------------------------------
+
+      ! ---------------------------------!
+      ! Flather boundary conditions     :!
+      ! ---------------------------------! 
+     
+      ! Fill temporary array with ssh data (here spgu):
+      jgrd = 1
+      DO jb = 1, nblenrim(jgrd)
+        ji = nbi(jb,jgrd)
+        jj = nbj(jb,jgrd)
+        spgu(ji, jj) = sshbdy(jb)
+#if defined key_bdy_tides
+        spgu(ji, jj) = spgu(ji, jj) + sshtide(jb)
+#endif
+      END DO
+
+      jgrd = 2  !: Flather bc on u-velocity; 
+                ! remember that flagu=-1 if normal velocity direction is outward
+                ! I think we should rather use after ssh ?
+
+      DO jb = 1, nblenrim(jgrd)
+        ji  = nbi(jb,jgrd)
+        jj  = nbj(jb,jgrd) 
+        jim1 = ji+MAX(0, INT(flagu(jb))) ! T pts i-indice inside the boundary
+        jip1 = ji-MIN(0, INT(flagu(jb))) ! T pts i-indice outside the boundary 
+     
+        corr = - flagu(jb) * sqrt (grav / (hu_e(ji, jj) + 1.e-20) )  &
+                                      * ( sshn_e(jim1, jj) - spgu(jip1,jj) )
+        ua_e(ji, jj) = ( ubtbdy(jb) + utide(jb) ) * hu_e(ji,jj)
+        if ( ln_bdy_fla ) then
+          ua_e(ji,jj) = ua_e(ji,jj) + corr * umask(ji,jj,1) * hu_e(ji,jj)
+        endif
+
+      END DO
+       
+      jgrd = 3  !: Flather bc on v-velocity
+                ! remember that flagv=-1 if normal velocity direction is outward
+                 
+      DO jb = 1, nblenrim(jgrd)
+        ji  = nbi(jb,jgrd)
+        jj  = nbj(jb,jgrd) 
+        jjm1 = jj+MAX(0, INT(flagv(jb))) ! T pts j-indice inside the boundary
+        jjp1 = jj-MIN(0, INT(flagv(jb))) ! T pts j-indice outside the boundary 
+     
+        corr = - flagv(jb) * sqrt (grav / (hv_e(ji, jj) + 1.e-20) )  &
+                                      * ( sshn_e(ji, jjm1) - spgu(ji,jjp1) )
+        va_e(ji, jj) = ( vbtbdy(jb) + vtide(jb) ) * hv_e(ji,jj)
+        if ( ln_bdy_fla ) then
+          va_e(ji,jj) = va_e(ji,jj) + corr * vmask(ji,jj,1) * hv_e(ji,jj)
+        endif
+
+      END DO
+
+      CALL lbc_lnk( ua_e, 'U', 1. ) ! Boundary points should be updated
+      CALL lbc_lnk( va_e, 'V', 1. ) !
+
+   END SUBROUTINE bdy_dyn_fla
+#endif
+
+#else
+   !!=================================================================================
+   !!                       ***  MODULE  bdydyn  ***
+   !! Ocean dynamics:   Flow relaxation scheme of velocities on unstruc. open boundary
+   !!=================================================================================
+CONTAINS
+
+   SUBROUTINE bdy_dyn_frs
+                              ! No Unstructured open boundaries ==> empty routine
+   END SUBROUTINE bdy_dyn_frs
+
+   SUBROUTINE bdy_dyn_fla
+                              ! No Unstructured open boundaries ==> empty routine
+   END SUBROUTINE bdy_dyn_fla
+
+#endif
+
+END MODULE bdydyn
Index: /trunk/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 911)
@@ -0,0 +1,533 @@
+
+ MODULE bdyini
+   !!=================================================================================
+   !!                       ***  MODULE  bdyini  ***
+   !! Initialization of unstructured open boundaries
+   !!=================================================================================
+#if defined key_bdy || defined key_bdy_tides
+   !!---------------------------------------------------------------------------------
+   !!   'key_bdy'                                Unstructured Open Boundary Conditions
+   !!---------------------------------------------------------------------------------
+   !!   bdy_init       : Initialization of unstructured open boundaries
+   !!---------------------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE bdy_oce         ! unstructured open boundary conditions
+   USE bdytides        ! tides at open boundaries initialization (tide_init routine)
+   USE in_out_manager  ! I/O units
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp         ! for mpp_sum  
+   USE iom             ! I/O
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC bdy_init        ! routine called by opa.F90
+
+   !! * Substitutions
+
+   !!---------------------------------------------------------------------------------
+   !!   OPA 9.0 , LODYC-IPSL  (2003)
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+   
+   SUBROUTINE bdy_init
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_init  ***
+      !!         
+      !! ** Purpose :   Initialization of the dynamics and tracer fields with 
+      !!      unstructured open boundaries.
+      !!
+      !! ** Method  :  Read initialization arrays (mask, indices) to identify 
+      !!               an unstructured open boundary
+      !!
+      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
+      !!
+      !! History :
+      !!  OPA 9.0 !  05-01  (J. Chanut, A. Sellar)  Original code
+      !!          !  07-01  (D. Storkey)  Update to use IOM module.
+      !!          !  07-01  (D. Storkey)  Tidal forcing.
+      !!----------------------------------------------------------------------      
+      !! * Local declarations
+      INTEGER ::   ji, jj, jk, jgrd,  & ! dummy loop indices
+                   jb, jr, icount,          &
+                   icountr, nb_rim, nb_len, nbr_max
+      INTEGER ::   iw, ie, is, in 
+      INTEGER ::   inum                 ! temporary logical unit
+      INTEGER ::   & ! temporary integers
+         dummy_id
+      INTEGER, DIMENSION (2)       ::   kdimsz
+      INTEGER, DIMENSION(jpbdta, jpbgrd) ::  &  !: Index arrays
+         nbidta, nbjdta, &              !: i and j indices of bdy dta
+         nbrdta		                !: Discrete distance from rim points
+      REAL(wp) :: &
+          efl, wfl, nfl, sfl            ! temporary scalars
+      REAL(wp) , DIMENSION(jpidta,jpjdta) ::   &
+         tmpmsk                         ! global domain mask
+      REAL(wp) , DIMENSION(jpbdta,1)      ::   &
+         ndta                           ! temporary array 
+      CHARACTER(LEN=80),DIMENSION(3)     :: bdyfile
+
+      NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V, &
+                      ln_bdy_clim, ln_bdy_vol, ln_bdy_fla, ln_bdy_mask,         &
+                      nbdy_dta, nb_rimwidth, volbdy
+
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of unstructured open boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
+
+      IF( jperio /= 0 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' E R R O R : Cyclic or symmetric,',   &
+            ' and unstructured open boundary condition are not compatible'
+         IF(lwp) WRITE(numout,*) ' ========== '
+         IF(lwp) WRITE(numout,*)
+         nstop = nstop + 1
+      END IF
+
+#if defined key_obc
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' E R R O R : Straight open boundaries,',   &
+            ' and unstructured open boundaries are not compatible'
+      IF(lwp) WRITE(numout,*) ' ========== '
+      IF(lwp) WRITE(numout,*)
+      nstop = nstop + 1
+#endif
+
+# if defined key_dynspg_rl
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' E R R O R : Rigid lid,',   &
+            ' and unstructured open boundaries are not compatible'
+      IF(lwp) WRITE(numout,*) ' ========== '
+      IF(lwp) WRITE(numout,*)
+      nstop = nstop + 1
+#endif
+
+      ! 0. Read namelist parameters
+      ! ---------------------------
+
+      REWIND( numnam )
+      READ  ( numnam, nambdy )
+
+      ! control prints
+      IF(lwp) WRITE(numout,*) '         nambdy'
+
+      IF ((nbdy_dta/=0).AND.(nbdy_dta/=1)) THEN ! Check nbdy_dta value
+        IF(lwp) WRITE(numout,*)
+        IF(lwp) WRITE(numout,*) ' E R R O R : nbdy_dta =',nbdy_dta,   &
+            'but it should have been 0 or 1'
+        IF(lwp) WRITE(numout,*) ' ========== '
+        IF(lwp) WRITE(numout,*)
+        nstop = nstop + 1
+      ELSE
+        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*) '         data in file (=1) or     nbdy_dta = ', nbdy_dta
+        IF(lwp) WRITE(numout,*) '         initial state used (=0)'
+      END IF 
+
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth
+
+      IF (ln_bdy_vol) THEN
+        IF    (volbdy==1) THEN ! Check volbdy value
+          IF(lwp) WRITE(numout,*) ' '
+          IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy
+          IF(lwp) WRITE(numout,*) '         The total volume will be constant'
+
+        ELSEIF (volbdy==0) THEN
+          IF(lwp) WRITE(numout,*) ' '
+          IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy
+          IF(lwp) WRITE(numout,*) '         The total volume will vary according to &
+                                          &the surface E-P flux'
+        ELSE
+          IF(lwp) WRITE(numout,*) ' '
+          IF(lwp) WRITE(numout,*) ' E R R O R : volbdy =',volbdy,   &
+            'but it should have been 0 or 1'
+          IF(lwp) WRITE(numout,*) ' ========== '
+          IF(lwp) WRITE(numout,*)
+          nstop = nstop + 1
+        END IF
+      ELSE
+        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries'
+        IF(lwp) WRITE(numout,*) ' '
+      ENDIF
+
+      IF (ln_bdy_fla) THEN
+        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*) 'Flather bc with unstructured open boundaries'
+        IF(lwp) WRITE(numout,*) ' '
+      ELSE
+        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*) 'NO Flather bc with unstructured open boundaries'
+        IF(lwp) WRITE(numout,*) ' '
+      ENDIF
+
+      ! 0.5 Read tides namelist 
+      ! ------------------------
+
+      IF ( lk_bdy_tides )  CALL tide_init
+
+      ! 1. Read arrays defining unstructured open boundaries
+      ! ----------------------------------------------------
+
+      ! 1.1 Read global 2D mask at T-points: bdytmask
+      ! *********************************************
+      ! bdytmask=1 on the computational domain AND on open boundaries
+      !         =0 elsewhere   
+ 
+      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
+        tmpmsk(: , : ) = 0.e0
+        tmpmsk(jpizoom+1:jpizoom+jpiglo-2,:  ) = 1.e0          
+      ELSE IF ( ln_bdy_mask ) THEN
+        CALL iom_open( filbdy_mask, inum )
+        CALL iom_get ( inum, jpdom_data, 'bdy_msk', tmpmsk(:,:) )
+        CALL iom_close( inum )
+      ELSE
+        tmpmsk(:,:) = 1.0
+      ENDIF
+
+      ! Save mask over local domain      
+      DO jj = 1, nlcj
+        DO ji = 1, nlci
+          bdytmask(ji,jj) = tmpmsk( mig(ji), mjg(jj))
+        END DO
+      END DO
+
+      ! Derive mask on U and V grid from mask on T grid
+      bdyumask(:,:)=0.e0
+      bdyvmask(:,:)=0.e0
+      DO jj=1, jpjm1
+        DO ji=1, jpim1
+          bdyumask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji+1, jj )
+          bdyvmask(ji,jj)=bdytmask(ji,jj)*bdytmask(ji  ,jj+1)  
+        END DO
+      END DO
+
+      ! Lateral boundary conditions
+      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
+      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
+
+      ! 1.2 Read discrete distance and mapping indices
+      ! **********************************************
+        nbidta(:,:)=0.
+        nbjdta(:,:)=0.
+        nbrdta(:,:)=0.
+
+      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
+
+        icount = 0
+      ! Define west boundary (from ji=2 to ji=1+nb_rimwidth):
+        DO jr=1,nb_rimwidth         
+          DO jj=3,jpjglo-2
+            icount=icount+1
+            nbidta(icount,:) = jr + 1 + (jpizoom-1)
+            nbjdta(icount,:) = jj + (jpjzoom-1) 
+            nbrdta(icount,:) = jr
+          END DO
+        END DO
+
+      ! Define east boundary (from ji=jpiglo-1 to ji=jpiglo-nb_rimwidth):
+        DO jr=1,nb_rimwidth         
+          DO jj=3,jpjglo-2
+            icount=icount+1
+            nbidta(icount,:) = jpiglo-jr + (jpizoom-1)
+            nbidta(icount,2) = jpiglo-jr-1 + (jpizoom-1) ! special case for u points
+            nbjdta(icount,:) = jj + (jpjzoom-1)
+            nbrdta(icount,:) = jr
+          END DO
+        END DO
+            
+      ELSE
+
+      ! Read indices and distances in unstructured boundary data files 
+
+        IF ( lk_bdy ) THEN 
+          bdyfile(1) = filbdy_data_T
+          bdyfile(2) = filbdy_data_U
+          bdyfile(3) = filbdy_data_V
+        ELSE 
+          ! In this case we have tides only at the boundaries
+          ! so read index arrays from tides files for first tidal component
+          bdyfile(1) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc'
+          bdyfile(2) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc'
+          bdyfile(3) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc'
+        ENDIF          
+
+        DO jgrd = 1,3
+           CALL iom_open( bdyfile(jgrd), inum )
+           dummy_id = iom_varid( inum, 'nbidta', kdimsz=kdimsz )  
+           WRITE(numout,*) 'kdimsz : ',kdimsz
+           nb_len = kdimsz(1)
+           IF (nb_len > jpbdta) THEN
+             IF(lwp) WRITE(numout,*)
+             IF(lwp) WRITE(numout,*) ' E R R O R : jpbdta is too small:'
+             IF(lwp) WRITE(numout,*) ' ==========  Boundary array length in file is ', nb_len
+             IF(lwp) WRITE(numout,*) '             But jpbdta is ', jpbdta
+             IF(lwp) WRITE(numout,*) '             File :        ', bdyfile(jgrd)
+             IF(lwp) WRITE(numout,*)
+             nstop = nstop + 1
+           ENDIF
+           CALL iom_get ( inum, jpdom_unknown, 'nbidta', ndta(1:nb_len,:) )
+           DO ji=1,nb_len
+             nbidta(ji,jgrd) = INT( ndta(ji,1) )
+           ENDDO
+           CALL iom_get ( inum, jpdom_unknown, 'nbjdta', ndta(1:nb_len,:) )
+           DO ji=1,nb_len
+             nbjdta(ji,jgrd) = INT( ndta(ji,1) )
+           ENDDO
+           CALL iom_get ( inum, jpdom_unknown, 'nbrdta', ndta(1:nb_len,:) )
+           DO ji=1,nb_len
+             nbrdta(ji,jgrd) = INT( ndta(ji,1) )
+           ENDDO
+           CALL iom_close( inum )
+
+           ! Check that rimwidth in file is big enough:
+           nbr_max = MAXVAL(nbrdta(:,jgrd))
+           IF (nbr_max < nb_rimwidth) THEN
+ 	        IF(lwp) WRITE(numout,*)
+ 	        IF(lwp) WRITE(numout,*) ' E R R O R : Maximum rimwidth in file is ', nbr_max
+ 	        IF(lwp) WRITE(numout,*) ' ==========  but nb_rimwidth is ', nb_rimwidth
+ 	        IF(lwp) WRITE(numout,*)
+ 	        nstop = nstop + 1
+           ELSE
+ 	        IF(lwp) WRITE(numout,*)
+ 	        IF(lwp) WRITE(numout,*) '             Maximum rimwidth in file is ', nbr_max
+ 	        IF(lwp) WRITE(numout,*)
+           END IF           
+
+        ENDDO
+
+      END IF 
+
+      ! 1.3 Dispatch mapping indices and discrete distances on each processor
+      ! *********************************************************************
+     
+      iw = mig(1)+1            ! if monotasking and no zoom, iw=2
+      ie = mig(1) + nlci-1-1   ! if monotasking and no zoom, ie=jpim1
+      is = mjg(1)+1            ! if monotasking and no zoom, is=2
+      in = mjg(1) + nlcj-1-1   ! if monotasking and no zoom, in=jpjm1
+
+      DO jgrd = 1, jpbgrd       
+        icount  = 0
+        icountr = 0
+        DO nb_rim=1, nb_rimwidth
+          DO jb = 1, jpbdta
+          ! check if point is in local domain and equals nb_rim
+            IF ( (nbidta(jb,jgrd) >= iw ).AND. &
+                 (nbidta(jb,jgrd) <= ie ).AND. &
+                 (nbjdta(jb,jgrd) >= is ).AND. &
+                 (nbjdta(jb,jgrd) <= in ).AND. &
+                 (nbrdta(jb,jgrd) == nb_rim ) ) THEN
+
+              icount = icount  + 1
+            
+              IF (nb_rim==1) icountr = icountr+1
+
+              IF (icount > jpbdim) THEN
+                IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'
+                nstop = nstop + 1
+              ELSE
+                nbi(icount, jgrd)  = nbidta(jb,jgrd)- mig(1)+1
+                nbj(icount, jgrd)  = nbjdta(jb,jgrd)- mjg(1)+1
+                nbr(icount, jgrd)  = nbrdta(jb,jgrd)
+                nbmap(icount,jgrd) = jb
+              END IF            
+            END IF
+          END DO
+        END DO
+        nblenrim(jgrd) = icountr !: length of rim boundary data on each proc
+        nblen   (jgrd) = icount  !: length of boundary data on each proc        
+      END DO 
+
+      ! 2. Compute rim weights
+      ! ----------------------
+      
+      DO  jgrd = 1, jpbgrd
+        DO jb = 1, nblen(jgrd)
+          ! tanh formulation
+          nbw(jb,jgrd) = 1.-TANH((FLOAT(nbr(jb,jgrd)-1))/2.)
+          ! quadratic
+!          nbw(jb,jgrd) = (FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth))**2
+          ! linear
+!          nbw(jb,jgrd) =  FLOAT(nb_rimwidth+1-nbr(jb,jgrd))/FLOAT(nb_rimwidth)
+        END DO
+      END DO 
+   
+      ! 3. Mask corrections
+      ! -------------------
+
+      DO jk=1, jpkm1
+        DO jj=1, jpj
+          DO ji=1, jpi
+            tmask(ji,jj,jk)=tmask(ji,jj,jk)*bdytmask(ji,jj)
+            umask(ji,jj,jk)=umask(ji,jj,jk)*bdyumask(ji,jj)
+            vmask(ji,jj,jk)=vmask(ji,jj,jk)*bdyvmask(ji,jj)
+            bmask(ji,jj)=bmask(ji,jj)*bdytmask(ji,jj)
+          END DO      
+        END DO
+      END DO
+
+      ! I am not sure that it is useful:  
+      DO jk=1, jpkm1
+        DO jj=2, jpjm1
+          DO ji=2, jpim1
+            fmask(ji,jj,jk) = fmask(ji,jj,jk) &
+                  &           * bdytmask(ji, jj ) * bdytmask(ji+1, jj )   &
+                  &           * bdytmask(ji,jj+1) * bdytmask(ji+1,jj+1)
+          END DO      
+        END DO
+      END DO
+
+      tmask_i(:,:) = tmask(:,:,1)*tmask_i(:,:)
+
+      bdytmask(:,:)=tmask(:,:,1)
+
+      ! bdy masks and bmask are now set to zero on boundary points:
+
+      jgrd=1 ! In the free surface case, bmask is at T-points
+      DO jb=1, nblenrim(jgrd)     
+        bmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0
+      END DO
+   
+      jgrd=1
+      DO jb=1, nblenrim(jgrd)      
+        bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0
+      END DO
+
+      jgrd=2
+      DO jb=1, nblenrim(jgrd)
+        bdyumask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0
+      END DO
+
+      jgrd=3
+      DO jb=1, nblenrim(jgrd)
+        bdyvmask(nbi(jb,jgrd), nbj(jb,jgrd)) = 0.e0
+      END DO
+
+       ! Lateral boundary conditions
+
+      CALL lbc_lnk( fmask, 'F', 1. )
+      CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
+      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
+      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
+
+      IF ((ln_bdy_vol).OR.(ln_bdy_fla)) THEN 
+
+      ! 4 Indices and directions of rim velocity components
+      ! ---------------------------------------------------
+        !flagu = -1 : u component is normal to the dynamical boundary 
+        !             but its direction is outward
+        !
+        !flagu =  0 : u is tangential
+        !
+        !flagu =  1 : u is normal to the boundary 
+        !             and is direction is inward
+ 
+        icount = 0 
+
+        flagu(:)=0.e0
+ 
+        jgrd=2 ! u-component 
+        DO jb=1, nblenrim(jgrd)  
+          efl=bdytmask(nbi(jb,jgrd)  , nbj(jb,jgrd))
+          wfl=bdytmask(nbi(jb,jgrd)+1, nbj(jb,jgrd))
+          IF ((efl+wfl)==2) THEN
+            icount = icount +1
+          ELSE
+            flagu(jb)=-efl+wfl
+          END IF
+        END DO
+
+        !flagv = -1 : u component is normal to the dynamical boundary 
+        !             but its direction is outward
+        !
+        !flagv =  0 : u is tangential
+        !
+        !flagv =  1 : u is normal to the boundary 
+        !             and is direction is inward
+
+        flagv(:)=0.e0
+
+        jgrd=3 ! v-component
+        DO jb=1, nblenrim(jgrd)  
+          nfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd))
+          sfl = bdytmask(nbi(jb,jgrd), nbj(jb,jgrd)+1)
+          IF ((nfl+sfl)==2) THEN
+            icount = icount +1
+          ELSE
+            flagv(jb)=-nfl+sfl
+          END IF
+        END DO
+ 
+        IF( icount /= 0 ) THEN
+           IF(lwp) WRITE(numout,*)
+           IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
+              ' are not boundary points. Check nbi, nbj, indices.'
+           IF(lwp) WRITE(numout,*) ' ========== '
+           IF(lwp) WRITE(numout,*)
+           nstop = nstop + 1
+        END IF 
+    
+      END IF
+
+      ! 5 Compute total lateral surface for volume correction:
+      ! ------------------------------------------------------
+ 
+      bdysurftot = 0.e0 
+    
+      IF (ln_bdy_vol) THEN  
+        jgrd=2 ! Lateral surface at U-points
+        DO jb=1, nblenrim(jgrd)
+          bdysurftot = bdysurftot + & 
+                    hu(nbi(jb,jgrd), nbj(jb,jgrd)) * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                   * ABS(flagu(jb))*tmask_i(nbi(jb,jgrd)  , nbj(jb,jgrd)) &
+                                   *tmask_i(nbi(jb,jgrd)+1, nbj(jb,jgrd))                   
+        END DO
+
+        jgrd=3 ! Add lateral surface at V-points
+        DO jb=1, nblenrim(jgrd)
+          bdysurftot = bdysurftot + & 
+                    hv(nbi(jb,jgrd), nbj(jb,jgrd)) * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                   * ABS(flagv(jb))*tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                                   *tmask_i(nbi(jb,jgrd), nbj(jb,jgrd)+1)
+        END DO
+
+        IF( lk_mpp ) CALL mpp_sum( bdysurftot ) ! sum over the global domain
+      END IF   
+
+
+      ! 6. Initialise bdy data arrays
+      ! -----------------------------
+
+      tbdy(:,:) = 0.e0
+      sbdy(:,:) = 0.e0
+      ubdy(:,:) = 0.e0
+      vbdy(:,:) = 0.e0
+      sshbdy(:) = 0.e0
+      ubtbdy(:) = 0.e0
+      vbtbdy(:) = 0.e0
+
+      ! 7. Read in tidal constituents and adjust for model start time
+      ! -------------------------------------------------------------
+
+      IF ( lk_bdy_tides )  CALL tide_data
+
+   END SUBROUTINE bdy_init
+
+#else
+   !!---------------------------------------------------------------------------------
+   !!   Dummy module                                   NO unstructured open boundaries
+   !!---------------------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_init      ! Dummy routine
+   END SUBROUTINE bdy_init
+#endif
+
+   !!=================================================================================
+END MODULE bdyini
Index: /trunk/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 911)
@@ -0,0 +1,862 @@
+MODULE bdytides
+   !!=================================================================================
+   !!                       ***  MODULE  bdytides  ***
+   !! Ocean dynamics:   Tidal forcing at open boundaries
+   !!=================================================================================
+#if defined key_bdy_tides
+   !!---------------------------------------------------------------------------------
+   !!   PUBLIC
+   !!   tide_init       : read of namelist
+   !!   tide_data       : read in and initialisation of tidal constituents at boundary
+   !!   tide_update     : calculation of tidal forcing at each timestep
+   !!   PRIVATE
+   !!   uvset           :\
+   !!   vday            : | Routines to correct tidal harmonics forcing for
+   !!   shpen           : | start time of integration
+   !!   ufset           : |
+   !!   vset            :/
+   !!---------------------------------------------------------------------------------
+
+   !!---------------------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers 
+   USE dom_oce         ! ocean space and time domain
+   USE iom
+   USE in_out_manager  ! I/O units
+   USE phycst          ! physical constants
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE bdy_par         ! Unstructured boundary parameters
+   USE bdy_oce         ! ocean open boundary conditions
+   USE daymod          ! calendar
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC tide_init     ! routine called in bdyini
+   PUBLIC tide_data     ! routine called in bdyini
+   PUBLIC tide_update   ! routine called in bdydyn
+
+   !! * Module variables
+   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy_tides = .TRUE. !: tidal forcing at boundaries. 
+
+   CHARACTER(len=80), PUBLIC :: &
+      filtide           !: Filename root for tidal input files
+
+   INTEGER, PUBLIC, PARAMETER ::  ntide_max = 15  ! Max number of tidal contituents
+   INTEGER, PUBLIC            ::  ntide           ! Actual number of tidal constituents
+
+   CHARACTER(len=4), PUBLIC, DIMENSION(ntide_max) :: &
+      tide_cpt         !: Names of tidal components used.
+
+   logical, PUBLIC      :: ln_tide_date ! if true correct tide phases and amplitude for model start date
+
+   REAL(wp), PUBLIC, DIMENSION(ntide_max) :: tide_speed ! Phase speed of tidal constituent (deg/hr)
+   INTEGER,          DIMENSION(ntide_max) :: indx
+   REAL(wp), DIMENSION(jpbdim,ntide_max) ::  &
+      ssh1, ssh2,	& !: Tidal constituents : SSH
+      u1  , u2	,       & !: Tidal constituents : U
+      v1  , v2	          !: Tidal constituents : V
+   
+   
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE tide_init
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE tide_init 
+      !!                     ***********************
+      !! ** Purpose : - Read in namelist for tides
+      !!
+      !! History :
+      !!   NEMO v2.0  !  07-01 (D.Storkey) Original
+      !!------------------------------------------------------------------------------
+      !! * Local declarations
+      INTEGER ::   jtide                  ! dummy loop index
+                                          ! different from nblendta!)
+ 
+      !!------------------------------------------------------------------------------
+      !!  OPA 9.0, LODYC-IPSL (2007)
+      !!------------------------------------------------------------------------------
+
+      NAMELIST/namtide/filtide, tide_cpt, tide_speed, ln_tide_date
+
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
+
+      ! 0. Read namelist parameters
+      ! ---------------------------
+
+      do jtide = 1, ntide_max
+        tide_cpt(jtide) = ''
+      enddo
+
+      REWIND( numnam )
+      READ  ( numnam, namtide )
+
+      ! control prints
+      IF(lwp) WRITE(numout,*) '         namtide'
+
+      ! Count number of components specified:
+      ntide=ntide_max
+      do jtide = 1, ntide_max
+        if ( tide_cpt(jtide) == '' ) then
+           ntide = jtide-1
+           exit
+        endif
+
+      enddo
+      ! find constients in standard list
+      do jtide = 1, ntide
+       indx(jtide) = 0
+       if ( TRIM(tide_cpt(jtide)) == 'Q1' ) indx(jtide) = 1
+       if ( TRIM(tide_cpt(jtide)) == 'O1' ) indx(jtide) = 2
+       if ( TRIM(tide_cpt(jtide)) == 'P1' ) indx(jtide) = 3
+       if ( TRIM(tide_cpt(jtide)) == 'S1' ) indx(jtide) = 4
+       if ( TRIM(tide_cpt(jtide)) == 'K1' ) indx(jtide) = 5
+       if ( TRIM(tide_cpt(jtide)) == '2N2' ) indx(jtide) = 6
+       if ( TRIM(tide_cpt(jtide)) == 'MU2' ) indx(jtide) = 7
+       if ( TRIM(tide_cpt(jtide)) == 'N2' ) indx(jtide) = 8
+       if ( TRIM(tide_cpt(jtide)) == 'NU2' ) indx(jtide) = 9
+       if ( TRIM(tide_cpt(jtide)) == 'M2' ) indx(jtide) = 10
+       if ( TRIM(tide_cpt(jtide)) == 'L2' ) indx(jtide) = 11
+       if ( TRIM(tide_cpt(jtide)) == 'T2' ) indx(jtide) = 12
+       if ( TRIM(tide_cpt(jtide)) == 'S2' ) indx(jtide) = 13
+       if ( TRIM(tide_cpt(jtide)) == 'K2' ) indx(jtide) = 14
+       if ( TRIM(tide_cpt(jtide)) == 'M4' ) indx(jtide) = 15
+      if (indx(jtide) == 0 ) then
+       if(lwp) write(numout,*) 'WARNING: constitunent', jtide,':', tide_cpt(jtide) &
+                                                         , 'not in standard list'
+      endif
+      enddo
+!
+      if ( ntide < 1 ) then 
+         if(lwp) write(numout,*)
+         if(lwp) write(numout,*) ' E R R O R : Did not find any tidal components in namelist.'
+         if(lwp) write(numout,*) ' ========== '
+         if(lwp) write(numout,*)
+         nstop = nstop + 1
+      else
+         if(lwp) write(numout,*)
+         if(lwp) write(numout,*) ntide,' tidal components specified : '
+         if(lwp) write(numout,*) tide_cpt(1:ntide)
+         if(lwp) write(numout,*) ntide,' phase speeds (deg/hr) : '
+         if(lwp) write(numout,*) tide_speed(1:ntide)
+      endif
+
+      ! Initialisation of tidal harmonics arrays
+
+      sshtide(:) = 0.e0
+      utide(:) = 0.e0
+      vtide(:) = 0.e0
+
+   END SUBROUTINE tide_init
+
+   SUBROUTINE tide_data
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE tide_data
+      !!                     ***********************
+      !! ** Purpose : - Read in tidal harmonics data and adjust for the start time of
+      !!                the model run. 
+      !!
+      !! History :
+      !!   NEMO v2.0  !  07-01 (D.Storkey) Original
+      !!              !  08-01 (J.Holt) Add date correction.
+      !!------------------------------------------------------------------------------
+      !! * Local declarations
+      CHARACTER(len=80) :: tidefile       ! full file name for tidal input file 
+      INTEGER ::   jtide, jgrd, jb        ! dummy loop indices
+      INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read)
+      INTEGER, DIMENSION(3) :: lendta     ! length of data in the file (note may be 
+                                          ! different from nblendta!)
+      INTEGER :: iday,imonth,iyear
+      REAL(wp) ::  arg, atde, btde,z1t,z2t           
+      REAL(wp), DIMENSION(jpbdta,1) ::   &
+         pdta                       ! temporary array for data fields
+      REAL(WP), DIMENSION(ntide_max) :: vplu, ftc
+ 
+      !!------------------------------------------------------------------------------
+      !!  OPA 9.0, LODYC-IPSL (2007)
+      !!------------------------------------------------------------------------------
+
+      ! 1. Open files and read in tidal forcing data
+      ! --------------------------------------------
+
+      ipj = 1
+
+      do jtide = 1 ,ntide
+
+         ! SSH fields
+         ! ----------
+
+         tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_T.nc'
+         if(lwp) write(numout,*) 'Reading data from file ',tidefile
+         CALL iom_open( tidefile, inum )
+
+         jgrd = 1
+         IF ( nblendta(jgrd) .le. 0 ) THEN 
+           idvar = iom_varid( inum,'z1' )
+           if(lwp) write(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar)
+           nblendta(jgrd) = iom_file(inum)%dimsz(1,idvar)
+           WRITE(numout,*) 'Dim size for z1 is ',nblendta(jgrd)
+         ENDIF
+         ipi=nblendta(jgrd)
+
+         CALL iom_get ( inum, jpdom_unknown, 'z1', pdta(1:ipi,1:ipj) )
+         DO jb=1, nblenrim(jgrd)
+            ssh1(jb,jtide) =  pdta(nbmap(jb,jgrd),1)
+         END DO
+
+         CALL iom_get ( inum, jpdom_unknown, 'z2', pdta(1:ipi,1:ipj) )
+         DO jb=1, nblenrim(jgrd)
+            ssh2(jb,jtide) =  pdta(nbmap(jb,jgrd),1)
+         END DO
+
+         CALL iom_close( inum )
+
+         ! U fields
+         ! --------
+
+         tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_U.nc'
+         if(lwp) write(numout,*) 'Reading data from file ',tidefile
+         CALL iom_open( tidefile, inum )
+
+         jgrd = 2
+         IF ( lendta(jgrd) .le. 0 ) THEN 
+           idvar = iom_varid( inum,'u1' )
+           lendta(jgrd) = iom_file(inum)%dimsz(1,idvar)
+           WRITE(numout,*) 'Dim size for u1 is ',lendta(jgrd)
+         ENDIF
+         ipi=lendta(jgrd)
+
+         CALL iom_get ( inum, jpdom_unknown, 'u1', pdta(1:ipi,1:ipj) )
+         DO jb=1, nblenrim(jgrd)
+            u1(jb,jtide) =  pdta(nbmap(jb,jgrd),1)
+         END DO
+
+         CALL iom_get ( inum, jpdom_unknown, 'u2', pdta(1:ipi,1:ipj) )
+         DO jb=1, nblenrim(jgrd)
+            u2(jb,jtide) =  pdta(nbmap(jb,jgrd),1)
+         END DO
+
+         CALL iom_close( inum )
+
+         ! V fields
+         ! --------
+
+         tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_V.nc'
+         if(lwp) write(numout,*) 'Reading data from file ',tidefile
+         CALL iom_open( tidefile, inum )
+
+         jgrd = 3
+         IF ( lendta(jgrd) .le. 0 ) THEN 
+           idvar = iom_varid( inum,'v1' )
+           lendta(jgrd) = iom_file(inum)%dimsz(1,idvar)
+           WRITE(numout,*) 'Dim size for v1 is ',lendta(jgrd)
+         ENDIF
+         ipi=lendta(jgrd)
+
+         CALL iom_get ( inum, jpdom_unknown, 'v1', pdta(1:ipi,1:ipj) )
+         DO jb=1, nblenrim(jgrd)
+            v1(jb,jtide) =  pdta(nbmap(jb,jgrd),1)
+         END DO
+
+         CALL iom_get ( inum, jpdom_unknown, 'v2', pdta(1:ipi,1:ipj) )
+         DO jb=1, nblenrim(jgrd)
+            v2(jb,jtide) =  pdta(nbmap(jb,jgrd),1)
+         END DO
+
+         CALL iom_close( inum )
+         enddo ! end loop on tidal components
+!
+!        correct for date factors
+!
+         if( ln_tide_date ) then
+
+
+! Calculate date corrects for 15 standard consituents
+          iyear  = int(ndate0 / 10000  )                           ! initial year
+          imonth = int((ndate0 - iyear * 10000 ) / 100 )          ! initial month
+          iday   = int(ndate0 - iyear * 10000 - imonth * 100)     ! initial day betwen 1 and 30
+
+            call uvset(0,iday,imonth,iyear,ftc,vplu)
+!
+            if(lwp) write(numout,*) 'Correcting tide for date:',iday,imonth,iyear
+
+            do jtide = 1, ntide
+!
+             if(indx(jtide) .ne. 0) then
+              arg=3.14159265d0*vplu(indx(jtide))/180.0d0
+              atde=ftc(indx(jtide))*cos(arg)
+              btde=ftc(indx(jtide))*sin(arg)
+              if(lwp) then
+                write(numout,'(2i5,8f10.6)') jtide,indx(jtide),tide_speed(jtide), &
+                                   ftc(indx(jtide)),vplu(indx(jtide))
+              endif            
+             else
+             atde = 1.0_wp
+             btde = 0.0_wp
+             endif
+!  elevation         
+           jgrd = 1
+           do jb=1, nblenrim(jgrd)                
+            z1t=atde*ssh1(jb,jtide)+btde*ssh2(jb,jtide)
+            z2t=atde*ssh2(jb,jtide)-btde*ssh1(jb,jtide)
+            ssh1(jb,jtide) = z1t
+            ssh2(jb,jtide) = z2t
+           end do
+!  u       
+           jgrd = 2
+           do jb=1, nblenrim(jgrd)                
+            z1t=atde*u1(jb,jtide)+btde*u2(jb,jtide)
+            z2t=atde*u2(jb,jtide)-btde*u1(jb,jtide)
+            u1(jb,jtide) = z1t
+            u2(jb,jtide) = z2t
+           end do
+!  v       
+           jgrd = 3
+           do jb=1, nblenrim(jgrd)                
+            z1t=atde*v1(jb,jtide)+btde*v2(jb,jtide)
+            z2t=atde*v2(jb,jtide)-btde*v1(jb,jtide)
+            v1(jb,jtide) = z1t
+            v2(jb,jtide) = z2t
+           end do
+
+      enddo ! end loop on tidal components
+      endif ! date correction
+
+   END SUBROUTINE tide_data
+
+   SUBROUTINE tide_update ( kt, jit )
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE tide_update
+      !!                     ************************
+      !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays. 
+      !!                
+      !!
+      !! History :
+      !!   NEMO v2.0  !  06-12 (D.Storkey) Original
+      !!------------------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter
+      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option)
+
+      !! * Local declarations
+      INTEGER ::   jtide, jgrd, jb       ! dummy loop indices
+      REAL(wp) ::  arg, sarg                      
+      REAL(wp), DIMENSION(ntide_max) :: sist, cost
+ 
+      !!------------------------------------------------------------------------------
+      !!  OPA 9.0, LODYC-IPSL (2003)
+      !!------------------------------------------------------------------------------
+
+      ! Note tide phase speeds are in deg/hour, so we need to convert the
+      ! elapsed time in seconds to hours by dividing by 3600.0
+      if ( jit .eq. 0 ) then  
+        arg = kt*rdt*rad/3600.0
+      else ! we are in a barotropic subcycle (for timesplitting option)
+        arg = ( (kt-1)*rdt + jit*rdtbt ) * rad/3600.0
+      endif
+
+      do jtide = 1,ntide
+        sarg = arg*tide_speed(jtide)
+        cost(jtide) = cos(sarg)
+        sist(jtide) = sin(sarg)
+      enddo
+
+      !! summing of tidal constituents into BDY arrays
+
+      sshtide(:) = 0.0
+      utide(:) = 0.0
+      vtide(:) = 0.0
+
+      do jtide = 1 ,ntide
+        jgrd=1 !: SSH on tracer grid.
+        do jb = 1, nblenrim(jgrd)
+          sshtide(jb) =sshtide(jb)+ ssh1(jb,jtide)*cost(jtide) + ssh2(jb,jtide)*sist(jtide)
+      !    if(lwp) write(numout,*) 'z',jb,jtide,sshtide(jb), ssh1(jb,jtide),ssh2(jb,jtide)
+        enddo
+        jgrd=2 !: U grid
+        do jb=1, nblenrim(jgrd)
+          utide(jb) = utide(jb)+ u1(jb,jtide)*cost(jtide) + u2(jb,jtide)*sist(jtide)
+      !    if(lwp) write(numout,*) 'u',jb,jtide,utide(jb), u1(jb,jtide),u2(jb,jtide)
+
+        enddo
+        jgrd=3 !: V grid
+        do jb=1, nblenrim(jgrd)
+          vtide(jb) = vtide(jb)+ v1(jb,jtide)*cost(jtide) + v2(jb,jtide)*sist(jtide)
+      !    if(lwp) write(numout,*) 'v',jb,jtide,vtide(jb), v1(jb,jtide),v2(jb,jtide)
+
+        enddo
+      enddo
+
+   END SUBROUTINE tide_update
+
+!
+!
+   SUBROUTINE uvset (ihs,iday,imnth,iyr,f,vplu)
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE uvset
+      !!                     ************************
+      !! ** Purpose : - adjust tidal forcing for date factors
+      !!                
+      !!
+      !! History :
+      !!
+      !! Origins POLCOMS v6.3 2007
+      !!   NEMO v2.3  !  Jason Holt
+      !!------------------------------------------------------------------------------
+      implicit none
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
+      INTEGER, INTENT( in ) ::   iday    ! start time days
+      INTEGER, INTENT( in ) ::   imnth	 ! start time month
+      INTEGER, INTENT( in ) ::   iyr	 ! start time year
+      INTEGER, PARAMETER ::  nc =15    ! maximum number of constituents
+     
+      REAL(WP)              ::   f(nc)    ! nodal correction
+      REAL(WP)              ::   vplu(nc) ! phase correction
+!    
+      INTEGER		    :: year,vd,ivdy,ndc,i,k
+      REAL(WP)   	    :: u(nc),v(nc),zig(nc),rtd
+      REAL(WP)   	    :: ss,h,p,en,p1
+      CHARACTER(len=8), DIMENSION(nc) :: cname
+!
+      !!------------------------------------------------------------------------------
+      !!  OPA 9.0, LODYC-IPSL (2007)
+      !!------------------------------------------------------------------------------
+
+      data cname/ 'q1', 'o1', 'p1', 's1', 'k1', &
+                '2n2','mu2', 'n2','nu2', 'm2', 'l2', 't2', 's2', 'k2', &
+                'm4'/
+      data zig/.2338507481,.2433518789,.2610826055,.2617993878 &
+     ,        .2625161701                                     &
+     ,        .4868657873,.4881373225,.4963669182,.4976384533 &
+     ,        .5058680490,.5153691799,.5228820265,.5235987756 &
+     ,        .5250323419                                     &
+     ,       1.011736098/
+!
+!     ihs             -  start time gmt on ...
+!     iday/imnth/iyr  -  date   e.g.   12/10/87
+!
+      call vday(iday,imnth,iyr,ivdy)
+      vd=ivdy
+!
+!rp   note change of year number for d. blackman shpen
+!rp   if(iyr.ge.1000) year=iyr-1900
+!rp   if(iyr.lt.1000) year=iyr
+      year = iyr
+!
+!.....year  =  year of required data
+!.....vd    =  day of required data..set up for 0000gmt day year
+!
+      ndc = nc
+!.....ndc   =  number of constituents allowed
+!
+      rtd = 360.0/6.2831852
+      do i = 1,ndc
+      zig(i) = zig(i)*rtd
+!      sigo(i)= zig(i)
+      enddo
+!
+      if(year == 0) then
+      return
+      endif
+      
+      call shpen (year,vd,ss,h,p,en,p1)
+      call ufset (p,en,u,f)
+      call vset (ss,h,p,en,p1,v)
+!
+      do  k=1,nc
+      
+      vplu(k)=v(k)+u(k)
+      vplu(k)=vplu(k)+dble(ihs)*zig(k)
+!
+      do while ( vplu(k) < 0 )
+        vplu(k) = vplu(k) + 360.0
+      enddo
+!
+      do while (vplu(k) > 360. )
+        vplu(k) = vplu(k) - 360.0
+      enddo
+!
+      enddo
+!
+      END SUBROUTINE uvset
+
+
+      SUBROUTINE vday(iday,imnth,iy,ivdy)
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE vday
+      !!                     ****************
+      !! ** Purpose : - adjust tidal forcing for date factors
+      !!                
+      !!
+      !! History :
+      !!
+      !! Origins POLCOMS v6.3 2007
+      !!   NEMO v2.3  !  Jason Holt
+      !!------------------------------------------------------------------------------
+      implicit none
+!
+!     subroutine arguments
+      integer :: iday,imnth,iy,ivdy
+!
+! local variables
+      integer iyr
+      !!------------------------------------------------------------------------------
+      !!  NEMO 2.3, LODYC-IPSL (2008)
+      !!------------------------------------------------------------------------------
+
+! =================================================================
+!
+!     calculate day number in year from day/month/year
+!
+! =================================================================
+       if(imnth.eq.1) ivdy=iday
+       if(imnth.eq.2) ivdy=iday+31
+       if(imnth.eq.3) ivdy=iday+59
+       if(imnth.eq.4) ivdy=iday+90
+       if(imnth.eq.5) ivdy=iday+120
+       if(imnth.eq.6) ivdy=iday+151
+       if(imnth.eq.7) ivdy=iday+181
+       if(imnth.eq.8) ivdy=iday+212
+       if(imnth.eq.9) ivdy=iday+243
+       if(imnth.eq.10) ivdy=iday+273
+       if(imnth.eq.11) ivdy=iday+304
+       if(imnth.eq.12) ivdy=iday+334
+      iyr=iy
+       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
+       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
+       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
+
+      RETURN
+      END SUBROUTINE vday
+
+      SUBROUTINE shpen (year,vd,s,h,p,en,p1)
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE shpen
+      !!                     *****************
+      !! ** Purpose : - calculate astronomical arguments for tides
+      !!                this version from d. blackman 30 nove 1990
+      !!                
+      !! History :
+      !!
+      !! Origins POLCOMS v6.3 2007
+      !!   NEMO v2.3  !  Jason Holt
+      !!------------------------------------------------------------------------------
+      implicit none
+
+!     subroutine arguments
+      integer  ::year,vd
+      real(wp) :: s,h,p,en,p1
+
+!     local variables
+
+      integer yr,ilc,icent,it,iday,ild,ipos,nn,iyd
+      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
+
+      !!------------------------------------------------------------------------------
+      !!  NEMO 2.3, LODYC-IPSL (2008)
+      !!------------------------------------------------------------------------------
+
+      data delt /-5.04,-3.90,-2.87,-0.58,0.71,1.80, &
+      3.08, 4.63, 5.86, 7.21, 8.58,10.50,12.10,    &
+     12.49,14.41,15.59,15.81,17.52,19.01,18.39,    &
+     19.55,20.36,21.01,21.81,21.76,22.35,22.68,    &
+     22.94,22.93,22.69,22.94,23.20,23.31,23.63,    &
+     23.47,23.68,23.62,23.53,23.59,23.99,23.80,    &
+     24.20,24.99,24.97,25.72,26.21,26.37,26.89,    &
+     27.68,28.13,28.94,29.42,29.66,30.29,30.96,    &
+     31.09,31.59,31.52,31.92,32.45,32.91,33.39,    &
+     33.80,34.23,34.73,35.40,36.14,36.99,37.87,    &
+     38.75,39.70,40.70,41.68,42.82,43.96,45.00,    &
+     45.98,47.00,48.03,49.10,50.10,50.97,51.81,    &
+     52.57/
+!
+      cycle = 360.0
+      ilc = 0
+      icent = year/100
+      yr = year - icent*100
+      t = icent - 20
+!
+!     for the following equations
+!     time origin is fixed at 00 hr of jan 1st,2000.
+!     see notes by cartwright
+!
+      it = icent - 20
+      if (it) 1,2,2
+    1 iday = it/4 -it
+      go to 3
+    2 iday = (it+3)/4 - it
+!
+!     t is in julian century
+!     correction in gegorian calander where only century year divisible
+!     by 4 is leap year.
+!
+    3 continue
+!
+      td = 0.0
+!
+      if (yr) 4,5,4
+!
+    4 iyd = 365*yr
+      ild = (yr-1)/4
+      if((icent - (icent/4)*4) .eq. 0) ilc = 1
+      td = iyd + ild + ilc
+!
+    5 td = td + iday + vd -1.0 - 0.5
+      t = t + (td/36525.0)
+!
+      ipos=year-1899
+      if (ipos .lt. 0) go to 7
+      if (ipos .gt. 83) go to 6
+!
+      delta = (delt(ipos+1)+delt(ipos))/2.0
+      go to 7
+!
+    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
+!
+    7 deltat = delta * 1.0e-6
+!
+      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
+      h   = 280.4661 + 36000.7698*t + 0.0003*t*t + 11.0*deltat
+      p   =  83.3535 + 4069.0139*t - 0.0103*t*t + deltat
+      en  = 234.9555 + 1934.1363*t - 0.0021*t*t + deltat
+      p1  = 282.9384 + 1.7195*t + 0.0005*t*t
+!
+      nn = s/cycle
+      s = s - nn*cycle
+      if ( s .lt. 0.0) s = s+cycle
+!
+      nn = h/cycle
+      h = h-cycle*nn
+      if (h .lt. 0.0) h = h+cycle
+!
+      nn = p/cycle
+      p = p- cycle*nn
+      if (p .lt. 0.0) p = p+cycle
+!
+      nn = en/cycle
+      en = en-cycle*nn
+      if(en .lt. 0.0) en = en + cycle
+      en = cycle - en
+!
+      nn = p1/cycle
+      p1 = p1 - nn*cycle
+!
+      RETURN
+!
+      END SUBROUTINE shpen
+
+
+      SUBROUTINE ufset (p,cn,b,a)
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE ufset
+      !!                     *****************
+      !! ** Purpose : - calculate nodal parameters for the tides
+      !!                
+      !!
+      !! History :
+      !!
+      !! Origins POLCOMS v6.3 2007
+      !!   NEMO v2.3  !  Jason Holt
+      !!------------------------------------------------------------------------------
+
+      implicit none
+      integer nc
+      parameter (nc=15)
+!     subroutine arguments
+      real(wp) p,cn
+!
+!
+!     local variables
+      real(wp) ::   w1,w2,w3,w4,w5,w6,w7,w8,nw,pw,rad
+      real(wp) ::   a(nc),b(nc)
+      integer ndc,k
+!
+      !!------------------------------------------------------------------------------
+      !!  NEMO 2.3, LODYC-IPSL (2008)
+      !!------------------------------------------------------------------------------
+
+      ndc=nc
+!
+!    a=f       ,  b =u
+!    t is  zero as compared to tifa.
+      rad = 6.2831852d0/360.0
+      pw = p*rad
+      nw = cn*rad
+      w1 = cos(nw)
+      w2 = cos(2*nw)
+      w3 = cos(3*nw)
+      w4 = sin(nw)
+      w5 = sin(2*nw)
+      w6 = sin(3*nw)
+      w7 = 1 -0.2505*cos(2*pw) -0.1102*cos(2*pw-nw) &
+            -0.156*cos(2*pw-2*nw) -0.037*cos(nw)
+      w8 = -0.2505*sin(2*pw) -0.1102*sin(2*pw-nw) &
+          -0.0156*sin(2*pw-2*nw) -0.037*sin(nw)
+!
+      a(1) = 1.0089+0.1871*w1-0.0147*w2+0.0014*w3
+      b(1) = 0.1885*w4 - 0.0234*w5+.0033*w6
+!   q1
+      a(2) = a(1)
+      b(2) = b(1)
+!   o1
+      a(3) = 1.0
+      b(3) = 0.0
+!   p1
+      a(4) = 1.0
+      b(4) = 0.0
+!   s1
+      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
+      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
+!   k1
+      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
+      b(6) = -0.0374*w4
+!  2n2
+      a(7) = a(6)
+      b(7) = b(6)
+!  mu2
+      a(8) = a(6)
+      b(8) = b(6)
+!   n2
+      a(9) = a(6)
+      b(9) = b(6)
+!  nu2
+      a(10) = a(6)
+      b(10) = b(6)
+!   m2
+      a(11) = sqrt(w7*w7+w8*w8)
+      b(11) = atan(w8/w7)
+      if(w7.lt.0) b(11) = b(11) + 3.141992
+!   l2
+      a(12) = 1.0
+      b(12) = 0.0
+!   t2
+      a(13)= a(12)
+      b(13)= b(12)
+!   s2
+      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
+      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
+!   k2
+      a(15) = a(6)*a(6)
+      b(15) = 2*b(6)
+!   m4
+!
+      do 40 k = 1,ndc
+      b(k) = b(k)/rad
+32    if (b(k)) 34,35,35
+34    b(k) = b(k) + 360.0
+      go to 32
+35    if (b(k)-360.0) 40,37,37
+37    b(k) = b(k)-360.0
+      go to 35
+40    continue
+      RETURN
+      END SUBROUTINE ufset
+
+      SUBROUTINE vset( s,h,p,en,p1,v)
+      !!------------------------------------------------------------------------------
+      !!                      SUBROUTINE vset
+      !!                     ****************
+      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
+      !!                
+      !!
+      !! History :
+      !!
+      !! Origins POLCOMS v6.3 2007
+      !!   NEMO v2.3  !  Jason Holt
+      !!------------------------------------------------------------------------------
+      implicit none
+      integer nc
+      parameter (nc=15)
+!     subroutine arguments
+      real(wp) s,h,p,en,p1
+      real(wp)   v(nc)
+!
+      integer ndc,k
+      !!------------------------------------------------------------------------------
+      !!  NEMO 2.3, LODYC-IPSL (2008)
+      !!------------------------------------------------------------------------------
+
+      ndc = nc
+!   v s  are computed here.
+      v(1) =-3*s +h +p +270
+!   q1
+      v(2) =-2*s +h +270.0
+!   o1
+      v(3) =-h +270
+!   p1
+      v(4) =180
+!   s1
+      v(5) =h +90.0
+!   k1
+      v(6) =-4*s +2*h +2*p
+!   2n2
+      v(7) =-4*(s-h)
+!   mu2
+      v(8) =-3*s +2*h +p
+!   n2
+      v(9) =-3*s +4*h -p
+!   mu2
+      v(10) =-2*s +2*h
+!   m2
+      v(11) =-s +2*h -p +180
+!   l2
+      v(12) =-h +p1
+!   t2
+      v(13) =0
+!   s2
+      v(14) =h+h
+!   k2
+      v(15) =2*v(10)
+!   m4
+!
+      do 72 k = 1, ndc
+69    if (v(k) ) 70,71,71
+70    v(k)=v(k)+360.0
+      go to 69
+71    if ( v(k) - 360.0) 72,73,73
+73    v(k)=v(k)-360.0
+      go to 71
+72    continue
+      RETURN
+      END SUBROUTINE vset
+
+#else
+   !!=================================================================================
+   !!                       ***  MODULE  bdytides ***
+   !!=================================================================================
+
+   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy_tides = .FALSE. !: tidal forcing at boundaries. 
+
+   CHARACTER(len=80), PUBLIC :: &
+      filtide           !: Filename root for tidal input files
+
+   CHARACTER(len=4), PUBLIC, DIMENSION(1) :: &
+      tide_cpt         !: Names of tidal components used.
+
+CONTAINS
+
+   SUBROUTINE tide_init
+                              ! No tidal forcing at boundaries ==> empty routine
+   END SUBROUTINE tide_init
+
+   SUBROUTINE tide_data
+                              ! No tidal forcing at boundaries ==> empty routine
+   END SUBROUTINE tide_data
+
+   SUBROUTINE tide_update( kt, jit )
+      INTEGER :: kt, jit
+      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit
+
+                              ! No tidal forcing at boundaries ==> empty routine
+   END SUBROUTINE tide_update
+
+#endif
+
+END MODULE bdytides
Index: /trunk/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 911)
@@ -0,0 +1,87 @@
+MODULE bdytra
+   !!=================================================================================
+   !!                       ***  MODULE  bdytra  ***
+   !! Ocean tracers:   Flow Relaxation Scheme of tracers on each open boundary
+   !!=================================================================================
+#if defined key_bdy
+   !!---------------------------------------------------------------------------------
+   !!   'key_bdy'      :                         Unstructured Open Boundary Conditions
+   !!---------------------------------------------------------------------------------
+   !!   bdy_tra        : Relaxation of tracers on unstructured open boundaries
+   !!---------------------------------------------------------------------------------
+   !! * Modules used
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables 
+   USE bdy_oce         ! ocean open boundary conditions
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC bdy_tra     ! routine called in tranxt.F90 
+
+   !! * Substitutions
+
+   !!---------------------------------------------------------------------------------
+   !!   OPA 9.0 , LODYC-IPSL  (2003)
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bdy_tra( kt )
+      !!------------------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_tra  ***
+      !!                    
+      !! ** Purpose : Apply the Flow Relaxation Scheme for tracers in the  
+      !!              case of unstructured open boundaries.
+      !! 
+      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in 
+      !!             a three-dimensional baroclinic ocean model with realistic
+      !!             topography. Tellus, 365-382.
+      !!  History :
+      !!   9.0  !  05-01 (J. Chanut, A. Sellar) Original
+      !!------------------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt
+
+      !! * Local declarations
+      REAL(wp) :: zwgt                       ! boundary weight
+      INTEGER ::   jb, jk, jgrd              ! dummy loop indices
+      INTEGER ::   ii, ij                    ! 2D addresses
+      !!------------------------------------------------------------------------------
+
+      jgrd=1 !: Everything is at T-points here
+ 
+      DO jb = 1, nblen(jgrd)
+        DO jk = 1, jpkm1
+          ii = nbi(jb,jgrd)
+          ij = nbj(jb,jgrd)
+          zwgt = nbw(jb,jgrd)
+
+          ! Temperature relaxation at the boundary   
+          ta(ii,ij,jk) = ( ta(ii,ij,jk)*(1.-zwgt)  +  tbdy(jb,jk)*zwgt ) &
+                                                          * tmask(ii,ij,jk)         
+
+          ! Salinity relaxation at the boundary   
+          sa(ii,ij,jk) = ( sa(ii,ij,jk)*(1.-zwgt)  +  sbdy(jb,jk)*zwgt ) &
+                                                          * tmask(ii,ij,jk)
+   
+        END DO
+      END DO 
+
+      CALL lbc_lnk( ta, 'T', 1. ) ! Boundary points should be updated
+      CALL lbc_lnk( sa, 'T', 1. ) !
+
+   END SUBROUTINE bdy_tra
+#else
+   !!---------------------------------------------------------------------------------
+   !!   Default option                                                    Empty module
+   !!---------------------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_tra      ! Empty routine
+   END SUBROUTINE bdy_tra
+#endif
+
+   !!=================================================================================
+END MODULE bdytra
Index: /trunk/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 911)
+++ /trunk/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 911)
@@ -0,0 +1,206 @@
+MODULE bdyvol
+   !!=================================================================================
+   !!                       ***  MODULE  bdyvol  ***
+   !! Ocean dynamic :  Volume constraint when unstructured boundary 
+   !!                  and Free surface are used
+   !!=================================================================================
+#if   defined key_bdy   &&   defined key_dynspg_flt
+   !!---------------------------------------------------------------------------------
+   !!   'key_bdy'               and              unstructured open boundary conditions
+   !!   'key_dynspg_flt'                                  constant volume free surface
+   !!---------------------------------------------------------------------------------
+   !! * Modules used
+   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 lib_mpp         ! for mppsum
+   USE in_out_manager  ! I/O manager
+   USE sbc_oce         ! ocean surface boundary conditions
+
+   IMPLICIT NONE
+   PRIVATE
+
+   !! * Accessibility
+   PUBLIC bdy_vol        ! routine called by dynspg_flt.h90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!---------------------------------------------------------------------------------
+   !!   OPA 9.0 , LODYC-IPSL  (2003)
+   !!---------------------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE bdy_vol ( kt )
+      !!------------------------------------------------------------------------------
+      !!                      ***  ROUTINE bdyvol  ***
+      !!
+      !! ** Purpose : 
+      !!      This routine is called in dynspg_flt to control 
+      !!      the volume of the system. A correction velocity is calculated
+      !!      to correct the total transport through the unstructured OBC. 
+      !!      The total depth used is constant (H0) to be consistent with the 
+      !!      linear free surface coded in OPA 8.2
+      !!
+      !! ** Method :  
+      !!      The correction velocity (zubtpecor here) is defined calculating
+      !!      the total transport through all open boundaries (trans_bdy) minus
+      !!      the cumulate E-P flux (zCflxemp) divided by the total lateral 
+      !!      surface (bdysurftot) of the unstructured boundary. 
+      !!
+      !!      zubtpecor = [trans_bdy - zCflxemp ]*(1./bdysurftot)
+      !!
+      !!      with zCflxemp => sum of (Evaporation minus Precipitation)
+      !!                       over all the domain in m3/s at each time step.
+      !!
+      !!      zCflxemp < 0 when precipitation dominate
+      !!      zCflxemp > 0 when evaporation dominate
+      !!
+      !!      There are 2 options (user's desiderata): 
+      !!
+      !!         1/ The volume changes according to E-P, this is the default
+      !!            option. In this case the cumulate E-P flux are setting to
+      !!            zero (zCflxemp=0) to calculate the correction velocity. So
+      !!            it will only balance the flux through open boundaries.
+      !!            (set volbdy to 0 in tne namelist for this option)
+      !!
+      !!         2/ The volume is constant even with E-P flux. In this case
+      !!            the correction velocity must balance both the flux 
+      !!            through open boundaries and the ones through the free
+      !!            surface. 
+      !!            (set volbdy to 1 in tne namelist for this option)
+      !!
+      !! History :
+      !!   8.5  !  05-01 (J. Chanut, A. Sellar) Original code
+      !!        !  06-01 (J. Chanut) Bug correction
+      !!----------------------------------------------------------------------------
+      !! * Arguments
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+
+      !! * Local declarations
+      INTEGER ::   ji,jj,jb, jgrd, jk
+      REAL(wp) ::   zubtpecor
+      REAL(wp) ::   zCflxemp
+      REAL(wp) ::   ztranst
+      !!-----------------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN 
+         IF(lwp) WRITE(numout,*)'        '
+         IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC'
+         IF(lwp) WRITE(numout,*)'~~~~~~~'
+         IF(lwp) WRITE(numout,*)'        '
+      END IF 
+
+      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain.
+      ! ---------------------------------------------------------------------------
+ 
+      zCflxemp = 0.e0
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zCflxemp = zCflxemp + ( (emp(ji,jj)*bdytmask(ji,jj)*tmask_i(ji,jj) )/rauw) &
+                                  *e1v(ji,jj)*e2u(ji,jj)
+         END DO
+      END DO
+      IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain
+
+      ! 2. Barotropic velocity through the unstructured open boundary
+      ! -------------------------------------------------------------
+
+      zubtpecor = 0.e0
+
+      jgrd = 2 ! cumulate u component contribution first 
+      DO jb = 1, nblenrim(jgrd)
+        DO jk = 1, jpkm1
+          zubtpecor = zubtpecor + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
+                                            * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                                            * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk)
+        END DO
+      END DO
+
+      jgrd = 3 ! then add v component contribution
+       DO jb = 1, nblenrim(jgrd)
+        DO jk = 1, jpkm1
+          zubtpecor = zubtpecor + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
+                                            * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                                            * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk) 
+        END DO
+      END DO
+
+      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
+
+
+      ! 3. The normal velocity correction
+      ! ---------------------------------
+
+      IF (volbdy==1) THEN
+        zubtpecor = (zubtpecor - zCflxemp)*(1./bdysurftot) 
+      ELSE
+        zubtpecor =  zubtpecor*(1./bdysurftot)
+      END IF
+
+      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
+         IF(lwp) WRITE(numout,*)'        '
+         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt
+         IF(lwp) WRITE(numout,*)'~~~~~~~ '
+         IF(lwp) WRITE(numout,*)'        '
+         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)'
+         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',bdysurftot,'(m2)'
+         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)'
+         IF(lwp) WRITE(numout,*)'        '
+      END IF 
+
+      ! 4. Correction of the total velocity on the unstructured 
+      !    boundary to respect the mass flux conservation
+      ! -------------------------------------------------------
+
+      ztranst = 0.e0
+
+      jgrd = 2 ! correct u component
+      DO jb = 1, nblenrim(jgrd)
+        DO jk = 1, jpkm1
+          ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) = ua(nbi(jb, jgrd), nbj(jb, jgrd), jk) &
+                            -flagu(jb) * zubtpecor * umask(nbi(jb, jgrd), nbj(jb, jgrd), jk)
+          ztranst = ztranst + flagu(jb) * ua (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
+                                        * e2u(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                                        * fse3u(nbi(jb,jgrd), nbj(jb,jgrd), jk)
+        END DO
+      END DO
+
+      jgrd = 3 ! correct v component
+       DO jb = 1, nblenrim(jgrd)
+        DO jk = 1, jpkm1
+          va(nbi(jb, jgrd), nbj(jb, jgrd), jk) = va(nbi(jb, jgrd), nbj(jb, jgrd), jk) &
+                            -flagv(jb) * zubtpecor * vmask(nbi(jb, jgrd), nbj(jb, jgrd), jk)
+          ztranst = ztranst + flagv(jb) * va (nbi(jb,jgrd), nbj(jb,jgrd), jk) &
+                                        * e1v(nbi(jb,jgrd), nbj(jb,jgrd)) &
+                                        * fse3v(nbi(jb,jgrd), nbj(jb,jgrd), jk)
+        END DO
+      END DO
+      
+      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
+ 
+      ! 5. Check the cumulate transport through unstructured OBC
+      !    once barotropic velocities corrected
+      ! --------------------------------------------------------
+
+      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
+         IF(lwp) WRITE(numout,*)'        '
+         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)'
+         IF(lwp) WRITE(numout,*)'        '
+      END IF 
+
+   END SUBROUTINE bdy_vol
+
+#else
+   !!---------------------------------------------------------------------------------
+   !!  Default option :                                                   Empty module
+   !!---------------------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_vol        ! Empty routine
+   END SUBROUTINE bdy_vol
+#endif
+
+   !!=================================================================================
+END MODULE bdyvol
Index: /trunk/NEMO/OPA_SRC/DOM/domvvl.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DOM/domvvl.F90	(revision 911)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History :   9.0  !  06-06  (B. Levier, L. Marie)  original code
+   !!              "   !  07-07  (D. Storkey) Bug fixes and code for BDY option.
    !!----------------------------------------------------------------------
 
@@ -24,4 +25,5 @@
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE bdy_oce         ! unstructured open boundary conditions
 
    IMPLICIT NONE
@@ -116,8 +118,15 @@
             zpk = mbathy(ji,jj) - 1
             DO jk = 1, zpk
+#if defined key_sigma_vvl
+               mut(ji,jj,jk) = 1./SUM( fsve3t(ji,jj,1:zpk) ) 
+               muu(ji,jj,jk) = 1./SUM( fsve3u(ji,jj,1:zpk) ) 
+               muv(ji,jj,jk) = 1./SUM( fsve3v(ji,jj,1:zpk) ) 
+               muf(ji,jj,jk) = 1./SUM( fsve3f(ji,jj,1:zpk) ) 
+#else
                mut(ji,jj,jk) = SUM( fsve3t(ji,jj,jk:zpk) ) / zmut(ji,jj)
                muu(ji,jj,jk) = SUM( fsve3u(ji,jj,jk:zpk) ) / zmuu(ji,jj)
                muv(ji,jj,jk) = SUM( fsve3v(ji,jj,jk:zpk) ) / zmuv(ji,jj)
                muf(ji,jj,jk) = SUM( fsve3f(ji,jj,jk:zpk) ) / zmuf(ji,jj)
+#endif
             END DO
          END DO
@@ -191,4 +200,9 @@
       END DO
 
+      IF( MINVAL(fse3t(:,:,:)) < 0.0 ) THEN
+         CALL ctl_stop('E R R O R : dom_vvl','Level thickness fse3t less than zero.')
+         nstop = nstop+1
+      ENDIF
+
       ! Local depth or Inverse of the local depth of the water column at u- and v-points
       ! ------------------------------
@@ -288,4 +302,12 @@
 #endif
 
+#if defined key_bdy || key_bdy_tides
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
+         END DO
+      END DO
+#endif
+
       CALL lbc_lnk( zhdiv, 'T', 1. )
 
Index: /trunk/NEMO/OPA_SRC/DYN/divcur.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/divcur.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DYN/divcur.F90	(revision 911)
@@ -14,4 +14,5 @@
    USE in_out_manager  ! I/O manager
    USE obc_oce         ! ocean lateral open boundary condition
+   USE bdy_oce	       ! Unstructured open boundaries variables
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
@@ -77,4 +78,5 @@
       !!   8.2  !  00-03  (G. Madec)  no slip accurate
       !!   9.0  !  03-08  (G. Madec)  merged of cur and div, free form, F90
+      !!        !  05-01  (J. Chanut, A. Sellar) unstructured open boundaries
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -134,4 +136,12 @@
          ENDIF
 #endif
+#endif         
+#if defined key_bdy || key_bdy_tides
+         ! unstructured open boundaries (div must be zero behind the open boundary)
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             hdivn(ji,jj,jk)=hdivn(ji,jj,jk)*bdytmask(ji,jj)
+           END DO
+         END DO
 #endif         
 #if defined key_agrif
@@ -291,4 +301,5 @@
       !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
       !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
+      !!        !  05-01  (J. Chanut) Unstructured open boundaries
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -344,4 +355,12 @@
 #endif
 #endif         
+#if defined key_bdy || key_bdy_tides
+         ! unstructured open boundaries (div must be zero behind the open boundary)
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             hdivn(ji,jj,jk)=hdivn(ji,jj,jk)*bdytmask(ji,jj)
+           END DO
+         END DO
+#endif        
 #if defined key_agrif
          if ( .NOT. AGRIF_Root() ) then
Index: /trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DYN/dynnxt.F90	(revision 911)
@@ -16,4 +16,7 @@
    USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
+   USE bdy_oce         ! unstructured open boundary conditions
+   USE bdydta          ! unstructured open boundary conditions
+   USE bdydyn          ! unstructured open boundary conditions
    USE dynspg_oce      ! type of surface pressure gradient
    USE lbclnk          ! lateral boundary condition (or mpp link)
@@ -66,4 +69,5 @@
       !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond.
       !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+      !!    "   !  07-07  (D. Storkey) Calls to BDY routines. 
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -186,4 +190,23 @@
             END DO
          END DO
+
+         IF( lk_vvl ) THEN
+            ! Unweight velocities prior to updating open boundaries.
+
+            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
+               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
+                  ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk)
+                  va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk)
+
+                  un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk)
+                  vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk)
+
+                  ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk)
+                  vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk)
+               END DO
+            END DO
+
+         ENDIF
+
 # if defined key_obc
          !                                             ! ===============
@@ -209,4 +232,48 @@
 
          IF ( ln_vol_cst ) CALL obc_vol( kt )
+
+      ENDIF
+
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+# elif defined key_bdy || key_bdy_tides
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+      ! Update (ua,va) along open boundaries (for exp or ts options).
+      IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN
+  
+         CALL bdy_dyn_frs( kt )
+
+         IF ( ln_bdy_fla ) THEN
+
+            ua_e(:,:)=0.0
+            va_e(:,:)=0.0
+
+            ! Set these variables for use in bdy_dyn_fla
+            hu_e(:,:) = hu(:,:)
+            hv_e(:,:) = hv(:,:)
+
+            DO jk = 1, jpkm1
+               !! Vertically integrated momentum trends
+               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+            END DO
+
+            DO jk = 1 , jpkm1
+               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
+               va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
+            END DO
+
+            CALL bdy_dta_bt( kt+1, 0)
+            CALL bdy_dyn_fla
+
+         ENDIF
+
+         DO jk = 1 , jpkm1
+            ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
+            va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
+         END DO
 
       ENDIF
Index: /trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90	(revision 911)
@@ -11,4 +11,5 @@
    !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
    !!            " "  !  06-07  (S. Masson)  distributed restart using iom
+   !!            " "  !  05-01  (J.Chanut, A.Sellar) Calls to BDY routines. 
    !!----------------------------------------------------------------------
 #if defined key_dynspg_flt   ||   defined key_esopa  
@@ -36,4 +37,7 @@
    USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
+   USE bdy_oce         ! Unstructured open boundaries condition
+   USE bdydyn          ! Unstructured open boundaries condition (bdy_dyn routine) 
+   USE bdyvol          ! Unstructured open boundaries condition (bdy_vol routine)
    USE cla_dynspg      ! cross land advection
    USE in_out_manager  ! I/O manager
@@ -249,4 +253,13 @@
       CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
 #endif
+#if defined key_bdy
+      ! Update velocities on unstructured boundary using the Flow Relaxation Scheme
+      CALL bdy_dyn_frs( kt )
+
+      IF (ln_bdy_vol) THEN
+      ! Correction of the barotropic component velocity to control the volume of the system
+        CALL bdy_vol( kt )
+      ENDIF
+#endif
 #if defined key_agrif
       CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces 
@@ -377,4 +390,10 @@
             spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
             spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
+#elif defined key_bdy
+            ! caution : grad D = 0 along open boundaries
+            ! Remark: The filtering force could be reduced here in the FRS zone
+            !         by multiplying spgu/spgv by (1-alpha) ??  
+            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
+            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)           
 #else
             spgu(ji,jj) = z2dt * ztdgu
Index: /trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90	(revision 911)
@@ -46,6 +46,7 @@
          un_b  , vn_b                     ! vertically integrated horizontal velocities (now)
       REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop
-         sshn_e, ssha_e,               &  ! sea surface heigth (now,after)
-         ua_e  , va_e                     ! vertically integrated horizontal velocities (after)
+         sshb_e, sshn_e, ssha_e,       &  ! sea surface heigth (before,now,after)
+         ua_e  , va_e,                 &  ! vertically integrated horizontal velocities (after)
+         hu_e  , hv_e                     ! depth arrays for the barotropic solution
 #endif
 
Index: /trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90	(revision 911)
@@ -7,5 +7,6 @@
    !!             " "  !  05-11  (V. Garnier, G. Madec)  optimization
    !!             9.0  !  06-08  (S. Masson)  distributed restart using iom
-   !!             9.0  !  08-01  (R. Benshila)  change averaging method
+   !!              "   !  08-01  (R. Benshila)  change averaging method
+   !!              "   !  07-07  (D. Storkey) calls to BDY routines
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -30,4 +31,9 @@
    USE obc_oce         ! Lateral open boundary condition
    USE obc_par         ! open boundary condition parameters
+   USE bdy_oce         ! unstructured open boundaries
+   USE bdy_par         ! unstructured open boundaries
+   USE bdydta          ! unstructured open boundaries
+   USE bdydyn          ! unstructured open boundaries
+   USE bdytides        ! tidal forcing at unstructured open boundaries.
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -102,9 +108,9 @@
          zua, zva, zub, zvb,                &  !     "        "
          zssha_b, zua_b, zva_b,             &  !     "        "
-         zsshb_e, zub_e, zvb_e,             &  !     "        "
+         zub_e, zvb_e,                      &  !     "        "
          zun_e, zvn_e                          !     "        "
       !! Variable volume
       REAL(wp), DIMENSION(jpi,jpj)     ::   &
-         zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e         ! 2D workspace
+         zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfse3un_e, zfse3vn_e  ! 3D workspace
       !!----------------------------------------------------------------------
@@ -130,4 +136,6 @@
          ua_e(:,:)   = un_b(:,:)
          va_e(:,:)   = vn_b(:,:)
+         hu_e(:,:)   = hu(:,:)
+         hv_e(:,:)   = hv(:,:)
 
          IF( ln_dynvor_een ) THEN
@@ -300,5 +308,5 @@
 
       ! variables for the barotropic equations
-      zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
+      sshb_e (:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
       sshn_e (:,:) = sshn_b(:,:)
       zub_e  (:,:) = un_b  (:,:)       ! barotropic transports issued from the barotropic equations (before and now)
@@ -309,9 +317,9 @@
       zua_b  (:,:) = 0.e0
       zva_b  (:,:) = 0.e0
+      hu_e   (:,:) = hu   (:,:)     ! (barotropic) ocean depth at u-point
+      hv_e   (:,:) = hv   (:,:)     ! (barotropic) ocean depth at v-point
       IF( lk_vvl ) THEN
          zsshun_e(:,:)    = sshu (:,:)     ! (barotropic) sea surface height (now) at u-point
          zsshvn_e(:,:)    = sshv (:,:)     ! (barotropic) sea surface height (now) at v-point
-         hu_e(:,:)        = hu   (:,:)     ! (barotropic) ocean depth at u-point
-         hv_e(:,:)        = hv   (:,:)     ! (barotropic) ocean depth at v-point
          zfse3un_e(:,:,:) = fse3u(:,:,:)   ! (barotropic) scale factors  at u-point
          zfse3un_e(:,:,:) = fse3v(:,:,:)   ! (barotropic) scale factors  at v-point
@@ -338,7 +346,9 @@
          ! Time interpolation of open boundary condition data
          IF( lk_obc )   CALL obc_dta_bt( kt, jit )
+         IF( lk_bdy .or. lk_bdy_tides)   CALL bdy_dta_bt( kt, jit+1 )
 
          ! Horizontal divergence of barotropic transports
          !--------------------------------------------------
+         zhdiv(:,:) = 0.e0
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -360,9 +370,17 @@
 #endif
 
+         IF( lk_bdy .or. lk_bdy_tides ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
+               END DO
+            END DO
+         ENDIF
+
          ! Sea surface height from the barotropic system
          !----------------------------------------------
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
+               ssha_e(ji,jj) = ( sshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
             &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
             END DO
@@ -457,5 +475,5 @@
          !         - Correct the barotropic transports
          IF( lk_obc )   CALL obc_fla_ts
-
+         IF( lk_bdy .or. lk_bdy_tides )   CALL bdy_dyn_fla
 
          ! ... Boundary conditions on ua_e, va_e, ssha_e
@@ -475,5 +493,5 @@
          ! ---------------------------------------
          IF( neuler == 0 .AND. jit == 1 ) THEN   ! Euler (forward) time stepping
-            zsshb_e(:,:) = sshn_e(:,:)
+            sshb_e (:,:) = sshn_e(:,:)
             zub_e  (:,:) = zun_e (:,:)
             zvb_e  (:,:) = zvn_e (:,:)
@@ -482,5 +500,5 @@
             zvn_e  (:,:) = va_e  (:,:)
          ELSE                                        ! Asselin filtering
-            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
+            sshb_e (:,:) = atfp * ( sshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
             zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:)
             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:)
@@ -571,4 +589,12 @@
       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
 #endif
+
+      IF ( lk_bdy .or. lk_bdy_tides ) THEN
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
+           END DO
+         END DO
+      ENDIF
 
       ! sea surface height
Index: /trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/DYN/wzvmod.F90	(revision 911)
@@ -7,4 +7,5 @@
    !!             7.0  !  96-01  (G. Madec)  Statement function for e3
    !!             8.5  !  02-07  (G. Madec)  Free form, F90
+   !!              "   !  07-07  (D. Storkey) Zero zhdiv at open boundary (BDY) 
    !!----------------------------------------------------------------------
    !!   wzv        : Compute the vertical velocity
@@ -18,4 +19,5 @@
    USE prtctl          ! Print control
    USE phycst
+   USE bdy_oce         ! unstructured open boundaries
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
 
@@ -54,4 +56,5 @@
 
       !! * Local declarations
+      INTEGER  ::     jgrd, jb           ! temporary scalars
       INTEGER  ::           jk           ! dummy loop indices
       !! Variable volume
@@ -88,4 +91,5 @@
          ! Horizontal divergence of barotropic transports
          !--------------------------------------------------
+         zhdiv(:,:) = 0.e0
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! vector opt.
@@ -105,4 +109,15 @@
          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
+#endif
+
+#if defined key_bdy || defined key_bdy_tides
+         jgrd=1 !: tracer grid.
+         DO jb = 1, nblenrim(jgrd)
+           ji = nbi(jb,jgrd)
+           jj = nbj(jb,jgrd)
+           zhdiv(ji, jj) = 0.e0
+           zhdiv(ji, jj) = 0.e0
+           zhdiv(ji, jj) = 0.e0
+         END DO
 #endif
 
Index: /trunk/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/IOM/iom.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/IOM/iom.F90	(revision 911)
@@ -6,4 +6,5 @@
    !! History :  9.0  ! 05 12  (J. Belier) Original code
    !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
+   !!             "   ! 07 07  (D. Storkey) Changes to iom_gettime
    !!--------------------------------------------------------------------
    !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes
@@ -646,5 +647,5 @@
 
 
-   SUBROUTINE iom_gettime( kiomid, cdvar, ptime )
+   SUBROUTINE iom_gettime( kiomid, ptime, cdvar, kntime, cdunits, cdcalendar )
       !!--------------------------------------------------------------------
       !!                   ***  SUBROUTINE iom_gettime  ***
@@ -652,15 +653,30 @@
       !! ** Purpose : read the time axis cdvar in the file 
       !!--------------------------------------------------------------------
-      INTEGER               , INTENT(in   ) ::   kiomid   ! file Identifier
-      CHARACTER(len=*)      , INTENT(in   ) ::   cdvar    ! time axis name
-      REAL(wp), DIMENSION(:), INTENT(  out) ::   ptime    ! the time axis
-      !
+      INTEGER                    , INTENT(in   ) ::   kiomid     ! file Identifier
+      REAL(wp), DIMENSION(:)     , INTENT(  out) ::   ptime      ! the time axis
+      CHARACTER(len=*), OPTIONAL , INTENT(in   ) ::   cdvar      ! time axis name
+      INTEGER         , OPTIONAL , INTENT(  out) ::   kntime     ! number of times in file
+      CHARACTER(len=*), OPTIONAL , INTENT(  out) ::   cdunits    ! units attribute of time coordinate
+      CHARACTER(len=*), OPTIONAL , INTENT(  out) ::   cdcalendar ! calendar attribute of 
+      !
+      INTEGER, DIMENSION(1) :: kdimsz
       INTEGER            ::   idvar    ! id of the variable
+      CHARACTER(LEN=32)  ::   tname    ! local name of time coordinate
       CHARACTER(LEN=100) ::   clinfo   ! info character
       !---------------------------------------------------------------------
       !
-      IF( kiomid > 0 ) THEN
-         clinfo = 'iom_gettime, file: '//trim(iom_file(kiomid)%name)//', var: '//trim(cdvar)
-         idvar = iom_varid( kiomid, cdvar )
+      IF ( PRESENT(cdvar) ) THEN
+         tname = cdvar
+      ELSE
+         tname = iom_file(kiomid)%uldname
+      ENDIF
+      IF( kiomid > 0 ) THEN
+         clinfo = 'iom_gettime, file: '//trim(iom_file(kiomid)%name)//', var: '//trim(tname)
+         IF ( PRESENT(kntime) ) THEN
+            idvar  = iom_varid( kiomid, tname, kdimsz = kdimsz )
+            kntime = kdimsz(1)
+         ELSE
+            idvar = iom_varid( kiomid, tname )
+         ENDIF
          !
          ptime(:) = 0. ! default definition
@@ -670,6 +686,6 @@
                   IF( iom_file(kiomid)%dimsz(1,idvar) == size(ptime) ) THEN
                      SELECT CASE (iom_file(kiomid)%iolib)
-                     CASE (jpioipsl )   ;   CALL iom_ioipsl_gettime( kiomid, idvar, ptime )
-                     CASE (jpnf90   )   ;   CALL iom_nf90_gettime(   kiomid, idvar, ptime )
+                     CASE (jpioipsl )   ;   CALL iom_ioipsl_gettime( kiomid, idvar, ptime, cdunits, cdcalendar )
+                     CASE (jpnf90   )   ;   CALL iom_nf90_gettime(   kiomid, idvar, ptime, cdunits, cdcalendar )
                      CASE (jprstdimg)   ;   CALL ctl_stop( TRIM(clinfo)//' case IO library == jprstdimg not coded...' )
                      CASE DEFAULT    
Index: /trunk/NEMO/OPA_SRC/IOM/iom_def.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 911)
@@ -5,4 +5,5 @@
    !!====================================================================
    !! History :  9.0  ! 06 09  (S. Masson) Original code
+   !!             "   ! 07 07  (D. Storkey) Add uldname
    !!--------------------------------------------------------------------
    !!---------------------------------------------------------------------------------
@@ -57,4 +58,5 @@
       INTEGER                                   ::   iduld    !: id of the unlimited dimension
       INTEGER                                   ::   irec     !: writing record position  
+      CHARACTER(LEN=32)                         ::   uldname  !: name of the unlimited dimension
       CHARACTER(LEN=32), DIMENSION(jpmax_vars)  ::   cn_var   !: names of the variables
       INTEGER, DIMENSION(jpmax_vars)            ::   nvid     !: id of the variables
Index: /trunk/NEMO/OPA_SRC/IOM/iom_ioipsl.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/IOM/iom_ioipsl.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/IOM/iom_ioipsl.F90	(revision 911)
@@ -6,4 +6,5 @@
    !! History :  9.0  ! 05 12  (J. Belier) Original code
    !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
+   !!             "   ! 07 07  (D. Storkey) Changes to iom_ioipsl_gettime
    !!--------------------------------------------------------------------
    !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes
@@ -254,5 +255,5 @@
 
 
-   SUBROUTINE iom_ioipsl_gettime( kiomid, kvid, ptime )
+   SUBROUTINE iom_ioipsl_gettime( kiomid, kvid, ptime, cdunits, cdcalendar )
       !!--------------------------------------------------------------------
       !!                   ***  SUBROUTINE iom_gettime  ***
@@ -260,11 +261,19 @@
       !! ** Purpose : read the time axis kvid in the file with IOIPSL (only fliocom module)
       !!--------------------------------------------------------------------
-      INTEGER               , INTENT(in   ) ::   kiomid   ! file Identifier
-      INTEGER               , INTENT(in   ) ::   kvid     ! variable id
-      REAL(wp), DIMENSION(:), INTENT(  out) ::   ptime    ! the time axis
+      INTEGER                   , INTENT(in   ) ::   kiomid     ! file Identifier
+      INTEGER                   , INTENT(in   ) ::   kvid       ! variable id
+      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   ptime      ! the time axis
+      CHARACTER(len=*), OPTIONAL, INTENT(  out) ::   cdunits    ! units attribute
+      CHARACTER(len=*), OPTIONAL, INTENT(  out) ::   cdcalendar ! calendar attribute
       !---------------------------------------------------------------------
       !
       CALL fliogetv( iom_file(kiomid)%nfid, TRIM(iom_file(kiomid)%cn_var(kvid)), ptime(:),   &
             &         start=(/ 1 /), count=(/ iom_file(kiomid)%dimsz(1, kvid) /) )
+      IF ( PRESENT(cdunits) ) THEN 
+         CALL fliogeta( iom_file(kiomid)%nfid, TRIM(iom_file(kiomid)%cn_var(kvid)), "units", cdunits )
+      ENDIF
+      IF ( PRESENT(cdcalendar) ) THEN 
+         CALL fliogeta( iom_file(kiomid)%nfid, TRIM(iom_file(kiomid)%cn_var(kvid)), "calendar", cdcalendar )
+      ENDIF
       !
    END SUBROUTINE iom_ioipsl_gettime
Index: /trunk/NEMO/OPA_SRC/IOM/iom_nf90.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 911)
@@ -6,4 +6,5 @@
    !! History :  9.0  ! 05 12  (J. Belier) Original code
    !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
+   !!             "   ! 07 07  (D. Storkey) Changes to iom_nf90_gettime
    !!--------------------------------------------------------------------
    !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes
@@ -121,4 +122,8 @@
          iom_file(kiomid)%irec   = -1   ! useless for NetCDF files, used to know if the file is in define mode 
          CALL iom_nf90_check(NF90_Inquire(if90id, unlimitedDimId = iom_file(kiomid)%iduld), clinfo)
+         IF ( iom_file(kiomid)%iduld .GE. 0 ) THEN
+           CALL iom_nf90_check(NF90_Inquire_Dimension(if90id, iom_file(kiomid)%iduld,   &
+        &                                               name = iom_file(kiomid)%uldname), clinfo)
+         ENDIF
          IF(lwp) WRITE(numout,*) '                   ---> '//TRIM(cdname)//' OK'
       ELSE
@@ -272,5 +277,5 @@
 
 
-   SUBROUTINE iom_nf90_gettime( kiomid, kvid, ptime )
+   SUBROUTINE iom_nf90_gettime( kiomid, kvid, ptime, cdunits, cdcalendar )
       !!--------------------------------------------------------------------
       !!                   ***  SUBROUTINE iom_gettime  ***
@@ -278,7 +283,9 @@
       !! ** Purpose : read the time axis kvid in the file with NF90
       !!--------------------------------------------------------------------
-      INTEGER               , INTENT(in   ) ::   kiomid   ! file Identifier
-      INTEGER               , INTENT(in   ) ::   kvid     ! variable id
-      REAL(wp), DIMENSION(:), INTENT(  out) ::   ptime    ! the time axis
+      INTEGER                   , INTENT(in   ) ::   kiomid     ! file Identifier
+      INTEGER                   , INTENT(in   ) ::   kvid       ! variable id
+      REAL(wp), DIMENSION(:)    , INTENT(  out) ::   ptime      ! the time axis
+      CHARACTER(len=*), OPTIONAL, INTENT(  out) ::   cdunits    ! units attribute
+      CHARACTER(len=*), OPTIONAL, INTENT(  out) ::   cdcalendar ! calendar attribute
       !
       CHARACTER(LEN=100) ::   clinfo     ! info character
@@ -287,4 +294,12 @@
       CALL iom_nf90_check(NF90_GET_VAR(iom_file(kiomid)%nfid, iom_file(kiomid)%nvid(kvid), ptime(:),   &
             &                           start=(/ 1 /), count=(/ iom_file(kiomid)%dimsz(1, kvid) /)), clinfo)
+      IF ( PRESENT(cdunits) ) THEN 
+         CALL iom_nf90_check(NF90_GET_ATT(iom_file(kiomid)%nfid, iom_file(kiomid)%nvid(kvid), "units", &
+            &                           values=cdunits), clinfo)
+      ENDIF
+      IF ( PRESENT(cdcalendar) ) THEN 
+         CALL iom_nf90_check(NF90_GET_ATT(iom_file(kiomid)%nfid, iom_file(kiomid)%nvid(kvid), "calendar", &
+            &                           values=cdcalendar), clinfo)
+      ENDIF
       !
    END SUBROUTINE iom_nf90_gettime
Index: /trunk/NEMO/OPA_SRC/OBC/obcdta.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/OBC/obcdta.F90	(revision 911)
@@ -8,4 +8,5 @@
    !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Free surface, F90
    !!   9.0  !  04-06 (F. Durand, A-M. Treguier) Netcdf BC files on input
+   !!    "   ! 07-07  (D. Storkey) Change to arguments for iom_gettime
    !!------------------------------------------------------------------------------
 #if defined key_obc
@@ -159,5 +160,5 @@
                IF ( lp_obc_east ) THEN
                   IF ( llnot_done ) THEN
-                     CALL iom_gettime ( id_e, TRIM(cl_vname), tcobc )
+                     CALL iom_gettime ( id_e, tcobc, TRIM(cl_vname) )
                      llnot_done = .FALSE.
                   ENDIF
@@ -166,5 +167,5 @@
                IF ( lp_obc_west ) THEN
                   IF ( llnot_done ) THEN
-                     CALL iom_gettime ( id_w, TRIM(cl_vname), tcobc )
+                     CALL iom_gettime ( id_w, tcobc, TRIM(cl_vname) )
                      llnot_done = .FALSE.
                  ENDIF
@@ -173,5 +174,5 @@
                IF ( lp_obc_north ) THEN
                   IF ( llnot_done ) THEN
-                     CALL iom_gettime ( id_n, TRIM(cl_vname), tcobc )
+                     CALL iom_gettime ( id_n, tcobc, TRIM(cl_vname) )
                      llnot_done = .FALSE.
                   ENDIF
@@ -180,5 +181,5 @@
                IF ( lp_obc_south ) THEN
                   IF ( llnot_done ) THEN
-                     CALL iom_gettime ( id_s, TRIM(cl_vname), tcobc )
+                     CALL iom_gettime ( id_s, tcobc, TRIM(cl_vname) )
                      llnot_done = .FALSE.
                   ENDIF
Index: /trunk/NEMO/OPA_SRC/TRA/tranxt.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/TRA/tranxt.F90	(revision 911)
@@ -27,4 +27,5 @@
    USE domvvl          ! variable volume
    USE obctra          ! open boundary condition (obc_tra routine)
+   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine)
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -159,4 +160,16 @@
       ! Update tracers on open boundaries.
       CALL obc_tra( kt )
+
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+#elif defined key_bdy
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+
+      ! Update tracers on open boundaries.
+      CALL bdy_tra( kt )
+
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
Index: /trunk/NEMO/OPA_SRC/opa.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/opa.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/opa.F90	(revision 911)
@@ -34,4 +34,5 @@
    !!    "   !  06-03  (L. Debreu, C. Mazauric)  Agrif implementation
    !!    "   !  06-04  (G. Madec, R. Benshila)  Step reorganization
+   !!    "   !  07-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -50,4 +51,7 @@
    USE obc_par         ! open boundary cond. parameters
    USE obcini          ! open boundary cond. initialization (obc_ini routine)
+   USE bdy_par         ! unstructured open boundary cond. parameters
+   USE bdyini          ! unstructured open boundary cond. initialization (bdy_init routine)
+   USE bdytides        ! tides at open boundaries initialization (lk_bdy_tides)
    USE istate          ! initial state setting          (istate_init routine)
    USE eosbn2          ! equation of state            (eos bn2 routine)
@@ -269,4 +273,7 @@
       IF( lk_obc    )   CALL obc_init       ! Open boundaries 
 
+      IF( lk_bdy .OR. lk_bdy_tides )   &
+         &              CALL bdy_init       ! Unstructured open boundaries
+
       CALL istate_init                      ! ocean initial state (Dynamics and tracers)
 
Index: /trunk/NEMO/OPA_SRC/step.F90
===================================================================
--- /trunk/NEMO/OPA_SRC/step.F90	(revision 910)
+++ /trunk/NEMO/OPA_SRC/step.F90	(revision 911)
@@ -20,4 +20,5 @@
    !!             " "  !  06-07  (S. Masson)  restart using iom
    !!             " "  !  06-08  (G. Madec)  surface module 
+   !!             " "  !  07-07  (J. Chanut, A. Sellar) Unstructured open boundaries (BDY)
    !!----------------------------------------------------------------------
 
@@ -74,4 +75,7 @@
    USE obcrad          ! open boundary cond. radiation    (obc_rad routine)
    USE obcspg          ! open boundary cond  spg          (obc_spg routine)
+
+   USE bdy_par         ! unstructured open boundary data variables
+   USE bdydta          ! unstructured open boundary data  (bdy_dta routine)
 
    USE divcur          ! hor. divergence and curl      (div & cur routines)
@@ -187,4 +191,6 @@
       IF( lk_obc     )   CALL obc_rad( kstp )         ! compute phase velocities at open boundaries
 
+      IF( lk_bdy     )   CALL bdy_dta( kstp )         ! update dynamic and tracer data at unstructured open boundary
+
       IF( ninist == 1 ) THEN                       ! Output the initial state and forcings
                         CALL dia_wri_state( 'output.init' )
