Index: NEMO/trunk/src/NST/agrif_oce_update.F90
===================================================================
--- NEMO/trunk/src/NST/agrif_oce_update.F90	(revision 14064)
+++ NEMO/trunk/src/NST/agrif_oce_update.F90	(revision 14086)
@@ -1,3 +1,3 @@
-#undef DECAL_FEEDBACK     /* SEPARATION of INTERFACES */
+#undef DECAL_FEEDBACK    /* SEPARATION of INTERFACES */
 #undef DECAL_FEEDBACK_2D  /* SEPARATION of INTERFACES (Barotropic mode) */
 #undef VOL_REFLUX         /* VOLUME REFLUXING*/
@@ -21,4 +21,5 @@
    USE zdf_oce        ! vertical physics: ocean variables 
    USE agrif_oce
+   USE dom_oce
    !
    USE in_out_manager ! I/O manager
@@ -34,5 +35,5 @@
 
    PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh
-   PUBLIC   Update_Scales
+   PUBLIC   Update_Scales, Agrif_Check_parent_bat
 
    !! * Substitutions
@@ -54,25 +55,20 @@
       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed()
 
-#if defined key_vertical
-! Effect of this has to be carrefully checked 
-! depending on what the nesting tools ensure for
-! volume conservation:
+      l_vremap                      = ln_vert_remap
+      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap
+      Agrif_SpecialValueFineGrid    = 0._wp
+      ! 
+# if ! defined DECAL_FEEDBACK
+      CALL Agrif_Update_Variable(ts_update_id, procname=updateTS)
+! near boundary update:
+!      CALL Agrif_Update_Variable(ts_update_id,locupdate=(/0,2/), procname=updateTS)
+# else
+      CALL Agrif_Update_Variable(ts_update_id, locupdate=(/1,0/),procname=updateTS)
+! near boundary update:
+!      CALL Agrif_Update_Variable(ts_update_id,locupdate=(/1,2/), procname=updateTS)
+# endif
+      !
       Agrif_UseSpecialValueInUpdate = .FALSE.
-#else
-      Agrif_UseSpecialValueInUpdate = .TRUE.
-#endif
-      Agrif_SpecialValueFineGrid    = 0._wp
-      ! 
-# if ! defined DECAL_FEEDBACK
-      CALL Agrif_Update_Variable(tsn_id, procname=updateTS)
-! near boundary update:
-!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS)
-# else
-      CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS)
-! near boundary update:
-!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS)
-# endif
-      !
-      Agrif_UseSpecialValueInUpdate = .FALSE.
+      l_vremap                      = .FALSE.
       !
       !
@@ -90,9 +86,27 @@
       Agrif_UseSpecialValueInUpdate = .FALSE.
       Agrif_SpecialValueFineGrid    = 0._wp
-
-      use_sign_north = .TRUE.
-      sign_north     = -1._wp
-
-      !     
+      l_vremap                      = ln_vert_remap
+      use_sign_north                = .TRUE.
+      sign_north                    = -1._wp     
+!
+# if ! defined DECAL_FEEDBACK_2D
+      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateU2d)
+      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d)
+# else
+      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateU2d)
+      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d)  
+# endif
+      ! 
+      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
+         ! Update time integrated transports
+#  if ! defined DECAL_FEEDBACK_2D
+         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateub2b)
+         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b)
+#  else
+         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateub2b)
+         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b)
+#  endif
+      END IF
+
 # if ! defined DECAL_FEEDBACK
       CALL Agrif_Update_Variable(un_update_id,procname = updateU)
@@ -108,33 +122,7 @@
 !      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV)
 # endif
-
-# if ! defined DECAL_FEEDBACK_2D
-      CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)
-      CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)  
-# else
-      CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)
-      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)  
-# endif
-      !
-# if ! defined DECAL_FEEDBACK_2D
-      ! Account for updated thicknesses at boundary edges
-      IF (.NOT.ln_linssh) THEN
-!         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)
-!         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)
-      ENDIF
-# endif
-      ! 
-      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN
-         ! Update time integrated transports
-#  if ! defined DECAL_FEEDBACK_2D
-         CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
-         CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
-#  else
-         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)
-         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)
-#  endif
-      END IF
       !
       use_sign_north = .FALSE.
