Index: NEMO/trunk/src/OCE/BDY/bdy_oce.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdy_oce.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdy_oce.F90	(revision 12377)
@@ -141,4 +141,6 @@
    LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lsend_bdyext   !: mark needed communication for given boundary, grid and neighbour
    LOGICAL, ALLOCATABLE, DIMENSION(:,:,:,:) ::   lrecv_bdyext   !:  when searching towards the exterior of the computational domain
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/trunk/src/OCE/BDY/bdydta.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdydta.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdydta.F90	(revision 12377)
@@ -23,5 +23,5 @@
    USE phycst         ! physical constants
    USE sbcapr         ! atmospheric pressure forcing
-   USE sbctide        ! Tidal forcing or not
+   USE tide_mod, ONLY: ln_tide ! tidal forcing
    USE bdy_oce        ! ocean open boundary conditions  
    USE bdytides       ! tidal forcing at boundaries
@@ -68,4 +68,6 @@
 !$AGRIF_END_DO_NOT_TREAT
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -75,5 +77,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_dta( kt, kit, kt_offset )
+   SUBROUTINE bdy_dta( kt, Kmm )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta  ***
@@ -85,11 +87,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
-      INTEGER, INTENT(in), OPTIONAL ::   kit          ! subcycle time-step index (for timesplitting option)
-      INTEGER, INTENT(in), OPTIONAL ::   kt_offset    ! time offset in units of timesteps. NB. if kit
-      !                                               ! is present then units = subcycle timesteps.
-      !                                               ! kt_offset = 0 => get data at "now" time level
-      !                                               ! kt_offset = -1 => get data at "before" time level
-      !                                               ! kt_offset = +1 => get data at "after" time level
-      !                                               ! etc.
+      INTEGER, INTENT(in)           ::   Kmm          ! ocean time level index
       !
       INTEGER ::  jbdy, jfld, jstart, jend, ib, jl    ! dummy loop indices
@@ -105,5 +101,5 @@
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt == nit000 .AND. .NOT.PRESENT(kit) ) THEN
+      IF( kt == nit000 ) THEN
 
          ! Calculate depth-mean currents
@@ -122,5 +118,5 @@
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%ssh(ib) = sshn(ii,ij) * tmask(ii,ij,1)         
+                     dta_bdy(jbdy)%ssh(ib) = ssh(ii,ij,Kmm) * tmask(ii,ij,1)         
                   END DO
                ENDIF
@@ -130,5 +126,5 @@
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%u2d(ib) = un_b(ii,ij) * umask(ii,ij,1)         
+                     dta_bdy(jbdy)%u2d(ib) = uu_b(ii,ij,Kmm) * umask(ii,ij,1)         
                   END DO
                   igrd = 3
@@ -136,5 +132,5 @@
                      ii = idx_bdy(jbdy)%nbi(ib,igrd)
                      ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                     dta_bdy(jbdy)%v2d(ib) = vn_b(ii,ij) * vmask(ii,ij,1)         
+                     dta_bdy(jbdy)%v2d(ib) = vv_b(ii,ij,Kmm) * vmask(ii,ij,1)         
                   END DO
                ENDIF
@@ -149,5 +145,5 @@
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%u3d(ib,ik) =  ( un(ii,ij,ik) - un_b(ii,ij) ) * umask(ii,ij,ik)         
+                        dta_bdy(jbdy)%u3d(ib,ik) =  ( uu(ii,ij,ik,Kmm) - uu_b(ii,ij,Kmm) ) * umask(ii,ij,ik)         
                      END DO
                   END DO
@@ -157,5 +153,5 @@
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vn(ii,ij,ik) - vn_b(ii,ij) ) * vmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%v3d(ib,ik) =  ( vv(ii,ij,ik,Kmm) - vv_b(ii,ij,Kmm) ) * vmask(ii,ij,ik)         
                      END DO
                   END DO
@@ -171,6 +167,6 @@
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%tem(ib,ik) = tsn(ii,ij,ik,jp_tem) * tmask(ii,ij,ik)         
-                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * tmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%tem(ib,ik) = ts(ii,ij,ik,jp_tem,Kmm) * tmask(ii,ij,ik)         
+                        dta_bdy(jbdy)%sal(ib,ik) = ts(ii,ij,ik,jp_sal,Kmm) * tmask(ii,ij,ik)         
                      END DO
                   END DO
@@ -216,6 +212,6 @@
          ! read/update all bdy data
          ! ------------------------
-         CALL fld_read( kt, 1, bf_alias, kit = kit, kt_offset = kt_offset )
-
+         ! BDY: use pt_offset=0.5 as applied at the end of the step and fldread is referenced at the middle of the step
+         CALL fld_read( kt, 1, bf_alias, pt_offset = 0.5_wp, Kmm = Kmm )
          ! apply some corrections in some specific cases...
          ! --------------------------------------------------
@@ -254,7 +250,7 @@
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                DO ik = 1, jpkm1
-                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
+                  dta_alias%u2d(ib) = dta_alias%u2d(ib) + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta_alias%u3d(ib,ik)
                END DO
-               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu_n(ii,ij)
+               dta_alias%u2d(ib) =  dta_alias%u2d(ib) * r1_hu(ii,ij,Kmm)
                DO ik = 1, jpkm1            ! compute barocline zonal velocity and put it in u3d
                   dta_alias%u3d(ib,ik) = dta_alias%u3d(ib,ik) - dta_alias%u2d(ib)
@@ -267,7 +263,7 @@
                ij   = idx_bdy(jbdy)%nbj(ib,igrd)
                DO ik = 1, jpkm1
-                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
+                  dta_alias%v2d(ib) = dta_alias%v2d(ib) + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta_alias%v3d(ib,ik)
                END DO
