Index: branches/2016/dev_NOC_2016/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 7342)
@@ -651,4 +651,5 @@
     ln_vol        = .false.               !  total volume correction (see nn_volctl parameter)
     nn_volctl     = 1                     !  = 0, the total water flux across open boundaries is zero
+    nb_jpk_bdy    = -1                    ! number of levels in the bdy data (set < 0 if consistent with planned run)
 /
 !-----------------------------------------------------------------------
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90	(revision 7342)
@@ -49,4 +49,5 @@
       LOGICAL                           ::  ll_tem
       LOGICAL                           ::  ll_sal
+      LOGICAL                           ::  ll_fvl
       REAL(wp), POINTER, DIMENSION(:)   ::  ssh
       REAL(wp), POINTER, DIMENSION(:)   ::  u2d
@@ -91,4 +92,5 @@
    !
    INTEGER                    ::   nb_bdy                   !: number of open boundary sets
+   INTEGER                    ::   nb_jpk_bdy               !: number of levels in the bdy data (set < 0 if consistent with planned run)
    INTEGER, DIMENSION(jp_bdy) ::   nn_rimwidth              !: boundary rim width for Flow Relaxation Scheme
    INTEGER                    ::   nn_volctl                !: = 0 the total volume will have the variability of the surface Flux E-P 
@@ -134,5 +136,9 @@
                                                                           !: =1 => some data to be read in from data files
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global        !: workspace for reading in global data arrays (unstr.  bdy)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_z      !: workspace for reading in global depth arrays (unstr.  bdy)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global_dz     !: workspace for reading in global depth arrays (unstr.  bdy)
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2       !: workspace for reading in global data arrays (struct. bdy)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_z     !: workspace for reading in global depth arrays (struct. bdy)
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), TARGET ::   dta_global2_dz    !: workspace for reading in global depth arrays (struct. bdy)
 !$AGRIF_DO_NOT_TREAT
    TYPE(OBC_INDEX), DIMENSION(jp_bdy), TARGET      ::   idx_bdy           !: bdy indices (local process)
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 7342)
@@ -267,6 +267,11 @@
 
                         jend = jstart + dta%nread(2) - 1
-                        CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  &
-                                     & kit=jit, kt_offset=time_offset )
+                        IF( ln_full_vel_array(ib_bdy) ) 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(ib_bdy)  )
+                        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
@@ -333,5 +338,5 @@
                      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 )
+                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, fvl=ln_full_vel_array(ib_bdy) )
                   ENDIF
                   ! If full velocities in boundary data then split into barotropic and baroclinic data
@@ -459,5 +464,5 @@
       NAMELIST/nambdy_dta/ bn_a_i, bn_ht_i, bn_ht_s
 #endif
-      NAMELIST/nambdy_dta/ ln_full_vel
+      NAMELIST/nambdy_dta/ ln_full_vel, nb_jpk_bdy
       !!---------------------------------------------------------------------------
       !
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90	(revision 7342)
@@ -57,4 +57,6 @@
          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 DEFAULT        ;   CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' )
          END SELECT
@@ -110,4 +112,53 @@
    END SUBROUTINE bdy_dyn3d_spe
 
+   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***
+      !!
+      !! ** Purpose : - Enforce a zero gradient of normal velocity
+      !!
+      !!----------------------------------------------------------------------
+      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  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
+      INTEGER  ::   fu, fv
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad')
+      !
+      igrd = 2                      ! Copying tangential velocity into bdy points
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            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 )
+         END DO
+      END DO
+      !
+      igrd = 3                      ! Copying tangential velocity into bdy points
+      DO jb = 1, idx%nblenrim(igrd)
+         DO jk = 1, jpkm1
+            ii   = idx%nbi(jb,igrd)
+            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( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated  
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )   
+      !
+      IF( kt .eq. nit000 ) CLOSE( unit = 102 )
+
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad')
+
+   END SUBROUTINE bdy_dyn3d_zgrad
 
    SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy )
@@ -296,4 +347,38 @@
    END SUBROUTINE bdy_dyn3d_dmp
 
+   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***
+      !!             
+      !!              - Apply Neumann condition to baroclinic velocities. 
+      !!              - Wrapper routine for bdy_nmn
+      !! 
+      !!
+      !!----------------------------------------------------------------------
+      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices
+      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index
+
+      INTEGER  ::   jb, igrd                               ! dummy loop indices
+      !!----------------------------------------------------------------------
+
+      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn')
+      !
+      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
+      !
+      igrd = 2      ! Neumann bc on u-velocity; 
+      !            
+      CALL bdy_nmn( idx, igrd, ua )
+
+      igrd = 3      ! Neumann bc on v-velocity
+      !  
+      CALL bdy_nmn( idx, igrd, va )
+      !
+      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated
+      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn')
+      !
+   END SUBROUTINE bdy_dyn3d_nmn
+
 #else
    !!----------------------------------------------------------------------
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 7342)
@@ -100,5 +100,5 @@
          &             cn_ice_lim, nn_ice_lim_dta,                           &
          &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
-         &             ln_vol, nn_volctl, nn_rimwidth
+         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
          !
       NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
@@ -209,4 +209,12 @@
              dta_bdy(ib_bdy)%ll_u3d = .true.
              dta_bdy(ib_bdy)%ll_v3d = .true.
+          CASE('neumann')
+             IF(lwp) WRITE(numout,*) '      Neumann conditions'
+             dta_bdy(ib_bdy)%ll_u3d = .false.
+             dta_bdy(ib_bdy)%ll_v3d = .false.
+          CASE('zerograd')
+             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities'
+             dta_bdy(ib_bdy)%ll_u3d = .false.
+             dta_bdy(ib_bdy)%ll_v3d = .false.
           CASE('zero')
              IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
@@ -377,4 +385,9 @@
           IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
           IF(lwp) WRITE(numout,*)
+        ENDIF
+        IF( nb_jpk_bdy > 0 ) THEN
+           IF(lwp) WRITE(numout,*) '*** open boundary will be interpolate in the vertical onto the native grid ***'
+        ELSE
+           IF(lwp) WRITE(numout,*) '*** open boundary will be read straight onto the native grid without vertical interpolation ***'
         ENDIF
      ENDIF
@@ -499,6 +512,25 @@
             &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
 
-         ALLOCATE( dta_global(jpbdtau, 1, jpk) )
-         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
+         IF( nb_jpk_bdy>0 ) THEN
+            ALLOCATE( dta_global(jpbdtau, 1, nb_jpk_bdy) )
+            ALLOCATE( dta_global_z(jpbdtau, 1, nb_jpk_bdy) )
+            ALLOCATE( dta_global_dz(jpbdtau, 1, nb_jpk_bdy) )
+         ELSE
+            ALLOCATE( dta_global(jpbdtau, 1, jpk) )
+            ALLOCATE( dta_global_z(jpbdtau, 1, jpk) ) ! needed ?? TODO
+            ALLOCATE( dta_global_dz(jpbdtau, 1, jpk) )! needed ?? TODO
+         ENDIF
+
+         IF ( icount>0 ) THEN
+            IF( nb_jpk_bdy>0 ) THEN
+               ALLOCATE( dta_global2(jpbdtas, nrimmax, nb_jpk_bdy) )
+               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, nb_jpk_bdy) )
+               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, nb_jpk_bdy) )
+            ELSE
+               ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
+               ALLOCATE( dta_global2_z(jpbdtas, nrimmax, jpk) ) ! needed ?? TODO
+               ALLOCATE( dta_global2_dz(jpbdtas, nrimmax, jpk) )! needed ?? TODO  
+            ENDIF
+         ENDIF
          ! 
       ENDIF
@@ -1092,4 +1124,6 @@
       !          = 0  elsewhere   
  
+      bdytmask(:,:) = ssmask(:,:)
+
       IF( ln_mask_file ) THEN
          CALL iom_open( cn_mask_file, inum )
@@ -1108,27 +1142,6 @@
          CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
 
-
-         ! Mask corrections
-         ! ----------------
-         DO ik = 1, jpkm1
-            DO ij = 1, jpj
-               DO ii = 1, jpi
-                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
-                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
-                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
-               END DO      
-            END DO
-            DO ij = 2, jpjm1
-               DO ii = 2, jpim1
-                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
-                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
-               END DO      
-            END DO
-         END DO
-         tmask_i (:,:) = ssmask(:,:) * tmask_i(:,:)
-         !
       ENDIF ! ln_mask_file=.TRUE.
       
-      bdytmask(:,:) = ssmask(:,:)
       IF( .NOT.ln_mask_file ) THEN
          ! If .not. ln_mask_file then we need to derive mask on U and V grid from mask on T grid here.
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90	(revision 7342)
@@ -27,4 +27,5 @@
    PUBLIC   bdy_orlanski_2d     ! routine called where?
    PUBLIC   bdy_orlanski_3d     ! routine called where?
+   PUBLIC   bdy_nmn             ! routine called where?
 
    !!----------------------------------------------------------------------
@@ -355,4 +356,76 @@
    END SUBROUTINE bdy_orlanski_3d
 
+   SUBROUTINE bdy_nmn( idx, igrd, phia )
+      !!----------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_nmn  ***
+      !!                    
+      !! ** Purpose : Duplicate the value at open boundaries, zero gradient.
+      !! 
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in)     ::   igrd     ! grid index
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)
+      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices
+      !! 
+      REAL(wp) ::   zcoef, zcoef1, zcoef2
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field
+      INTEGER  ::   ib, ik   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn')
+      !
+      SELECT CASE(igrd)
+         CASE(1)
+            pmask => tmask(:,:,:)
+            bdypmask => bdytmask(:,:)
+         CASE(2)
+            pmask => umask(:,:,:)
+            bdypmask => bdyumask(:,:)
+         CASE(3)
+            pmask => vmask(:,:,:)
+            bdypmask => bdyvmask(:,:)
+         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' )
+      END SELECT
+      DO ib = 1, idx%nblenrim(igrd)
+         ii = idx%nbi(ib,igrd)
+         ij = idx%nbj(ib,igrd)
+         DO ik = 1, jpkm1
+            ! search the sense of the gradient
+            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik)
+            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik)
+            IF ( nint(zcoef1+zcoef2) == 0) THEN
+               ! corner **** we probably only want to set the tangentail component for the dynamics here
+               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik)
+               IF (zcoef > .5_wp) THEN ! Only set none isolated points.
+                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + &
+                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + &
+                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + &
+                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)
+                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik)
+               ELSE
+                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik)
+               ENDIF
+            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN
+               ! oblique corner **** we probably only want to set the normal component for the dynamics here
+               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + &
+                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  )
+               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + &
+                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + &
+                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + &
+                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  )
+ 
+               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik)
+            ELSE
+               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik))
+               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik))
+               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik)
+            ENDIF
+         END DO
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn')
+      !
+   END SUBROUTINE bdy_nmn
 
 #else
@@ -367,4 +440,7 @@
       WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_orlanski_3d
+   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine
+      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_nmn
 #endif
 
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 7342)
@@ -155,6 +155,7 @@
       !
       REAL(wp) ::   zwgt           ! boundary weight
-      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
-      INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses
+      REAL(wp) ::   zcoef, zcoef1,zcoef2
+      INTEGER  ::   ib, ik, igrd   ! dummy loop indices
+      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses
       !!----------------------------------------------------------------------
       !
@@ -169,5 +170,5 @@
             zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  )
             zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1)
-            IF ( zcoef1+zcoef2 == 0) THEN
+            IF ( NINT(zcoef1+zcoef2) == 0) THEN
                ! corner
                zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik)
@@ -176,13 +177,13 @@
                  &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + &
                  &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik)
-               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)
                tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + &
                  &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + &
                  &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + &
                  &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik)
-               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik)
+               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)
             ELSE
-               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )
-               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)
+               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ))
+               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1))
                tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik)
                tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik)
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90	(revision 7342)
@@ -24,8 +24,9 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   !
+   USE bdy_oce      
    USE in_out_manager  ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE lib_mpp         !
+   USE iom
    USE wrk_nemo        ! Memory allocation
    USE timing          ! Timing
@@ -102,5 +103,5 @@
       INTEGER  ::   iif, iil, ii0, ii1, ii   ! local integers
       INTEGER  ::   ijf, ijl, ij0, ij1       !   -       -
-      INTEGER  ::   ios
+      INTEGER  ::   ios, inum
       INTEGER  ::   isrow                    ! index for ORCA1 starting row
       INTEGER , POINTER, DIMENSION(:,:) ::  imsk
@@ -108,4 +109,13 @@
       !!
       NAMELIST/namlbc/ rn_shlat, ln_vorlat
+#if defined key_bdy 
+      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 &
+         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
+         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
+         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
+         &             cn_ice_lim, nn_ice_lim_dta,                           &
+         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
+         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
+#endif
       !!---------------------------------------------------------------------
       !
@@ -155,4 +165,30 @@
       END DO  
       
+#if defined key_bdy 
+      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist :Unstructured open boundaries  
+      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist :Unstructured open boundaries
+      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nambdy )
+
+     IF( ln_mask_file ) THEN ! correct for bdy mask
+         CALL iom_open( cn_mask_file, inum )
+         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
+         CALL iom_close( inum )
+
+         ! Mask corrections
+         ! ----------------
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  tmask(ji,jj,jk) = tmask(ji,jj,jk) * bdytmask(ji,jj)
+               END DO
+            END DO
+         END DO
+      ENDIF
+#endif
       ! (ISF) define barotropic mask and mask the ice shelf point
       ssmask(:,:)=tmask(:,:,1) ! at this stage ice shelf is not masked
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90	(revision 7342)
@@ -24,4 +24,7 @@
    USE wrk_nemo        ! Memory Allocation
    USE timing          ! Timing
+#if defined key_bdy 
+   USE bdy_oce         ! ocean open boundary conditions
+#endif
 
    IMPLICIT NONE
@@ -78,4 +81,9 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke
       REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 
+#if defined key_bdy
+      INTEGER  ::   jb                 ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers
+      INTEGER  ::   fu, fv
+#endif
       !!----------------------------------------------------------------------
       !
@@ -98,4 +106,32 @@
       zhke(:,:,jpk) = 0._wp
       
+#if defined key_bdy
+      ! Maria Luneva & Fred Wobus: July-2016
+      ! compensate for lack of turbulent kinetic energy on liquid bdy points
+      DO ib_bdy = 1, nb_bdy
+         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
+            igrd = 2           ! Copying normal velocity into points outside bdy
+            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               DO jk = 1, jpkm1
+                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )
+                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)
+               END DO
+            END DO
+            !
+            igrd = 3           ! Copying normal velocity into points outside bdy
+            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               DO jk = 1, jpkm1
+                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)
+                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)
+                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )
+                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)
+               END DO
+            END DO
+         ENDIF
+      ENDDO  
+#endif 
+
       SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==!
       !
@@ -133,4 +169,12 @@
          !
       END SELECT
+
+#if defined key_bdy
+      ! restore velocity masks at points outside boundary
+      un(:,:,:) = un(:,:,:) * umask(:,:,:)
+      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)
+#endif      
+
+
       !
       DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==!
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 7342)
@@ -9,4 +9,5 @@
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie and G. Reffray)  add C1D case  
    !!            3.6  ! 2014-15  DIMG format removed
+   !!            3.6  ! 2015-15  (J. Harle) Added procedure to read REAL attributes
    !!--------------------------------------------------------------------
 
@@ -67,5 +68,5 @@
    END INTERFACE
    INTERFACE iom_getatt
-      MODULE PROCEDURE iom_g0d_intatt
+      MODULE PROCEDURE iom_g0d_intatt, iom_g0d_ratt
    END INTERFACE
    INTERFACE iom_rstput
@@ -979,5 +980,5 @@
          IF( iom_file(kiomid)%nfid > 0 ) THEN
             SELECT CASE (iom_file(kiomid)%iolib)
-            CASE (jpnf90   )   ;   CALL iom_nf90_getatt( kiomid, cdatt, pvar )
+            CASE (jpnf90   )   ;   CALL iom_nf90_getatt( kiomid, cdatt, pv_i0d=pvar )
             CASE DEFAULT
                CALL ctl_stop( 'iom_g0d_att: accepted IO library is only jpnf90' )
@@ -987,4 +988,24 @@
    END SUBROUTINE iom_g0d_intatt
 
+   SUBROUTINE iom_g0d_ratt( kiomid, cdatt, pvar, cdvar )
+      INTEGER         , INTENT(in   )                 ::   kiomid    ! Identifier of the file
+      CHARACTER(len=*), INTENT(in   )                 ::   cdatt     ! Name of the attribute
+      REAL(wp)        , INTENT(  out)                 ::   pvar      ! written field
+      CHARACTER(len=*), INTENT(in   ), OPTIONAL       ::   cdvar     ! Name of the variable
+      !
+      IF( kiomid > 0 ) THEN
+         IF( iom_file(kiomid)%nfid > 0 ) THEN
+            SELECT CASE (iom_file(kiomid)%iolib)
+            CASE (jpnf90   )   ;   IF( PRESENT(cdvar) ) THEN
+                                      CALL iom_nf90_getatt( kiomid, cdatt, pv_r0d=pvar, cdvar=cdvar )
+                                   ELSE
+                                      CALL iom_nf90_getatt( kiomid, cdatt, pv_r0d=pvar )
+                                   ENDIF
+            CASE DEFAULT    
+               CALL ctl_stop( 'iom_g0d_att: accepted IO library is only jpnf90' )
+            END SELECT
+         ENDIF
+      ENDIF
+   END SUBROUTINE iom_g0d_ratt
 
    !!----------------------------------------------------------------------
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/IOM/iom_nf90.F90	(revision 7342)
@@ -7,4 +7,5 @@
    !!            9.0  ! 06 02  (S. Masson) Adaptation to NEMO
    !!             "   ! 07 07  (D. Storkey) Changes to iom_nf90_gettime
+   !!            3.6  ! 2015-15  (J. Harle) Added procedure to read REAL attributes
    !!--------------------------------------------------------------------
    !!gm  caution add !DIR nec: improved performance to be checked as well as no result changes
@@ -35,5 +36,5 @@
    END INTERFACE
    INTERFACE iom_nf90_getatt
-      MODULE PROCEDURE iom_nf90_intatt
+      MODULE PROCEDURE iom_nf90_att
    END INTERFACE
    INTERFACE iom_nf90_rstput
@@ -313,7 +314,7 @@
 
 
-   SUBROUTINE iom_nf90_intatt( kiomid, cdatt, pvar )
-      !!-----------------------------------------------------------------------
-      !!                  ***  ROUTINE  iom_nf90_intatt  ***
+   SUBROUTINE iom_nf90_att( kiomid, cdatt, pv_i0d, pv_r0d, cdvar)
+      !!-----------------------------------------------------------------------
+      !!                  ***  ROUTINE  iom_nf90_att  ***
       !!
       !! ** Purpose : read an integer attribute with NF90
@@ -321,7 +322,10 @@
       INTEGER         , INTENT(in   ) ::   kiomid   ! Identifier of the file
       CHARACTER(len=*), INTENT(in   ) ::   cdatt    ! attribute name
-      INTEGER         , INTENT(  out) ::   pvar     ! read field
+      INTEGER         , INTENT(  out), OPTIONAL       ::   pv_i0d    ! read field
+      REAL(wp),         INTENT(  out), OPTIONAL       ::   pv_r0d    ! read field 
+      CHARACTER(len=*), INTENT(in   ), OPTIONAL       ::   cdvar     ! name of the variable
       !
       INTEGER                         ::   if90id   ! temporary integer
+      INTEGER                         ::   ivarid           ! NetCDF  variable Id
       LOGICAL                         ::   llok     ! temporary logical
       CHARACTER(LEN=100)              ::   clinfo   ! info character
@@ -329,14 +333,40 @@
       ! 
       if90id = iom_file(kiomid)%nfid
-      llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr
+      IF( PRESENT(cdvar) ) THEN
+         llok = NF90_INQ_VARID( if90id, TRIM(cdvar), ivarid ) == nf90_noerr   ! does the variable exist in the file
+         IF( llok ) THEN
+            llok = NF90_Inquire_attribute(if90id, ivarid, cdatt) == nf90_noerr
+         ELSE
+            CALL ctl_warn('iom_nf90_getatt: no variable '//cdvar//' found')
+         ENDIF
+      ELSE
+         llok = NF90_Inquire_attribute(if90id, NF90_GLOBAL, cdatt) == nf90_noerr
+      ENDIF 
+!
       IF( llok) THEN
          clinfo = 'iom_nf90_getatt, file: '//TRIM(iom_file(kiomid)%name)//', att: '//TRIM(cdatt)
-         CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pvar), clinfo)
+         IF(     PRESENT(pv_r0d) ) THEN
+            IF( PRESENT(cdvar) ) THEN
+               CALL iom_nf90_check(NF90_GET_ATT(if90id, ivarid, cdatt, values=pv_r0d), clinfo)
+            ELSE
+               CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pv_r0d), clinfo)
+            ENDIF
+         ELSE
+            IF( PRESENT(cdvar) ) THEN
+               CALL iom_nf90_check(NF90_GET_ATT(if90id, ivarid, cdatt, values=pv_i0d), clinfo)
+            ELSE
+               CALL iom_nf90_check(NF90_GET_ATT(if90id, NF90_GLOBAL, cdatt, values=pv_i0d), clinfo)
+            ENDIF
+         ENDIF
       ELSE
          CALL ctl_warn('iom_nf90_getatt: no attribute '//cdatt//' found')
-         pvar = -999
+         IF(     PRESENT(pv_r0d) ) THEN
+            pv_r0d = -999._wp
+         ELSE
+            pv_i0d = -999
+         ENDIF
       ENDIF
       ! 
-   END SUBROUTINE iom_nf90_intatt
+   END SUBROUTINE iom_nf90_att
 
 
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/LBC/mppini_2.h90	(revision 7342)
@@ -41,4 +41,7 @@
       USE in_out_manager  ! I/O Manager
       USE iom
+#if defined key_bdy
+      USE bdy_oce
+#endif
       !! 
       INTEGER :: ji, jj, jn, jproc, jarea     ! dummy loop indices
@@ -73,5 +76,13 @@
       ! read namelist for ln_zco
       NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
-
+#if defined key_bdy
+      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,                 &
+         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
+         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &  
+         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
+         &             cn_ice_lim, nn_ice_lim_dta,                           &
+         &             rn_ice_tem, rn_ice_sal, rn_ice_age,                 &
+         &             ln_vol, nn_volctl, nn_rimwidth, nb_jpk_bdy
+#endif
       !!----------------------------------------------------------------------
       !!  OPA 9.0 , LOCEAN-IPSL (2005) 
@@ -137,4 +148,24 @@
       imask(:,:)=1
       WHERE ( zdta(:,:) - zdtaisf(:,:) <= rn_isfhmin ) imask = 0
+
+#if defined key_bdy
+      ! Adjust imask with bdy_msk if exists
+
+      REWIND( numnam_ref )              ! Namelist nambdy in reference namelist : BDY
+      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in reference namelist (mppini_2)', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nambdy in configuration namelist : BDY
+      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 904 )
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy in configuration namelist (mppini_2)', lwp )
+      IF(lwm) WRITE ( numond, namzgr )
+
+      IF( ln_mask_file ) THEN
+         CALL iom_open( cn_mask_file, inum )
+         CALL iom_get ( inum, jpdom_unknown, 'bdy_msk', zdta(:,:), kstart=(/jpizoom,jpjzoom/), kcount=(/jpiglo,jpjglo/) )
+         CALL iom_close( inum )
+         WHERE ( zdta(:,:) <= 0. ) imask = 0
+      ENDIF
+#endif
 
       !  1. Dimension arrays for subdomains
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 7342)
@@ -7,4 +7,5 @@
    !!                 !  05-2008  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid
    !!                 !  10-2013  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation
+   !!                 !  12-2015  (J. Harle) Adding BDY on-the-fly interpolation
    !!----------------------------------------------------------------------
 
@@ -67,4 +68,6 @@
       INTEGER                         ::   nreclast     ! last record to be read in the current file
       CHARACTER(len = 256)            ::   lsmname      ! current name of the NetCDF mask file acting as a key
+      INTEGER                         ::   igrd         ! grid type for bdy data
+      INTEGER                         ::   ibdy         ! bdy set id number
    END TYPE FLD
 
@@ -114,5 +117,5 @@
 CONTAINS
 
-   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset )
+   SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_read  ***
@@ -135,4 +138,7 @@
       !                                                     !   kt_offset = +1 => fields at "after"  time level
       !                                                     !   etc.
+      INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data
+      LOGICAL  , INTENT(in   ), OPTIONAL     ::   fvl   ! number of vertical levels in the BDY data
+      !!
       INTEGER  ::   itmp         ! local variable
       INTEGER  ::   imf          ! size of the structure sd
@@ -171,5 +177,9 @@
          DO jf = 1, imf 
             IF( PRESENT(map) ) imap = map(jf)
-            CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped)
+               IF( PRESENT(jpk_bdy) ) THEN
+                  CALL fld_init( kn_fsbc, sd(jf), imap, jpk_bdy, fvl )  ! read each before field (put them in after as they will be swapped)
+               ELSE
+                  CALL fld_init( kn_fsbc, sd(jf), imap )  ! read each before field (put them in after as they will be swapped)
+               ENDIF
          END DO
          IF( lwp ) CALL wgt_print()                ! control print
@@ -263,6 +273,9 @@
 
                ! read after data
-               CALL fld_get( sd(jf), imap )
-
+               IF( PRESENT(jpk_bdy) ) THEN
+                  CALL fld_get( sd(jf), imap, jpk_bdy, fvl)
+               ELSE
+                  CALL fld_get( sd(jf), imap )
+               ENDIF
             ENDIF   ! read new data?
          END DO                                    ! --- end loop over field --- !
@@ -302,5 +315,5 @@
 
 
-   SUBROUTINE fld_init( kn_fsbc, sdjf, map )
+   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_init  ***
@@ -309,7 +322,9 @@
       !!               - if time interpolation, read before data 
       !!----------------------------------------------------------------------
-      INTEGER  , INTENT(in   ) ::   kn_fsbc   ! sbc computation period (in time step) 
-      TYPE(FLD), INTENT(inout) ::   sdjf      ! input field related variables
-      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
+      INTEGER  , INTENT(in   ) ::   kn_fsbc      ! sbc computation period (in time step) 
+      TYPE(FLD), INTENT(inout) ::   sdjf         ! input field related variables
+      TYPE(MAP_POINTER),INTENT(in) ::   map      ! global-to-local mapping indices
+      INTEGER  , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data
+      LOGICAL  , INTENT(in), OPTIONAL :: fvl     ! number of vertical levels in the BDY data
       !!
       LOGICAL :: llprevyr              ! are we reading previous year  file?
@@ -405,5 +420,10 @@
          ENDIF
          !
-         CALL fld_get( sdjf, map )         ! read before data in after arrays(as we will swap it later)
+         ! read before data in after arrays(as we will swap it later)
+         IF( PRESENT(jpk_bdy) ) THEN
+            CALL fld_get( sdjf, map, jpk_bdy, fvl )
+         ELSE
+            CALL fld_get( sdjf, map )
+         ENDIF
          !
          clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
@@ -581,5 +601,5 @@
 
 
-   SUBROUTINE fld_get( sdjf, map )
+   SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_get  ***
@@ -589,4 +609,6 @@
       TYPE(FLD)        , INTENT(inout) ::   sdjf   ! input field related variables
       TYPE(MAP_POINTER), INTENT(in   ) ::   map    ! global-to-local mapping indices
+      INTEGER  , INTENT(in), OPTIONAL  ::   jpk_bdy ! number of vertical levels in the bdy data
+      LOGICAL  , INTENT(in), OPTIONAL  ::   fvl     ! number of vertical levels in the bdy data
       !
       INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
@@ -601,7 +623,15 @@
       !
       IF( ASSOCIATED(map%ptr) ) THEN
-         IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
-         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map )
-         ENDIF
+         IF( PRESENT(jpk_bdy) ) THEN
+            IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2),                &
+                                                        sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )
+            ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ),                &
+                                                        sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )
+            ENDIF
+         ELSE
+            IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map )
+            ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map )
+            ENDIF
+         ENDIF        
       ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
          CALL wgt_list( sdjf, iw )
@@ -658,6 +688,5 @@
    END SUBROUTINE fld_get
 
-
-   SUBROUTINE fld_map( num, clvar, dta, nrec, map )
+   SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_map  ***
@@ -667,11 +696,14 @@
       !!----------------------------------------------------------------------
 #if defined key_bdy
-      USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays
+      USE bdy_oce, ONLY:  idx_bdy, dta_global, dta_global_z, dta_global_dz, dta_global2, dta_global2_z, dta_global2_dz                 ! workspace to read in global data arrays
 #endif 
       INTEGER                   , INTENT(in ) ::   num     ! stream number
       CHARACTER(LEN=*)          , INTENT(in ) ::   clvar   ! variable name
-      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta   ! output field on model grid (2 dimensional)
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta     ! output field on model grid (2 dimensional)
       INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
       TYPE(MAP_POINTER)         , INTENT(in ) ::   map     ! global-to-local mapping indices
+      INTEGER  , INTENT(in), OPTIONAL         ::   igrd, ibdy, jpk_bdy  ! grid type, set number and number of vertical levels in the bdy data
+      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl     ! grid type, set number and number of vertical levels in the bdy data
+      INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1
       !!
       INTEGER                                 ::   ipi      ! length of boundary data on local process
@@ -682,5 +714,8 @@
       INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
       INTEGER                                 ::   ierr
-      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data
+      REAL(wp)                                ::   fv          ! fillvalue 
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_dz ! work space for global data
       !!---------------------------------------------------------------------
       !
@@ -696,6 +731,20 @@
       IF( map%ll_unstruc) THEN   ! unstructured open boundary data file
          dta_read => dta_global
-      ELSE                       ! structured open boundary data file
+         IF( PRESENT(jpk_bdy) ) THEN
+            IF( jpk_bdy>0 ) THEN
+               dta_read_z => dta_global_z
+               dta_read_dz => dta_global_dz
+               jpkm1_bdy = jpk_bdy-1
+            ENDIF
+         ENDIF
+      ELSE                       ! structured open boundary file
          dta_read => dta_global2
+         IF( PRESENT(jpk_bdy) ) THEN
+            IF( jpk_bdy>0 ) THEN
+               dta_read_z => dta_global2_z
+               dta_read_dz => dta_global2_dz
+               jpkm1_bdy = jpk_bdy-1
+            ENDIF
+         ENDIF
       ENDIF
 #endif
@@ -705,26 +754,373 @@
       !
       SELECT CASE( ipk )
-      CASE(1)        ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
-      CASE DEFAULT   ;   CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )
+      CASE(1)        ;   
+      CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1    ), nrec )
+         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
+            DO ib = 1, ipi
+              DO ik = 1, ipk
+                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
+              END DO
+            END DO
+         ELSE ! we assume that this is a structured open boundary file
+            DO ib = 1, ipi
+               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+               ji=map%ptr(ib)-(jj-1)*ilendta
+               DO ik = 1, ipk
+                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
+               END DO
+            END DO
+         ENDIF
+
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      ! Do we include something here to adjust barotropic velocities !
+      ! in case of a depth difference between bdy files and          !
+      ! bathymetry in the case ln_full_vel = .false. and jpk_bdy>0?  !
+      ! [as the enveloping and parital cells could change H]         !
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      CASE DEFAULT   ;
+
+      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation
+         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar )
+         dta_read(:,:,:) = -ABS(fv)
+         dta_read_z(:,:,:) = 0._wp
+         dta_read_dz(:,:,:) = 0._wp
+         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:jpk_bdy), nrec )
+         SELECT CASE( igrd )                  
+            CASE(1)
+               CALL iom_get ( num, jpdom_unknown, 'gdept', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
+               CALL iom_get ( num, jpdom_unknown, 'e3t',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
+            CASE(2)  
+               CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
+               CALL iom_get ( num, jpdom_unknown, 'e3u',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
+            CASE(3)
+               CALL iom_get ( num, jpdom_unknown, 'gdepv', dta_read_z(1:ilendta,1:ipj,1:jpk_bdy) )
+               CALL iom_get ( num, jpdom_unknown, 'e3v',  dta_read_dz(1:ilendta,1:ipj,1:jpk_bdy) )
+         END SELECT
+
+#if defined key_bdy
+         CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
+#endif
+      ELSE ! boundary data assumed to be on model grid
+         
+         CALL iom_get ( num, jpdom_unknown, clvar, dta_read(1:ilendta,1:ipj,1:ipk), nrec )                    
+         IF ( map%ll_unstruc) THEN ! unstructured open boundary data file
+            DO ib = 1, ipi
+              DO ik = 1, ipk
+                dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
+              END DO
+            END DO
+         ELSE ! we assume that this is a structured open boundary file
+            DO ib = 1, ipi
+               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+               ji=map%ptr(ib)-(jj-1)*ilendta
+               DO ik = 1, ipk
+                  dta(ib,1,ik) =  dta_read(ji,jj,ik)
+               END DO
+            END DO
+         ENDIF
+      ENDIF ! PRESENT(jpk_bdy)
       END SELECT
-      !
-      IF( map%ll_unstruc ) THEN  ! unstructured open boundary data file
+
+   END SUBROUTINE fld_map
+   
+#if defined key_bdy
+   SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
+
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE fld_bdy_interp  ***
+      !!
+      !! ** Purpose :   on the fly vertical interpolation to allow the use of
+      !!                boundary data from non-native vertical grid
+      !!----------------------------------------------------------------------
+      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation
+
+      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read      ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_z    ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:), INTENT(in )     ::   dta_read_dz   ! work space for global data
+      REAL(wp) , INTENT(in)                                ::   fv            ! fillvalue and alternative -ABS(fv)
+      REAL(wp), DIMENSION(:,:,:), INTENT(out) ::   dta                        ! output field on model grid (2 dimensional)
+      TYPE(MAP_POINTER)         , INTENT(in ) ::   map                        ! global-to-local mapping indices
+      LOGICAL  , INTENT(in), OPTIONAL         ::   fvl                        ! grid type, set number and number of vertical levels in the bdy data
+      INTEGER  , INTENT(in)                   ::   igrd, ibdy, jpk_bdy        ! number of levels in bdy data
+      INTEGER  , INTENT(in)                   ::   ilendta                    ! length of data in file
+      !!
+      INTEGER                                 ::   ipi                        ! length of boundary data on local process
+      INTEGER                                 ::   ipj                        ! length of dummy dimension ( = 1 )
+      INTEGER                                 ::   ipk                        ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
+      INTEGER                                 ::   jpkm1_bdy                  ! number of levels in bdy data minus 1
+      INTEGER                                 ::   ib, ik, ikk                ! loop counters
+      INTEGER                                 ::   ji, jj, zij, zjj           ! temporary indices
+      REAL(wp)                                ::   zl, zi, zh                 ! tmp variable for current depth and interpolation factor
+      REAL(wp)                                ::   fv_alt, ztrans, ztrans_new ! fillvalue and alternative -ABS(fv)
+      CHARACTER (LEN=10)                      ::   ibstr
+      !!---------------------------------------------------------------------
+     
+
+      ipi       = SIZE( dta, 1 )
+      ipj       = SIZE( dta_read, 2 )
+      ipk       = SIZE( dta, 3 )
+      jpkm1_bdy = jpk_bdy-1
+      
+      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
+      DO ib = 1, ipi
+            zij = idx_bdy(ibdy)%nbi(ib,igrd)
+            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+         IF(narea==2) WRITE(*,*) 'MAPI', ib, igrd, map%ptr(ib), narea-1, zij, zjj
+      ENDDO
+      !
+      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file
+
          DO ib = 1, ipi
-            DO ik = 1, ipk
-               dta(ib,1,ik) =  dta_read(map%ptr(ib),1,ik)
+            DO ik = 1, jpk_bdy
+               IF( ( dta_read(map%ptr(ib),1,ik) == fv ) ) THEN
+                  dta_read_z(map%ptr(ib),1,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data
+                  dta_read_dz(map%ptr(ib),1,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct
+               ENDIF
+            ENDDO
+         ENDDO 
+
+         DO ib = 1, ipi
+            zij = idx_bdy(ibdy)%nbi(ib,igrd)
+            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+            zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
+            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
+            SELECT CASE( igrd )                         
+               CASE(1)
+                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN
+                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
+                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
+                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
+                  ENDIF
+               CASE(2)
+                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN
+                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
+                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
+                     IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u_n(zij,zjj,:), mask=umask(zij,zjj,:)==1),  sum(umask(zij,zjj,:)), &
+                       &                hu_n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1  , &
+                        &                dta_read(map%ptr(ib),1,:)
+                  ENDIF
+               CASE(3)
+                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN
+                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
+                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
+                  ENDIF
+            END SELECT
+            DO ik = 1, ipk                      
+               SELECT CASE( igrd )                       
+                  CASE(1)
+                     zl =  gdept_n(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_n?
+                  CASE(2)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                     ELSE
+                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) ) 
+                     ENDIF
+                  CASE(3)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                     ELSE
+                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
+                     ENDIF
+               END SELECT
+               IF( zl < dta_read_z(map%ptr(ib),1,1) ) THEN                                         ! above the first level of external data
+                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,1)
+               ELSEIF( zl > MAXVAL(dta_read_z(map%ptr(ib),1,:),1) ) THEN                           ! below the last level of external data 
+                  dta(ib,1,ik) =  dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1))
+               ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1
+                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
+                     IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) &
+                    &    .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN
+                        zi           = ( zl - dta_read_z(map%ptr(ib),1,ikk) ) / &
+                       &               ( dta_read_z(map%ptr(ib),1,ikk+1) - dta_read_z(map%ptr(ib),1,ikk) )
+                        dta(ib,1,ik) = dta_read(map%ptr(ib),1,ikk) + &
+                       &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
          END DO
-      ELSE                       ! structured open boundary data file
+
+         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
+            DO ib = 1, ipi
+              zij = idx_bdy(ibdy)%nbi(ib,igrd)
+              zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
+              ztrans = 0._wp
+              ztrans_new = 0._wp
+              DO ik = 1, jpk_bdy                            ! calculate transport on input grid
+                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
+              ENDDO
+              DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
+              ENDDO
+              DO ik = 1, ipk                                ! make transport correction
+                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
+                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                    IF( ABS(ztrans * r1_hu_n(zij,zjj)) > 0.01_wp ) &
+                   &   CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at')
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu_n(zij,zjj) * umask(zij,zjj,ik)
+                 ENDIF
+              ENDDO
+            ENDDO
+         ENDIF
+
+         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
+            DO ib = 1, ipi
+              zij = idx_bdy(ibdy)%nbi(ib,igrd)
+              zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+              zh  = SUM(dta_read_dz(map%ptr(ib),1,:) )
+              ztrans = 0._wp
+              ztrans_new = 0._wp
+              DO ik = 1, jpk_bdy                            ! calculate transport on input grid
+                  ztrans     = ztrans     + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
+              ENDDO
+              DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
+              ENDDO
+              DO ik = 1, ipk                                ! make transport correction
+                 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
+                 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv_n(zij,zjj) * vmask(zij,zjj,ik)
+                 ENDIF
+              ENDDO
+            ENDDO
+         ENDIF
+  
+      ELSE ! structured open boundary file
+
          DO ib = 1, ipi
             jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
             ji=map%ptr(ib)-(jj-1)*ilendta
-            DO ik = 1, ipk
-               dta(ib,1,ik) =  dta_read(ji,jj,ik)
+            DO ik = 1, jpk_bdy                      
+               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN
+                  dta_read_z(ji,jj,ik)  = fv_alt ! safety: put fillvalue into external depth field so consistent with data
+                  dta_read_dz(ji,jj,ik) = 0._wp  ! safety: put 0._wp into external thickness factors to ensure transport is correct
+               ENDIF
+            ENDDO
+         ENDDO 
+       
+
+         DO ib = 1, ipi
+            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+            ji=map%ptr(ib)-(jj-1)*ilendta
+            zij = idx_bdy(ibdy)%nbi(ib,igrd)
+            zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+            zh  = SUM(dta_read_dz(ji,jj,:) )
+            ! Warnings to flag differences in the input and model topgraphy - is this useful/necessary?
+            SELECT CASE( igrd )                         
+               CASE(1)
+                  IF( ABS( (zh - ht_n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN
+                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
+                     CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
+                 !   IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t_n(zij,zjj,:), mask=tmask(zij,zjj,:)==1),  ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj
+                  ENDIF
+               CASE(2)
+                  IF( ABS( (zh - hu_n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN
+                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
+                     CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
+                  ENDIF
+               CASE(3)
+                  IF( ABS( (zh - hv_n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN
+                     WRITE(ibstr,"(I10.10)") map%ptr(ib) 
+                     CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%')
+                  ENDIF
+            END SELECT
+            DO ik = 1, ipk                      
+               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min
+                  CASE(1)
+                     zl =  gdept_n(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_n?
+                  CASE(2)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                     ELSE
+                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij+1,zjj,ik) )
+                     ENDIF
+                  CASE(3)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_n(zij,zjj,ik) + gdept_n(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_n?
+                     ELSE
+                        zl =  MIN( gdept_n(zij,zjj,ik), gdept_n(zij,zjj+1,ik) )
+                     ENDIF
+               END SELECT
+               IF( zl < dta_read_z(ji,jj,1) ) THEN                                      ! above the first level of external data
+                  dta(ib,1,ik) =  dta_read(ji,jj,1)
+               ELSEIF( zl > MAXVAL(dta_read_z(ji,jj,:),1) ) THEN                        ! below the last level of external data 
+                  dta(ib,1,ik) =  dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1))
+               ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1
+                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_n(ikk) < zl < gdept_n(ikk+1)
+                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) &
+                    &    .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN
+                        zi           = ( zl - dta_read_z(ji,jj,ikk) ) / &
+                       &               ( dta_read_z(ji,jj,ikk+1) - dta_read_z(ji,jj,ikk) )
+                        dta(ib,1,ik) = dta_read(ji,jj,ikk) + &
+                       &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
          END DO
-      ENDIF
-      !
-   END SUBROUTINE fld_map
-
+
+         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
+            DO ib = 1, ipi
+               jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+               ji=map%ptr(ib)-(jj-1)*ilendta
+               zij = idx_bdy(ibdy)%nbi(ib,igrd)
+               zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+               zh = SUM(dta_read_dz(ji,jj,:) )
+               ztrans = 0._wp
+               ztrans_new = 0._wp
+               DO ik = 1, jpk_bdy                            ! calculate transport on input grid
+                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
+               ENDDO
+               DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_n(zij,zjj,ik) * umask(zij,zjj,ik)
+               ENDDO
+               DO ik = 1, ipk                                ! make transport correction
+                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
+                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hu_n(zij,zjj) ) * umask(zij,zjj,ik)
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDIF
+
+         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
+            DO ib = 1, ipi
+               jj  = 1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
+               ji  = map%ptr(ib)-(jj-1)*ilendta
+               zij = idx_bdy(ibdy)%nbi(ib,igrd)
+               zjj = idx_bdy(ibdy)%nbj(ib,igrd)
+               zh  = SUM(dta_read_dz(ji,jj,:) )
+               ztrans = 0._wp
+               ztrans_new = 0._wp
+               DO ik = 1, jpk_bdy                            ! calculate transport on input grid
+                  ztrans     = ztrans     + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
+               ENDDO
+               DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_n(zij,zjj,ik) * vmask(zij,zjj,ik)
+               ENDDO
+               DO ik = 1, ipk                                ! make transport correction
+                  IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data
+                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
+                  ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero
+                     dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * r1_hv_n(zij,zjj) ) * vmask(zij,zjj,ik)
+                  ENDIF
+               ENDDO
+            ENDDO
+         ENDIF
+
+      ENDIF ! endif unstructured or structured
+
+   END SUBROUTINE fld_bdy_interp
+
+!  SUBROUTINE fld_bdy_conserve(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
+
+!  END SUBROUTINE fld_bdy_conserve
+
+#endif
 
    SUBROUTINE fld_rot( kt, sd )
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90	(revision 7342)
@@ -162,5 +162,8 @@
                   &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
                   &          + avtb(jk) * tmask(ji,jj,jk)
-               !                                            ! Add the background coefficient on eddy viscosity
+            END DO
+         END DO
+         DO jj = 2, jpjm1                                   ! Add the background coefficient on eddy viscosity
+            DO ji = fs_2, fs_jpim1
                avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
                avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/TOP_SRC/trcrst.F90	(revision 7342)
@@ -295,5 +295,5 @@
          IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift
       END DO
-      WRITE(numout,*) 
+      IF(lwp) WRITE(numout,*) 
 9000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, &
       &      '    max :',e18.10,'    drift :',e18.10, ' %')
Index: branches/2016/dev_NOC_2016/NEMOGCM/NEMO/TOP_SRC/trcsub.F90
===================================================================
--- branches/2016/dev_NOC_2016/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 7341)
+++ branches/2016/dev_NOC_2016/NEMOGCM/NEMO/TOP_SRC/trcsub.F90	(revision 7342)
@@ -16,10 +16,4 @@
    USE in_out_manager
    USE lbclnk
-#if defined key_zdftke
-   USE zdftke          ! twice TKE (en)
-#endif
-#if defined key_zdfgls
-   USE zdfgls, ONLY: en
-#endif
    USE trabbl
    USE zdf_oce