+      l_vremap = .FALSE.
       !
    END SUBROUTINE Agrif_Update_Dyn
@@ -150,7 +138,7 @@
       Agrif_SpecialValueFineGrid = 0._wp
 # if ! defined DECAL_FEEDBACK_2D
-      CALL Agrif_Update_Variable(sshn_id,procname = updateSSH)
+      CALL Agrif_Update_Variable(sshn_id,locupdate=(/  nn_shift_bar,-2/), procname = updateSSH) 
 # else
-      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
+      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1+nn_shift_bar,-2/),procname = updateSSH)
 # endif
       !
@@ -158,14 +146,14 @@
       !
 #  if defined VOL_REFLUX
-      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
+      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
          use_sign_north = .TRUE.
          sign_north = -1._wp
          ! Refluxing on ssh:
 #  if defined DECAL_FEEDBACK_2D
-         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)
-         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv)
+         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/nn_shift_bar,nn_shift_bar/),locupdate2=(/1+nn_shift_bar,1+nn_shift_bar/),procname = reflux_sshu)
+         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,1+nn_shift_bar/),locupdate2=(/nn_shift_bar,nn_shift_bar/),procname = reflux_sshv)
 #  else
-         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu)
-         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv)
+         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1+nn_shift_bar,-1+nn_shift_bar/),locupdate2=(/nn_shift_bar, nn_shift_bar/),procname = reflux_sshu)
+         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ nn_shift_bar, nn_shift_bar/),locupdate2=(/-1+nn_shift_bar,-1+nn_shift_bar/),procname = reflux_sshv)
 #  endif
          use_sign_north = .FALSE.
@@ -320,5 +308,4 @@
 #endif
 
-#if defined key_vertical
 
    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
@@ -340,166 +327,131 @@
       !
       IF (before) THEN
-!jc_alt
-!         AGRIF_SpecialValue = -999._wp
-         DO jn = n1,n2-1
+         IF ( l_vremap ) THEN
+            DO jn = n1,n2-1
+               DO jk=k1,k2
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
+                     END DO
+                  END DO
+               END DO
+            END DO
             DO jk=k1,k2
                DO jj=j1,j2
                   DO ji=i1,i2
-!jc_alt
-!                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) &
-!                                         &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp
-                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)
+                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
                   END DO
                END DO
             END DO
-         END DO
-         DO jk=k1,k2
+         ELSE
+            DO jn = 1,jpts
+               DO jk=k1,k2
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
+                     END DO
+                  END DO
+               END DO
+            END DO
+
+         ENDIF
+      ELSE
+         IF ( l_vremap ) THEN
+            tabres_child(:,:,:,:) = 0._wp
+            AGRIF_SpecialValue = 0._wp
             DO jj=j1,j2
                DO ji=i1,i2
-!jc_alt
-!                  tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) &
-!                                      &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp
-                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
-               END DO
-            END DO
-         END DO
-      ELSE
-         tabres_child(:,:,:,:) = 0._wp
-         AGRIF_SpecialValue = 0._wp
-         DO jj=j1,j2
-            DO ji=i1,i2
-               N_in = 0
-               DO jk=k1,k2 !k2 = jpk of child grid
-! jc_alt
-!                  IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT
-                  IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT
-                  N_in = N_in + 1
-                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
-                  h_in(N_in) = tabres(ji,jj,jk,n2)
+                  N_in = 0
+                  DO jk=k1,k2 !k2 = jpk of child grid
+                     IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp  ) EXIT
+                     N_in = N_in + 1
+                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
+                     h_in(N_in) = tabres(ji,jj,jk,n2)
+                  ENDDO
+                  N_out = 0
+                  DO jk=1,jpk ! jpk of parent grid
+                     IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
+                     N_out = N_out + 1
+                     h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 
+                  ENDDO
+                  IF (N_in*N_out > 0) THEN !Remove this?
+                     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
-               N_out = 0
-               DO jk=1,jpk ! jpk of parent grid
-                  IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF
-                  N_out = N_out + 1
-                  h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 
-               ENDDO
-               IF (N_in*N_out > 0) THEN !Remove this?
-                  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
-         ENDDO
-
-         IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
-            ! Add asselin part
+
+            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
+               ! Add asselin part
+               DO jn = 1,jpts
+                  DO jk = 1, jpkm1
+                     DO jj = j1, j2
+                        DO ji = i1, i2
+                           IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
+                              ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
+                              ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
+                              ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
+                              ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
+                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
+                           ENDIF
+                        END DO
+                     END DO
+                  END DO
+               END DO
+            ENDIF
             DO jn = 1,jpts
                DO jk = 1, jpkm1
                   DO jj = j1, j2
                      DO ji = i1, i2
