Index: /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8821)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8822)
@@ -665,6 +665,6 @@
              IF (ptab(ji,jj,jk,n2) == 0) EXIT
              N_in = N_in + 1
-             tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2)
-             h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj))
+             tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2))
+             h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy)
            END DO
            N_out = 0
@@ -716,9 +716,4 @@
 
          !
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         !
          zrhox = Agrif_Rhox()
          ! 
@@ -939,5 +934,5 @@
                   N_in = N_in + 1
                   tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
-                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)) 
+                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 
               ENDDO
          
@@ -968,5 +963,4 @@
               ENDIF
               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
-         
             ENDDO
          ENDDO
@@ -1050,5 +1044,5 @@
                   N_in = N_in + 1
                   tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
-                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj))
+                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
                enddo
                IF (N_in == 0) THEN
Index: /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8821)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8822)
@@ -206,6 +206,6 @@
                   DO jj=j1,j2
                      DO ji=i1,i2
-                        tabres(ji,jj,jk,jn) = (zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &
-                                              * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp
+                        tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &
+                                              * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp
                      END DO
                   END DO
@@ -215,5 +215,5 @@
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  &
+                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy  &
                                               + (tmask(ji,jj,jk)-1)*999._wp
                   END DO
@@ -240,5 +240,5 @@
            N_in = 0
            DO jk=k1,k2 !k2 = jpk of child grid
-             IF (tabres(ji,jj,jk,n2) > -900) EXIT
+             IF (tabres(ji,jj,jk,n2) == 0  ) EXIT
              N_in = N_in + 1
              tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
@@ -247,23 +247,15 @@
            N_out = 0
            DO jk=1,jpk ! jpk of parent grid
-             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF
+             IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF
              N_out = N_out + 1
              h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above
            ENDDO
-!           IF(ji.EQ.i1 .AND. jj.EQ.j1) print *,'1st parent point',sum(h_in(1:N_in)), sum(h_out(1:N_out))
            IF (N_in > 0) THEN
              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
-!             IF abs(h_diff > 1.e-8) THEN
-!               N_in = N_in + 1
-!               h_in(N_in) = h_diff
-!               tabin(N_in,:) = tabin(N_in-1,:)
              IF (h_diff < -1.e-4) THEN
-             print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,e1e2t(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out
-             print *,h_in(1:N_in)
-             print *,h_out(1:N_out)
-             STOP
-!               N_out = N_out + 1
-!               h_out(N_out) = - h_diff
+                print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,e1e2t(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out
+                print *,h_in(1:N_in)
+                print *,h_out(1:N_out)
+                STOP
              ENDIF
              DO jn=n1,n2-1
@@ -276,10 +268,8 @@
             tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts)
 #endif
-         
 ! WARNING :
 ! tabres replaced by tabres_child in the following
 ! k1 k2 replaced by 1 jpk
 ! VERTICAL REFINEMENT END
-
          IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
             ! Add asselin part
@@ -329,5 +319,5 @@
       REAL(wp) :: h_out(1:jpk)
       INTEGER :: N_in, N_out
-      REAL(wp) :: h_diff
+      REAL(wp) :: h_diff, excess, thick
       REAL(wp) :: tabin(k1:k2)
 ! VERTICAL REFINEMENT END
@@ -343,7 +333,7 @@
                DO ji=i1,i2
                   tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  &
-                                       - (umask(ji,jj,jk)-1)*999._wp
-                  tabres(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
-                                       - (umask(ji,jj,jk)-1)*999._wp
+                                       + (umask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  &
+                                       + (umask(ji,jj,jk)-1)*999._wp
                END DO
             END DO
@@ -357,10 +347,12 @@
          tabres_child(:,:,:) = 0.
 # if defined key_vertical
-         AGRIF_SpecialValue = -999._wp
+         AGRIF_SpecialValue = 0._wp
          DO jj=j1,j2
          DO ji=i1,i2
            N_in = 0
+           h_in(:) = 0._wp
+           tabin(:) = 0._wp
            DO jk=k1,k2 !k2=jpk of child grid
-             IF (tabres(ji,jj,jk,2) > -900) EXIT
+             IF( tabres(ji,jj,jk,2) < -900) EXIT
              N_in = N_in + 1
              tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
@@ -375,17 +367,25 @@
            IF (N_in * N_out > 0) THEN
              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
-             if (h_diff < -1.e-4) then
-             print *,'CHECK YOUR bathy U points ...',ji,jj,h_diff,e2u(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out
-             stop
-!             else ! Extends with 0
-!             N_in = N_in + 1
-!             tabin(N_in) = 0.
-!             h_in(N_in) = h_diff
-             endif
+             IF (h_diff < -1.e-4) then
+!Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
+!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
+                excess = 0._wp
+                DO jk=N_in,1,-1
+                   thick = MIN(-1*h_diff, h_in(jk))
+                   excess = excess + tabin(jk)*thick*e2u(ji,jj)
+                   tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
+                   h_diff = h_diff + thick
+                   IF ( h_diff == 0) THEN
+                      N_in = jk
+                      h_in(jk) = h_in(jk) - thick
+                      EXIT
+                   ENDIF
+                ENDDO
+             ENDIF
              CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
-          ENDIF
+             tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
+           ENDIF
          ENDDO
-         ENDDO
+      ENDDO
 #else
             DO jk=1,jpk
@@ -397,17 +397,8 @@
             END DO
 #endif
-         
-! WARNING :
-! tabres replaced by tabres_child in the following
-! k1 k2 replaced by 1 jpk
-! remove division by fs e3u (already done)
-! VERTICAL REFINEMENT END
 
          DO jk=1,jpk
             DO jj=j1,j2
                DO ji=i1,i2
-!Following line now replaced by division higher up in vertical refinement case
-!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)
-                  !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      ub(ji,jj,jk) = ub(ji,jj,jk) & 
@@ -439,5 +430,5 @@
       REAL(wp) :: h_out(1:jpk)
       INTEGER :: N_in, N_out
-      REAL(wp) :: h_diff
+      REAL(wp) :: h_diff, excess, thick
       REAL(wp) :: tabin(k1:k2)
 ! VERTICAL REFINEMENT END
@@ -447,9 +438,12 @@
          zrhox = Agrif_Rhox()
 #if defined key_vertical
+         AGRIF_SpecialValue = -999._wp
          DO jk=k1,k2
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)
-                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) 
+                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) &
+                                       + (vmask(ji,jj,jk)-1)*999._wp
+                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) &
+                                       + (vmask(ji,jj,jk)-1)*999._wp
                END DO
             END DO
