Index: /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
===================================================================
--- /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8134)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90	(revision 8135)
@@ -17,4 +17,5 @@
 
    PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90
+   PUBLIC reconstructandremap
 
    !                                              !!* Namelist namagrif: AGRIF parameters
@@ -104,4 +105,225 @@
    END FUNCTION agrif_oce_alloc
 
+      subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
+      implicit none
+      integer N, Nout
+      real tabin(N), tabout(Nout)
+      real hin(N), hout(Nout)
+      real coeffremap(N,3),zwork(N,3)
+      real zwork2(N+1,3)
+      integer k
+      double precision, parameter :: dsmll=1.0d-8  
+      real q,q01,q02,q001,q002,q0
+      real z_win(1:N+1), z_wout(1:Nout+1)
+      real,parameter :: dpthin = 1.D-3
+      integer :: k1, kbox, ktop, ka, kbot
+      real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+
+      z_win(1)=0.; z_wout(1)= 0.
+      do k=1,N
+       z_win(k+1)=z_win(k)+hin(k)
+      enddo 
+      
+      do k=1,Nout
+       z_wout(k+1)=z_wout(k)+hout(k)       
+      enddo       
+
+        do k=2,N
+          zwork(k,1)=1./(hin(k-1)+hin(k))
+        enddo
+        
+        do k=2,N-1
+          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
+          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
+          zwork(k,3)=q0
+        enddo       
+      
+        do k= 2,N
+        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
+        enddo
+        
+        coeffremap(:,1) = tabin(:)
+ 
+         do k=2,N-1
+        q001 = hin(k)*zwork2(k+1,1)
+        q002 = hin(k)*zwork2(k,1)        
+        if (q001*q002 < 0) then
+          q001 = 0.
+          q002 = 0.
+        endif
+        q=zwork(k,2)
+        q01=q*zwork2(k+1,1)
+        q02=q*zwork2(k,1)
+        if (abs(q001) > abs(q02)) q001 = q02
+        if (abs(q002) > abs(q01)) q002 = q01
+        
+        q=(q001-q002)*zwork(k,3)
+        q001=q001-q*hin(k+1)
+        q002=q002+q*hin(k-1)
+        
+        coeffremap(k,3)=coeffremap(k,1)+q001
+        coeffremap(k,2)=coeffremap(k,1)-q002
+        
+        zwork2(k,1)=(2.*q001-q002)**2
+        zwork2(k,2)=(2.*q002-q001)**2
+        enddo
+        
+        do k=1,N
+        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
+        coeffremap(k,3) = coeffremap(k,1)
+        coeffremap(k,2) = coeffremap(k,1)
+        zwork2(k,1) = 0.
+        zwork2(k,2) = 0.
+        endif
+        enddo
+        
+        do k=2,N
+        q002=max(zwork2(k-1,2),dsmll)
+        q001=max(zwork2(k,1),dsmll)
+        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002)
+        enddo
+        
+        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
+        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
+ 
+        do k=1,N
+        q01=zwork2(k+1,3)-coeffremap(k,1)
+        q02=coeffremap(k,1)-zwork2(k,3)
+        q001=2.*q01
+        q002=2.*q02
+        if (q01*q02<0) then
+          q01=0.
+          q02=0.
+        elseif (abs(q01)>abs(q002)) then
+          q01=q002
+        elseif (abs(q02)>abs(q001)) then
+          q02=q001
+        endif
+        coeffremap(k,2)=coeffremap(k,1)-q02
+        coeffremap(k,3)=coeffremap(k,1)+q01
+        enddo
+
+      zbot=0.0
+      kbot=1
+      do k=1,Nout
+        ztop=zbot  !top is bottom of previous layer
+        ktop=kbot
+        if     (ztop.ge.z_win(ktop+1)) then
+          ktop=ktop+1
+        endif
+        
+        zbot=z_wout(k+1)
+        zthk=zbot-ztop
+
+        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
+
+          kbot=ktop
+          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
+            kbot=kbot+1
+          enddo
+          zbox=zbot
+          do k1= k+1,Nout
+            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
+              exit !thick layer
+            else
+              zbox=z_wout(k1+1)  !include thin adjacent layers
+              if     (zbox.eq.z_wout(Nout+1)) then
+                exit !at bottom
+              endif
+            endif
+          enddo
+          zthk=zbox-ztop
+
+          kbox=ktop
+          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
+            kbox=kbox+1
+          enddo
+          
+          if     (ktop.eq.kbox) then
+
+
+            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then
+
+              if     (hin(kbox).gt.dpthin) then
+                q001 = (zbox-z_win(kbox))/hin(kbox)
+                q002 = (ztop-z_win(kbox))/hin(kbox)
+                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
+                q02=q01-1.+(q001+q002)
+                q0=1.-q01-q02
+              else
+                q0 = 1.0
+                q01 = 0.
+                q02 = 0.
+              endif
+          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+              
+            else
+            tabout(k) = tabin(kbox)
+              
+            endif 
+
+          else
+
+            if     (ktop.le.k .and. kbox.ge.k) then
+              ka = k
+            elseif (kbox-ktop.ge.3) then
+              ka = (kbox+ktop)/2
+            elseif (hin(ktop).ge.hin(kbox)) then
+              ka = ktop
+            else
+              ka = kbox
+            endif !choose ka
+
+            offset=coeffremap(ka,1)
+
+            qtop = z_win(ktop+1)-ztop !partial layer thickness
+            if     (hin(ktop).gt.dpthin) then
+              q=(ztop-z_win(ktop))/hin(ktop)
+              q01=q*(q-1.)
+              q02=q01+q
+              q0=1-q01-q02            
+            else
+              q0 = 1.
+              q01 = 0.
+              q02 = 0.
+            endif
+            
+            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
+
+            do k1= ktop+1,kbox-1
+              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
+            enddo !k1
+
+            
+            qbot = zbox-z_win(kbox) !partial layer thickness
+            if     (hin(kbox).gt.dpthin) then
+              q=qbot/hin(kbox)
+              q01=(q-1.)**2
+              q02=q01-1.+q
+              q0=1-q01-q02                            
+            else
+              q0 = 1.0
+              q01 = 0.
+              q02 = 0.
+            endif
+           
+            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
+            
+            rpsum=1.0d0/zthk
+              tabout(k)=offset+tsum*rpsum
+              
+          endif !single or multiple layers
+        else
+        if (k==1) then
+        write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1)
+        endif
+         tabout(k) = tabout(k-1)
+          
+        endif !normal:thin layer
+      enddo !k
+            
+      return
+      end subroutine reconstructandremap
+      
 #endif
    !!======================================================================
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 8134)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8135)
@@ -576,4 +576,6 @@
       !!----------------------------------------------------------------------  
       !