-               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv_n(ii,ij)
+               dta_alias%v2d(ib) =  dta_alias%v2d(ib) * r1_hv(ii,ij,Kmm)
                DO ik = 1, jpkm1            ! compute barocline meridional velocity and put it in v3d
                   dta_alias%v3d(ib,ik) = dta_alias%v3d(ib,ik) - dta_alias%v2d(ib)
@@ -275,10 +271,4 @@
             END DO
          ENDIF   ! ltotvel
-
-         ! update tidal harmonic forcing
-         IF( PRESENT(kit) .AND. nn_dyn2d_dta(jbdy) .GE. 2 ) THEN
-            CALL bdytide_update( kt = kt, idx = idx_bdy(jbdy), dta = dta_alias, td = tides(jbdy),   & 
-               &                 kit = kit, kt_offset = kt_offset )
-         ENDIF
 
          !  atm surface pressure : add inverted barometer effect to ssh if it was read
@@ -343,19 +333,22 @@
                   nblen => idx_bdy(jbdy)%nblen
                   nblenrim => idx_bdy(jbdy)%nblenrim
-                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN ; ilen1(:)=nblen(:) ; ELSE ; ilen1(:)=nblenrim(:) ; ENDIF 
-                     IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
-                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
-                     IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
-                  ENDIF
-               END DO
-            ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
-               !
-               CALL bdy_dta_tides( kt=kt, kt_offset=kt_offset )
-            ENDIF
+                  IF( cn_dyn2d(jbdy) == 'frs' ) THEN   ;   ilen1(:)=nblen(:)
+                  ELSE                                 ;   ilen1(:)=nblenrim(:)
+                  ENDIF
+                  IF ( dta_bdy(jbdy)%lneed_ssh   ) dta_bdy_s(jbdy)%ssh(1:ilen1(1)) = dta_bdy(jbdy)%ssh(1:ilen1(1))
+                  IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%u2d(1:ilen1(2)) = dta_bdy(jbdy)%u2d(1:ilen1(2))
+                  IF ( dta_bdy(jbdy)%lneed_dyn2d ) dta_bdy_s(jbdy)%v2d(1:ilen1(3)) = dta_bdy(jbdy)%v2d(1:ilen1(3))
+               ENDIF
+            END DO
+         ELSE ! Add tides if not split-explicit free surface else this is done in ts loop
+            !
+            ! BDY: use pt_offset=1.0 as applied at the end of the step and bdy_dta_tides is referenced at the middle of the step
+            CALL bdy_dta_tides( kt=kt, pt_offset = 1._wp )
          ENDIF
-         !
-         IF( ln_timing )   CALL timing_stop('bdy_dta')
-         !
-      END SUBROUTINE bdy_dta
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop('bdy_dta')
+      !
+   END SUBROUTINE bdy_dta
 
 
@@ -373,4 +366,6 @@
       INTEGER ::   ierror, ios     ! 
       !
+      INTEGER ::   nbdy_rdstart, nbdy_loc
+      CHARACTER(LEN=50)                      ::   cerrmsg       ! error string
       CHARACTER(len=3)                       ::   cl3           ! 
       CHARACTER(len=100)                     ::   cn_dir        ! Root directory for location of data files
@@ -415,5 +410,5 @@
       ! Read namelists
       ! --------------
-      REWIND(numnam_cfg)
+      nbdy_rdstart = 1
       DO jbdy = 1, nb_bdy
 
@@ -421,6 +416,5 @@
          WRITE(ctmp2, '(a,i2)') 'block nambdy_dta number ', jbdy
 
-         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we do a rewind 
-         REWIND(numnam_ref)
+         ! There is only one nambdy_dta block in namelist_ref -> use it for each bdy so we read from the beginning
          READ  ( numnam_ref, nambdy_dta, IOSTAT = ios, ERR = 901)
 901      IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_dta in reference namelist' )
@@ -431,6 +425,19 @@
             & .OR. ( dta_bdy(jbdy)%lneed_tra   .AND.       nn_tra_dta(jbdy)    == 1 )   &
             & .OR. ( dta_bdy(jbdy)%lneed_ice   .AND.       nn_ice_dta(jbdy)    == 1 )   )   THEN
-            ! WARNING: we don't do a rewind here, each bdy reads its own nambdy_dta block one after another
-            READ  ( numnam_cfg, nambdy_dta, IOSTAT = ios, ERR = 902 )
+            !
+            ! Need to support possibility of reading more than one
+            ! nambdy_dta from the namelist_cfg internal file.
+            ! Do this by finding the jbdy'th occurence of nambdy_dta in the
+            ! character buffer as the starting point.
+            !
+            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_dta' )
+            IF( nbdy_loc .GT. 0 ) THEN
+               nbdy_rdstart = nbdy_rdstart + nbdy_loc
+            ELSE
+               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',jbdy,' of nambdy_dta not found'
+               ios = -1
+               CALL ctl_nam ( ios , cerrmsg )
+            ENDIF
+            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_dta, IOSTAT = ios, ERR = 902)
 902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_dta in configuration namelist' )
             IF(lwm) WRITE( numond, nambdy_dta )           
