Index: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 4622)
+++ branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 4694)
@@ -67,4 +67,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
 
@@ -113,5 +115,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 )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_read  ***
@@ -134,4 +136,6 @@
                                                             !   kt_offset = +1 => fields at "after"  time level
                                                             !   etc.
+      !!
+      INTEGER  , INTENT(in   ), OPTIONAL     ::   jpk_bdy   ! number of vertical levels in the BDY data
       !!
       INTEGER  ::   itmp       ! temporary variable
@@ -169,5 +173,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 )  ! 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
@@ -261,6 +269,9 @@
 
                ! read after data
-               CALL fld_get( sd(jf), imap )
-
+               IF( PRESENT(jpk_bdy) ) THEN
+                  CALL fld_get( sd(jf), imap, jpk_bdy)
+               ELSE
+                  CALL fld_get( sd(jf), imap )
+               ENDIF
             ENDIF   ! read new data?
          END DO                                    ! --- end loop over field --- !
@@ -303,5 +314,5 @@
 
 
-   SUBROUTINE fld_init( kn_fsbc, sdjf, map )
+   SUBROUTINE fld_init( kn_fsbc, sdjf, map , jpk_bdy)
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_init  ***
@@ -310,7 +321,8 @@
       !!               - 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 :: llprevyr              ! are we reading previous year  file?
@@ -407,5 +419,9 @@
 
          ! read before data in after arrays(as we will swap it later)
-         CALL fld_get( sdjf, map )
+         IF( PRESENT(jpk_bdy) ) THEN
+            CALL fld_get( sdjf, map, jpk_bdy )
+         ELSE
+            CALL fld_get( sdjf, map )
+         ENDIF
 
          clfmt = "('fld_init : time-interpolation for ', a, ' read previous record = ', i6, ' at time = ', f7.2, ' days')"
@@ -577,5 +593,5 @@
 
 
-   SUBROUTINE fld_get( sdjf, map )
+   SUBROUTINE fld_get( sdjf, map, jpk_bdy )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_get  ***
@@ -583,6 +599,7 @@
       !! ** Purpose :   read the data
       !!----------------------------------------------------------------------
-      TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables
-      TYPE(MAP_POINTER),INTENT(in) ::   map   ! global-to-local mapping indices
+      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
       !!
       INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
@@ -597,13 +614,23 @@
       !
       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%ptr )
-         ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map%ptr )
+         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%ptr, sdjf%igrd, sdjf%ibdy, jpk_bdy )
+            ELSE                      ;   
+               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map%ptr, sdjf%igrd, sdjf%ibdy, jpk_bdy )
+            ENDIF
+         ELSE
+            IF( sdjf%ln_tint ) THEN   ;   
+               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), map%ptr )
+            ELSE                      ;   
+               CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), map%ptr )
+            ENDIF
          ENDIF
       ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN
          CALL wgt_list( sdjf, iw )
-         IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), & 
-              & sdjf%nrec_a(1), sdjf%lsmname )
-         ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), &
-              & sdjf%nrec_a(1), sdjf%lsmname )
+         IF( sdjf%ln_tint ) THEN   ;   
+            CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fdta(:,:,:,2), sdjf%nrec_a(1), sdjf%lsmname )
+         ELSE                      ;   
+            CALL fld_interp( sdjf%num, sdjf%clvar, iw , ipk  , sdjf%fnow(:,:,:  ), sdjf%nrec_a(1), sdjf%lsmname )
          ENDIF
       ELSE
@@ -654,5 +681,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 )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_map  ***
@@ -669,4 +696,6 @@
       INTEGER                   , INTENT(in ) ::   nrec    ! record number to read (ie time slice)
       INTEGER,  DIMENSION(:)    , 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
+      INTEGER                                 ::   jpkm1_bdy! number of vertical levels in the bdy data minus 1
       !!
       INTEGER                                 ::   ipi      ! length of boundary data on local process
@@ -677,5 +706,6 @@
       INTEGER                                 ::   ib, ik, ji, jj   ! loop counters
       INTEGER                                 ::   ierr