+      return
+      
       zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
       IF( zalpha > 1. )   zalpha = 1.
@@ -603,10 +605,80 @@
       REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
       REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
-      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------
-      !
+      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
+      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER :: N_in, N_out
+      REAL(wp) :: h_diff
+      REAL(wp) :: zrhoxy
+! VERTICAL REFINEMENT END
+
+      zrhoxy = Agrif_rhox()*Agrif_rhoy()
       IF (before) THEN         
-         ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
-      ELSE
+         DO jn = n1,n2-1
+            DO jk=k1,k2
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
+                  END DO
+               END DO
+            END DO
+         END DO
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+
+      ELSE 
+! VERTICAL REFINEMENT BEGIN
+
+         ptab_child(:,:,:,:) = 0.
+         do jj=j1,j2
+         do ji=i1,i2
+           N_in = 0
+           DO jk=k1,k2 !k2 = jpk of parent grid
+             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)*zrhoxy)
+           END DO
+           N_out = 0
+           DO jk=1,jpk ! jpk of child grid
+             IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model
+             N_out = N_out + 1
+             h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
+           ENDDO
+           IF (N_in > 0) THEN
+             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+!             if (h_diff > 0) then
+!               h_in(N_in+1) = h_diff
+!               N_in = N_in + 1
+!             else
+!               h_out(N_out+1) = -h_diff
+!               N_out = N_out + 1
+!             endif 
+           ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for?????
+           do jn=1,jpts
+             call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
+           enddo
+!           if (abs(h_diff) > 1000.) then
+!           do jn=1,jpts
+!             do jk=1,N_out
+!             print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn)
+!             enddo
+!           enddo
+!         endif
+           ENDIF
+         enddo
+         enddo
+   
+
+! VERTICAL REFINEMENT END
+
          !
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -640,5 +712,5 @@
          IF( eastern_side ) THEN
             DO jn = 1, jpts
-               tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
+               tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
@@ -660,5 +732,5 @@
          IF( northern_side ) THEN            
             DO jn = 1, jpts