-                        IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN
-                           ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
-                           ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a)
-                           ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
-                        ENDIF
+                        IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 
+                           ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
+                        END IF
                      END DO
                   END DO
                END DO
             END DO
-         ENDIF
-         DO jn = 1,jpts
-            DO jk = 1, jpkm1
-               DO jj = j1, j2
-                  DO ji = i1, i2
-                     IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 
-                        ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn)
-                     END IF
+         ELSE
+            DO jn = 1,jpts
+               tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
+                                            & * tmask(i1:i2,j1:j2,k1:k2)
+            ENDDO
+ 
+            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
+               ! Add asselin part
+               DO jn = 1,jpts
+                  DO jk = k1, k2
+                     DO jj = j1, j2
+                        DO ji = i1, i2
+                           IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
+                              ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
+                              ztnu = tabres(ji,jj,jk,jn)
+                              ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
+                              ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
+                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
+                           ENDIF
+                        END DO
+                     END DO
                   END DO
                END DO
-            END DO
-         END DO
-         !
+            ENDIF
+            DO jn = 1,jpts
+               DO jk=k1,k2
+                  DO jj=j1,j2
+                     DO ji=i1,i2
+                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
+                           ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
+                        END IF
+                     END DO
+                  END DO
+               END DO
+            END DO
+            !
+         ENDIF
          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
             ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a)
-         ENDIF
+         ENDIF         
       ENDIF
       ! 
    END SUBROUTINE updateTS
 
-# else
-
-   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updateT ***
-      !!---------------------------------------------
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL, INTENT(in) :: before
-      !!
-      INTEGER :: ji,jj,jk,jn
-      REAL(wp) :: ztb, ztnu, ztno
-      !!---------------------------------------------
-      !
-      IF (before) THEN
-         DO jn = 1,jpts
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
-!> jc tmp
-                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)
-!                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a)
-!< jc tmp
-                  END DO
-               END DO
-            END DO
-         END DO
-      ELSE
-!> jc tmp
-         DO jn = 1,jpts
-            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &
-                                         & * tmask(i1:i2,j1:j2,k1:k2)
-         ENDDO
-!< jc tmp
-         IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN
-            ! Add asselin part
-            DO jn = 1,jpts
-               DO jk = k1, k2
-                  DO jj = j1, j2
-                     DO ji = i1, i2
-                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN
-                           ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
-                           ztnu = tabres(ji,jj,jk,jn)
-                           ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)
-                           ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  & 
-                                     &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)
-                        ENDIF
-                     END DO
-                  END DO
-               END DO
-            END DO
-         ENDIF
-         DO jn = 1,jpts
-            DO jk=k1,k2
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
-                        ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)
-                     END IF
-                  END DO
-               END DO
-            END DO
-         END DO
-         !
-         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
-            ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)
-         ENDIF
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE updateTS
-
-# endif
-
-# if defined key_vertical
 
    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
@@ -513,10 +465,11 @@
       INTEGER ::   ji, jj, jk
       REAL(wp)::   zrhoy, zub, zunu, zuno
+      REAL(wp), DIMENSION(jpi,jpj) ::   zpgu   ! 2D workspace
 ! VERTICAL REFINEMENT BEGIN
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
       REAL(wp) :: h_in(k1:k2)
       REAL(wp) :: h_out(1:jpk)
