Index: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3868)
+++ branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 3970)
@@ -29,4 +29,8 @@
    USE iom             ! IOM library
    USE in_out_manager  ! I/O logical units
+   ! bg jchanut tschanges
+   USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag
+   ! end jchanut tschanges
+
 #if defined key_lim2
    USE ice_2
@@ -314,8 +318,10 @@
                      END DO
                   ENDIF
-                  IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
-                     CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  &
-                                        & td=tides(ib_bdy), time_offset=time_offset )
-                  ENDIF
+                  ! bg jchanut tschanges
+                  !IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                  !   CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  &
+                  !                      & td=tides(ib_bdy), time_offset=time_offset )
+                  !ENDIF
+                  ! end jchanut tschanges
                ENDIF
             ENDIF
@@ -323,4 +329,11 @@
          END IF ! nn_dta(ib_bdy) = 1
       END DO  ! ib_bdy
+
+      ! bg jchanut tschanges
+#if defined key_tide
+      ! Add tides if not split-explicit free surface else this is done in ts loop
+      IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset )
+#endif
+      ! end jchanut tschanges
 
       IF ( ln_apr_obc ) THEN
@@ -476,5 +489,5 @@
             IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 
 
-               IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
+               IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading
                   jfld = jfld + 1
                   blf_i(jfld) = bn_ssh
@@ -572,6 +585,6 @@
             ! Recalculate field counts
             !-------------------------
-            nb_bdy_fld_sum = 0
             IF( ib_bdy .eq. 1 ) THEN 
+               nb_bdy_fld_sum = 0
                nb_bdy_fld(ib_bdy) = jfld
                nb_bdy_fld_sum     = jfld              
@@ -616,5 +629,5 @@
                ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) )
                ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) )
-               IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN
+               IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN
                   jfld = jfld + 1
                   dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1)
Index: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 3868)
+++ branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 3970)
@@ -30,4 +30,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  !
+   USE domvvl          ! variable volume
 
    IMPLICIT NONE
@@ -84,13 +85,25 @@
       pu2d(:,:) = 0.e0
       pv2d(:,:) = 0.e0
-      DO jk = 1, jpkm1   !! Vertically integrated momentum trends
-          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
-          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
-      END DO
-      pu2d(:,:) = pu2d(:,:) * phur(:,:)
-      pv2d(:,:) = pv2d(:,:) * phvr(:,:)
+      ! bg jchanut tschanges (not specifically related to ts; this is a bug)
+      IF (lk_vvl) THEN
+         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+            pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+            pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+         END DO
+         pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) )
+         pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) )
+      ! end jchanut tschanges
+      ELSE
+         DO jk = 1, jpkm1   !! Vertically integrated momentum trends
+            pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+            pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+         END DO
+         pu2d(:,:) = pu2d(:,:) * phur(:,:)
+         pv2d(:,:) = pv2d(:,:) * phvr(:,:)
+      ENDIF
+
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:)
-         va(:,:,jk) = va(:,:,jk) - pv2d(:,:)
+         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk)
+         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk)
       END DO
 
Index: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90
===================================================================
--- branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 3868)
+++ branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90	(revision 3970)
@@ -6,4 +6,5 @@
    !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite
    !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!----------------------------------------------------------------------
 #if defined key_bdy 
@@ -27,4 +28,5 @@
 
    PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn
+   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv
 
    !!----------------------------------------------------------------------
@@ -135,7 +137,10 @@
       REAL(wp) ::   zcorr                            ! Flather correction
       REAL(wp) ::   zforc                            ! temporary scalar
+      REAL(wp) ::   zflag, z1_2                      !    "        "
       !!----------------------------------------------------------------------
 
       IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla')
+
+      z1_2 = 0.5_wp
 
       ! ---------------------------------!
@@ -164,6 +169,11 @@
          !
          zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
-         zforc = dta%u2d(jb)
-         pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
+         ! bg jchanut tschanges: Set zflag to 0 below to revert to Flather scheme
+!!         zforc = dta%u2d(jb)
+         zflag = ABS(idx%flagu(jb))
+         iim1 = ii + idx%flagu(jb)
+         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pu2d(iim1,ij)
+         pu2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1) 
+         ! end jchanut tschanges
       END DO
       !
@@ -177,6 +187,11 @@
          !
          zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
-         zforc = dta%v2d(jb)
-         pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
+         ! bg jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme
+!!         zforc = dta%v2d(jb)
+         zflag = ABS(idx%flagv(jb))
+         ijm1 = ij + idx%flagv(jb)
+         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pv2d(ii,ijm1)
+         pv2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1)
+         ! end jchanut tschanges
       END DO
       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated
@@ -186,4 +201,51 @@
       !
    END SUBROUTINE bdy_dyn2d_fla
+
+   SUBROUTINE bdy_ssh( zssh )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_ssh  ***
+      !!
+      !! ** Purpose : Duplicate sea level across open boundaries
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level
+      !!
+      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers
+      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       "
+
+      igrd = 1                       ! Everything is at T-points here
+
+      DO ib_bdy = 1, nb_bdy
+         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+            ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+            ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+            ! Set gradient direction:
+            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
+            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
+            IF ( zcoef1+zcoef2 == 0 ) THEN
+               ! corner
+!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1)
+!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + &
+!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + &
+!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + &
+!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1)
+               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1)
+               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + &
+                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + &
+                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + &
+                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1)
+               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1)
+            ELSE
+               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
+               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1)
+            ENDIF
+         END DO
+
+         ! Boundary points should be updated
+         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy )
+      END DO
+
+   END SUBROUTINE bdy_ssh
 #else
    !!----------------------------------------------------------------------
@@ -192,8 +254,11 @@
 CONTAINS
    SUBROUTINE bdy_dyn2d( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt
+      INTEGER, intent(in) :: kt
+      WRITE(*,*) 'bdy_dyn2: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_dyn2d
+
 #endif
 
    !!======================================================================
 END MODULE bdydyn2d
+
Index: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3868)
+++ branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 3970)
@@ -161,4 +161,5 @@
 
       DO ib_bdy = 1,nb_bdy
+        icount = 0 ! flag to set max rimwidth to 1 if no relaxation
         IF(lwp) WRITE(numout,*) ' ' 
         IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
@@ -175,4 +176,5 @@
           CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
@@ -196,4 +198,5 @@
           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
@@ -215,4 +218,5 @@
               CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
+              icount = icount + 1
               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
@@ -228,4 +232,5 @@
           CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
           CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
@@ -248,4 +253,5 @@
               CALL ctl_stop( 'Use FRS OR relaxation' )
            ELSE
+              icount = icount + 1
               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
@@ -262,4 +268,5 @@
           CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'        
           CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
+          icount = icount + 1
           CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
         END SELECT
@@ -273,7 +280,10 @@
         IF(lwp) WRITE(numout,*)
 #endif
-
-        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
-        IF(lwp) WRITE(numout,*)
+        IF ( icount>0 ) THEN
+           IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
+           IF(lwp) WRITE(numout,*)
+        ELSE
+           nn_rimwidth(ib_bdy) = 1 ! no relaxation
+        ENDIF
 
       ENDDO
@@ -391,5 +401,4 @@
             ENDDO
             CALL iom_close( inum )
-
          ENDIF 
 
@@ -398,5 +407,4 @@
       IF (nb_bdy>0) THEN
          jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
-
          ! Allocate arrays
          !---------------
@@ -446,5 +454,5 @@
          ENDIF 
 
-      ENDDO      
+      ENDDO     
     
       ! 2. Now fill indices corresponding to straight open boundary arrays:
@@ -752,5 +760,6 @@
                ! check if point is in local domain
                IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
-                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
+                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
+                  & nbrdta(ib,igrd,ib_bdy) <= nn_rimwidth(ib_bdy)     ) THEN      
                   !
                   icount = icount  + 1
@@ -765,5 +774,7 @@
          ! Allocate index arrays for this boundary set
          !--------------------------------------------
-         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
+
+         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(1:jpbgrd))
+         ilen1 = MAX(1,ilen1)
          ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
@@ -773,5 +784,5 @@
          ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
          ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
-         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
+         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )      
 
          ! Dispatch mapping indices and discrete distances on each processor
@@ -1008,5 +1019,8 @@
       ! bdytmask = 1  on the computational domain AND on open boundaries
       !          = 0  elsewhere   
- 
+      bdytmask(:,:) = 1.e0
+      bdyumask(:,:) = 1.e0
+      bdyvmask(:,:) = 1.e0
+
       IF( ln_mask_file ) THEN
          CALL iom_open( cn_mask_file, inum )
Index: branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3868)
+++ branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 3970)
@@ -9,4 +9,5 @@
    !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
    !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
+   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
    !!----------------------------------------------------------------------
 #if defined key_bdy
@@ -32,4 +33,5 @@
 !   USE tide_mod       ! Useless ??
    USE fldread, ONLY: fld_map
+   USE dynspg_oce, ONLY: lk_dynspg_ts
 
    IMPLICIT NONE
@@ -38,4 +40,5 @@
    PUBLIC   bdytide_init     ! routine called in bdy_init
    PUBLIC   bdytide_update   ! routine called in bdy_dta