-               tsa(i1:i2,nlcj,k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)
+               tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO ji = imin,imax
@@ -680,5 +752,5 @@
          IF( western_side ) THEN            
             DO jn = 1, jpts
-               tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
+               tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
@@ -699,5 +771,5 @@
          IF( southern_side ) THEN           
             DO jn = 1, jpts
-               tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
+               tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn)
                DO jk = 1, jpk      
                   DO ji=imin,imax
@@ -720,17 +792,17 @@
          ! East south
          IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)
+            tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts)
          ENDIF
          ! East north
          IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)
+            tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts)
          ENDIF
          ! West south
          IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
-            tsa(2,2,:,:) = ptab(2,2,:,:)
+            tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts)
          ENDIF
          ! West north
          IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
-            tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)
+            tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts)
          ENDIF
          !
@@ -771,22 +843,104 @@
       !!----------------------------------------------------------------------
       !!   *** ROUTINE interpun ***
-      !!----------------------------------------------------------------------
-      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
-      LOGICAL                               , INTENT(in   ) ::   before
-      !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhoy	
-      !!----------------------------------------------------------------------
-      !
-      IF( before ) THEN 
-         DO jk = k1, jpk
-            ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
+      !!---------------------------------------------    
+      !!
+      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
+      LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb , ndir
+      !!
+      INTEGER :: ji,jj,jk
+      REAL(wp) :: zrhoy
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
+      REAL(wp), DIMENSION(k1:k2) :: tabin
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER :: N_in, N_out
+      REAL(wp) :: h_diff
+      LOGICAL :: western_side, eastern_side
+      INTEGER :: iref
+
+! VERTICAL REFINEMENT END
+      !!---------------------------------------------    
+      !
+      zrhoy = Agrif_rhoy()
+      IF (before) THEN 
+         !We can't use zero as the special value because we need to include zeros
+         !when interpolating the scale factors
+         IF(Agrif_UseSpecialValue) THEN 
+             Agrif_SpecialValue = -999._wp
+         ELSE
+             Agrif_SpecialValue = 0._wp
+         ENDIF
+         DO jk=1,jpk
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - &
+                                   & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue)
+                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
+               END DO
+            END DO
          END DO
       ELSE
-         zrhoy = Agrif_Rhoy()
+! VERTICAL REFINEMENT BEGIN
+         western_side  = (nb == 1).AND.(ndir == 1)
+         eastern_side  = (nb == 1).AND.(ndir == 2)
+
+         Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed
+
+         ptab_child(:,:,:) = 0.
+         DO jj=j1,j2
+            DO ji=i1,i2
+               iref = ji
+               IF (western_side) iref = 2
+               IF (eastern_side) iref = nlci-2
+
+               N_in = 0
+               DO jk=k1,k2
+                  IF (ptab(ji,jj,jk,2) == 0) EXIT
+                  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) * zrhoy) 
+              ENDDO
+         
+              IF (N_in == 0) THEN
+                 ptab_child(ji,jj,:) = 0.
+                 CYCLE
+              ENDIF
+         
+              N_out = 0
+              DO jk=1,jpk
+                 if (umask(ji,jj,jk) == 0) EXIT
+                 N_out = N_out + 1
+                 h_out(N_out) = e3u_n(ji,jj,jk)
+              ENDDO
+         
+              IF (N_out == 0) THEN
+                 ptab_child(ji,jj,:) = 0.
+                 CYCLE
+              ENDIF
+         
+!              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 < 0.) then
+!                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+!                    stop
+!                 endif
+!              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
+
+! in the following
+! remove division of ua by fs e3u (already done) and also zrhoy and e2u
+! VERTICAL REFINEMENT END
+
          DO jk = 1, jpkm1
             DO jj=j1,j2
-               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) )
+               ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
+!/(zrhoy*e2u(i1:i2,jj)))
             END DO
          END DO
@@ -800,20 +954,96 @@
       !!   *** ROUTINE interpvn ***
       !!----------------------------------------------------------------------
-      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
-      LOGICAL                               , INTENT(in   ) ::   before
-      !
-      INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zrhox	
-      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
+      LOGICAL, INTENT(in) :: before
+      INTEGER, INTENT(in) :: nb , ndir
+      !
+      INTEGER :: ji,jj,jk
+      REAL(wp) :: zrhox
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
+      REAL(wp), DIMENSION(k1:k2) :: tabin
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER :: N_in, N_out
+      REAL(wp) :: h_diff
+      LOGICAL :: northern_side,southern_side
+      INTEGER :: jref
+
+! VERTICAL REFINEMENT END
+      !!---------------------------------------------    
       !      