-      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read  ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read    ! work space for global data
+      REAL(wp), POINTER, DIMENSION(:,:,:)     ::   dta_read_z  ! work space for global data
       !!---------------------------------------------------------------------
             
@@ -689,8 +719,20 @@
 #if defined key_bdy
       ipj = iom_file(num)%dimsz(2,idvar)
-      IF (ipj == 1) THEN ! we assume that this is a structured open boundary file
+      IF (ipj == 1) THEN 
          dta_read => dta_global
-      ELSE
+         IF( PRESENT(jpk_bdy) ) THEN
+            IF( jpk_bdy>0 ) THEN
+               dta_read_z => dta_global_z
+               jpkm1_bdy = jpk_bdy-1
+            ENDIF
+         ENDIF
+      ELSE ! we assume that this is a structured open boundary file
          dta_read => dta_global2
+         IF( PRESENT(jpk_bdy) ) THEN
+            IF( jpk_bdy>0 ) THEN
+               dta_read_z => dta_global2_z
+               jpkm1_bdy = jpk_bdy-1
+            ENDIF
+         ENDIF
       ENDIF
 #endif
@@ -701,11 +743,100 @@
       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 DEFAULT   ;   
+  
+      IF( PRESENT(jpk_bdy) .AND. jpk_bdy>0 ) THEN       ! boundary data not on model grid: vertical interpolation
+         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) )
+            CASE(2); CALL iom_get ( num, jpdom_unknown, 'gdepu', dta_read_z(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) )
+         END SELECT
+         CALL iom_getatt(num, '_FillValue', fv, cdvar=clvar )
+         CALL fld_bdy_interp(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
+      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 (ipj==1) THEN ! we assume that this is an un-structured open boundary file
+            DO ib = 1, ipi
+              DO ik = 1, ipk
+                dta(ib,1,ik) =  dta_read(map(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(ib)-1)/REAL(ilendta))
+               ji=map(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 (ipj==1) THEN
+ 
+   END SUBROUTINE fld_map
+   
+   SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, map, jpk_bdy, igrd, ibdy, fv, dta)
+
+      !!---------------------------------------------------------------------
+      !!                    ***  ROUTINE fld_bdy_interp  ***
+      !!
+      !! ** Purpose :   on the fly vertical interpolation to allow the use of
+      !!                boundary data from non-native vertical grid
+      !!----------------------------------------------------------------------
+#if defined key_bdy
+      USE bdy_oce, ONLY:  idx_bdy         ! indexing for map <-> ij transformation
+#endif 
+
+      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), DIMENSION(:,:,:), INTENT(out) ::   dta        ! output field on model grid (2 dimensional)
+      INTEGER,  DIMENSION(:)    , INTENT(in ) ::   map        ! global-to-local mapping indices
+      INTEGER  , INTENT(in)                   ::   igrd, ib_bdy, jpk_bdy      ! number of levels in bdy data
+      INTEGER                                 ::   jpkm1_bdy    ! number of levels in bdy data minus 1
+      !!
+      INTEGER                                 ::   ipi        ! length of boundary data on local process
+      INTEGER                                 ::   ipj        ! length of dummy dimension ( = 1 )
+      INTEGER                                 ::   ipk, ipkm1 ! number of vertical levels of dta ( 2D: ipk=1 ; 3D: ipk=jpk )
+      INTEGER                                 ::   ilendta    ! length of data in file
+      INTEGER                                 ::   ib, ik, ikk! loop counters
+      REAL(wp)                                ::   zl, zi     ! tmp variable for current depth and interpolation factor
+      REAL(wp)                                ::   fv, fv_alt ! fillvalue and alternative -ABS(fv)
+      !!---------------------------------------------------------------------
+
+      ipi       = SIZE( dta, 1 )
+      ipj       = SIZE( dta_read, 2 )
+      ipk       = SIZE( dta, 3 )
+      ilendta   = SIZE( dta_read, 1 )
+      jpkm1_bdy = jpk_bdy-1
+ 
+      fv_alt = -ABS(fv)  ! set _FillValue < 0 as we make use of MAXVAL and MAXLOC later
+      !
+      IF (ipj==1) THEN ! we assume that this is an un-structured open boundary file
          DO ib = 1, ipi
