Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydta.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydta.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydta.F90	(revision 10957)
@@ -65,5 +65,5 @@
 CONTAINS
 
-      SUBROUTINE bdy_dta( kt, jit, time_offset )
+      SUBROUTINE bdy_dta( kt, Kmm, time_offset )
       !!----------------------------------------------------------------------
       !!                   ***  SUBROUTINE bdy_dta  ***
@@ -75,5 +75,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in)           ::   kt           ! ocean time-step index 
-      INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option)
+      INTEGER, INTENT(in)           ::   Kmm          ! ocean time level index
       INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit
       !                                               ! is present then units = subcycle timesteps.
@@ -94,5 +94,5 @@
       ! Initialise data arrays once for all from initial conditions where required
       !---------------------------------------------------------------------------
-      IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN
+      IF( kt == nit000 ) THEN
 
          ! Calculate depth-mean currents
@@ -112,5 +112,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
@@ -120,5 +120,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 
                ENDIF
@@ -128,5 +128,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
@@ -141,5 +141,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 
@@ -151,5 +151,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 
@@ -165,5 +165,5 @@
                         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)%tem(ib,ik) = ts(ii,ij,ik,jp_tem,Kmm) * tmask(ii,ij,ik)         
                      END DO
                   END DO 
@@ -175,5 +175,5 @@
                         ii = idx_bdy(jbdy)%nbi(ib,igrd)
                         ij = idx_bdy(jbdy)%nbj(ib,igrd)
-                        dta_bdy(jbdy)%sal(ib,ik) = tsn(ii,ij,ik,jp_sal) * 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 
@@ -227,144 +227,84 @@
          dta => dta_bdy(jbdy)
          IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required
-      
-            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
+            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
             ELSE
-               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 )
-                  !
+               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)
-                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) )
+                     DO ik = 1, jpkm1
+                        dta%u2d(ib) = dta%u2d(ib) &
+                             &                       + e3u(ii,ij,ik,Kmm) * 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
-                  !
                   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)
-                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) )
+                     DO ik = 1, jpkm1
+                        dta%v2d(ib) = dta%v2d(ib) &
+                             &                       + e3v(ii,ij,ik,Kmm) * 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
-               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
-                     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
-                  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
+
+            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
             jstart = jstart + dta%nread(1)
          ENDIF    ! nn_dta(jbdy) = 1
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn.F90	(revision 10957)
@@ -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_a(:,:)
+      zva2d(:,:) = zva2d(:,:) * r1_hv_a(:,:)
 
       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_a(:,:), r1_hv_a(:,:), 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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn3d.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn3d.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn3d.F90	(revision 10957)
@@ -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   ! loop index
@@ -49,11 +51,11 @@
          SELECT CASE( cn_dyn3d(ib_bdy) )
          CASE('none')        ;   CYCLE
-         CASE('frs' )        ;   CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE('specified')   ;   CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE('zero')        ;   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, ll_npo=.false. )
-         CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
-         CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE('neumann')     ;   CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy )
+         CASE('frs' )        ;   CALL bdy_dyn3d_frs(           puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+         CASE('specified')   ;   CALL bdy_dyn3d_spe(           puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+         CASE('zero')        ;   CALL bdy_dyn3d_zro(           puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+         CASE('orlanski' )   ;   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. )
+         CASE('orlanski_npo');   CALL bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. )
+         CASE('zerograd')    ;   CALL bdy_dyn3d_zgrad(         puu, pvv, Kaa, idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy )
+         CASE('neumann')     ;   CALL bdy_dyn3d_nmn(           puu, pvv, Kaa, idx_bdy(ib_bdy), ib_bdy )
          CASE DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
          END SELECT
@@ -63,5 +65,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_spe( idx, dta, kt , ib_bdy )
+   SUBROUTINE bdy_dyn3d_spe( puu, pvv, Kaa, idx, dta, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_spe  ***
@@ -71,8 +73,9 @@
       !!
       !!----------------------------------------------------------------------
-      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    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   jb, jk         ! dummy loop indices
@@ -86,5 +89,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
@@ -95,16 +98,14 @@
             ii   = idx%nbi(jb,igrd)
             ij   = idx%nbj(jb,igrd)
-            va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
-         END DO
-      END DO
-      CALL lbc_bdy_lnk( 'bdydyn3d', ua, 'U', -1., ib_bdy )   ! Boundary points should be updated  
-      CALL lbc_bdy_lnk( 'bdydyn3d', va, 'V', -1., ib_bdy )   
-      !
-      IF( kt == nit000 )   CLOSE( unit = 102 )
+            pvv(ii,ij,jk,Kaa) = dta%v3d(jb,jk) * vmask(ii,ij,jk)
+         END DO
+      END DO
+      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )   ! Boundary points should be updated  
+      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
       !
    END SUBROUTINE bdy_dyn3d_spe
 
 
-   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy )
+   SUBROUTINE bdy_dyn3d_zgrad( puu, pvv, Kaa, idx, dta, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
@@ -113,8 +114,9 @@
       !!
       !!----------------------------------------------------------------------
-      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
+      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
@@ -130,6 +132,6 @@
             ij   = idx%nbj(jb,igrd)
             fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 )
