Index: /branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- /branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 7823)
+++ /branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 7824)
@@ -714,12 +714,31 @@
 ! 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)
+         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
@@ -728,19 +747,29 @@
          do jj=j1,j2
          do ji=i1,i2
-           h_in(k1:k2) = interp_scales_t(ji,jj,k1:k2)
-           h_out(1:jpk) = e3t_n(ji,jj,1:jpk)
-           h_diff = sum(h_out(1:jpk-1))-sum(h_in(k1:k2-1))
-           N_in = k2-1
+           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 = jpk-1
-           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,:)
+           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. Is it just a GYRE issue???????
+!             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(ptab(ji,jj,1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
+             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
@@ -751,6 +780,8 @@
 !           enddo
 !         endif
+           ENDIF
          enddo
          enddo
+   
 
 ! VERTICAL REFINEMENT END
@@ -950,6 +981,6 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = e2u(ji,jj) * un(ji,jj,jk) * e3u_n(ji,jj,jk)/zrhoy
-                  ptab(ji,jj,jk,2) = e2u(ji,jj) * e3u_n(ji,jj,jk)/zrhoy
+                  ptab(ji,jj,jk,1) = e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)
+                  ptab(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)
                END DO
             END DO
@@ -972,5 +1003,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
          
@@ -992,8 +1023,12 @@
               ENDIF
          
-              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-              IF (abs(h_diff) > 0.) THEN
-                 CALL ctl_stop
-              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)
          
@@ -1009,5 +1044,4 @@
                ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
 !/(zrhoy*e2u(i1:i2,jj)))
-            write(numout,*) 'ua=',ua(i1:i2,jj,jk)
             END DO
          END DO
@@ -1043,12 +1077,11 @@
       !!---------------------------------------------    
       !      
+      zrhox = Agrif_rhox()
       IF (before) THEN          
-         !interpv entre 1 et k2 et interpv2d en jpkp1
-         DO jk=k1,jpk
+         DO jk=k1,k2
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = e1v(ji,jj) * vn(ji,jj,jk)
-                  ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * e3v_n(ji,jj,jk)
-                  ptab(ji,jj,jk,2) = e3v_n(ji,jj,jk)
+                  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
@@ -1060,54 +1093,50 @@
          northern_side = (nb == 2).AND.(ndir == 2)
          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 (interp_scales_v(ji,jj,jk) == 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)
-         enddo
-         IF (N_in == 0) THEN
-           ptab_child(ji,jj,:) = 0.
-           CYCLE
-         ENDIF
+            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
+               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)
          
-         h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-         IF (h_diff > 0.) THEN
-           N_in = N_in + 1
-           h_in(N_in) = h_diff
-           tabin(N_in) = 0.
-         ELSE
-           h_out(N_out) = h_out(N_out) - h_diff
-         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)
-         
-         ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / e3v_n(ji,jj,N_out)
-
-         enddo
+            enddo
          enddo
 ! in the following
-! remove division of va by fs e3v (already done)
+! remove division of va by fs e3v, zrhox and e1v (already done)
 ! VERTICAL REFINEMENT END
-         zrhox= Agrif_Rhox()
          DO jk=1,jpkm1
             DO jj=j1,j2
-               va(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj)))
+               va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
             END DO
          END DO
Index: /branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 7823)
+++ /branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 7824)
@@ -1,3 +1,3 @@
-#undef TWO_WAY        /* TWO WAY NESTING */
+#define TWO_WAY        /* TWO WAY NESTING */
 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
  
@@ -121,6 +121,4 @@
       ! 
       IF (MOD(nbcline,nbclineupdate) == 0) THEN
-         WRITE(numout,*) 'TG print 1'
-         CALL FLUSH(numout)
 # if ! defined DECAL_FEEDBACK
          CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
@@ -128,9 +126,5 @@
          CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
 # endif
-         WRITE(numout,*) 'TG print 2'
-         CALL FLUSH(numout)
-      ELSE
-         WRITE(numout,*) 'TG print 3'
-         CALL FLUSH(numout)
+      ELSE
 # if ! defined DECAL_FEEDBACK
          CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
@@ -138,6 +132,4 @@
          CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
 # endif
-         WRITE(numout,*) 'TG print 4'
-         CALL FLUSH(numout)
       ENDIF
       !
@@ -380,6 +372,4 @@
       ENDIF
       ! 
-      WRITE(numout,*) 'I got to end of updateTS before=',before
-      CALL FLUSH(numout)
    END SUBROUTINE updateTS
 
