Index: /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdydta.F90	(revision 12058)
+++ /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdydta.F90	(revision 12059)
@@ -65,5 +65,5 @@
 CONTAINS
 
-      SUBROUTINE bdy_dta( kt, time_offset )
+      SUBROUTINE bdy_dta( kt, jit, time_offset )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta  ***
@@ -75,5 +75,7 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
-      INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps.
+      INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option)
+      INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
+      !                                               ! is present then units = subcycle timesteps.
       !                                               ! time_offset = 0 => get data at "now" time level
       !                                               ! time_offset = -1 => get data at "before" time level
@@ -92,5 +94,5 @@
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt == nit000 ) THEN
+      IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN
 
          ! Calculate depth-mean currents
@@ -226,84 +228,143 @@
          IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required
       
-            IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition
-               jend = nb_bdy_fld(jbdy)
-               CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
-                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
-               !
-               igrd = 2                      ! zonal velocity
-               DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                  dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
-               END DO
-               !
-               igrd = 3                      ! meridional velocity
-               DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
-                  ii   = idx_bdy(jbdy)%nbi(ib,igrd)
-                  ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                  dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
-               END DO
+            IF( PRESENT(jit) ) THEN
+               ! Update barotropic boundary conditions only
+               ! jit is optional argument for fld_read and bdytide_update
+               IF( cn_dyn2d(jbdy) /= 'none' ) THEN
+                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
+                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
+                     IF( dta%ll_u3d ) dta%v2d(:) = 0._wp
+                  ENDIF
+                  IF (cn_tra(jbdy) /= 'runoff') THEN
+                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN
+
+                        jend = jstart + dta%nread(2) - 1
+                        IF( ln_full_vel_array(jbdy) ) THEN
+                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
+                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy,   &
+                                     & fvl=ln_full_vel_array(jbdy)  )
+                        ELSE
+                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
+                                     & kit=jit, kt_offset=time_offset  )
+                        ENDIF
+
+                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields
+                        IF( ln_full_vel_array(jbdy) .AND.                                             &
+                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  &
+                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN
+
+                           igrd = 2                      ! zonal velocity
+                           dta%u2d(:) = 0._wp
+                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                              DO ik = 1, jpkm1
+                                 dta%u2d(ib) = dta%u2d(ib) &
+                       &                          + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
+                              END DO
+                              dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
+                           END DO
+                           igrd = 3                      ! meridional velocity
+                           dta%v2d(:) = 0._wp
+                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                              ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                              ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                              DO ik = 1, jpkm1
+                                 dta%v2d(ib) = dta%v2d(ib) &
+                       &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
+                              END DO
+                              dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
+                           END DO
+                        ENDIF                    
+                     ENDIF
+                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing
+                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   & 
+                          &                 jit=jit, time_offset=time_offset )
+                     ENDIF
+                  ENDIF
+               ENDIF
             ELSE
-               IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
-                  IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
-                  IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
-                  IF( dta%ll_v2d ) dta%v2d(:) = 0._wp
-               ENDIF
-               IF( dta%nread(1) .gt. 0 ) THEN ! update external data
-                  jend = jstart + dta%nread(1) - 1
-                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
-                     & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   &
-                     & fvl=ln_full_vel_array(jbdy) )
-               ENDIF
-               ! If full velocities in boundary data then split into barotropic and baroclinic data
-               IF( ln_full_vel_array(jbdy) .and.                                             &
-                  & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. &
-                  &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN
+               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition
+                  jend = nb_bdy_fld(jbdy)
+                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  &
+                               & map=nbmap_ptr(jstart:jend), kt_offset=time_offset )
+                  !
                   igrd = 2                      ! zonal velocity
-                  dta%u2d(:) = 0._wp
                   DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      ii   = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                     DO ik = 1, jpkm1
-                        dta%u2d(ib) = dta%u2d(ib) &
-              &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
-                     END DO
-                     dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
-                     DO ik = 1, jpkm1
-                        dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib)
-                     END DO
+                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
                   END DO
+                  !
                   igrd = 3                      ! meridional velocity
-                  dta%v2d(:) = 0._wp
                   DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
                      ii   = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij   = idx_bdy(jbdy)%nbj(ib,igrd)
-                     DO ik = 1, jpkm1
-                        dta%v2d(ib) = dta%v2d(ib) &
-              &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
+                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+                  END DO
+               ELSE
+                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays
+                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp
+                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp
+                     IF( dta%ll_v2d ) dta%v2d(:) = 0._wp
+                  ENDIF
+                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data
+                     jend = jstart + dta%nread(1) - 1
+                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), &
+                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   &
+                                  & fvl=ln_full_vel_array(jbdy) )
+                  ENDIF
+                  ! If full velocities in boundary data then split into barotropic and baroclinic data
+                  IF( ln_full_vel_array(jbdy) .and.                                             &
+                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. &
+                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN
+                     igrd = 2                      ! zonal velocity
+                     dta%u2d(:) = 0._wp
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                        DO ik = 1, jpkm1
+                           dta%u2d(ib) = dta%u2d(ib) &
+                 &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik)
+                        END DO
+                        dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij)
+                        DO ik = 1, jpkm1
+                           dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib)
+                        END DO
                      END DO