-            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) &
-                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )
+            puu(ii,ij,jk,Kaa) = puu(ii,ij,jk,Kaa) * REAL( 1 - fu) + ( puu(ii,ij+fu,jk,Kaa) * umask(ii,ij+fu,jk) &
+                        &+ puu(ii,ij-fu,jk,Kaa) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )
          END DO
       END DO
@@ -141,17 +143,15 @@
             ij   = idx%nbj(jb,igrd)
             fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 )
-            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) &
-                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )
-         END DO
-      END DO
-      CALL lbc_bdy_lnk( 'bdydyn3d', ua, 'U', -1., ib_bdy )   ! Boundary points should be updated  
-      CALL lbc_bdy_lnk( 'bdydyn3d', va, 'V', -1., ib_bdy )   
-      !
-      IF( kt == nit000 )   CLOSE( unit = 102 )
+            pvv(ii,ij,jk,Kaa) = pvv(ii,ij,jk,Kaa) * REAL( 1 - fv ) + ( pvv(ii+fv,ij,jk,Kaa) * vmask(ii+fv,ij,jk) &
+                        &+ pvv(ii-fv,ij,jk,Kaa) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )
+         END DO
+      END DO
+      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )   ! Boundary points should be updated  
+      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
       !
    END SUBROUTINE bdy_dyn3d_zgrad
 
 
-   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
+   SUBROUTINE bdy_dyn3d_zro( puu, pvv, Kaa, idx, dta, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_zro  ***
@@ -160,8 +160,9 @@
       !!
       !!----------------------------------------------------------------------
-      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    ) ::   ib_bdy    ! BDY set index
       !
       INTEGER  ::   ib, ik         ! dummy loop indices
@@ -175,5 +176,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
@@ -184,16 +185,14 @@
          ij = idx%nbj(ib,igrd)
          DO ik = 1, jpkm1
-            va(ii,ij,ik) = 0._wp
-         END DO
-      END DO
-      !
-      CALL lbc_bdy_lnk( 'bdydyn3d', ua, 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( 'bdydyn3d', va, 'V', -1.,ib_bdy )   ! Boundary points should be updated
-      !
-      IF( kt == nit000 )   CLOSE( unit = 102 )
+            pvv(ii,ij,ik,Kaa) = 0._wp
+         END DO
+      END DO
+      !
+      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )   ;   CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1.,ib_bdy )   ! Boundary points should be updated
       !
    END SUBROUTINE bdy_dyn3d_zro
 
 
-   SUBROUTINE bdy_dyn3d_frs( idx, dta, kt, ib_bdy )
+   SUBROUTINE bdy_dyn3d_frs( puu, pvv, Kaa, idx, dta, ib_bdy )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_frs  ***
@@ -206,8 +205,9 @@
       !!               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    ) ::   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
@@ -222,5 +222,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
@@ -232,16 +232,14 @@
             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 
-      CALL lbc_bdy_lnk( 'bdydyn3d', ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
-      CALL lbc_bdy_lnk( 'bdydyn3d', va, 'V', -1., ib_bdy )   
-      !
-      IF( kt == nit000 )   CLOSE( unit = 102 )
+      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )    ! Boundary points should be updated
+      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
       !
    END SUBROUTINE bdy_dyn3d_frs
 
 