@@ -466,11 +460,11 @@
       ELSE
          tabres_child(:,:,:) = 0.
-! VERTICAL REFINEMENT BEGIN
 #if defined key_vertical
+         AGRIF_SpecialValue = 0._wp
          DO jj=j1,j2
          DO ji=i1,i2
            N_in = 0
            DO jk=k1,k2
-             IF (tabres(ji,jj,jk,2) == 0) EXIT
+             IF (tabres(ji,jj,jk,2) < -900) EXIT
              N_in = N_in + 1
              tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
@@ -485,14 +479,22 @@
            IF (N_in * N_out > 0) THEN
              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
-             if (h_diff < -1.e-4) then
-             print *,'CHECK YOUR BATHY ...'
-             stop
-!             else ! Extends with 0
-!             N_in = N_in + 1
-!             tabin(N_in) = 0.
-!             h_in(N_in) = h_diff
-             endif
+             IF (h_diff < -1.e-4) then
+!Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
+!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
+                excess = 0._wp
+                DO jk=N_in,1,-1
+                   thick = MIN(-1*h_diff, h_in(jk))
+                   excess = excess + tabin(jk)*thick*e2u(ji,jj)
+                   tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
+                   h_diff = h_diff + thick
+                   IF ( h_diff == 0) THEN
+                      N_in = jk
+                      h_in(jk) = h_in(jk) - thick
+                      EXIT
+                   ENDIF
+                ENDDO
+             ENDIF
              CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+             tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
           ENDIF
          ENDDO
@@ -507,16 +509,8 @@
             END DO
 #endif
-         
-! WARNING :
-! tabres replaced by tabres_child in the following
-! k1 k2 replaced by 1 jpk
-! remove division by fs e3v (already done)
-! VERTICAL REFINEMENT END
 
          DO jk=k1,jpk
             DO jj=j1,j2
                DO ji=i1,i2
-!Following line now replaced by division higher up in vertical refinement case
-!                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)
                   !
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
@@ -761,5 +755,4 @@
          END DO
       ENDIF
-      !
    END SUBROUTINE updatevb2b
 