-      INTEGER  :: N_in, N_out
-      REAL(wp) :: h_diff, excess, thick
+      INTEGER  :: N_in, N_out, N_in_save, N_out_save
+      REAL(wp) :: zhmin, zd
       REAL(wp) :: tabin(k1:k2)
 ! VERTICAL REFINEMENT END
@@ -525,65 +478,78 @@
       IF( before ) THEN
          zrhoy = Agrif_Rhoy()
-!jc_alt
-!         AGRIF_SpecialValue = -999._wp
          DO jk=k1,k2
+            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) & 
+                                     &   * umask(i1:i2,j1:j2,jk) * uu(i1:i2,j1:j2,jk,Kmm_a)  
+         END DO
+
+         IF ( l_vremap ) THEN
+            DO jk=k1,k2
+               tabres(i1:i2,j1:j2,jk,2) = zrhoy * umask(i1:i2,j1:j2,jk) * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a)
+            END DO
+         ENDIF
+
+      ELSE
+
+         tabres_child(:,:,:) = 0._wp
+
+         IF ( l_vremap ) THEN
+
             DO jj=j1,j2
                DO ji=i1,i2
-!jc_alt
-!                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)  &
-!                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp
-                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)  
-!jc_alt
-!                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  &
-!                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp
-                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
-               END DO
-            END DO
-         END DO
-      ELSE
-         tabres_child(:,:,:) = 0.
-         AGRIF_SpecialValue = 0._wp
-         DO jj=j1,j2
-            DO ji=i1,i2
-               N_in = 0
-               h_in(:) = 0._wp
-               tabin(:) = 0._wp
-               DO jk=k1,k2 !k2=jpk of child grid
-!jc_alt
-!                  IF( tabres(ji,jj,jk,2) < -900._wp) EXIT
-                  IF( tabres(ji,jj,jk,2) == 0.) EXIT
-                  N_in = N_in + 1
-                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
-                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)
+                  N_in = 0
+                  h_in(:) = 0._wp
+                  tabin(:) = 0._wp
+                  DO jk=k1,k2 !k2=jpk of child grid
+                     IF( tabres(ji,jj,jk,2)*r1_e2u(ji,jj) <= 1.e-6_wp ) EXIT
+                     N_in = N_in + 1
+                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+                     h_in(N_in) = tabres(ji,jj,jk,2) * r1_e2u(ji,jj)
+                  ENDDO
+                  N_out = 0
+                  DO jk=1,jpk
+                     IF (umask(ji,jj,jk) == 0._wp) EXIT
+                     N_out = N_out + 1
+                     h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
+                  ENDDO
+                  IF (N_in * N_out > 0) THEN
+                     ! Deal with potentially different depths at velocity points:
+                     N_in_save  = N_in
+                     N_out_save = N_out
+                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN
+                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in)))
+                        zd = 0._wp
+                        DO jk=1, N_in_save
+                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN
+                              N_in = jk
+                              h_in(jk) = zhmin - zd
+                              EXIT 
+                           ENDIF
+                           zd = zd + h_in(jk)
+                        END DO
+                        zd = 0._wp
+                        DO jk=1, N_out_save
+                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN
+                              N_out = jk
+                              h_out(jk) = zhmin - zd
+                              EXIT 
+                           ENDIF
+                           zd = zd + h_out(jk)
+                        END DO
+                     END IF
+                     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)
+                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out)
+                  ENDIF
                ENDDO
-               N_out = 0
-               DO jk=1,jpk
-                  IF (umask(ji,jj,jk) == 0) EXIT
-                  N_out = N_out + 1
-                  h_out(N_out) = e3u(ji,jj,jk,Kmm_a)
-               ENDDO
-               IF (N_in * N_out > 0) THEN
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-                  excess = 0._wp
-                  IF (h_diff < -1.e-4) THEN
-!Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 
-!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
-                     DO jk=N_in,1,-1
-                        thick = MIN(-1*h_diff, h_in(jk))
-                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
-                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
-                        h_diff = h_diff + thick
-                        IF ( h_diff == 0) THEN
-                           N_in = jk
-                           h_in(jk) = h_in(jk) - thick
-                           EXIT
-                        ENDIF
-                     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,1)
-                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))
-               ENDIF
             ENDDO