-   SUBROUTINE bdy_dyn3d_orlanski( idx, dta, ib_bdy, ll_npo )
+   SUBROUTINE bdy_dyn3d_orlanski( Kbb, puu, pvv, Kaa, idx, dta, ib_bdy, ll_npo )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn3d_orlanski  ***
@@ -253,29 +251,31 @@
       !! 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) ::   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    ) ::   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 )
+      CALL bdy_orlanski_3d( idx, igrd, puu(:,:,:,Kbb), puu(:,:,:,Kaa), dta%u3d, ll_npo )
 
       igrd = 3      ! Orlanski bc on v-velocity
       !  
-      CALL bdy_orlanski_3d( idx, igrd, vb, va, dta%v3d, ll_npo )
-      !
-      CALL lbc_bdy_lnk( 'bdydyn3d', ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
-      CALL lbc_bdy_lnk( 'bdydyn3d', va, 'V', -1., ib_bdy )   
+      CALL bdy_orlanski_3d( idx, igrd, pvv(:,:,:,Kbb), pvv(:,:,:,Kaa), dta%v3d, ll_npo )
+      !
+      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )    ! Boundary points should be updated
+      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )   
       !
    END SUBROUTINE bdy_dyn3d_orlanski
 
 
-   SUBROUTINE bdy_dyn3d_dmp( kt )
+   SUBROUTINE bdy_dyn3d_dmp( kt, Kbb, puu, pvv, Krhs )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_dyn3d_dmp  ***
@@ -284,5 +284,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
@@ -302,6 +304,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
@@ -313,6 +315,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
@@ -320,5 +322,5 @@
       END DO
       !
-      CALL lbc_lnk_multi( 'bdydyn3d', ua, 'U', -1.,  va, 'V', -1. )   ! Boundary points should be updated
+      CALL lbc_lnk_multi( 'bdydyn3d', puu(:,:,:,Krhs), 'U', -1.,  pvv(:,:,:,Krhs), 'V', -1. )   ! Boundary points should be updated
       !
       IF( ln_timing )   CALL timing_stop('bdy_dyn3d_dmp')
@@ -327,5 +329,5 @@
 
 
-   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy )
+   SUBROUTINE bdy_dyn3d_nmn( puu, pvv, Kaa, idx, ib_bdy )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
@@ -336,22 +338,24 @@
       !!
       !!----------------------------------------------------------------------
-      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
-      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
+      INTEGER                             , INTENT( in    ) ::   ib_bdy    ! BDY set index
 
       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      ! Neumann bc on u-velocity; 
       !            
-      CALL bdy_nmn( idx, igrd, ua )
+      CALL bdy_nmn( idx, igrd, puu(:,:,:,Kaa) )
 
       igrd = 3      ! Neumann bc on v-velocity
       !  