-      IF( before ) THEN       !interpv entre 1 et k2 et interpv2d en jpkp1
-         DO jk = k1, jpk
-            ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk)
-         END DO
-      ELSE          
-         zrhox= Agrif_Rhox()
-         DO jk = 1, jpkm1
-            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) )
+      zrhox = Agrif_rhox()
+      IF (before) THEN          
+         IF(Agrif_UseSpecialValue) THEN 
+             Agrif_SpecialValue = -999._wp
+         ELSE
+             Agrif_SpecialValue = 0._wp
+         ENDIF
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE        
+! VERTICAL REFINEMENT BEGIN
+         ptab_child(:,:,:) = 0.
+         southern_side = (nb == 2).AND.(ndir == 1)
+         northern_side = (nb == 2).AND.(ndir == 2)
+
+         Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done
+
+         do jj=j1,j2
+            jref = jj
+            IF (southern_side) jref = 2
+            IF (northern_side) jref = nlcj-2
+            do ji=i1,i2
+
+               N_in = 0
+               do jk=k1,k2
+                  if (ptab(ji,jj,jk,2) == 0) EXIT
+                  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)*zrhox)
+               enddo
+               IF (N_in == 0) THEN
+                  ptab_child(ji,jj,:) = 0.
+                  CYCLE
+               ENDIF
+         
+               N_out = 0
+               do jk=1,jpk
+                  if (vmask(ji,jref,jk) == 0) EXIT
+                  N_out = N_out + 1
+                  h_out(N_out) = e3v_n(ji,jj,jk)
+               enddo
+               IF (N_out == 0) THEN
+                 ptab_child(ji,jj,:) = 0.
+                 CYCLE
+               ENDIF
+
+!               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 < 0.) then
+!                     print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
+!                     stop
+!                  endif
+!               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
+! in the following
+! remove division of va by fs e3v, zrhox and e1v (already done)
+! VERTICAL REFINEMENT END
+         DO jk=1,jpkm1
+            DO jj=j1,j2
+               va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
+            END DO
          END DO
       ENDIF
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 8134)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8135)
@@ -158,4 +158,6 @@
       INTEGER, INTENT(in) :: kt
       !       
+      return
+      
       IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN
 #  if defined TWO_WAY
@@ -185,26 +187,89 @@
       !!
       INTEGER :: ji,jj,jk,jn
-      !!---------------------------------------------
-      !
-      IF (before) THEN
-         DO jn = n1,n2
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER  :: N_in, N_out
+      REAL(wp) :: h_diff
+      REAL(wp) :: zrho_xy
+      REAL(wp) :: tabin(k1:k2,n1:n2)
+! VERTICAL REFINEMENT END
+      !!---------------------------------------------
+      !
+      IF (before) THEN
+         zrho_xy = Agrif_rhox() * Agrif_rhoy() 
+         DO jn = n1,n2-1
             DO jk=k1,k2
                DO jj=j1,j2
                   DO ji=i1,i2
-                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)
+                     tabres(ji,jj,jk,jn) = zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
                   END DO
                END DO
             END DO
          END DO
-      ELSE
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * zrho_xy * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+         
+      ELSE
+! VERTICAL REFINEMENT BEGIN
+         tabres_child(:,:,:,:) = 0.
+
+         DO jj=j1,j2
+         DO ji=i1,i2
+           N_in = 0
+           DO jk=k1,k2 !k2 = jpk of child grid
+             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)
+             h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj)
+           ENDDO
+           N_out = 0
+           DO jk=1,jpk ! jpk of parent grid
+             IF (tmask(ji,jj,jk) == 0) 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 (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 < 0) THEN
+             print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out))
+             print *,'Nval = ',N_out,mbathy(ji,jj)
+             print *,'BATHY = ',gdepw_0(ji,jj,mbathy(ji,jj)+1),sum(e3t_0(ji,jj,1:mbathy(ji,jj)))
+             STOP
+!               N_out = N_out + 1
+!               h_out(N_out) = - h_diff
+             ENDIF
+             DO jn=n1,n2-1
+               CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
+             ENDDO
+          ENDIF
+         ENDDO
+         ENDDO
+         
+! 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
-            DO jn = n1,n2
-               DO jk=k1,k2
+            DO jn = n1,n2-1
+               DO jk=1,jpk
                   DO jj=j1,j2
                      DO ji=i1,i2