-            DO ik = 1, ipk
-               dta(ib,1,ik) =  dta_read(map(ib),1,ik)
+            DO ik = 1, ipk                      
+               IF( ( dta_read(map(ib),1,ik) == fv ) ) THEN
+                  dta_read_z(map(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data
+               ENDIF
+                  dta(ib,1,ik) = fv_alt    ! put fillvalue into new field as if all goes well all wet points will be replaced
+            ENDDO
+         ENDDO 
+      !
+         DO ib = 1, ipi
+            DO ik = 1, ipk                      
+               zl =  gdept_1(idx_bdy(ib_bdy)%nbi(ib,igrd),idx_bdy(ib_bdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_1?
+               IF( zl < dta_read_z(map(ib),1,1) ) THEN                                         ! above the first level of external data
+                  dta(ib,1,ik) =  dta_read(map(ib),1,1)
+               ELSEIF( zl > MAXVAL(dta_read_z(map(ib),1,:),1) ) THEN                           ! below the last level of external data 
+                  dta(ib,1,ik) =  dta_read(map(ib),1,MAXLOC(dta_read_z(map(ib),1,:),1))
+               ELSE                                                                          ! inbetween : vertical interpolation between ikk & ikk+1
+                  DO ikk = 1, jpkm1_bdy                                                          ! when  gdept_1(ikk) < zl < gdept_1(ikk+1)
+                     IF( ( (zl-dta_read_z(map(ib),1,ikk)) * (zl-dta_read_z(map(ib),1,ikk+1)) <= 0._wp)   &
+                    &    .AND. (dta_read_z(map(ib),1,ikk+1) /= fv_alt)) THEN
+                        zi = ( zl - dta_read_z(map(ib),1,ikk) ) / (dta_read_z(map(ib),1,ikk+1)-dta_read_z(map(ib),1,ikk))
+                        dta(ib,1,ik) = dta_read(map(ib),1,ikk) + &
+                      &                ( dta_read(map(ib),1,ikk+1) -  dta_read(map(ib),1,ikk) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
          END DO
@@ -714,12 +845,36 @@
             jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
             ji=map(ib)-(jj-1)*ilendta
-            DO ik = 1, ipk
-               dta(ib,1,ik) =  dta_read(ji,jj,ik)
+            DO ik = 1, ipk                      
+               IF( ( dta_read(ji,jj,ik) == fv ) ) THEN
+                  dta_read_z(map(ib),1,ik) = fv_alt ! safety: put fillvalue into external depth field so consistent with data
+               ENDIF
+                  dta(ib,1,ik) = fv_alt    ! put fillvalue into new field as if all goes well all wet points will be replaced
+            ENDDO
+         ENDDO 
+      !
+         DO ib = 1, ipi
+            jj=1+floor(REAL(map(ib)-1)/REAL(ilendta))
+            ji=map(ib)-(jj-1)*ilendta
+            DO ik = 1, ipk                      
+               zl =  gdept_1(idx_bdy(ib_bdy)%nbi(ib,igrd),idx_bdy(ib_bdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_1?
+               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,1)
+               ELSEIF( zl > MAXVAL(dta_read_z(ji,ji,:),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_1(ikk) < zl < gdept_1(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,1,ikk+1) -  dta_read(ji,jj,ikk) ) * zi
+                     ENDIF
+                  END DO
+               ENDIF
             END DO
          END DO
       ENDIF
 
-   END SUBROUTINE fld_map
-
+   END SUBROUTINE fld_bdy_interp
 
    SUBROUTINE fld_rot( kt, sd )
@@ -905,4 +1060,6 @@
             &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.')
          sdf(jf)%nreclast = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn
+         sdf(jf)%igrd     = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init
+         sdf(jf)%ibdy     = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during bdy_dta_init
       END DO
 