-                     dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
-                     DO ik = 1, jpkm1
-                        dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib)
+                     igrd = 3                      ! meridional velocity
+                     dta%v2d(:) = 0._wp
+                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd)
+                        ii   = idx_bdy(jbdy)%nbi(ib,igrd)
+                        ij   = idx_bdy(jbdy)%nbj(ib,igrd)
+                        DO ik = 1, jpkm1
+                           dta%v2d(ib) = dta%v2d(ib) &
+                 &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik)
+                        END DO
+                        dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij)
+                        DO ik = 1, jpkm1
+                           dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib)
+                        END DO
                      END DO
-                  END DO
-               ENDIF
-
-            ENDIF
-#if defined key_si3
-            ! convert N-cat fields (input) into jpl-cat (output)
-            IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN
-               jfld_hti = jfld_htit(jbdy)
-               jfld_hts = jfld_htst(jbdy)
-               jfld_ai  = jfld_ait(jbdy)
-               IF    ( jpl /= 1 .AND. nice_cat == 1 ) THEN                       ! case input cat = 1
-                  CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
-                     &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
-               ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl
-                  CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), &
-                     &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
-               ENDIF
-            ENDIF
-#endif
+                  ENDIF
+
+               ENDIF
+#if defined key_si3
+               ! convert N-cat fields (input) into jpl-cat (output)
+               IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN
+                  jfld_hti = jfld_htit(jbdy)
+                  jfld_hts = jfld_htst(jbdy)
+                  jfld_ai  = jfld_ait(jbdy)
+                  IF    ( jpl /= 1 .AND. nice_cat == 1 ) THEN                       ! case input cat = 1
+                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
+                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
+                  ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl
+                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), &
+                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    )
+                  ENDIF
+               ENDIF
+#endif
+            ENDIF
             jstart = jstart + dta%nread(1)
          ENDIF    ! nn_dta(jbdy) = 1
Index: /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdytides.F90
===================================================================
--- /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdytides.F90	(revision 12058)
+++ /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdytides.F90	(revision 12059)
@@ -30,4 +30,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
 
@@ -262,4 +263,93 @@
 
 
+   SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdytide_update  ***
+      !!                
+      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 
+      !!                
+      !!----------------------------------------------------------------------
+      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
+      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
+      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
+      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
+      INTEGER, OPTIONAL, INTENT(in   ) ::   jit         ! Barotropic timestep counter (for timesplitting option)
+      INTEGER, OPTIONAL, INTENT(in   ) ::   time_offset ! time offset in units of timesteps. NB. if jit
+      !                                                 ! is present then units = subcycle timesteps.
+      !                                                 ! time_offset = 0  => get data at "now"    time level
+      !                                                 ! time_offset = -1 => get data at "before" time level
+      !                                                 ! time_offset = +1 => get data at "after"  time level
+      !                                                 ! etc.
+      !
+      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
+      INTEGER  ::   time_add              ! time offset in units of timesteps
+      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
+      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars    
+      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
+      !!----------------------------------------------------------------------
+      !
+      ilen0(1) =  SIZE(td%ssh(:,1,1))
+      ilen0(2) =  SIZE(td%u(:,1,1))
+      ilen0(3) =  SIZE(td%v(:,1,1))
+
+      zflag=1
+      IF ( PRESENT(jit) ) THEN
+        IF ( jit /= 1 ) zflag=0
+      ENDIF
+
+      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
+        !
+        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
+        !
+        IF(lwp) THEN
+           WRITE(numout,*)
+           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
+           WRITE(numout,*) '~~~~~~~~~~~~~~ '
+        ENDIF
+        !
+        CALL tide_init_elevation ( idx, td )
+        CALL tide_init_velocities( idx, td )
+        !
+      ENDIF 
+
+      time_add = 0
+      IF( PRESENT(time_offset) ) THEN
+         time_add = time_offset
+      ENDIF
+         
+      IF( PRESENT(jit) ) THEN  
+         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
+      ELSE                              
+         z_arg = ((kt-kt_tide)+time_add) * rdt
+      ENDIF
+
+      ! Linear ramp on tidal component at open boundaries 
+      zramp = 1._wp
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rn_tide_ramp_dt*rday),0._wp),1._wp)
+
+      DO itide = 1, nb_harmo
+         z_sarg = z_arg * tide_harmonics(itide)%omega
+         z_cost(itide) = COS( z_sarg )
+         z_sist(itide) = SIN( z_sarg )
+      END DO
+
+      DO itide = 1, nb_harmo
+         igrd=1                              ! SSH on tracer grid
+         DO ib = 1, ilen0(igrd)
+            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
+         END DO
+         igrd=2                              ! U grid
+         DO ib = 1, ilen0(igrd)
+            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
+         END DO
+         igrd=3                              ! V grid
+         DO ib = 1, ilen0(igrd) 
+            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
+         END DO
+      END DO
+      !
+   END SUBROUTINE bdytide_update
+
+
    SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
       !!----------------------------------------------------------------------