-      CALL bdy_nmn( idx, igrd, va )
-      !
-      CALL lbc_bdy_lnk( 'bdydyn3d', ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
-      CALL lbc_bdy_lnk( 'bdydyn3d', va, 'V', -1., ib_bdy )
+      CALL bdy_nmn( idx, igrd, pvv(:,:,:,Kaa) )
+      !
+      CALL lbc_bdy_lnk( 'bdydyn3d', puu(:,:,:,Kaa), 'U', -1., ib_bdy )    ! Boundary points should be updated
+      CALL lbc_bdy_lnk( 'bdydyn3d', pvv(:,:,:,Kaa), 'V', -1., ib_bdy )
       !
    END SUBROUTINE bdy_dyn3d_nmn
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdylib.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdylib.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdylib.F90	(revision 10957)
@@ -34,5 +34,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_frs( idx, pta, dta )
+   SUBROUTINE bdy_frs( idx, phia, dta )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_frs  ***
@@ -44,5 +44,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
@@ -57,5 +57,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
@@ -64,5 +64,5 @@
 
 
-   SUBROUTINE bdy_spe( idx, pta, dta )
+   SUBROUTINE bdy_spe( idx, phia, dta )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_spe  ***
@@ -73,5 +73,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
@@ -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, ll_npo )
+   SUBROUTINE bdy_orl( idx, phib, phia, dta, 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,                             INTENT(in) ::   ll_npo  ! switch for NPO version
       !!
@@ -111,5 +111,5 @@
       igrd = 1                       ! Everything is at T-points here
       !
-      CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, ll_npo )
+      CALL bdy_orlanski_3d( idx, igrd, phib(:,:,:), phia(:,:,:), dta, ll_npo )
       !
    END SUBROUTINE bdy_orl
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdytra.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdytra.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdytra.F90	(revision 10957)
@@ -40,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE bdy_tra( kt )
+   SUBROUTINE bdy_tra( kt, Kbb, pts, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  SUBROUTINE bdy_tra  ***
@@ -47,7 +47,9 @@
       !!
       !!----------------------------------------------------------------------
-      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   ! Loop indeces
+      INTEGER                        :: ib_bdy, jn, igrd   ! Loop indices
       TYPE(ztrabdy), DIMENSION(jpts) :: zdta               ! Temporary data structure
       !!----------------------------------------------------------------------
@@ -63,14 +65,14 @@
             SELECT CASE( TRIM(cn_tra(ib_bdy)) )
             CASE('none'        )   ;   CYCLE
-            CASE('frs'         )   ;   CALL bdy_frs ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra )
-            CASE('specified'   )   ;   CALL bdy_spe ( idx_bdy(ib_bdy),                tsa(:,:,:,jn), zdta(jn)%tra )
-            CASE('neumann'     )   ;   CALL bdy_nmn ( idx_bdy(ib_bdy), igrd         , tsa(:,:,:,jn)               )
-            CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), zdta(jn)%tra, ll_npo=.false. )
-            CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), tsb(:,:,:,jn), tsa(:,:,:,jn), zdta(jn)%tra, ll_npo=.true. )
-            CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                tsa(:,:,:,jn),               jn )
+            CASE('frs'         )   ;   CALL bdy_frs ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa), zdta(jn)%tra )
+            CASE('specified'   )   ;   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)               )
+            CASE('orlanski'    )   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), zdta(jn)%tra, ll_npo=.false. )
+            CASE('orlanski_npo')   ;   CALL bdy_orl ( idx_bdy(ib_bdy), pts(:,:,:,jn,Kbb), pts(:,:,:,jn,Kaa), zdta(jn)%tra, ll_npo=.true. )
+            CASE('runoff'      )   ;   CALL bdy_rnf ( idx_bdy(ib_bdy),                pts(:,:,:,jn,Kaa),               jn )
             CASE DEFAULT           ;   CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' )
             END SELECT
             ! Boundary points should be updated
-            CALL lbc_bdy_lnk( 'bdytra', tsa(:,:,:,jn), 'T', 1., ib_bdy )
+            CALL lbc_bdy_lnk( 'bdytra', pts(:,:,:,jn,Kaa), 'T', 1., ib_bdy )
             ! 
          END DO
@@ -80,5 +82,5 @@
 
 
-   SUBROUTINE bdy_rnf( idx, pta, jpa )
+   SUBROUTINE bdy_rnf( idx, pt, jpa )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_rnf  ***
@@ -90,5 +92,5 @@
       !!----------------------------------------------------------------------
       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
       !
@@ -105,6 +107,6 @@
             ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
             jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
