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 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90	(revision 11603)
@@ -45,4 +45,6 @@
    PUBLIC   interpe3t, interpumsk, interpvmsk
 
+   INTEGER ::   bdy_tinterp = 0
+
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
@@ -501,9 +503,11 @@
          CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
          !
+         bdy_tinterp = 1
          CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
-         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )
+         CALL Agrif_Bc_variable( vnb_id        , calledweight=1._wp, procname=interpvnb  )  
          !
+         bdy_tinterp = 2
          CALL Agrif_Bc_variable( unb_id        , calledweight=0._wp, procname=interpunb  ) ! Before
-         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )         
+         CALL Agrif_Bc_variable( vnb_id        , calledweight=0._wp, procname=interpvnb  )   
       ELSE ! Linear interpolation
          !
@@ -703,7 +707,5 @@
                ENDDO
                IF (N_in > 0) THEN
-                  DO jn=1,jpts
-                     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
+                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),ptab_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
                ENDIF
             ENDDO
@@ -800,13 +802,5 @@
               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 < -1.e4) 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),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+              CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
             ENDDO
          ENDDO
@@ -881,5 +875,5 @@
                  CYCLE
                ENDIF
-               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
             END DO
          END DO
@@ -916,22 +910,23 @@
          DO ji = i1, i2
             DO jj = j1, j2
-               IF    ( utint_stage(ji,jj) == 1  ) THEN
-                  ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-                     &               - zt0**2._wp * (       zt0 - 1._wp)        )
-               ELSEIF( utint_stage(ji,jj) == 2  ) THEN
-                  ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-                     &               - zt0        * (       zt0 - 1._wp)**2._wp )
-               ELSEIF( utint_stage(ji,jj) == 0  ) THEN                
-                  ztcoeff = 1._wp
-               ELSE
-                  ztcoeff = 0._wp
-               ENDIF
-               !   
-               ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
-               !            
-               IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
-                  ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
-                  utint_stage(ji,jj) = 3
-               ELSE
+               IF ( utint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
+                  IF    ( utint_stage(ji,jj) == 1  ) THEN
+                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
+                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
+                  ELSEIF( utint_stage(ji,jj) == 2  ) THEN
+                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
+                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
+                  ELSEIF( utint_stage(ji,jj) == 0  ) THEN                
+                     ztcoeff = 1._wp
+                  ELSE
+                     ztcoeff = 0._wp
+                  ENDIF
+                  !   
+                  ubdy(ji,jj) = ubdy(ji,jj) + ztcoeff * ptab(ji,jj)
+                  !            
+                  IF (( utint_stage(ji,jj) == 2 ).OR.( utint_stage(ji,jj) == 0 )) THEN
+                     ubdy(ji,jj) = ubdy(ji,jj) / (zrhoy*e2u(ji,jj)) * umask(ji,jj,1)
+                  ENDIF
+                  !
                   utint_stage(ji,jj) = utint_stage(ji,jj) + 1
                ENDIF
@@ -966,22 +961,23 @@
          DO ji = i1, i2
             DO jj = j1, j2
-               IF    ( vtint_stage(ji,jj) == 1  ) THEN
-                  ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-                     &               - zt0**2._wp * (       zt0 - 1._wp)        )
-               ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
-                  ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-                     &               - zt0        * (       zt0 - 1._wp)**2._wp )
-               ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                
-                  ztcoeff = 1._wp
-               ELSE
-                  ztcoeff = 0._wp
-               ENDIF
-               !   
-               vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
-               !            
-               IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
-                  vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
-                  vtint_stage(ji,jj) = 3
-               ELSE
+               IF ( vtint_stage(ji,jj) < (bdy_tinterp + 1) ) THEN
+                  IF    ( vtint_stage(ji,jj) == 1  ) THEN
+                     ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
+                        &               - zt0**2._wp * (       zt0 - 1._wp)        )
+                  ELSEIF( vtint_stage(ji,jj) == 2  ) THEN
+                     ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
+                        &               - zt0        * (       zt0 - 1._wp)**2._wp )
+                  ELSEIF( vtint_stage(ji,jj) == 0  ) THEN                
+                     ztcoeff = 1._wp
+                  ELSE
+                     ztcoeff = 0._wp
+                  ENDIF
+                  !   
+                  vbdy(ji,jj) = vbdy(ji,jj) + ztcoeff * ptab(ji,jj)
+                  !            
+                  IF (( vtint_stage(ji,jj) == 2 ).OR.( vtint_stage(ji,jj) == 0 )) THEN
+                     vbdy(ji,jj) = vbdy(ji,jj) / (zrhox*e1v(ji,jj)) * vmask(ji,jj,1)
+                  ENDIF
+                  !
                   vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
                ENDIF