@@ -442,9 +449,11 @@
             IF( nn_ice_dta(jbdy) == 1 ) THEN   ! if we get ice bdy data from netcdf file
                CALL fld_fill(  bf(jp_bdya_i,jbdy:jbdy), bn_a_i, cn_dir, 'bdy_dta', 'a_i'//' '//ctmp1, ctmp2 )   ! use namelist info
-               CALL fld_clopn( bf(jp_bdya_i,jbdy), nyear, nmonth, nday )   ! not a problem when we call it again after
+               CALL fld_def( bf(jp_bdya_i,jbdy) )
+               CALL iom_open( bf(jp_bdya_i,jbdy)%clname, bf(jp_bdya_i,jbdy)%num )
                idvar = iom_varid( bf(jp_bdya_i,jbdy)%num, bf(jp_bdya_i,jbdy)%clvar, kndims=indims, kdimsz=i4dimsz, lduld=lluld )
                IF( indims == 4 .OR. ( indims == 3 .AND. .NOT. lluld ) ) THEN   ;   ipl = i4dimsz(3)   ! xylt or xyl
                ELSE                                                            ;   ipl = 1            ! xy or xyt
                ENDIF
+               CALL iom_close( bf(jp_bdya_i,jbdy)%num )
                bf(jp_bdya_i,jbdy)%clrootname = 'NOT USED'   ! reset to default value as this subdomain may not need to read this bdy
             ENDIF
Index: NEMO/trunk/src/OCE/BDY/bdydyn.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdydyn.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdydyn.F90	(revision 12377)
@@ -37,5 +37,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn( kt, dyn3d_only )
+   SUBROUTINE bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn  ***
@@ -44,10 +44,12 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)           ::   kt           ! Main time step counter
-      LOGICAL, INTENT(in), OPTIONAL ::   dyn3d_only   ! T => only update baroclinic velocities
+      INTEGER                             , INTENT(in)    ::   kt           ! Main time step counter
+      INTEGER                             , INTENT(in)    ::   Kbb, Kaa     ! Ocean time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::   puu, pvv     ! Ocean velocities (to be updated at open boundaries)
+      LOGICAL, OPTIONAL                   , INTENT(in)    ::   dyn3d_only   ! T => only update baroclinic velocities
       !
       INTEGER ::   jk, ii, ij, ib_bdy, ib, igrd     ! Loop counter
       LOGICAL ::   ll_dyn2d, ll_dyn3d, ll_orlanski
-      REAL(wp), DIMENSION(jpi,jpj) :: pua2d, pva2d     ! after barotropic velocities
+      REAL(wp), DIMENSION(jpi,jpj) :: zua2d, zva2d     ! after barotropic velocities
       !!----------------------------------------------------------------------
       !
@@ -70,16 +72,16 @@
 
       !                          ! "After" velocities: 
-      pua2d(:,:) = 0._wp
-      pva2d(:,:) = 0._wp     
+      zua2d(:,:) = 0._wp
+      zva2d(:,:) = 0._wp     
       DO jk = 1, jpkm1
-         pua2d(:,:) = pua2d(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
-         pva2d(:,:) = pva2d(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
+         zua2d(:,:) = zua2d(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
+         zva2d(:,:) = zva2d(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
       END DO
-      pua2d(:,:) = pua2d(:,:) * r1_hu_a(:,:)
-      pva2d(:,:) = pva2d(:,:) * r1_hv_a(:,:)
+      zua2d(:,:) = zua2d(:,:) * r1_hu(:,:,Kaa)
+      zva2d(:,:) = zva2d(:,:) * r1_hv(:,:,Kaa)
 
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) - pua2d(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) - pva2d(:,:) ) * vmask(:,:,jk)
+         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zua2d(:,:) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zva2d(:,:) ) * vmask(:,:,jk)
       END DO
 
@@ -87,6 +89,6 @@
       IF( ll_orlanski ) THEN     ! "Before" velocities (Orlanski condition only) 
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = ( ub(:,:,jk) - ub_b(:,:) ) * umask(:,:,jk)
-            vb(:,:,jk) = ( vb(:,:,jk) - vb_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) - uu_b(:,:,Kbb) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) - vv_b(:,:,Kbb) ) * vmask(:,:,jk)
          END DO
       ENDIF
@@ -97,7 +99,7 @@
       !-------------------------------------------------------
 
-      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, pua2d, pva2d, ub_b, vb_b, r1_hu_a(:,:), r1_hv_a(:,:), ssha )
+      IF( ll_dyn2d )   CALL bdy_dyn2d( kt, zua2d, zva2d, uu_b(:,:,Kbb), vv_b(:,:,Kbb), r1_hu(:,:,Kaa), r1_hv(:,:,Kaa), ssh(:,:,Kaa) )
 
-      IF( ll_dyn3d )   CALL bdy_dyn3d( kt )
+      IF( ll_dyn3d )   CALL bdy_dyn3d( kt, Kbb, puu, pvv, Kaa )
 
       !-------------------------------------------------------
@@ -106,12 +108,12 @@
       !
       DO jk = 1 , jpkm1
-         ua(:,:,jk) = ( ua(:,:,jk) + pua2d(:,:) ) * umask(:,:,jk)
-         va(:,:,jk) = ( va(:,:,jk) + pva2d(:,:) ) * vmask(:,:,jk)
+         puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) + zua2d(:,:) ) * umask(:,:,jk)
+         pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) + zva2d(:,:) ) * vmask(:,:,jk)
       END DO
       !
       IF ( ll_orlanski ) THEN
          DO jk = 1 , jpkm1
-            ub(:,:,jk) = ( ub(:,:,jk) + ub_b(:,:) ) * umask(:,:,jk)
-            vb(:,:,jk) = ( vb(:,:,jk) + vb_b(:,:) ) * vmask(:,:,jk)
+            puu(:,:,jk,Kbb) = ( puu(:,:,jk,Kbb) + uu_b(:,:,Kbb) ) * umask(:,:,jk)
+            pvv(:,:,jk,Kbb) = ( pvv(:,:,jk,Kbb) + vv_b(:,:,Kbb) ) * vmask(:,:,jk)
          END DO
       END IF
