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 5704)
+++ /branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 5705)
@@ -106,4 +106,5 @@
    INTEGER                            ::   nxt_wgt = 1  ! point to next available space in ref_wgts array
    REAL(wp), PARAMETER                ::   undeff_lsm = -999.00_wp
+#  include "domzgr_substitute.h90"
 
 !$AGRIF_END_DO_NOT_TREAT
@@ -839,22 +840,23 @@
       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), 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                                 ::   jpkm1_bdy    ! number of levels in bdy data minus 1
-      REAL(wp) , INTENT(in)                                ::   fv ! fillvalue and alternative -ABS(fv)
-      !!
-      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, INTENT(in)                                 ::   ilendta    ! length of data in file
-      INTEGER                                 ::   ib, ik, ikk! loop counters
-      INTEGER                                 ::   ji, jj ! loop counters
-      REAL(wp)                                ::   zl, zi, zh, zz, zdz    ! tmp variable for current depth and interpolation factor
+      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
       !!---------------------------------------------------------------------
      
@@ -868,11 +870,11 @@
       
       !
-      IF ( map%ll_unstruc ) THEN ! unstructured open boundary data file
+      IF ( map%ll_unstruc ) THEN                            ! unstructured open boundary data file
 
          DO ib = 1, ipi
             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
+                  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
@@ -880,17 +882,36 @@
 
          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,:) )
             DO ik = 1, ipk                      
-               zl =  gdept_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_0?
+               SELECT CASE( igrd )                       
+                  CASE(1)
+                     zl =  gdept_0(zij,zjj,ik)                                          ! if using in step could use fsdept instead of gdept_0?
+                  CASE(2)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0?
+                     ELSE
+                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(zij+1,zjj,ik) ) 
+                     ENDIF
+                  CASE(3)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0?
+                     ELSE
+                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(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_0(ikk) < zl < gdept_0(ikk+1)
-                     IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp)   &
+               ELSE                                                                                ! inbetween : vertical interpolation between ikk & ikk+1
+                  DO ikk = 1, jpkm1_bdy                                                            ! when  gdept_0(ikk) < zl < gdept_0(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))
+                        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
+                       &               ( dta_read(map%ptr(ib),1,ikk+1) - dta_read(map%ptr(ib),1,ikk) ) * zi
                      ENDIF
                   END DO
@@ -899,23 +920,24 @@
          END DO
 
-         IF(igrd == 2) THEN ! do we need to adjust the transport term?
+         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
             DO ib = 1, ipi
-              zh = SUM(dta_read_dz(map%ptr(ib),1,:) )
+              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
-                  ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
+              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
-                  zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik) 
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
+              DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_0(zij,zjj,ik) * umask(zij,zjj,ik)
               ENDDO
-              DO ik = 1, ipk
-                 zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   
-                 zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))       ! Is array HUR the correct one to use????
+              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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) / hu_0(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 ) * ( zdz * zz )
+                    IF( ABS(ztrans / hu_0(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 ) / hu_0(zij,zjj) * umask(zij,zjj,ik)
                  ENDIF
               ENDDO
@@ -923,23 +945,22 @@
          ENDIF
 
-         IF(igrd == 3) THEN ! do we need to adjust the transport term?
+         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
             DO ib = 1, ipi
-              zh = SUM(dta_read_dz(map%ptr(ib),1,:) )
+              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
-                  ztrans = ztrans + dta_read(map%ptr(ib),1,ik) * dta_read_dz(map%ptr(ib),1,ik)
+              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
-                  zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik) 
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
+              DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_0(zij,zjj,ik) * vmask(zij,zjj,ik)
               ENDDO
-              DO ik = 1, ipk
-                 zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   
-                 zz  =  hvr(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))   
+              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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) / hv_0(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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) / hv_0(zij,zjj) * vmask(zij,zjj,ik)
                  ENDIF
               ENDDO
@@ -947,8 +968,4 @@
          ENDIF
   