@@ -1249,5 +1245,5 @@
                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)
+                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out,1)
                ENDIF
             ENDDO
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 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90	(revision 11603)
@@ -346,7 +346,5 @@
                ENDDO
                IF (N_in > 0) THEN
-                  DO jn=1,jpts
-                     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
+                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
                ENDIF
             ENDDO
@@ -488,5 +486,5 @@
                  endif
               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)
+              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,1)
          
             ENDDO
@@ -635,5 +633,5 @@
                  endif
               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)
+              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,1)
             ENDDO
          ENDDO
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 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90	(revision 11603)
@@ -339,7 +339,5 @@
                      STOP
                   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
+                  CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
                ENDIF
             ENDDO
@@ -524,5 +522,5 @@
                      ENDDO
                   ENDIF
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+                  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,1)
                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
                ENDIF
@@ -709,5 +707,5 @@
                      ENDDO
                   ENDIF
-                  CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
+                  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,1)
                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
                ENDIF
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 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_interp.F90	(revision 11603)
@@ -106,7 +106,5 @@
                ENDDO
                IF (N_in > 0) THEN
-                  DO jn=1,jptra
-                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
-                  ENDDO
+                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,ptab_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra)
                ENDIF
             ENDDO
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 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_sponge.F90	(revision 11603)
@@ -117,7 +117,5 @@
                ENDDO
                IF (N_in > 0) THEN
-                  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)
-                  ENDDO
+                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in,tabres_child(ji,jj,1:N_out,1:jptra),h_out,N_in,N_out,jptra)
                ENDIF
             ENDDO
Index: /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_update.F90
===================================================================
--- /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_update.F90	(revision 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_update.F90	(revision 11603)
@@ -120,7 +120,5 @@
                      STOP
                   ENDIF
-                  DO jn=1,jptra
-                     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
+                  CALL reconstructandremap(tabin(1:N_in,1:jptra),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
                ENDIF
             ENDDO
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 11602)
+++ /NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90	(revision 11603)
@@ -30,5 +30,5 @@
 
 #if ! defined PPR_LIB
-   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)      
+   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)      
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE  reconstructandremap ***
@@ -47,212 +47,218 @@
       !!                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
+      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
+      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
+      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
+      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   phin       ! Input thicknesses
+      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   phout      ! Output thicknesses
+      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
+      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Remapped data
+      !
+      INTEGER             :: jk, jn, k1, kbox, ktop, ka, kbot
+      REAL(wp), PARAMETER :: dpthin = 1.D-3, dsmll = 1.0D-8
+      REAL(wp)            :: q, q01, q02, q001, q002, q0
+      REAL(wp)            :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
+      REAL(wp)            :: coeffremap(kjpk_in,3), zwork(kjpk_in,3), zwork2(kjpk_in+1,3)
+      REAL(wp)            :: z_win(1:kjpk_in+1), z_wout(1:kjpk_out+1)
       !!-----------------------------------------------------------------------
 
-      z_win(1)=0.; z_wout(1)= 0.
-      DO jk=1,N
-         z_win(jk+1)=z_win(jk)+hin(jk)
-      ENDDO 
+      z_win(1)=0._wp ; z_wout(1)= 0._wp
+      DO jk = 1, kjpk_in
+         z_win(jk+1)=z_win(jk)+phin(jk)
+      END DO 
       
-      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 = 1, kjpk_out
+         z_wout(jk+1)=z_wout(jk)+phout(jk)       
+      END DO       
+
+      DO jk = 2, kjpk_in
+         zwork(jk,1)=1./(phin(jk-1)+phin(jk))
+      END DO
+        
+      DO jk = 2, kjpk_in-1
+         q0 = 1._wp / (phin(jk-1)+phin(jk)+phin(jk+1))
+         zwork(jk,2) = phin(jk-1) + 2._wp*phin(jk) + phin(jk+1)
+         zwork(jk,3) = q0
+      END DO       
      