Index: NEMO/trunk/src/OCE/BDY/bdydyn3d.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdydyn3d.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdydyn3d.F90	(revision 12377)
@@ -33,5 +33,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_dyn3d( kt )
+   SUBROUTINE bdy_dyn3d( kt, Kbb, puu, pvv, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d  ***
@@ -40,5 +40,7 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! Main time step counter
+      INTEGER                             , INTENT( in    ) ::   kt        ! Main time step counter
+      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
       !
       INTEGER  ::   ib_bdy, ir     ! BDY set index, rim index
@@ -58,13 +60,13 @@
             CASE('none')        ;   CYCLE
             CASE('frs' )        ! treat the whole boundary at once
-               IF( ir == 0) CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+                       IF( ir == 0) CALL bdy_dyn3d_frs( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
             CASE('specified')   ! treat the whole rim      at once
-               IF( ir == 0) CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+                       IF( ir == 0) CALL bdy_dyn3d_spe( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
             CASE('zero')        ! treat the whole rim      at once
-               IF( ir == 0) CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-            CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.false. )
-            CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.true.  )
-            CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy, llrim0 )
-            CASE('neumann')     ;   CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy, llrim0 )
+                       IF( ir == 0) CALL bdy_dyn3d_zro( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+            CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.false. )
+            CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, llrim0, ll_npo=.true.  )
+            CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy, llrim0 )
+            CASE('neumann')     ;   CALL bdy_dyn3d_nmn( puu, pvv, Kaa, idx_bdy(ib_bdy), ib_bdy, llrim0 )
             CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
             END SELECT
@@ -97,8 +99,8 @@
          !
          IF( ANY(llsend2) .OR. ANY(llrecv2) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', ua, 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
+            CALL lbc_lnk( 'bdydyn2d', puu(:,:,:,Kaa), 'U', -1., kfillmode=jpfillnothing ,lsend=llsend2, lrecv=llrecv2 )
          END IF
          IF( ANY(llsend3) .OR. ANY(llrecv3) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdydyn2d', va, 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
+            CALL lbc_lnk( 'bdydyn2d', pvv(:,:,:,Kaa), 'V', -1., kfillmode=jpfillnothing ,lsend=llsend3, lrecv=llrecv3 )
          END IF
       END DO   ! ir
@@ -107,5 +109,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
+   SUBROUTINE bdy_dyn3d_spe( puu, pvv, Kaa, idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
@@ -115,8 +117,10 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER        , INTENT(in) ::   kt      ! time step index
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
-      INTEGER        , INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   kt        ! Time step
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -129,5 +133,5 @@
             ii   = idx%nbi(jb,igrd)
             ij   = idx%nbj(jb,igrd)
-            ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk)
+            puu(ii,ij,jk,Kaa) = dta%u3d(jb,jk) * umask(ii,ij,jk)
          END DO
       END DO
@@ -138,5 +142,5 @@
             ii   = idx%nbi(jb,igrd)
             ij   = idx%nbj(jb,igrd)
-            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
+            pvv(ii,ij,jk,Kaa) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
          END DO
       END DO
@@ -145,5 +149,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt, ib_bdy, llrim0 )
+   SUBROUTINE bdy_dyn3d_zgrad( puu, pvv, Kaa, idx, dta, kt, ib_bdy, llrim0 )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
@@ -152,9 +156,11 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER                     ::   kt
-      TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta      ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy   ! BDY set index
-      LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   kt
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
+      LOGICAL                             , INTENT( in    ) ::   llrim0   ! indicate if rim 0 is treated
       !!
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -178,5 +184,5 @@
             !
             DO jk = 1, jpkm1
-               ua(ii,ij,jk) = ua(ii,ij+flagv,jk) * umask(ii,ij+flagv,jk)
+               puu(ii,ij,jk,Kaa) = puu(ii,ij+flagv,jk,Kaa) * umask(ii,ij+flagv,jk)
             END DO
             !
@@ -198,5 +204,5 @@
             !
             DO jk = 1, jpkm1
-               va(ii,ij,jk) = va(ii+flagu,ij,jk) * vmask(ii+flagu,ij,jk)
+               pvv(ii,ij,jk,Kaa) = pvv(ii+flagu,ij,jk,Kaa) * vmask(ii+flagu,ij,jk)
             END DO
             !
@@ -207,5 +213,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
+   SUBROUTINE bdy_dyn3d_zro( puu, pvv, Kaa, idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
@@ -214,8 +220,10 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER        , INTENT(in) ::   kt      ! time step index
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER                             , INTENT( in    ) ::   kt        ! time step index
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   ib, ik         ! dummy loop indices
@@ -228,5 +236,5 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            ua(ii,ij,ik) = 0._wp
+            puu(ii,ij,ik,Kaa) = 0._wp
          END DO
       END DO
@@ -237,5 +245,5 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            va(ii,ij,ik) = 0._wp
+            pvv(ii,ij,ik,Kaa) = 0._wp
          END DO
       END DO
@@ -244,5 +252,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
+   SUBROUTINE bdy_dyn3d_frs( puu, pvv, Kaa, idx, dta, kt, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
@@ -255,8 +263,10 @@
       !!               topography. Tellus, 365-382.
       !!----------------------------------------------------------------------
-      INTEGER        , INTENT(in) ::   kt      ! time step index
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA) , INTENT(in) ::   dta     ! OBC external data
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      INTEGER                             , INTENT( in    ) ::   kt        ! time step index
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -271,5 +281,5 @@
             ij   = idx%nbj(jb,igrd)
             zwgt = idx%nbw(jb,igrd)
-            ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta%u3d(jb,jk) - ua(ii,ij,jk) ) ) * umask(ii,ij,jk)
+            puu(ii,ij,jk,Kaa) = ( puu(ii,ij,jk,Kaa) + zwgt * ( dta%u3d(jb,jk) - puu(ii,ij,jk,Kaa) ) ) * umask(ii,ij,jk)
          END DO
       END DO