@@ -396,5 +386,5 @@
       REAL(wp) ::   zrhoy
 ! VERTICAL REFINEMENT BEGIN
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,2) :: ptab_child
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
       REAL(wp) :: h_in(k1:k2)
       REAL(wp) :: h_out(1:jpk)
@@ -405,6 +395,4 @@
       !!---------------------------------------------
       ! 
-      WRITE(numout,*) 'TG print 5: Start of updateu before = ',before
-      CALL FLUSH(numout)
       IF( before ) THEN
          zrhoy = Agrif_Rhoy()
@@ -412,5 +400,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)
+                  ptab(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)
                   ptab(ji,jj,jk,2) = umask(ji,jj,jk) * zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) 
                END DO
@@ -419,5 +407,5 @@
       ELSE
 ! VERTICAL REFINEMENT BEGIN
-         ptab_child(:,:,:,:) = 0.
+         ptab_child(:,:,:) = 0.
          
          DO jj=j1,j2
@@ -427,5 +415,5 @@
              IF (ptab(ji,jj,jk,2) == 0) EXIT
              N_in = N_in + 1
-             tabin(jk) = ptab(ji,jj,jk)/ptab(ji,jj,jk,2)
+             tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
              h_in(N_in) = ptab(ji,jj,jk,2)/e2u(ji,jj)
            ENDDO
@@ -447,5 +435,5 @@
 !             h_in(N_in) = h_diff
              endif
-             CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out,1),h_out(1:N_out),N_in,N_out)
+             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)
           ENDIF
          ENDDO
@@ -466,8 +454,8 @@
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      ub(ji,jj,jk) = ub(ji,jj,jk) & 
-                           & + atfp * ( ptab_child(ji,jj,jk,1) - un(ji,jj,jk) ) * umask(ji,jj,jk)
+                           & + atfp * ( ptab_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)
                   ENDIF
                   !
-                  un(ji,jj,jk) = ptab_child(ji,jj,jk,1) * umask(ji,jj,jk)
+                  un(ji,jj,jk) = ptab_child(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
@@ -475,6 +463,4 @@
       ENDIF
       ! 
-      WRITE(numout,*) 'TG print 6: End of updateu before = ',before
-      CALL FLUSH(numout)
    END SUBROUTINE updateu
 
@@ -491,5 +477,5 @@
       REAL(wp) :: zrhox
 ! VERTICAL REFINEMENT BEGIN
-      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,2) :: ptab_child
+      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
       REAL(wp) :: h_in(k1:k2)
       REAL(wp) :: h_out(1:jpk)
@@ -500,6 +486,4 @@
       !!---------------------------------------------      
       !
-      WRITE(numout,*) 'TG print 7: Start of updatev before = ',before
-      CALL FLUSH(numout)
       IF (before) THEN
          zrhox = Agrif_Rhox()
@@ -507,5 +491,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  ptab(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
+                  ptab(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)
                   ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) 
                END DO
@@ -520,7 +504,7 @@
            N_in = 0
            DO jk=k1,k2
-             IF (update_scales_v(ji,jj,jk) == 0) EXIT
+             IF (ptab(ji,jj,jk,2) == 0) EXIT
              N_in = N_in + 1
-             tabin(jk) = ptab(ji,jj,1)/ptab(ji,jj,jk,2)
+             tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
              h_in(N_in) = ptab(ji,jj,jk,2)/e1v(ji,jj)
            ENDDO
@@ -542,5 +526,5 @@
 !             h_in(N_in) = h_diff
              endif
-             CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out,1),h_out(1:N_out),N_in,N_out)
+             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)
           ENDIF
          ENDDO
@@ -553,5 +537,5 @@
 ! VERTICAL REFINEMENT END
 
-         DO jk=k1,k2
+         DO jk=k1,jpk
             DO jj=j1,j2
                DO ji=i1,i2
@@ -561,14 +545,12 @@
                   IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part
                      vb(ji,jj,jk) = vb(ji,jj,jk) & 
-                           & + atfp * ( ptab_child(ji,jj,jk,1) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
+                           & + atfp * ( ptab_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)
                   ENDIF
                   !
-                  vn(ji,jj,jk) = ptab_child(ji,jj,jk,1) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ENDIF
-      WRITE(numout,*) 'TG print 8: End of updatev before = ',before
-      ! 
+                  vn(ji,jj,jk) = ptab_child(ji,jj,jk) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
    END SUBROUTINE updatev
 