-      DO jk= 2,N
-         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
-      ENDDO
-        
-      coeffremap(:,1) = tabin(:)
+      DO jn = 1, kn_var
+
+         DO jk = 2,kjpk_in
+            zwork2(jk,1) = zwork(jk,1)*(ptin(jk,jn)-ptin(jk-1,jn))
+         END DO
+        
+         coeffremap(:,1) = ptin(:,jn)
  
-      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 = 2, kjpk_in-1
+            q001 = phin(jk)*zwork2(jk+1,1)
+            q002 = phin(jk)*zwork2(jk,1)        
+            IF (q001*q002 < 0) then
+               q001 = 0._wp
+               q002 = 0._wp
+            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*phin(jk+1)
+            q002=q002+q*phin(jk-1)
+        
+            coeffremap(jk,3)=coeffremap(jk,1)+q001
+            coeffremap(jk,2)=coeffremap(jk,1)-q002
+        
+            zwork2(jk,1)=(2._wp*q001-q002)**2
+            zwork2(jk,2)=(2._wp*q002-q001)**2
+         ENDDO
+        
+         DO jk = 1, kjpk_in
+            IF(jk.EQ.1 .OR. jk.EQ.kjpk_in .OR. phin(jk).LE.dpthin) THEN
+               coeffremap(jk,3) = coeffremap(jk,1)
+               coeffremap(jk,2) = coeffremap(jk,1)
+               zwork2(jk,1) = 0._wp
+               zwork2(jk,2) = 0._wp
+            ENDIF
+         END DO
+        
+         DO jk = 2, kjpk_in
+            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)
+         END DO
+        
+         zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
+         zwork2(kjpk_in+1,3)=2*coeffremap(kjpk_in,1)-zwork2(kjpk_in,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
+         DO jk = 1, kjpk_in
+            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._wp
+               q02=0._wp
+            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._wp
+         kbot=1
+         DO jk=1,kjpk_out
+            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(kjpk_out+1)) THEN
+ 
+               kbot=ktop
+               DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.kjpk_in)
+                  kbot=kbot+1
+               ENDDO
+               zbox=zbot
+               DO k1= jk+1,kjpk_out
+                  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(kjpk_out+1)) THEN
+                        exit !at bottom
+                     ENDIF
+                  ENDIF
+               ENDDO
+               zthk=zbox-ztop
+
+               kbox=ktop
+               DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.kjpk_in)
+                  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(phin(kbox).GT.dpthin) THEN
+                        q001 = (zbox-z_win(kbox))/phin(kbox)
+                        q002 = (ztop-z_win(kbox))/phin(kbox)
+                        q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
+                        q02=q01-1.+(q001+q002)
+                        q0=1._wp-q01-q02
+                     ELSE
+                        q0  = 1._wp
+                        q01 = 0._wp
+                        q02 = 0._wp
+                     ENDIF
+                     ptout(jk,jn)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
+                  ELSE
+                     ptout(jk,jn) = ptin(kbox,jn)
+                  ENDIF 
                ELSE
-                  zbox=z_wout(k1+1)  !include thin adjacent layers
-                  IF(zbox.EQ.z_wout(Nout+1)) THEN
-                     exit !at bottom
+                  IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
+                     ka = jk
+                  ELSEIF (kbox-ktop.GE.3) THEN
+                     ka = (kbox+ktop)/2
+                  ELSEIF (phin(ktop).GE.phin(kbox)) THEN
+                     ka = ktop
+                  ELSE
+                     ka = kbox
+                  ENDIF !choose ka
+    
+                  offset=coeffremap(ka,1)
+    
+                  qtop = z_win(ktop+1)-ztop !partial layer thickness
+                  IF(phin(ktop).GT.dpthin) THEN
+                     q=(ztop-z_win(ktop))/phin(ktop)
+                     q01=q*(q-1._wp)
+                     q02=q01+q
+                     q0=1._wp-q01-q02            
+                  ELSE
+                     q0  = 1._wp
+                     q01 = 0._wp
+                     q02 = 0._wp
                   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)*phin(k1)