-         ENDDO
+
+         ELSE
+            DO jk=1,jpk
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) /  e3u(ji,jj,jk,Kmm_a)
+                  END DO
+               END DO
+            END DO
+         ENDIF
          !
          DO jk=1,jpk
@@ -603,4 +569,30 @@
          END DO
          !
+         ! Correct now and before transports:
+         DO jj=j1,j2
+            DO ji=i1,i2
+               zpgu(ji,jj) = 0._wp
+               DO jk=1,jpkm1
+                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
+               END DO
+               !
+               DO jk=1,jpkm1              
+                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + &
+                      &  (uu_b(ji,jj,Kmm_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kmm_a)) * umask(ji,jj,jk)           
+               END DO
+               !
+               zpgu(ji,jj) = 0._wp
+               DO jk=1,jpkm1
+                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
+               END DO
+               !
+               DO jk=1,jpkm1              
+                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + &
+                      &  (uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)) * umask(ji,jj,jk)           
+               END DO
+               !
+            END DO
+         END DO
+         !
          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
             uu(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a)
@@ -611,114 +603,22 @@
    END SUBROUTINE updateu
 
-#else
-
-   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updateu ***
-      !!---------------------------------------------
-      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL                                     , INTENT(in   ) :: before
-      !
-      INTEGER  :: ji, jj, jk
-      REAL(wp) :: zrhoy, zub, zunu, zuno
-      !!---------------------------------------------
-      ! 
-      IF( before ) THEN
-         zrhoy = Agrif_Rhoy()
-         DO jk = k1, k2
-            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a)
-         END DO
-      ELSE
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 
-                  !
-                  IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
-                     zub  = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a)  ! fse3t_b prior update should be used
-                     zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a)
-                     zunu = tabres(ji,jj,jk,1)
-                     uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &      
-                                    & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)
-                  ENDIF
-                  !
-                  uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a)
-               END DO
-            END DO
-         END DO
-         !
-         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
-            uu(i1:i2,j1:j2,k1:k2,Kbb_a)  = uu(i1:i2,j1:j2,k1:k2,Kmm_a)
-         ENDIF
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE updateu
-
-# endif
-
-   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
-      !!---------------------------------------------
-      !!           *** ROUTINE correct_u_bdy ***
+
+   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
+      !!---------------------------------------------
+      !!           *** ROUTINE updatev ***
       !!---------------------------------------------
       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
       LOGICAL                                     , INTENT(in   ) :: before
-      INTEGER                                     , INTENT(in)    :: nb, ndir
-      !!
-      LOGICAL :: western_side, eastern_side 
-      !
-      INTEGER  :: jj, jk
-      REAL(wp) :: zcor
-      !!---------------------------------------------
-      ! 
-      IF( .NOT.before ) THEN
-         !
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         !
-         IF (western_side) THEN
-            DO jj=j1,j2
-               zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a)
-               uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor
-               DO jk=1,jpkm1
-                  uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk)
-               END DO 
-            END DO
-         ENDIF
-         !
-         IF (eastern_side) THEN
-            DO jj=j1,j2
-               zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a)
-               uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor
-               DO jk=1,jpkm1
-                  uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk)
-               END DO 
-            END DO
-         ENDIF
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE correct_u_bdy
-
-# if  defined key_vertical
-
-   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
-      !!---------------------------------------------
-      !!           *** ROUTINE updatev ***
-      !!---------------------------------------------
-      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL                                     , INTENT(in   ) :: before
       !
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zrhox, zvb, zvnu, zvno
+      REAL(wp), DIMENSION(jpi,jpj) ::   zpgv   ! 2D workspace
 ! VERTICAL REFINEMENT BEGIN
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
       REAL(wp) :: h_in(k1:k2)
       REAL(wp) :: h_out(1:jpk)