@@ -281,5 +291,5 @@
             ij   = idx%nbj(jb,igrd)
             zwgt = idx%nbw(jb,igrd)
-            va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta%v3d(jb,jk) - va(ii,ij,jk) ) ) * vmask(ii,ij,jk)
+            pvv(ii,ij,jk,Kaa) = ( pvv(ii,ij,jk,Kaa) + zwgt * ( dta%v3d(jb,jk) - pvv(ii,ij,jk,Kaa) ) ) * vmask(ii,ij,jk)
          END DO
       END DO   
@@ -288,5 +298,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, llrim0, ll_npo )
+   SUBROUTINE bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx, dta, ib_bdy, llrim0, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
@@ -298,27 +308,29 @@
       !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)    
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
-      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data
-      INTEGER,                      INTENT(in) ::   ib_bdy   ! BDY set index
-      LOGICAL,                      INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
-      LOGICAL,                      INTENT(in) ::   ll_npo   ! switch for NPO version
+      INTEGER                             , INTENT( in    ) ::   Kbb, Kaa  ! Time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      TYPE(OBC_DATA)                      , INTENT( in    ) ::   dta       ! OBC external data
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
+      LOGICAL                             , INTENT( in    ) ::   llrim0    ! indicate if rim 0 is treated
+      LOGICAL                             , INTENT( in    ) ::   ll_npo    ! switch for NPO version
 
       INTEGER  ::   jb, igrd                               ! dummy loop indices
       !!----------------------------------------------------------------------
       !
-      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
+      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities. 
       !
       igrd = 2      ! Orlanski bc on u-velocity; 
       !            
-      CALL bdy_orlanski_3d( idx, igrd, ub, ua, dta%u3d, ll_npo, llrim0 )
+      CALL bdy_orlanski_3d( idx, igrd, puu(:,:,:,Kbb), puu(:,:,:,Kaa), dta%u3d, ll_npo, llrim0 )
 
       igrd = 3      ! Orlanski bc on v-velocity
       !  
-      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo, llrim0 )
+      CALL bdy_orlanski_3d( idx, igrd, pvv(:,:,:,Kbb), pvv(:,:,:,Kaa), dta%v3d, ll_npo, llrim0 )
       !
    END SUBROUTINE bdy_dyn3d_orlanski
 
 
-   SUBROUTINE bdy_dyn3d_dmp( kt )
+   SUBROUTINE bdy_dyn3d_dmp( kt, Kbb, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
@@ -327,5 +339,7 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! time step index
+      INTEGER                             , INTENT( in    ) ::   kt        ! time step
+      INTEGER                             , INTENT( in    ) ::   Kbb, Krhs ! Time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities and trends (to be updated at open boundaries)
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -345,6 +359,6 @@
                zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
                DO jk = 1, jpkm1
-                  ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
-                                   ub(ii,ij,jk) + ub_b(ii,ij)) ) * umask(ii,ij,jk)
+                  puu(ii,ij,jk,Krhs) = ( puu(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - &
+                                   puu(ii,ij,jk,Kbb) + uu_b(ii,ij,Kbb)) ) * umask(ii,ij,jk)
                END DO
             END DO
@@ -356,6 +370,6 @@
                zwgt = idx_bdy(ib_bdy)%nbd(jb,igrd)
                DO jk = 1, jpkm1
-                  va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
-                                   vb(ii,ij,jk) + vb_b(ii,ij)) ) * vmask(ii,ij,jk)
+                  pvv(ii,ij,jk,Krhs) = ( pvv(ii,ij,jk,Krhs) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) -  &
+                                   pvv(ii,ij,jk,Kbb) + vv_b(ii,ij,Kbb)) ) * vmask(ii,ij,jk)
                END DO
             END DO
@@ -368,5 +382,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy, llrim0 )
+   SUBROUTINE bdy_dyn3d_nmn( puu, pvv, Kaa, idx, ib_bdy, llrim0 )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
@@ -377,19 +391,21 @@
       !!
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX), INTENT(in) ::   idx      ! OBC indices
-      INTEGER,         INTENT(in) ::   ib_bdy   ! BDY set index
-      LOGICAL,         INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
+      INTEGER                             , INTENT( in    ) ::   Kaa       ! Time level index
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT( inout ) ::   puu, pvv  ! Ocean velocities (to be updated at open boundaries)
+      TYPE(OBC_INDEX)                     , INTENT( in    ) ::   idx       ! OBC indices
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
+      LOGICAL                             , INTENT( in    ) ::   llrim0    ! indicate if rim 0 is treated
       INTEGER  ::   igrd                        ! dummy indice
       !!----------------------------------------------------------------------
       !
-      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
+      !! Note that at this stage the puu(:,:,:,Kbb) and puu(:,:,:,Kaa) arrays contain the baroclinic velocities. 
       !
       igrd = 2      ! Neumann bc on u-velocity; 
       !            
-      CALL bdy_nmn( idx, igrd, ua, llrim0 )   ! ua is masked
+      CALL bdy_nmn( idx, igrd, puu(:,:,:,Kaa), llrim0 )
 
       igrd = 3      ! Neumann bc on v-velocity
       !  
-      CALL bdy_nmn( idx, igrd, va, llrim0 )   ! va is masked
+      CALL bdy_nmn( idx, igrd, pvv(:,:,:,Kaa), llrim0 )
       !
    END SUBROUTINE bdy_dyn3d_nmn
Index: NEMO/trunk/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdyini.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdyini.F90	(revision 12377)
@@ -22,5 +22,5 @@
    USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
    USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
-   USE sbctide        ! Tidal forcing or not
+   USE tide_mod, ONLY: ln_tide ! tidal forcing
    USE phycst   , ONLY: rday
    !
