Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90	(revision 11590)
@@ -17,7 +17,5 @@
 
    PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90
-#if defined key_vertical
-   PUBLIC reconstructandremap ! remapping routine
-#endif   
+  
    !                                              !!* Namelist namagrif: AGRIF parameters
    LOGICAL , PUBLIC ::   ln_agrif_2way = .TRUE.    !: activate two way nesting 
@@ -95,222 +93,4 @@
    END FUNCTION agrif_oce_alloc
 
-#if defined key_vertical
-   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
-      !!----------------------------------------------------------------------
-      !!                ***  FUNCTION reconstructandremap  ***
-      !!----------------------------------------------------------------------
-      IMPLICIT NONE
-      INTEGER N, Nout
-      REAL(wp) tabin(N), tabout(Nout)
-      REAL(wp) hin(N), hout(Nout)
-      REAL(wp) coeffremap(N,3),zwork(N,3)
-      REAL(wp) zwork2(N+1,3)
-      INTEGER jk
-      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8  
-      REAL(wp) q,q01,q02,q001,q002,q0
-      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
-      REAL(wp),PARAMETER :: dpthin = 1.D-3
-      INTEGER :: k1, kbox, ktop, ka, kbot
-      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
-
-      z_win(1)=0.; z_wout(1)= 0.
-      DO jk=1,N
-         z_win(jk+1)=z_win(jk)+hin(jk)
-      ENDDO 
-      
-      DO jk=1,Nout
-         z_wout(jk+1)=z_wout(jk)+hout(jk)       
-      ENDDO       
-
-      DO jk=2,N
-         zwork(jk,1)=1./(hin(jk-1)+hin(jk))
-      ENDDO
-        
-      DO jk=2,N-1
-         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
-         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
-         zwork(jk,3)=q0
-      ENDDO       
-     
-      DO jk= 2,N
-         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
-      ENDDO
-        
-      coeffremap(:,1) = tabin(:)
- 
-      DO jk=2,N-1
-         q001 = hin(jk)*zwork2(jk+1,1)
-         q002 = hin(jk)*zwork2(jk,1)        
-         IF (q001*q002 < 0) then
-            q001 = 0.
-            q002 = 0.
-         ENDIF
-         q=zwork(jk,2)
-         q01=q*zwork2(jk+1,1)
-         q02=q*zwork2(jk,1)
-         IF (abs(q001) > abs(q02)) q001 = q02
-         IF (abs(q002) > abs(q01)) q002 = q01
-        
-         q=(q001-q002)*zwork(jk,3)
-         q001=q001-q*hin(jk+1)
-         q002=q002+q*hin(jk-1)
-        
-         coeffremap(jk,3)=coeffremap(jk,1)+q001
-         coeffremap(jk,2)=coeffremap(jk,1)-q002
-        
-         zwork2(jk,1)=(2.*q001-q002)**2
-         zwork2(jk,2)=(2.*q002-q001)**2
-      ENDDO
-        
-      DO jk=1,N
-         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
-            coeffremap(jk,3) = coeffremap(jk,1)
-            coeffremap(jk,2) = coeffremap(jk,1)
-            zwork2(jk,1) = 0.
-            zwork2(jk,2) = 0.
-         ENDIF
-      ENDDO
-        
-      DO jk=2,N
-         q002=max(zwork2(jk-1,2),dsmll)
-         q001=max(zwork2(jk,1),dsmll)
-         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,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 jk=1,N
-         q01=zwork2(jk+1,3)-coeffremap(jk,1)
-         q02=coeffremap(jk,1)-zwork2(jk,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(jk,2)=coeffremap(jk,1)-q02
-         coeffremap(jk,3)=coeffremap(jk,1)+q01
-      ENDDO
-
-      zbot=0.0
-      kbot=1
-      DO jk=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(jk+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= jk+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(jk).NE.z_win(kbox).OR.z_wout(jk+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(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
-               ELSE
-                  tabout(jk) = tabin(kbox)
-               ENDIF 
-            ELSE
-               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
-                  ka = jk
-               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(jk)=offset+tsum*rpsum
-                 
-            ENDIF !single or multiple layers
-         ELSE
-            IF (jk==1) THEN
-               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
-            ENDIF
-            tabout(jk) = tabout(jk-1)
-             
-         ENDIF !normal:thin layer
-      ENDDO !jk
-            
-      return
-      end subroutine reconstructandremap
-#endif
-
 #endif
    !!======================================================================
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90	(revision 11590)
@@ -33,4 +33,5 @@
    USE agrif_oce_sponge
    USE lib_mpp
+   USE vremap
  
    IMPLICIT NONE
@@ -661,4 +662,5 @@
       REAL(wp), DIMENSION(k1:k2) :: h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
+      !!----------------------------------------------------------------------
 
       IF( before ) THEN         
@@ -698,9 +700,9 @@
                   IF (tmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(ji,jj,jk)
+                  h_out(jk) = e3t_a(ji,jj,jk)
                ENDDO
                IF (N_in > 0) THEN
                   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)
+                     call reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),ptab_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)
                   ENDDO
                ENDIF
@@ -737,5 +739,5 @@
    END SUBROUTINE interpsshn
 
-   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
+   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
       !!----------------------------------------------------------------------
       !!                  *** ROUTINE interpun ***
@@ -745,5 +747,4 @@
       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
@@ -752,7 +753,6 @@
       REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER  :: N_in, N_out, iref
+      INTEGER  :: N_in, N_out
       REAL(wp) :: h_diff
-      LOGICAL  :: western_side, eastern_side
       !!---------------------------------------------    
       !
@@ -772,11 +772,6 @@
 # if defined key_vertical
 ! VERTICAL REFINEMENT BEGIN
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
 
          DO ji=i1,i2
-            iref = ji
-            IF (western_side) iref = MAX(2,ji)
-            IF (eastern_side) iref = MIN(nlci-2,ji)
             DO jj=j1,j2
                N_in = 0
@@ -795,7 +790,7 @@
               N_out = 0
               DO jk=1,jpk
-                 if (umask(iref,jj,jk) == 0) EXIT
+                 if (umask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3u_a(iref,jj,jk)
+                 h_out(N_out) = e3u_a(ji,jj,jk)
               ENDDO
          
@@ -829,5 +824,5 @@
    END SUBROUTINE interpun
 
-   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
+   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
       !!----------------------------------------------------------------------
       !!                  *** ROUTINE interpvn ***
@@ -837,5 +832,4 @@
       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
@@ -844,7 +838,6 @@
       REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
-      INTEGER  :: N_in, N_out, jref
+      INTEGER  :: N_in, N_out
       REAL(wp) :: h_diff
-      LOGICAL  :: northern_side,southern_side
       !!---------------------------------------------    
       !      
@@ -864,11 +857,5 @@
 # if defined key_vertical
 
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-
          DO jj=j1,j2
-            jref = jj
-            IF (southern_side) jref = MAX(2,jj)
-            IF (northern_side) jref = MIN(nlcj-2,jj)
             DO ji=i1,i2
                N_in = 0
@@ -886,7 +873,7 @@
                N_out = 0
                DO jk=1,jpk
-                  if (vmask(ji,jref,jk) == 0) EXIT
+                  if (vmask(ji,jj,jk) == 0) EXIT
                   N_out = N_out + 1
-                  h_out(N_out) = e3v_a(ji,jref,jk)
+                  h_out(N_out) = e3v_a(ji,jj,jk)
                END DO
                IF (N_out == 0) THEN
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90	(revision 11590)
@@ -23,4 +23,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE iom
+   USE vremap
 
    IMPLICIT NONE
@@ -319,5 +320,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 
                END DO
             END DO
@@ -342,11 +343,9 @@
                   IF (tmask(ji,jj,jk) == 0) EXIT 
                   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
+                  h_out(jk) = e3t_b(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))
-                  tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
                   DO jn=1,jpts
-                     call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
+                     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
@@ -446,5 +445,5 @@
                   tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
 # if defined key_vertical
-                  tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
+                  tabres(ji,jj,jk,m2) = e3u_b(ji,jj,jk)*umask(ji,jj,jk)
 # endif
                END DO
@@ -475,5 +474,5 @@
                  if (umask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3u_n(ji,jj,jk)
+                 h_out(N_out) = e3u_b(ji,jj,jk)
               ENDDO
          
@@ -507,7 +506,7 @@
             DO jj = j1,j2
                DO ji = i1+1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
-                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
-                                     &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
+                  zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_b(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
+                                     &   -e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
                END DO
             END DO
@@ -599,5 +598,5 @@
                   tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
 # if defined key_vertical
-                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
+                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_b(ji,jj,jk)
 # endif
                END DO
@@ -627,5 +626,5 @@
                  if (vmask(ji,jj,jk) == 0) EXIT
                  N_out = N_out + 1
-                 h_out(N_out) = e3v_n(ji,jj,jk)
+                 h_out(N_out) = e3v_b(ji,jj,jk)
               ENDDO
          
@@ -653,7 +652,7 @@
             DO jj = j1+1,j2
                DO ji = i1,i2   ! vector opt.
-                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
-                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
-                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
+                  zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * fsahm_spt(ji,jj)
+                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_b(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
+                                     &  -e1v(ji,jj-1) * e3v_b(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
                END DO
             END DO
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90	(revision 11590)
@@ -24,4 +24,5 @@
    USE lib_mpp        ! MPP library
    USE domvvl         ! Need interpolation routines 
+   USE vremap         ! Vertical remapping
 
    IMPLICIT NONE
@@ -319,5 +320,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) < -900._wp  ) EXIT
                   N_in = N_in + 1
                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_interp.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_interp.F90	(revision 11590)
@@ -48,5 +48,5 @@
    END SUBROUTINE Agrif_trc
 
-   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
+   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!----------------------------------------------------------------------
       !!                  *** ROUTINE interptrn ***
@@ -55,10 +55,9 @@
       INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
       LOGICAL                                     , INTENT(in   ) ::   before
-      INTEGER                                     , INTENT(in   ) ::   nb , ndir
       !
-      INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jn, ibdy, jbdy   ! dummy loop indices
       INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
       REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+
       ! vertical interpolation:
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
@@ -66,5 +65,5 @@
       REAL(wp), DIMENSION(k1:k2) :: h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
-      REAL(wp) :: h_diff
+      !!----------------------------------------------------------------------
 
       IF( before ) THEN         
@@ -91,15 +90,6 @@
 
 # if defined key_vertical
-         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)
-
          DO jj=j1,j2
             DO ji=i1,i2
-               iref = ji
-               jref = jj
-               if(western_side) iref=MAX(2,ji)
-               if(eastern_side) iref=MIN(nlci-1,ji)
-               if(southern_side) jref=MAX(2,jj)
-               if(northern_side) jref=MIN(nlcj-1,jj)
                N_in = 0
                DO jk=k1,k2 !k2 = jpk of parent grid
@@ -111,7 +101,7 @@
                N_out = 0
                DO jk=1,jpk ! jpk of child grid
-                  IF (tmask(iref,jref,jk) == 0) EXIT 
+                  IF (tmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(iref,jref,jk)
+                  h_out(jk) = e3t_a(ji,jj,jk)
                ENDDO
                IF (N_in > 0) THEN
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_sponge.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_sponge.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_sponge.F90	(revision 11590)
@@ -93,5 +93,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 
                END DO
             END DO
@@ -114,9 +114,7 @@
                   IF (tmask(ji,jj,jk) == 0) EXIT 
                   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
+                  h_out(jk) = e3t_b(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))
-                  tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
                   DO jn=1,jptra
                      call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90	(revision 11574)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90	(revision 11590)
@@ -324,7 +324,7 @@
       ! 3. Location of interpolation
       !-----------------------------
-      CALL Agrif_Set_bc(       tsn_id, (/0,ind1/) )
-      CALL Agrif_Set_bc( un_interp_id, (/0,ind1/) )
-      CALL Agrif_Set_bc( vn_interp_id, (/0,ind1/) )
+      CALL Agrif_Set_bc(       tsn_id, (/0,ind1-1/) )
+      CALL Agrif_Set_bc( un_interp_id, (/0,ind1-1/) )
+      CALL Agrif_Set_bc( vn_interp_id, (/0,ind1-1/) )
 
       CALL Agrif_Set_bc( tsn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) )  ! if west and rhox=3 and sponge=2 and ghost=1: columns 2 to 9 
@@ -594,5 +594,5 @@
       ! 3. Location of interpolation
       !-----------------------------
-      CALL Agrif_Set_bc(trn_id,(/0,ind1/))
+      CALL Agrif_Set_bc(trn_id,(/0,ind1-1/))
       CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/))
 