+   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
 
    TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
@@ -49,4 +52,5 @@
 
    TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
+   TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
 
    !!----------------------------------------------------------------------
@@ -252,4 +256,14 @@
             ENDIF
             !
+            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during
+                                     ! time splitting integration
+               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
+               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
+               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
+               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0
+               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0
+               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0
+            ENDIF
+            !
          ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
          !
@@ -318,5 +332,5 @@
          
       IF( PRESENT(jit) ) THEN  
-         z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) )
+         z_arg = ((kt-kt_tide) * rdt + (jit+0.5*(time_add-1)) * rdt / REAL(nn_baro,wp) )
       ELSE                              
          z_arg = ((kt-kt_tide)+time_add) * rdt
@@ -324,6 +338,6 @@
 
       ! Linear ramp on tidal component at open boundaries 
-      zramp = 1.
-      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.)
+      zramp = 1._wp
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
 
       DO itide = 1, nb_harmo
@@ -351,4 +365,130 @@
       !
    END SUBROUTINE bdytide_update
+
+   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_dta_tides  ***
+      !!                
+      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
+      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit
+                                                        ! is present then units = subcycle timesteps.
+                                                        ! time_offset = 0  => get data at "now"    time level
+                                                        ! time_offset = -1 => get data at "before" time level
+                                                        ! time_offset = +1 => get data at "after"  time level
+                                                        ! etc.
+      !!
+      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step
+      INTEGER,          DIMENSION(jpbgrd) :: ilen0 
+      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts
+      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices
+      INTEGER  :: time_add                                    ! time offset in units of timesteps
+      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist      
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides')
+
+      lk_first_btstp=.TRUE.
+      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
+
+      time_add = 0
+      IF( PRESENT(time_offset) ) THEN
+         time_add = time_offset
+      ENDIF
+      
+      ! Absolute time from model initialization:   
+      IF( PRESENT(kit) ) THEN  
+         z_arg = ( kt + (kit+0.5*(time_add-1)) / REAL(nn_baro,wp) ) * rdt
+      ELSE                              
+         z_arg = ( kt + time_add ) * rdt
+      ENDIF
+
+      ! Linear ramp on tidal component at open boundaries 
+      zramp = 1.
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
+
+      DO ib_bdy = 1,nb_bdy
+
+         ! line below should be simplified (runoff case)
+         IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN
+
+            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
+            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
+
+            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN
+               ilen0(:)=nblen(:)
+            ELSE
+               ilen0(:)=nblenrim(:)
+            ENDIF     
+
+            ! We refresh nodal factors every day below
+            ! This should be done somewhere else
+            IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. lk_first_btstp ) THEN
+               !
+               kt_tide = kt               
+               !
+               IF(lwp) THEN
+               WRITE(numout,*)
+               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
+               WRITE(numout,*) '~~~~~~~~~~~~~~ '
+               ENDIF
+               !
+               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
+               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
+               !
+            ENDIF
+            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
+            !
+            ! If time splitting, save data at first barotropic iteration
+            IF ( PRESENT(kit) ) THEN
+               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data:
+                  dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1))
+                  dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2))
+                  dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3))
+
+               ELSE ! Initialize arrays from slow varying open boundary data:            
+                  dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
+                  dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
+                  dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
+               ENDIF
+            ENDIF
+            !
+            ! Update open boundary data arrays:
+            DO itide = 1, nb_harmo
+               !
+               z_sarg = (z_arg + zoff) * omega_tide(itide)
+               z_cost = zramp * COS( z_sarg )
+               z_sist = zramp * SIN( z_sarg )
+               !
+               igrd=1                              ! SSH on tracer grid
+               DO ib = 1, ilen0(igrd)
+                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
+                     &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
+                     &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
+               END DO
+               !
+               igrd=2                              ! U grid
+               DO ib = 1, ilen0(igrd)
+                  dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
+                     &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
+                     &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
+               END DO
+               !
+               igrd=3                              ! V grid
+               DO ib = 1, ilen0(igrd) 
+                  dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
+                     &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
+                     &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
+               END DO
+            END DO
+         END IF
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides')
+      !
+   END SUBROUTINE bdy_dta_tides
 
    SUBROUTINE tide_init_elevation( idx, td )
@@ -457,6 +597,10 @@
       WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit
    END SUBROUTINE bdytide_update
+   SUBROUTINE bdy_dta_tides( kt, jit )   ! Empty routine
+      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit
+   END SUBROUTINE bdy_dta_tides
 #endif
 
    !!======================================================================
 END MODULE bdytides
+