-         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-         ! At this point write out a single velocity profile/dz/H for model and input data             !
-         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
       ELSE ! structured open boundary file
 
@@ -958,6 +975,6 @@
             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
+                  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
@@ -968,17 +985,36 @@
             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,:) )
             DO ik = 1, ipk                      
-               zl =  gdept_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)   ! if using in step could use fsdept instead of gdept_0?
-               IF( zl < dta_read_z(ji,jj,1) ) THEN                                         ! above the first level of external data
+               SELECT CASE( igrd )                          ! coded for sco - need zco and zps option using min
+                  CASE(1)
+                     zl =  gdept_0(zij,zjj,ik)              ! if using in step could use fsdept instead of gdept_0?
+                  CASE(2)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij+1,zjj,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0?
+                     ELSE
+                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(zij+1,zjj,ik) )
+                     ENDIF
+                  CASE(3)
+                     IF(ln_sco) THEN
+                        zl =  ( gdept_0(zij,zjj,ik) + gdept_0(zij,zjj+1,ik) ) * 0.5_wp  ! if using in step could use fsdept instead of gdept_0?
+                     ELSE
+                        zl =  MIN( gdept_0(zij,zjj,ik), gdept_0(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 
+               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_0(ikk) < zl < gdept_0(ikk+1)
-                     IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp)   &
+               ELSE                                                                     ! inbetween : vertical interpolation between ikk & ikk+1
+                  DO ikk = 1, jpkm1_bdy                                                 ! when  gdept_0(ikk) < zl < gdept_0(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))
+                        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
+                       &               ( dta_read(ji,jj,ikk+1) - dta_read(ji,jj,ikk) ) * zi
                      ENDIF
                   END DO
@@ -987,25 +1023,24 @@
          END DO
 
-         IF(igrd == 2) THEN ! do we need to adjust the transport term?
+         IF(igrd == 2) THEN                                 ! do we need to adjust the transport term?
             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)
             DO ib = 1, ipi
               zh = SUM(dta_read_dz(ji,jj,:) )
               ztrans = 0._wp
               ztrans_new = 0._wp
-              DO ik = 1, jpk_bdy
+              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
-                  zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
+              DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3u_0(zij,zjj,ik) * umask(zij,zjj,ik)
               ENDDO
-              DO ik = 1, ipk
-                 zdz =  e3u_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
-                 zz  =  hur(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))
+              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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hur(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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hur(zij,zjj) ) * umask(zij,zjj,ik)
                  ENDIF
               ENDDO
@@ -1013,25 +1048,24 @@
          ENDIF
 
-         IF(igrd == 3) THEN ! do we need to adjust the transport term?
-            jj=1+floor(REAL(map%ptr(ib)-1)/REAL(ilendta))
-            ji=map%ptr(ib)-(jj-1)*ilendta
+         IF(igrd == 3) THEN                                 ! do we need to adjust the transport term?
            DO ib = 1, ipi
-              zh = SUM(dta_read_dz(ji,jj,:) )
+              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
-                  ztrans = ztrans + dta_read(ji,jj,ik) * dta_read_dz(ji,jj,ik)
+              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
-                  zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
-                  ztrans_new = ztrans_new + dta(ib,1,ik) * zdz
+              DO ik = 1, ipk                                ! calculate transport on model grid
+                  ztrans_new = ztrans_new + dta(ib,1,ik) * e3v_0(zij,zjj,ik) * vmask(zij,zjj,ik)
               ENDDO
-              DO ik = 1, ipk
-                 zdz =  e3v_0(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd),ik)
-                 zz  =  hvr(idx_bdy(ibdy)%nbi(ib,igrd),idx_bdy(ibdy)%nbj(ib,igrd))
+              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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * hvr(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 ) * ( zdz * zz )
+                    dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp  - ztrans_new ) * hvr(zij,zjj) ) * vmask(zij,zjj,ik)
                  ENDIF
               ENDDO