Index: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90
===================================================================
--- NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90	(revision 11590)
+++ NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90	(revision 11590)
@@ -0,0 +1,340 @@
+#undef PPR_LIB      /* USE PPR library */
+MODULE vremap
+!$AGRIF_DO_NOT_TREAT
+   !!======================================================================
+   !!                       ***  MODULE  vremap  ***
+   !! Ocean physics:  Vertical remapping routines
+   !!
+   !!======================================================================
+   !! History : 4.0  !  2019-09  (Jérôme Chanut)  Original code
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!
+   !!----------------------------------------------------------------------
+   USE par_oce
+#if defined PPR_LIB
+   USE ppr_1d   ! D. Engwirda piecewise polynomial reconstruction library
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   reconstructandremap
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+   !! $Id: vremap 11573 2019-09-19 09:18:03Z jchanut $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+#if ! defined PPR_LIB
+   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE  reconstructandremap ***
+      !!
+      !! ** Purpose :   Brief description of the routine
+      !!
+      !! ** Method  :   description of the methodoloy used to achieve the
+      !!                objectives of the routine. Be as clear as possible!
+      !!
+      !! ** Action  : - first action (share memory array/varible modified
+      !!                in this routine
+      !!              - second action .....
+      !!              - .....
+      !!
+      !! References :   Author et al., Short_name_review, Year
+      !!                Give references if exist otherwise suppress these lines
+      !!-----------------------------------------------------------------------
+      INTEGER N, Nout
+      REAL(wp) tabin(N), tabout(Nout)
+      REAL(wp) hin(N), hout(Nout)
+      REAL(wp) coeffremap(N,3),zwork(N,3)
+      REAL(wp) zwork2(N+1,3)
+      INTEGER jk
+      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8  
+      REAL(wp) q,q01,q02,q001,q002,q0
+      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
+      REAL(wp),PARAMETER :: dpthin = 1.D-3
+      INTEGER :: k1, kbox, ktop, ka, kbot
+      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+      !!-----------------------------------------------------------------------
+
+      z_win(1)=0.; z_wout(1)= 0.
+      DO jk=1,N
+         z_win(jk+1)=z_win(jk)+hin(jk)
+      ENDDO 
+      
+      DO jk=1,Nout
+         z_wout(jk+1)=z_wout(jk)+hout(jk)       
+      ENDDO       
+
+      DO jk=2,N
+         zwork(jk,1)=1./(hin(jk-1)+hin(jk))
+      ENDDO
+        
+      DO jk=2,N-1
+         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
+         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
+         zwork(jk,3)=q0
+      ENDDO       
+     
+      DO jk= 2,N
+         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
+      ENDDO
+        
+      coeffremap(:,1) = tabin(:)
+ 
+      DO jk=2,N-1
+         q001 = hin(jk)*zwork2(jk+1,1)
+         q002 = hin(jk)*zwork2(jk,1)        
+         IF (q001*q002 < 0) then
+            q001 = 0.
+            q002 = 0.
+         ENDIF
+         q=zwork(jk,2)
+         q01=q*zwork2(jk+1,1)
+         q02=q*zwork2(jk,1)
+         IF (abs(q001) > abs(q02)) q001 = q02
+         IF (abs(q002) > abs(q01)) q002 = q01
+        
+         q=(q001-q002)*zwork(jk,3)
+         q001=q001-q*hin(jk+1)
+         q002=q002+q*hin(jk-1)
+        
+         coeffremap(jk,3)=coeffremap(jk,1)+q001
+         coeffremap(jk,2)=coeffremap(jk,1)-q002
+        
+         zwork2(jk,1)=(2.*q001-q002)**2
+         zwork2(jk,2)=(2.*q002-q001)**2
+      ENDDO
+        
+      DO jk=1,N
+         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
+            coeffremap(jk,3) = coeffremap(jk,1)
+            coeffremap(jk,2) = coeffremap(jk,1)
+            zwork2(jk,1) = 0.
+            zwork2(jk,2) = 0.
+         ENDIF
+      ENDDO
+        
+      DO jk=2,N
+         q002=max(zwork2(jk-1,2),dsmll)
+         q001=max(zwork2(jk,1),dsmll)
+         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,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 jk=1,N
+         q01=zwork2(jk+1,3)-coeffremap(jk,1)
+         q02=coeffremap(jk,1)-zwork2(jk,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(jk,2)=coeffremap(jk,1)-q02
+         coeffremap(jk,3)=coeffremap(jk,1)+q01
+      ENDDO
+
+      zbot=0.0
+      kbot=1
+      DO jk=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(jk+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= jk+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(jk).NE.z_win(kbox).OR.z_wout(jk+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(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+               ELSE
+                  tabout(jk) = tabin(kbox)
+               ENDIF 
+            ELSE
+               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
+                  ka = jk
+               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(jk)=offset+tsum*rpsum
+                 
+            ENDIF !single or multiple layers
+         ELSE
+            IF (jk==1) THEN
+               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
+            ENDIF
+            tabout(jk) = tabout(jk-1)
+             
+         ENDIF !normal:thin layer
+      ENDDO !jk
+            
+   END SUBROUTINE reconstructandremap
+
+#else
+
+   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, jkin, jkout)
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE  reconstructandremap ***
+      !!
+      !! ** Purpose :   Brief description of the routine
+      !!
+      !! ** Method  :   description of the methodoloy used to achieve the
+      !!                objectives of the routine. Be as clear as possible!
+      !!
+      !! ** Action  : - first action (share memory array/varible modified
+      !!                in this routine
+      !!              - second action .....
+      !!              - .....
+      !!
+      !! References :   Author et al., Short_name_review, Year
+      !!                Give references if exist otherwise suppress these lines
+      !!-----------------------------------------------------------------------
+      INTEGER , INTENT(in   )                     ::   jkin, jkout       ! 
+      REAL(wp), INTENT(in   ), DIMENSION(jkin)    ::   phin, ptin        !
+      REAL(wp), INTENT(in   ), DIMENSION(jkout)   ::   phout             !
+      REAL(wp), INTENT(inout), DIMENSION(jkout)   ::   ptout             !
+      !
+      INTEGER, PARAMETER :: ndof = 1, nvar = 1
+      INTEGER  :: jk
+      REAL(wp) ::  zwin(jkin+1),   ztin(ndof, nvar, jkin)
+      REAL(wp) :: zwout(jkout+1), ztout(ndof, nvar, jkout)
+      TYPE(rmap_work) :: work
+      TYPE(rmap_opts) :: opts
+      TYPE(rcon_ends) :: bc_l(nvar)
+      TYPE(rcon_ends) :: bc_r(nvar)
+      !!--------------------------------------------------------------------
+     
+      ! Set interfaces and input data:
+      zwin(1) = 0._wp
+      DO jk = 2, jkin + 1
+         zwin(jk) = zwin(jk-1) + phin(jk-1) 
+         ztin(ndof, nvar,jk-1) =  ptin(jk-1)
+      END DO
+
+      zwout(1) = 0._wp
+      DO jk = 2, jkout + 1
+         zwout(jk) = zwout(jk-1) + phout(jk-1) 
+      END DO
+
+      ! specify methods
+!      opts%edge_meth = p3e_method     ! 3rd-order edge interp.
+!      opts%cell_meth = ppm_method     ! PPM method in cells
+!      opts%cell_lims = null_limit     ! no lim.     
+      opts%edge_meth = p5e_method     ! 5th-order edge interp.
+      opts%cell_meth = pqm_method     ! PPM method in cells
+      opts%cell_lims = mono_limit     ! monotone limiter   
+ 
+      ! set boundary conditions
+      bc_l%bcopt = bcon_loose         ! "loose" = extrapolate
+      bc_r%bcopt = bcon_loose
+!      bc_l%bcopt = bcon_slope        
+!      bc_r%bcopt = bcon_slope
+
+      ! init. method workspace
+      CALL work%init(jkin+1, nvar, opts)
+
+      ! remap
+      CALL rmap1d(jkin+1, jkout+1, nvar, ndof, &
+      &           zwin, zwout, ztin, ztout,    &
+      &           bc_l, bc_r, work, opts)
+
+      ! clear method workspace
+      CALL work%free()
+
+      DO jk =1, jkout
+         ptout(jk) = ztout(1,1,jk)
+      END DO
+            
+   END SUBROUTINE reconstructandremap
+#endif
+
+   !!======================================================================
+!$AGRIF_END_DO_NOT_TREAT
+END MODULE vremap