-                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN
+                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN
                            tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 
-                                 & + atfp * ( tabres(ji,jj,jk,jn) &
+                                 & + atfp * ( tabres_child(ji,jj,jk,jn) &
                                  &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk)
                         ENDIF
@@ -214,10 +279,10 @@
             ENDDO
          ENDIF
-         DO jn = n1,n2
-            DO jk=k1,k2
+         DO jn = n1,n2-1
+            DO jk=1,jpk
                DO jj=j1,j2
                   DO ji=i1,i2
-                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
-                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk)
+                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
+                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)
                      END IF
                   END DO
@@ -230,33 +295,87 @@
 
 
-   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )
+   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!---------------------------------------------
       !!           *** ROUTINE updateu ***
       !!---------------------------------------------
-      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
-      LOGICAL                               , INTENT(in   ) :: before
+      INTEGER                                 , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: tabres
+      LOGICAL                                 , INTENT(in   ) :: before
       !
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zrhoy
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER :: N_in, N_out
+      REAL(wp) :: h_diff
+      REAL(wp) :: tabin(k1:k2)
+! VERTICAL REFINEMENT END
       !!---------------------------------------------
       ! 
       IF( before ) THEN
          zrhoy = Agrif_Rhoy()
-         DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
-         END DO
-      ELSE
          DO jk=k1,k2
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)
+                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)
+                  tabres(ji,jj,jk,2) = umask(ji,jj,jk) * zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) 
+               END DO
+            END DO
+         END DO
+      ELSE
+! VERTICAL REFINEMENT BEGIN
+         tabres_child(:,:,:) = 0.
+         
+         DO jj=j1,j2
+         DO ji=i1,i2
+           N_in = 0
+           DO jk=k1,k2 !k2=jpk of child grid
+             IF (tabres(ji,jj,jk,2) == 0) EXIT
+             N_in = N_in + 1
+             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+             h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
+           ENDDO
+           N_out = 0
+           DO jk=1,jpk
+             IF (umask(ji,jj,jk) == 0) EXIT
+             N_out = N_out + 1
+             h_out(N_out) = e3u_n(ji,jj,jk)
+           ENDDO
+           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 < 0.) 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
+             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
+         ENDDO
+         ENDDO
+         
+! 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) & 
-                           & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
+                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
                   ENDIF
                   !
-                  un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk)
+                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
@@ -271,10 +390,18 @@
       !!           *** ROUTINE updatev ***
       !!---------------------------------------------
-      INTEGER :: i1,i2,j1,j2,k1,k2
+      INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2
       INTEGER :: ji,jj,jk
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: tabres
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: tabres
       LOGICAL :: before
       !!
       REAL(wp) :: zrhox
+! VERTICAL REFINEMENT BEGIN
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      INTEGER :: N_in, N_out
+      REAL(wp) :: h_diff
+      REAL(wp) :: tabin(k1:k2)
+! VERTICAL REFINEMENT END
       !!---------------------------------------------      
       !
@@ -284,20 +411,62 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ELSE
-         DO jk=k1,k2
+                  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) 
+               END DO
+            END DO
+         END DO
+      ELSE
+! VERTICAL REFINEMENT BEGIN
+         tabres_child(:,:,:) = 0.
+         
+         DO jj=j1,j2
+         DO ji=i1,i2
+           N_in = 0
+           DO jk=k1,k2
+             IF (tabres(ji,jj,jk,2) == 0) EXIT
+             N_in = N_in + 1
+             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+             h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
+           ENDDO
+           N_out = 0
+           DO jk=1,jpk
+             IF (vmask(ji,jj,jk) == 0) EXIT
+             N_out = N_out + 1
+             h_out(N_out) = e3v_n(ji,jj,jk)
+           ENDDO
+           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 < 0.) 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
+             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
+         ENDDO
+         ENDDO
+         
+! 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
-                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)
+!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
                      vb(ji,jj,jk) = vb(ji,jj,jk) & 
-                           & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
+                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk)
+                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
Index: /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
===================================================================
--- /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8134)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8135)
@@ -346,11 +346,11 @@
    ! 1. Declaration of the type of variable which have to be interpolated
    !---------------------------------------------------------------------
-   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id)
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id)
    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id)
 
-   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id)
-   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id)
-   CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id)
-   CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id)
+   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id)
+   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id)
    CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id)
    CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id)