-      INTEGER :: N_in, N_out
-      REAL(wp) :: h_diff, excess, thick
+      INTEGER  :: N_in, N_out, N_in_save, N_out_save
+      REAL(wp) :: zhmin, zd
       REAL(wp) :: tabin(k1:k2)
 ! VERTICAL REFINEMENT END
@@ -727,63 +627,76 @@
       IF( before ) THEN
          zrhox = Agrif_Rhox()
-!jc_alt
-!         AGRIF_SpecialValue = -999._wp
          DO jk=k1,k2
+            tabres(i1:i2,j1:j2,jk,1) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) & 
+                                     &   * vmask(i1:i2,j1:j2,jk) * vv(i1:i2,j1:j2,jk,Kmm_a)  
+         END DO
+
+         IF ( l_vremap ) THEN
+            DO jk=k1,k2
+               tabres(i1:i2,j1:j2,jk,2) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk)
+            END DO
+         ENDIF
+
+      ELSE
+
+         tabres_child(:,:,:) = 0._wp
+
+         IF ( l_vremap ) THEN
+
             DO jj=j1,j2
                DO ji=i1,i2
-!jc_alt
-!                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) &
-!                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) 
-!jc_alt
-!                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) &
-!                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp
-                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
-      ELSE
-         tabres_child(:,:,:) = 0.
-         AGRIF_SpecialValue = 0._wp
-         DO jj=j1,j2
-            DO ji=i1,i2
-               N_in = 0
-               DO jk=k1,k2
-!jc_alt
-!                  IF (tabres(ji,jj,jk,2) < -900._wp) EXIT
-                  IF (tabres(ji,jj,jk,2) == 0) EXIT
-                  N_in = N_in + 1
-                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
-                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)
+                  N_in = 0
+                  DO jk=k1,k2
+                     IF (tabres(ji,jj,jk,2)* r1_e1v(ji,jj) <= 1.e-6_wp) EXIT
+                     N_in = N_in + 1
+                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)
+                     h_in(N_in) = tabres(ji,jj,jk,2) * r1_e1v(ji,jj)
+                  ENDDO
+                  N_out = 0
+                  DO jk=1,jpk
+                     IF (vmask(ji,jj,jk) == 0) EXIT
+                     N_out = N_out + 1
+                     h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
+                  ENDDO
+                  IF (N_in * N_out > 0) THEN
+                     ! Deal with potentially different depths at velocity points:
+                     N_in_save  = N_in
+                     N_out_save = N_out
+                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN
+                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in)))
+                        zd = 0._wp
+                        DO jk=1, N_in_save
+                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN
+                              N_in = jk
+                              h_in(jk) = zhmin - zd
+                              EXIT 
+                           ENDIF
+                           zd = zd + h_in(jk)
+                        END DO
+                        zd = 0._wp
+                        DO jk=1, N_out_save
+                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN
+                              N_out = jk
+                              h_out(jk) = zhmin - zd
+                              EXIT 
+                           ENDIF
+                           zd = zd + h_out(jk)
+                        END DO
+                     END IF
+                     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)
+                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out)
+                  ENDIF
                ENDDO
-               N_out = 0
-               DO jk=1,jpk
-                  IF (vmask(ji,jj,jk) == 0) EXIT
-                  N_out = N_out + 1
-                  h_out(N_out) = e3v(ji,jj,jk,Kmm_a)
-               ENDDO
-               IF (N_in * N_out > 0) THEN
-                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
-                  excess = 0._wp
-                  IF (h_diff < -1.e-4) then
-!Even if bathy at T points match it's possible for the V points to be deeper in the child grid. 
-!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.
-                     DO jk=N_in,1,-1
-                        thick = MIN(-1*h_diff, h_in(jk))
-                        excess = excess + tabin(jk)*thick*e2u(ji,jj)
-                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))
-                        h_diff = h_diff + thick
-                        IF ( h_diff == 0) THEN
-                           N_in = jk
-                           h_in(jk) = h_in(jk) - thick
-                           EXIT
-                        ENDIF
-                     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,1)
-                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))
-               ENDIF
             ENDDO