@@ -75,5 +75,4 @@
       ! Read namelist parameters
       ! ------------------------
-      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries
       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' )
@@ -93,5 +92,4 @@
       cn_ice         (2:jp_bdy) = cn_ice         (1)
       nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
-      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
       READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
@@ -364,5 +362,4 @@
       ! -------------------------------------------------
 
-      REWIND( numnam_cfg )     
       nblendta(:,:) = 0
       nbdysege = 0
@@ -1080,12 +1077,26 @@
       INTEGER          ::   ios                 ! Local integer output status for namelist read
       INTEGER          ::   nbdyind, nbdybeg, nbdyend
+      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc
       CHARACTER(LEN=1) ::   ctypebdy   !     -        - 
+      CHARACTER(LEN=50)::   cerrmsg    !     -        - 
       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
       !!----------------------------------------------------------------------
-
-      ! No REWIND here because may need to read more than one nambdy_index namelist.
-      ! Read only namelist_cfg to avoid unseccessfull overwrite 
-      ! keep full control of the configuration namelist
-      READ  ( numnam_cfg, nambdy_index, IOSTAT = ios, ERR = 904 )
+      ! Need to support possibility of reading more than one nambdy_index from
+      ! the namelist_cfg internal file.
+      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the
+      ! character buffer as the starting point.
+      nbdy_rdstart = 1
+      DO nbdy_count = 1, kb_bdy
+       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
+       IF( nbdy_loc .GT. 0 ) THEN
+          nbdy_rdstart = nbdy_rdstart + nbdy_loc
+       ELSE
+          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
+          ios = -1
+          CALL ctl_nam ( ios , cerrmsg )
+       ENDIF
+      END DO
+      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
+      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
 904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
       IF(lwm) WRITE ( numond, nambdy_index )
Index: NEMO/trunk/src/OCE/BDY/bdylib.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdylib.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdylib.F90	(revision 12377)
@@ -35,5 +35,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_frs( idx, pta, dta )
+   SUBROUTINE bdy_frs( idx, phia, dta )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_frs  ***
@@ -45,5 +45,5 @@
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
       REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       !!
       REAL(wp) ::   zwgt           ! boundary weight
@@ -58,5 +58,5 @@
             ij = idx%nbj(ib,igrd)
             zwgt = idx%nbw(ib,igrd)
-            pta(ii,ij,ik) = ( pta(ii,ij,ik) + zwgt * (dta(ib,ik) - pta(ii,ij,ik) ) ) * tmask(ii,ij,ik)
+            phia(ii,ij,ik) = ( phia(ii,ij,ik) + zwgt * (dta(ib,ik) - phia(ii,ij,ik) ) ) * tmask(ii,ij,ik)
          END DO
       END DO
@@ -65,5 +65,5 @@
 
 
-   SUBROUTINE bdy_spe( idx, pta, dta )
+   SUBROUTINE bdy_spe( idx, phia, dta )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_spe  ***
@@ -74,5 +74,5 @@
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
       REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       !!
       INTEGER  ::   ib, ik, igrd   ! dummy loop indices
@@ -85,5 +85,5 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            pta(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik)
+            phia(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik)
          END DO
       END DO
@@ -92,5 +92,5 @@
 
 
-   SUBROUTINE bdy_orl( idx, ptb, pta, dta, lrim0, ll_npo )
+   SUBROUTINE bdy_orl( idx, phib, phia, dta, lrim0, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_orl  ***
@@ -102,6 +102,6 @@
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx  ! OBC indices
       REAL(wp), DIMENSION(:,:),            INTENT(in) ::   dta  ! OBC external data
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptb  ! before tracer field
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta  ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phib  ! before tracer field
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   phia  ! tracer trend
       LOGICAL                 , OPTIONAL,  INTENT(in) ::   lrim0   ! indicate if rim 0 is treated
       LOGICAL,                             INTENT(in) ::   ll_npo  ! switch for NPO version
@@ -112,5 +112,5 @@
       igrd = 1                       ! Everything is at T-points here
       !
-      CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, lrim0, ll_npo )
+      CALL bdy_orlanski_3d( idx, igrd, phib(:,:,:), phia(:,:,:), dta, lrim0, ll_npo )
       !
    END SUBROUTINE bdy_orl
Index: NEMO/trunk/src/OCE/BDY/bdytides.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdytides.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdytides.F90	(revision 12377)
@@ -18,5 +18,5 @@
    USE phycst         ! physical constants
    USE bdy_oce        ! ocean open boundary conditions
-   USE tideini        ! 
+   USE tide_mod       ! 
    USE daymod         ! calendar
    !
@@ -30,5 +30,4 @@
 
    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
 
@@ -45,4 +44,8 @@
    TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
 
+   INTEGER ::   kt_tide
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -64,5 +67,4 @@
       CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
       LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
-      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
       !!
       INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
@@ -71,4 +73,6 @@
       INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
       INTEGER                                   ::   ios                 ! Local integer output status for namelist read
+      INTEGER                                   ::   nbdy_rdstart, nbdy_loc
+      CHARACTER(LEN=50)                         ::   cerrmsg             !: error string
       CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file 
       REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
@@ -77,5 +81,5 @@
       TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
       !!
-      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
+      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
       !!----------------------------------------------------------------------
       !
@@ -84,6 +88,6 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
 
-      REWIND(numnam_cfg)
-
+
+      nbdy_rdstart = 1
       DO ib_bdy = 1, nb_bdy
          IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
@@ -94,9 +98,21 @@
             filtide(:) = ''
 
-            REWIND( numnam_ref )
             READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
 901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