+                  ENDDO !k1
+               
+                  qbot = zbox-z_win(kbox) !partial layer thickness
+                  IF(phin(kbox).GT.dpthin) THEN
+                     q=qbot/phin(kbox)
+                     q01=(q-1._wp)**2
+                     q02=q01-1._wp+q
+                     q0=1_wp-q01-q02                            
+                  ELSE
+                     q0  = 1._wp
+                     q01 = 0._wp
+                     q02 = 0._wp
+                  ENDIF
+              
+                  tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
+               
+                  rpsum=1._wp / zthk
+                  ptout(jk,jn)=offset+tsum*rpsum
+                 
+               ENDIF !single or multiple layers
+            ELSE
+               IF (jk==1) THEN
+                  write(*,'(a7,i4,i4,3f12.5)')'problem = ',kjpk_in,kjpk_out,zthk,z_wout(jk+1),phout(1)
                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)
+               ptout(jk,jn) = ptout(jk-1,jn)
              
-         ENDIF !normal:thin layer
-      ENDDO !jk
+            ENDIF !normal:thin layer
+         ENDDO !jk
+
+      END DO ! loop over variables
             
    END SUBROUTINE reconstructandremap
@@ -260,55 +266,61 @@
 #else
 
-   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, jkin, jkout)
+   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)
       !!----------------------------------------------------------------------
-      !!                    ***  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
+      !!                    *** ROUTINE  reconstructandremap ***
+      !!
+      !! ** Purpose :   Conservative remapping of a vertical column 
+      !!                from one set of layers to an other one.
+      !!
+      !! ** Method  :   Uses D. Engwirda Piecewise Polynomial Reconstruction library.
+      !!                https://github.com/dengwirda/PPR
+      !!                
+      !!
+      !! References :   Engwirda, Darren & Kelley, Maxwell. (2015). A WENO-type 
+      !!                slope-limiter for a family of piecewise polynomial methods. 
+      !!                https://arxiv.org/abs/1606.08188
       !!-----------------------------------------------------------------------
-      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 , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
+      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
+      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
+      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   phin       ! Input thicknesses
+      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   phout      ! Output thicknesses
+      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
+      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Remapped data
       !
-      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)
+      INTEGER, PARAMETER :: ndof = 1
+      INTEGER  :: jk, jn
+      REAL(wp) ::  zwin(kjpk_in+1) ,  ztin(ndof, kn_var, kjpk_in)
+      REAL(wp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out)
       TYPE(rmap_work) :: work
       TYPE(rmap_opts) :: opts
-      TYPE(rcon_ends) :: bc_l(nvar)
-      TYPE(rcon_ends) :: bc_r(nvar)
+      TYPE(rcon_ends) :: bc_l(kn_var)
+      TYPE(rcon_ends) :: bc_r(kn_var)
       !!--------------------------------------------------------------------
      
       ! Set interfaces and input data:
       zwin(1) = 0._wp
-      DO jk = 2, jkin + 1
+      DO jk = 2, kjpk_in + 1
          zwin(jk) = zwin(jk-1) + phin(jk-1) 
-         ztin(ndof, nvar,jk-1) =  ptin(jk-1)
+      END DO
+      
+      DO jn = 1, kn_var 
+         DO jk = 1, kjpk_in
+            ztin(ndof, jn, jk) =  ptin(jk, jn)
+         END DO
       END DO
 
       zwout(1) = 0._wp
-      DO jk = 2, jkout + 1
+      DO jk = 2, kjpk_out + 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   
+      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     ! PQM method in cells
+!      opts%cell_lims = mono_limit     ! monotone limiter   
  
       ! set boundary conditions
@@ -319,9 +331,9 @@
 
       ! init. method workspace
-      CALL work%init(jkin+1, nvar, opts)
+      CALL work%init(kjpk_in+1, kn_var, opts)
 
       ! remap
-      CALL rmap1d(jkin+1, jkout+1, nvar, ndof, &
-      &           zwin, zwout, ztin, ztout,    &
+      CALL rmap1d(kjpk_in+1, kjpk_out+1, kn_var, ndof, &
+      &           zwin, zwout, ztin, ztout,            &
       &           bc_l, bc_r, work, opts)
 
@@ -329,6 +341,8 @@
       CALL work%free()
 
-      DO jk =1, jkout
-         ptout(jk) = ztout(1,1,jk)
+      DO jn = 1, kn_var 
+         DO jk = 1, kjpk_out
+            ptout(jk, jn) = ztout(1, jn, jk)
+         END DO
       END DO
             