-         ENDDO
+
+         ELSE
+            DO jk=1,jpk
+               DO jj=j1,j2
+                  DO ji=i1,i2
+                     tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) /  e3v(ji,jj,jk,Kmm_a)
+                  END DO
+               END DO
+            END DO
+         ENDIF
          !
          DO jk=1,jpkm1
@@ -803,4 +716,30 @@
          END DO
          !
+         ! Correct now and before transports:
+         DO jj=j1,j2
+            DO ji=i1,i2
+               zpgv(ji,jj) = 0._wp
+               DO jk=1,jpkm1
+                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
+               END DO
+               !
+               DO jk=1,jpkm1              
+                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + &
+                      &  (vv_b(ji,jj,Kmm_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kmm_a)) * vmask(ji,jj,jk)           
+               END DO
+               !
+               zpgv(ji,jj) = 0._wp
+               DO jk=1,jpkm1
+                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
+               END DO
+               !
+               DO jk=1,jpkm1              
+                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + &
+                      &  (vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)) * vmask(ji,jj,jk)           
+               END DO
+               !
+            END DO
+         END DO
+         !
          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
             vv(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a)
@@ -810,100 +749,4 @@
       ! 
    END SUBROUTINE updatev
-
-# else
-
-   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
-      !!---------------------------------------------
-      !!           *** ROUTINE updatev ***
-      !!---------------------------------------------
-      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL                                     , INTENT(in   ) :: before
-      !
-      INTEGER  :: ji, jj, jk
-      REAL(wp) :: zrhox, zvb, zvnu, zvno
-      !!---------------------------------------------      
-      !
-      IF (before) THEN
-         zrhox = Agrif_Rhox()
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
-               END DO
-            END DO
-         END DO
-      ELSE
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)
-                  !
-                  IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part
-                     zvb  = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used
-                     zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)
-                     zvnu = tabres(ji,jj,jk,1)
-                     vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &      
-                                    & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)
-                  ENDIF
-                  !
-                  vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a)
-               END DO
-            END DO
-         END DO
-         !
-         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN
-            vv(i1:i2,j1:j2,k1:k2,Kbb_a)  = vv(i1:i2,j1:j2,k1:k2,Kmm_a)
-         ENDIF
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE updatev
-
-# endif
-
-   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
-      !!---------------------------------------------
-      !!           *** ROUTINE correct_v_bdy ***
-      !!---------------------------------------------
-      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
-      LOGICAL                                     , INTENT(in   ) :: before
-      INTEGER                                     , INTENT(in)    :: nb, ndir
-      !!
-      LOGICAL :: southern_side, northern_side 
-      !
-      INTEGER  :: ji, jk
-      REAL(wp) :: zcor
-      !!---------------------------------------------
-      ! 
-      IF( .NOT.before ) THEN
-         !
-         southern_side = (nb == 2).AND.(ndir == 1)
-         northern_side = (nb == 2).AND.(ndir == 2)
-         !
-         IF (southern_side) THEN
-            DO ji=i1,i2
-               zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a)
-               vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor
-               DO jk=1,jpkm1
-                  vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk)
-               END DO 
-            END DO
-         ENDIF
-         !
-         IF (northern_side) THEN
-            DO ji=i1,i2
-               zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a)
-               vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor
-               DO jk=1,jpkm1
-                  vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk)
-               END DO 
-            END DO
-         ENDIF
-         !
-      ENDIF
-      ! 
-   END SUBROUTINE correct_v_bdy
 
 
@@ -935,15 +778,4 @@
                tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)  
                !    
-               ! Update "now" 3d velocities:
-               zpgu(ji,jj) = 0._wp
-               DO jk=1,jpkm1
-                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)
-               END DO
-               !
-               zcorr = (tabres(ji,jj) - zpgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)
-               DO jk=1,jpkm1              
-                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)           
-               END DO
-               !
                ! Update barotropic velocities:
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
@@ -955,15 +787,4 @@
                uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1)
                !       