-            ! Don't REWIND here - may need to read more than one of these namelists. 
-            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
+            !
+            ! Need to support possibility of reading more than one
+            ! nambdy_tide from the namelist_cfg internal file.
+            ! Do this by finding the ib_bdy'th occurence of nambdy_tide in the
+            ! character buffer as the starting point.
+            !
+            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' )
+            IF( nbdy_loc .GT. 0 ) THEN
+               nbdy_rdstart = nbdy_rdstart + nbdy_loc
+            ELSE
+               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found'
+               ios = -1
+               CALL ctl_nam ( ios , cerrmsg )
+            ENDIF
+            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_tide, IOSTAT = ios, ERR = 902)
 902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
             IF(lwm) WRITE ( numond, nambdy_tide )
@@ -105,10 +121,9 @@
             IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
             IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
-            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
             IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
             IF(lwp) THEN 
                     WRITE(numout,*) '             Tidal components: ' 
                DO itide = 1, nb_harmo
-                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
+                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide 
                END DO
             ENDIF 
@@ -151,6 +166,6 @@
                igrd = 1                       ! Everything is at T-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
-                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
+                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
+                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -168,6 +183,6 @@
                igrd = 2                       ! Everything is at U-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -185,6 +200,6 @@
                igrd = 3                       ! Everything is at V-points here
                DO itide = 1, nb_harmo
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
-                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
+                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
                   DO ib = 1, ilen0(igrd)
                      ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -210,5 +225,5 @@
                DO itide = 1, nb_harmo
                   !                                                              ! SSH fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
+                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
                   CALL iom_open( clfile, inum )
                   CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
@@ -218,5 +233,5 @@
                   CALL iom_close( inum )
                   !                                                              ! U fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
+                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
                   CALL iom_open( clfile, inum )
                   CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
@@ -226,5 +241,5 @@
                   CALL iom_close( inum )
                   !                                                              ! V fields
-                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
+                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
                   CALL iom_open( clfile, inum )
                   CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
@@ -240,10 +255,4 @@
             ENDIF ! ln_bdytide_2ddta=.true.
             !
-            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
-               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
-               td%u0  (:,:,2) = - td%u0  (:,:,2)
-               td%v0  (:,:,2) = - td%v0  (:,:,2)
-            ENDIF
-            !
             ! Allocate slow varying data in the case of time splitting:
             ! Do it anyway because at this stage knowledge of free surface scheme is unknown
@@ -262,94 +271,5 @@
 
 
-   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_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   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
-      !                                                 ! is present then units = subcycle timesteps.
-      !                                                 ! kt_offset = 0  => get data at "now"    time level
-      !                                                 ! kt_offset = -1 => get data at "before" time level
-      !                                                 ! kt_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(kit) ) THEN
-        IF ( kit /= 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(kt_offset) ) THEN
-         time_add = kt_offset
-      ENDIF
-         
-      IF( PRESENT(kit) ) THEN  
-         z_arg = ((kt-kt_tide) * rdt + (kit+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)/(rdttideramp*rday),0._wp),1._wp)
-
-      DO itide = 1, nb_harmo
-         z_sarg = z_arg * omega_tide(itide)
-         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, kt_offset )
+   SUBROUTINE bdy_dta_tides( kt, kit, pt_offset )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dta_tides  ***
@@ -360,17 +280,11 @@
       INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
       INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
-      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
-      !                                              ! is present then units = subcycle timesteps.
-      !                                              ! kt_offset = 0  => get data at "now"    time level
-      !                                              ! kt_offset = -1 => get data at "before" time level
-      !                                              ! kt_offset = +1 => get data at "after"  time level
-      !                                              ! etc.
+      REAL(wp),OPTIONAL, INTENT(in) ::   pt_offset   ! time offset in units of timesteps
       !
       LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
       INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
-      INTEGER  ::   time_add                  ! time offset in units of timesteps
       INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
       INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
-      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist      
+      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset   
       !!----------------------------------------------------------------------
       !
@@ -378,19 +292,17 @@
       IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
 
-      time_add = 0
-      IF( PRESENT(kt_offset) ) THEN
-         time_add = kt_offset
-      ENDIF
+      zt_offset = 0._wp
+      IF( PRESENT(pt_offset) )   zt_offset = pt_offset
       
       ! Absolute time from model initialization:   
       IF( PRESENT(kit) ) THEN  
-         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
+         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_baro, wp) ) * rdt
       ELSE                              
-         z_arg = ( kt + time_add ) * rdt
+         z_arg = ( REAL(kt, wp) + zt_offset ) * 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.)
+      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - REAL(nit000,wp)*rdt)/(rn_tide_ramp_dt*rday),0.),1.)
 
       DO ib_bdy = 1,nb_bdy
@@ -409,5 +321,5 @@
             IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
                !
-               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
+               kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rdt)/rdt)
                !
                IF(lwp) THEN
@@ -421,5 +333,5 @@
                !
             ENDIF
-            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
+            zoff = REAL(-kt_tide,wp) * rdt ! time offset relative to nodal factor computation time
             !
             ! If time splitting, initialize arrays from slow varying open boundary data:
@@ -433,5 +345,5 @@
             DO itide = 1, nb_harmo
                !
-               z_sarg = (z_arg + zoff) * omega_tide(itide)
+               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega
                z_cost = zramp * COS( z_sarg )
                z_sist = zramp * SIN( z_sarg )
@@ -491,6 +403,6 @@
          END DO
          DO ib = 1 , ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
          ENDDO
          DO ib = 1 , ilen0(igrd)
@@ -530,6 +442,6 @@
          END DO
          DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
          ENDDO
          DO ib = 1, ilen0(igrd)
@@ -551,6 +463,6 @@
          END DO
          DO ib = 1, ilen0(igrd)
