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 8860)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 8861)
@@ -576,5 +576,4 @@
       !!----------------------------------------------------------------------  
       !
-      return
       
       zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
@@ -602,46 +601,32 @@
       !
       INTEGER  ::   ji, jj, jk, jn, iref, jref   ! dummy loop indices
-      INTEGER  ::   imin, imax, jmin, jmax
+      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
       REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
       REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
       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
+      REAL(wp) :: h_diff, zrhoxy
 
       zrhoxy = Agrif_rhox()*Agrif_rhoy()
-
-
-
       IF (before) THEN         
-            IF(Agrif_UseSpecialValue) THEN 
-!               Agrif_SpecialValue = -999._wp
-               Agrif_SpecialValue = 0._wp
-            ELSE
-               Agrif_SpecialValue = 0._wp
-            ENDIF
-            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)! * tmask(ji,jj,jk) - &
-!                                          & (tmask(ji,jj,jk)-1) * Agrif_SpecialValue
-                     END DO
-                  END DO
-               END DO
-            END DO
+         DO jn = n1,n2-1
             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
+                 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 
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -649,70 +634,44 @@
          southern_side = (nb == 2).AND.(ndir == 1)
          northern_side = (nb == 2).AND.(ndir == 2)
-         Agrif_SpecialValue = 0._wp !reset now interpolation is done
-! VERTICAL REFINEMENT BEGIN
 #ifdef key_vertical         
          ptab_child(:,:,:,:) = 0.
-         do jj=j1,j2
-         do ji=i1,i2
-           iref = ji
-           jref = jj
-           if(western_side) iref=2
-           if(eastern_side) iref=nlci-1
-           if(southern_side) jref=2
-           if(northern_side) jref=nlcj-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 = 0
-           DO jk=1,jpk ! jpk of child grid
-             IF (tmask(iref,jref,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(iref,jref,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
+         DO jj=j1,j2
+            DO ji=i1,i2
+               iref = ji
+               jref = jj
+               if(western_side) iref=2
+               if(eastern_side) iref=nlci-1
+               if(southern_side) jref=2
+               if(northern_side) jref=nlcj-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 = 0
+               DO jk=1,jpk ! jpk of child grid
+                  IF (tmask(iref,jref,jk) == 0) EXIT 
+                  N_out = N_out + 1
+                  h_out(jk) = e3t_n(iref,jref,jk)
+               ENDDO
+               IF (N_in > 0) THEN
+                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
+                  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
+               ENDIF
+            ENDDO
+         ENDDO
 #else
-         do jk=k1,k2
-            do jj=j1,j2
-               do ji=i1,i2
-                 ! This will be slow - Is there a way to speed it up and avoid divide by zero?
-                 IF (ptab(ji,jj,jk,n2) .NE. 0) THEN
-                    ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2)
-                 ELSE
-                    ptab_child(ji,jj,jk,n1:n2-1) = 0._wp
-                 ENDIF
-               enddo
-            enddo
-         enddo
+         DO jk=k1,k2
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/(zrhoxy*e1e2t(ji,jj)*e3t_n(ji,jj,jk)
+               ENDDO
+            ENDDO
+         ENDDO
 #endif
-   
-
-! VERTICAL REFINEMENT END
-
          !
          zrhox = Agrif_Rhox()
@@ -996,5 +955,4 @@
       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
@@ -1006,5 +964,4 @@
       INTEGER :: jref
 
-! VERTICAL REFINEMENT END
       !!---------------------------------------------    
       !      
@@ -1029,5 +986,4 @@
          ptab_child(:,:,:) = 0.
 #ifdef key_vertical
-! VERTICAL REFINEMENT BEGIN
          southern_side = (nb == 2).AND.(ndir == 1)
          northern_side = (nb == 2).AND.(ndir == 2)
@@ -1061,19 +1017,7 @@
                  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)
             DO jk=1,jpkm1
                DO jj=j1,j2
@@ -1081,5 +1025,4 @@
                END DO
             END DO
-! VERTICAL REFINEMENT END
 #else 
             DO jk=1,jpkm1
@@ -1448,17 +1391,64 @@
 # if defined key_zdftke
 
-   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
+   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interavm  ***
       !!----------------------------------------------------------------------  
-      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab
+      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
+      REAL(wp), DIMENSION(k1:k2) :: tabin
+      REAL(wp) :: h_in(k1:k2)
+      REAL(wp) :: h_out(1:jpk)
+      REAL(wp) :: zrhoxy
+      INTEGER  :: N_in, N_out, ji, jj, jk
       !!----------------------------------------------------------------------  
       !      
-      IF( before ) THEN
-         ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
-      ELSE
-         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
+      zrhoxy = Agrif_rhox()*Agrif_rhoy()
+      IF (before) THEN         
+         DO jk=k1,k2
+            DO jj=j1,j2
+              DO ji=i1,i2
+                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
+              END DO
+           END DO
+        END DO
+#ifdef key_vertical         
+        DO jk=k1,k2
+           DO jj=j1,j2
+              DO ji=i1,i2
+                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e1e2t(ji,jj) * e3w_n(ji,jj,jk) 
+              END DO
+           END DO
+        END DO
+#else
+      ptab(i1:i2,j1:j2,k1:k2,2) = 0._wp
+#endif
+      ELSE 
+#ifdef key_vertical         
+         avm_k(i1:i2,j1:j2,1:jpk) = 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,2) == 0) EXIT
+                  N_in = N_in + 1
+                  tabin(jk) = ptab(ji,jj,jk,1)
+                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1e2t(ji,jj)*zrhoxy)
+               END DO
+               N_out = 0
+               DO jk=1,jpk ! jpk of child grid
+                  IF (wmask(ji,jj,jk) == 0) EXIT 
+                  N_out = N_out + 1
+                  h_out(jk) = e3t_n(ji,jj,jk)
+               ENDDO
+               IF (N_in > 0) THEN
+                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
+               ENDIF
+            ENDDO
+         ENDDO
+#else
+         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
+#endif
       ENDIF
       !
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 8860)
+++ /branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90	(revision 8861)
@@ -374,5 +374,5 @@
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id)
    CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id)
-   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id)
+   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,2/),avm_id)
 # endif
 
@@ -757,4 +757,6 @@
    ENDIF
 
+   AGRIF_Parent(ln_vertical)=ln_vertical ! Will this work?
+
    CALL Agrif_Update_trc(0)
    !
@@ -832,5 +834,5 @@
    INTEGER  ::   ios                 ! Local integer output status for namelist read
    INTEGER  ::   iminspon
-   NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy
+   NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy, ln_vertical
    !!--------------------------------------------------------------------------------------
    !
@@ -854,4 +856,5 @@
       WRITE(numout,*) '      use special values for dynamics   ln_spc_dyn    = ', ln_spc_dyn
       WRITE(numout,*) '      check bathymetry                  ln_chk_bathy  = ', ln_chk_bathy
+      WRITE(numout,*) '      use vertical interpolation        ln_vertical   = ', ln_vertical
       WRITE(numout,*) 
    ENDIF