-               ! Correct "before" velocities to hold correct bt component:
-               zpgu(ji,jj) = 0.e0
-               DO jk=1,jpkm1
-                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)
-               END DO
-               !
-               zcorr = uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)
-               DO jk=1,jpkm1              
-                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)           
-               END DO
-               !
             END DO
          END DO
@@ -1002,16 +823,4 @@
             DO ji=i1,i2
                tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)  
-               !    
-               ! Update "now" 3d velocities:
-               zpgv(ji,jj) = 0.e0
-               DO jk=1,jpkm1
-                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)
-               END DO
-               !
-               zcorr = (tabres(ji,jj) - zpgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)
-               DO jk=1,jpkm1              
-                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)           
-               END DO
-               !
                ! Update barotropic velocities:
                IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN
@@ -1023,15 +832,4 @@
                vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1)
                !       
-               ! Correct "before" velocities to hold correct bt component:
-               zpgv(ji,jj) = 0.e0
-               DO jk=1,jpkm1
-                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)
-               END DO
-               !
-               zcorr = vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)
-               DO jk=1,jpkm1              
-                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)           
-               END DO
-               !
             END DO
          END DO
@@ -1120,5 +918,5 @@
                ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
                ! Update corrective fluxes:
-               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) un_bf(ji,jj)  = un_bf(ji,jj) + zcor
                ! Update half step back fluxes:
                ub2_b(ji,jj) = tabres(ji,jj)
@@ -1208,5 +1006,5 @@
                vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
                ! Update corrective fluxes:
-               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler)))  vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
                ! Update half step back fluxes:
                vb2_b(ji,jj) = tabres(ji,jj)
@@ -1479,4 +1277,63 @@
 #endif
 
+   SUBROUTINE Agrif_Check_parent_bat( )
+      !!----------------------------------------------------------------------
+      !!                   *** ROUTINE Agrif_Check_parent_bat ***
+      !!----------------------------------------------------------------------
+      ! 
+      IF (( .NOT.ln_agrif_2way ).OR.(.NOT.ln_chk_bathy).OR.(Agrif_Root())) RETURN
+      !
+      Agrif_UseSpecialValueInUpdate = .FALSE.
+      !
+      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*) 'AGRIF: Check parent volume at Level:', Agrif_Level()
+      !
+# if ! defined DECAL_FEEDBACK
+      CALL Agrif_Update_Variable(batupd_id,procname = update_bat)
+# else
+      CALL Agrif_Update_Variable(batupd_id,locupdate=(/1,0/),procname = update_bat)
+# endif
+      !
+      kindic_agr = Agrif_Parent(kindic_agr)
+      CALL mpp_sum( 'Agrif_Check_parent_bat', kindic_agr )
+
+      IF( kindic_agr /= 0 ) THEN
+         CALL ctl_stop('==> Averaged Bathymetry does not match parent volume') 
+      ELSE
+         IF(lwp) WRITE(numout,*) '==> Averaged Bathymetry matches parent ' 
+         IF(lwp) WRITE(numout,*) ''
+      ENDIF
+      !
+   END SUBROUTINE Agrif_Check_parent_bat
+
+   SUBROUTINE update_bat(ptab, i1, i2, j1, j2, before )
+      !!---------------------------------------------
+      !!           *** ROUTINE update_bat ***
+      !!---------------------------------------------
+      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ptab
+      INTEGER, INTENT(in) :: i1, i2, j1, j2
+      LOGICAL, INTENT(in) :: before
+      INTEGER :: ji, jj
+      !
+      !!---------------------------------------------
+      !
+      IF( before ) THEN
+         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
+      ELSE
+         kindic_agr = 0
+         !
+         DO jj=j1,j2
+            DO ji=i1,i2
+               IF ( (ssmask(ji,jj).NE.0._wp).AND.&
+               & (ABS(ptab(ji,jj)-ht_0(ji,jj)).GE.1.e-6) ) THEN 
+                  kindic_agr = kindic_agr + 1 
+               ENDIF
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE update_bat
+
 #else
    !!----------------------------------------------------------------------