-            mod_tide(ib)=mod_tide(ib)*ftide(itide)
-            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
+            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
+            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
          ENDDO
          DO ib = 1, ilen0(igrd)
Index: NEMO/trunk/src/OCE/BDY/bdytra.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdytra.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdytra.F90	(revision 12377)
@@ -40,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_tra( kt )
+   SUBROUTINE bdy_tra( kt, Kbb, pts, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_tra  ***
@@ -47,5 +47,7 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   ! Main time step counter
+      INTEGER                                  , INTENT(in)    :: kt        ! Main time step counter
+      INTEGER                                  , INTENT(in)    :: Kbb, Kaa  ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! tracer fields
       !
       INTEGER                        :: ib_bdy, jn, igrd, ir   ! Loop indeces
@@ -70,13 +72,13 @@
                CASE('none'        )   ;   CYCLE
                CASE('frs'         )   ! treat the whole boundary at once
-                  IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra )
+                  IF( ir == 0 ) CALL bdy_frs ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
                CASE('specified'   )   ! treat the whole rim      at once
-                  IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra )
-               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , tsa(:,:,:,jn), llrim0 )   ! tsa masked
-               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), &
+                  IF( ir == 0 ) CALL bdy_spe ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
+               CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , pts(:,:,:,jn,Kaa), llrim0 )   ! tsa masked
+               CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), &
                     & zdta(jn)%tra, llrim0, ll_npo=.false. )
-               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), &
+               CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), &
                     & zdta(jn)%tra, llrim0, ll_npo=.true.  )
-               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), jn, llrim0 )
+               CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), jn, llrim0 )
                CASE DEFAULT           ;   CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
                END SELECT
@@ -98,5 +100,5 @@
          END DO
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdytra', tsa, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+            CALL lbc_lnk( 'bdytra', pts(:,:,:,jn,Kaa), 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
          END IF
          !
@@ -106,5 +108,5 @@
 
 
-   SUBROUTINE bdy_rnf( idx, pta, jpa, llrim0 )
+   SUBROUTINE bdy_rnf( idx, pt, jpa, llrim0 )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_rnf  ***
@@ -116,15 +118,14 @@
       !!----------------------------------------------------------------------
       TYPE(OBC_INDEX),                     INTENT(in) ::   idx      ! OBC indices
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pta      ! tracer trend
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pt       ! tracer trend
       INTEGER,                             INTENT(in) ::   jpa      ! TRA index
       LOGICAL,                             INTENT(in) ::   llrim0   ! indicate if rim 0 is treated
       !
       INTEGER  ::   ib, ii, ij, igrd   ! dummy loop indices
-      INTEGER  ::   ik, ip, jp ! 2D addresses
       !!----------------------------------------------------------------------
       !
       igrd = 1                       ! Everything is at T-points here
       IF(      jpa == jp_tem ) THEN
-         CALL bdy_nmn( idx, igrd, pta, llrim0 )
+         CALL bdy_nmn( idx, igrd, pt, llrim0 )
       ELSE IF( jpa == jp_sal ) THEN
          IF( .NOT. llrim0 )   RETURN
@@ -132,5 +133,5 @@
             ii = idx%nbi(ib,igrd)
             ij = idx%nbj(ib,igrd)
-            pta(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1)
+            pt(ii,ij,1:jpkm1) = 0.1 * tmask(ii,ij,1:jpkm1)
          END DO
       END IF
@@ -139,5 +140,5 @@
 
 
-   SUBROUTINE bdy_tra_dmp( kt )
+   SUBROUTINE bdy_tra_dmp( kt, Kbb, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_tra_dmp  ***
@@ -146,5 +147,7 @@
       !! 
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt   !
+      INTEGER                                  , INTENT(in)    :: kt        ! time step
+      INTEGER                                  , INTENT(in)    :: Kbb, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation
       !
       REAL(wp) ::   zwgt           ! boundary weight
@@ -165,8 +168,8 @@
                zwgt = idx_bdy(ib_bdy)%nbd(ib,igrd)
                DO ik = 1, jpkm1
-                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) * tmask(ii,ij,ik)
-                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) * tmask(ii,ij,ik)
-                  tsa(ii,ij,ik,jp_tem) = tsa(ii,ij,ik,jp_tem) + zta
-                  tsa(ii,ij,ik,jp_sal) = tsa(ii,ij,ik,jp_sal) + zsa
+                  zta = zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - pts(ii,ij,ik,jp_tem,Kbb) ) * tmask(ii,ij,ik)
+                  zsa = zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - pts(ii,ij,ik,jp_sal,Kbb) ) * tmask(ii,ij,ik)
+                  pts(ii,ij,ik,jp_tem,Krhs) = pts(ii,ij,ik,jp_tem,Krhs) + zta
+                  pts(ii,ij,ik,jp_sal,Krhs) = pts(ii,ij,ik,jp_sal,Krhs) + zsa
                END DO
             END DO
Index: NEMO/trunk/src/OCE/BDY/bdyvol.F90
===================================================================
--- NEMO/trunk/src/OCE/BDY/bdyvol.F90	(revision 12148)
+++ NEMO/trunk/src/OCE/BDY/bdyvol.F90	(revision 12377)
@@ -14,7 +14,7 @@
    USE bdy_oce        ! ocean open boundary conditions
    USE sbc_oce        ! ocean surface boundary conditions
+   USE isf_oce, ONLY : fwfisf_cav, fwfisf_par  ! ice shelf
    USE dom_oce        ! ocean space and time domain 
    USE phycst         ! physical constants
-   USE sbcisf         ! ice shelf
    !
    USE in_out_manager ! I/O manager
@@ -77,5 +77,5 @@
       ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
       ! -----------------------------------------------------------------------
-      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
+      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0
 
       ! Compute bdy surface each cycle if non linear free surface