-            if (jpa == jp_tem) pta(ii,ij,ik) = pta(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
-            if (jpa == jp_sal) pta(ii,ij,ik) =                 0.1 * tmask(ii,ij,ik)
+            if (jpa == jp_tem) pt(ii,ij,ik) = pt(ii+ip,ij+jp,ik) * tmask(ii,ij,ik)
+            if (jpa == jp_sal) pt(ii,ij,ik) =                 0.1 * tmask(ii,ij,ik)
          END DO
       END DO
@@ -113,5 +115,5 @@
 
 
-   SUBROUTINE bdy_tra_dmp( kt )
+   SUBROUTINE bdy_tra_dmp( kt, Kbb, pts, Krhs )
       !!----------------------------------------------------------------------
       !!                 ***  SUBROUTINE bdy_tra_dmp  ***
@@ -120,5 +122,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
@@ -139,8 +143,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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynnxt.F90	(revision 10957)
@@ -50,5 +50,5 @@
 #if defined key_agrif
    USE agrif_oce_interp
-#endif
+#endif   
 
    IMPLICIT NONE
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE dyn_nxt ( kt, Kmm )
+   SUBROUTINE dyn_nxt ( kt, Kbb, Kmm, Kaa )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE dyn_nxt  ***
@@ -92,6 +92,6 @@
       !!               un,vn   now horizontal velocity of next time-step
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-      INTEGER, INTENT( in ) ::   Kmm     ! time level index
+      INTEGER, INTENT( in ) ::   kt             ! ocean time-step index
+      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa  ! time level indices
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
@@ -148,6 +148,10 @@
       !
       !                                !* BDY open boundaries
-      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
-      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )
+      !! IMMERSE development : Following the general pattern for the new code we want to pass in the 
+      !!                       velocities to bdy_dyn as arguments so here we use "uu" and "vv" even 
+      !!                       though we haven't converted the velocity names in the rest of dynnxt.F90
+      !!                       because it will be completely rewritten. DS.
+      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa )
+      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, uu, vv, Kaa, dyn3d_only=.true. )
 
 !!$   Do we need a call to bdy_vol here??
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranxt.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranxt.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/tranxt.F90	(revision 10957)
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE tra_nxt( kt, Kbb, Kmm, Krhs )
+   SUBROUTINE tra_nxt( kt, Kbb, Kmm, Krhs, Kaa )
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE tranxt  ***
@@ -86,6 +86,6 @@
       !! ** Action  : - ts(Kbb) & ts(Kmm) ready for the next time step
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
-      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! time level indices
+      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
+      INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs, Kaa  ! time level indices
       !!
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
@@ -110,5 +110,9 @@
       CALL lbc_lnk_multi( 'tranxt', ts(:,:,:,jp_tem,Krhs), 'T', 1., ts(:,:,:,jp_sal,Krhs), 'T', 1. )
       !
-      IF( ln_bdy )   CALL bdy_tra( kt )  ! BDY open boundaries
+      !! IMMERSE development : Following the general pattern for the new code we want to pass in the 
+      !!                       velocities to bdy_dyn as arguments so here we use "ts" even 
+      !!                       though we haven't converted the tracer names in the rest of tranxt.F90
+      !!                       because it will be completely rewritten. DS.
+      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, ts, Kaa )  ! BDY open boundaries
  
       ! set time step size (Euler/Leapfrog)
Index: /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
===================================================================
--- /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10956)
+++ /NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90	(revision 10957)
@@ -118,6 +118,6 @@
       IF( ln_tide    )   CALL sbc_tide( kstp )                   ! update tide potential
       IF( ln_apr_dyn )   CALL sbc_apr ( kstp )                   ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib) 
-      IF( ln_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-                         CALL sbc     ( kstp, Nbb, Nnn )         ! Sea Boundary Condition (including sea-ice)
+      IF( ln_bdy     )   CALL bdy_dta ( kstp,      Nnn, time_offset=+1 )  ! update dynamic & tracer data at open boundaries
+                         CALL sbc     ( kstp, Nbb, Nnn )                  ! Sea Boundary Condition (including sea-ice)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -188,5 +188,5 @@
       IF(  lk_asminc .AND. ln_asmiau .AND. ln_dyninc )   &
                &         CALL dyn_asm_inc   ( kstp )  ! apply dynamics assimilation increment
-      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp )  ! bdy damping trends
+      IF( ln_bdy     )   CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs )  ! bdy damping trends
 #if defined key_agrif
       IF(.NOT. Agrif_Root())  & 
@@ -247,5 +247,5 @@
       IF( ln_trabbl  )   CALL tra_bbl       ( kstp, Nbb, Nnn, Nrhs )  ! advective (and/or diffusive) bottom boundary layer scheme
       IF( ln_tradmp  )   CALL tra_dmp       ( kstp, Nbb, Nnn, Nrhs )  ! internal damping trends
-      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp )  ! bdy damping trends
+      IF( ln_bdy     )   CALL bdy_tra_dmp   ( kstp, Nbb, ts, Nrhs )  ! bdy damping trends
 #if defined key_agrif
       IF(.NOT. Agrif_Root())  & 
@@ -281,6 +281,6 @@
 !! 
 !!jc2: dynnxt must be the latest call. e3t_b are indeed updated in that routine
-                         CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs )  ! finalize (bcs) tracer fields at next time step and swap
-                         CALL dyn_nxt       ( kstp,      Nnn )        ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
+                         CALL tra_nxt       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! finalize (bcs) tracer fields at next time step and swap
+                         CALL dyn_nxt       ( kstp, Nbb, Nnn, Naa  )  ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt)
                          CALL ssh_swp       ( kstp )  ! swap of sea surface height
       IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors
