Index: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 6140)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 7646)
@@ -4,12 +4,26 @@
    !! Ocean forcing:  read input field for surface boundary condition
    !!=====================================================================
-   !! History :  2.0  !  06-2006  (S. Masson, G. Madec)  Original code
-   !!                 !  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
+   !! History :  2.0  !  2006-06  (S. Masson, G. Madec)  Original code
+   !!            3.0  !  2008-05  (S. Alderson)  Modified for Interpolation in memory from input grid to model grid
+   !!            3.4  !  2013-10  (D. Delrosso, P. Oddo)  suppression of land point prior to interpolation
+   !!                 !  12-2015  (J. Harle) Adding BDY on-the-fly interpolation
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   fld_read      : read input fields used for the computation of the
-   !!                   surface boundary condition
+   !!   fld_read      : read input fields used for the computation of the surface boundary condition
+   !!   fld_init      : initialization of field read
+   !!   fld_rec       : determined the record(s) to be read
+   !!   fld_get       : read the data
+   !!   fld_map       : read global data from file and map onto local data using a general mapping (use for open boundaries)
+   !!   fld_rot       : rotate the vector fields onto the local grid direction
+   !!   fld_clopn     : update the data file name and close/open the files
+   !!   fld_fill      : fill the data structure with the associated information read in namelist
+   !!   wgt_list      : manage the weights used for interpolation
+   !!   wgt_print     : print the list of known weights
+   !!   fld_weight    : create a WGT structure and fill in data from file, restructuring as required
+   !!   apply_seaoverland : fill land with ocean values
+   !!   seaoverland   : create shifted matrices for seaoverland application
+   !!   fld_interp    : apply weights to input gridded data to create data on model grid
+   !!   ksec_week     : function returning the first 3 letters of the first day of the weekly file
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and tracers
@@ -67,4 +81,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 +130,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 +151,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 +190,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 +286,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 --- !
@@ -274,9 +300,9 @@
             IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation
                IF(lwp .AND. kt - nit000 <= 100 ) THEN 
-                  clfmt = "('fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
+                  clfmt = "('   fld_read: var ', a, ' kt = ', i8, ' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')"
                   WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &            
                      & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday
-                  WRITE(numout, *) 'it_offset is : ',it_offset
+                  WRITE(numout, *) '      it_offset is : ',it_offset
                ENDIF
                ! temporal interpolation weights
@@ -286,5 +312,5 @@
             ELSE   ! nothing to do...
                IF(lwp .AND. kt - nit000 <= 100 ) THEN
-                  clfmt = "('fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
+                  clfmt = "('   fld_read: var ', a, ' kt = ', i8,' (', f9.4,' days), Y/M/D = ', i4.4,'/', i2.2,'/', i2.2," //   &
                      &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')"
                   WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    &
@@ -302,5 +328,5 @@
 
 
-   SUBROUTINE fld_init( kn_fsbc, sdjf, map )
+   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy, fvl)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_init  ***
@@ -309,7 +335,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,7 +433,12 @@
          ENDIF
          !
-         CALL fld_get( sdjf, map )         ! read before data in after arrays(as we will swap it later)
-         !
-         clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
+         ! 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')"
          IF(lwp) WRITE(numout, clfmt) TRIM(sdjf%clvar), sdjf%nrec_a(1), REAL(sdjf%nrec_a(2),wp)/rday
          !
@@ -581,5 +614,5 @@
 
 
-   SUBROUTINE fld_get( sdjf, map )
+   SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_get  ***
@@ -589,4 +622,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 +636,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 +701,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  ***
@@ -666,12 +708,15 @@
       !!                using a general mapping (for open boundaries)
       !!----------------------------------------------------------------------
-#if defined key_bdy
-      USE bdy_oce, ONLY:  dta_global, dta_global2         ! workspace to read in global data arrays
-#endif 
+
+      USE bdy_oce, ONLY: ln_bdy, 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
+
       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 +727,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
       !!---------------------------------------------------------------------
       !
@@ -692,12 +740,26 @@
       ilendta = iom_file(num)%dimsz(1,idvar)
 
-#if defined key_bdy
-      ipj = iom_file(num)%dimsz(2,idvar)
-      IF( map%ll_unstruc) THEN   ! unstructured open boundary data file
-         dta_read => dta_global
-      ELSE                       ! structured open boundary data file
-         dta_read => dta_global2
+      IF ( ln_bdy ) THEN
+         ipj = iom_file(num)%dimsz(2,idvar)
+         IF( map%ll_unstruc) THEN   ! unstructured open boundary data file
+            dta_read => dta_global
+            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
-#endif
 
       IF(lwp) WRITE(numout,*) 'Dim size for ',        TRIM(clvar),' is ', ilendta
@@ -705,25 +767,366 @@
       !
       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 ( ln_bdy ) & 
+         CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta)
+
+      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
+   
+   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
 
 
@@ -791,5 +1194,5 @@
       !!                    ***  ROUTINE fld_clopn  ***
       !!
-      !! ** Purpose :   update the file name and open the file
+      !! ** Purpose :   update the file name and close/open the files
       !!----------------------------------------------------------------------
       TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
@@ -882,18 +1285,20 @@
 
 
-   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam )
+   SUBROUTINE fld_fill( sdf, sdf_n, cdir, cdcaller, cdtitle, cdnam, knoprint )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_fill  ***
       !!
-      !! ** Purpose :   fill sdf with sdf_n and control print
+      !! ** Purpose :   fill the data structure (sdf) with the associated information 
+      !!              read in namelist (sdf_n) and control print
       !!----------------------------------------------------------------------
-      TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
-      TYPE(FLD_N), DIMENSION(:), INTENT(in   ) ::   sdf_n      ! array of namelist information structures
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdcaller   ! 
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdtitle    ! 
-      CHARACTER(len=*)         , INTENT(in   ) ::   cdnam      ! 
-      !
-      INTEGER  ::   jf       ! dummy indices
+      TYPE(FLD)  , DIMENSION(:)          , INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
+      TYPE(FLD_N), DIMENSION(:)          , INTENT(in   ) ::   sdf_n      ! array of namelist information structures
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdir       ! Root directory for location of flx files
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdcaller   ! name of the calling routine
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdtitle    ! description of the calling routine 
+      CHARACTER(len=*)                   , INTENT(in   ) ::   cdnam      ! name of the namelist from which sdf_n comes
+      INTEGER                  , OPTIONAL, INTENT(in   ) ::   knoprint   ! no calling routine information printed
+      !
+      INTEGER  ::   jf   ! dummy indices
       !!---------------------------------------------------------------------
       !
@@ -922,18 +1327,20 @@
       IF(lwp) THEN      ! control print
          WRITE(numout,*)
-         WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
-         WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
-         WRITE(numout,*) '          '//TRIM( cdnam )//' Namelist'
-         WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
+         IF( .NOT.PRESENT( knoprint) ) THEN
+            WRITE(numout,*) TRIM( cdcaller )//' : '//TRIM( cdtitle )
+            WRITE(numout,*) (/ ('~', jf = 1, LEN_TRIM( cdcaller ) ) /)
+         ENDIF
+         WRITE(numout,*) '   fld_fill : fill data structure with information from namelist '//TRIM( cdnam )
+         WRITE(numout,*) '   ~~~~~~~~'
+         WRITE(numout,*) '      list of files and frequency (>0: in hours ; <0 in months)'
          DO jf = 1, SIZE(sdf)
-            WRITE(numout,*) '               root filename: '  , TRIM( sdf(jf)%clrootname ),   &
-               &                          ' variable name: '  , TRIM( sdf(jf)%clvar      )
-            WRITE(numout,*) '               frequency: '      ,       sdf(jf)%nfreqh      ,   &
-               &                          ' time interp: '    ,       sdf(jf)%ln_tint     ,   &
-               &                          ' climatology: '    ,       sdf(jf)%ln_clim     ,   &
-               &                          ' weights    : '    , TRIM( sdf(jf)%wgtname    ),   &
-               &                          ' pairing    : '    , TRIM( sdf(jf)%vcomp      ),   &
-               &                          ' data type: '      ,       sdf(jf)%cltype      ,   &
-               &                          ' land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
+            WRITE(numout,*) '      root filename: '  , TRIM( sdf(jf)%clrootname ), '   variable name: ', TRIM( sdf(jf)%clvar )
+            WRITE(numout,*) '         frequency: '      ,       sdf(jf)%nfreqh      ,   &
+               &                  '   time interp: '    ,       sdf(jf)%ln_tint     ,   &
+               &                  '   climatology: '    ,       sdf(jf)%ln_clim
+            WRITE(numout,*) '         weights: '        , TRIM( sdf(jf)%wgtname    ),   &
+               &                  '   pairing: '        , TRIM( sdf(jf)%vcomp      ),   &
+               &                  '   data type: '      ,       sdf(jf)%cltype      ,   &
+               &                  '   land/sea mask:'   , TRIM( sdf(jf)%lsmname    )
             call flush(numout)
          END DO
@@ -947,8 +1354,7 @@
       !!                    ***  ROUTINE wgt_list  ***
       !!
-      !! ** Purpose :   search array of WGTs and find a weights file
-      !!                entry, or return a new one adding it to the end
-      !!                if it is a new entry, the weights data is read in and
-      !!                restructured (fld_weight)
+      !! ** Purpose :   search array of WGTs and find a weights file entry,
+      !!                or return a new one adding it to the end if new entry.
+      !!                the weights data is read in and restructured (fld_weight)
       !!----------------------------------------------------------------------
       TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
@@ -1019,6 +1425,6 @@
       !!                    ***  ROUTINE fld_weight  ***
       !!
-      !! ** Purpose :   create a new WGT structure and fill in data from  
-      !!                file, restructuring as required
+      !! ** Purpose :   create a new WGT structure and fill in data from file,
+      !!              restructuring as required
       !!----------------------------------------------------------------------
       TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
@@ -1163,5 +1569,5 @@
 
    SUBROUTINE apply_seaoverland( clmaskfile, zfieldo, jpi1_lsm, jpi2_lsm, jpj1_lsm,   &
-                          &      jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
+      &                          jpj2_lsm, itmpi, itmpj, itmpz, rec1_lsm, recn_lsm )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE apply_seaoverland  ***
@@ -1492,11 +1898,10 @@
       !!                    ***  FUNCTION kshift_week *** 
       !!
-      !! ** Purpose :  
-      !!---------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
-      !!
-      INTEGER                        ::   ksec_week  ! output variable
-      INTEGER                        ::   ijul       !temp variable
-      INTEGER                        ::   ishift     !temp variable
+      !! ** Purpose :   return the first 3 letters of the first day of the weekly file
+      !!---------------------------------------------------------------------
+      CHARACTER(len=*), INTENT(in)   ::   cdday   ! first 3 letters of the first day of the weekly file
+      !!
+      INTEGER                        ::   ksec_week      ! output variable
+      INTEGER                        ::   ijul, ishift   ! local integer
       CHARACTER(len=3),DIMENSION(7)  ::   cl_week 
       !!----------------------------------------------------------------------
