Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8802)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90	(revision 8803)
@@ -1,4 +1,5 @@
 #define TWO_WAY        /* TWO WAY NESTING */
-#undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/
+#undef DECAL_FEEDBACK  /* SEPARATION of INTERFACES*/
+#undef VOL_REFLUX      /* VOLUME REFLUXING*/
  
 MODULE agrif_opa_update
@@ -17,5 +18,5 @@
    PRIVATE
 
-   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl
+   PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl, Agrif_Update_ssh
 
 # if defined key_zdftke
@@ -100,27 +101,29 @@
       CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)  
 # endif
-
+      !
+      nbcline = nbcline + 1
+      ! 
       IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
          ! Update time integrated transports
-         IF (mod(nbcline,nbclineupdate) == 0) THEN
 #  if ! defined DECAL_FEEDBACK
-            CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)
-            CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)
+         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)
+         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
-         ELSE
-#  if ! defined DECAL_FEEDBACK
-            CALL Agrif_Update_Variable(ub2b_update_id,locupdate=(/0,1/),procname = updateub2b)
-            CALL Agrif_Update_Variable(vb2b_update_id,locupdate=(/0,1/),procname = updatevb2b)
-#  else
-            CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateub2b)
-            CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updatevb2b)
-#  endif
-         ENDIF
       END IF
-      !
-      nbcline = nbcline + 1
+#endif
+      !
+   END SUBROUTINE Agrif_Update_Dyn
+
+   SUBROUTINE Agrif_Update_ssh( )
+      !!---------------------------------------------
+      !!   *** ROUTINE Agrif_Update_ssh ***
+      !!---------------------------------------------
+      ! 
+      IF (Agrif_Root()) RETURN
+      !
+#if defined TWO_WAY
       !
       Agrif_UseSpecialValueInUpdate = .TRUE.
@@ -131,9 +134,18 @@
       CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH)
 # endif
+      !
       Agrif_UseSpecialValueInUpdate = .FALSE.
-      ! 
+      !
+#  if defined DECAL_FEEDBACK && defined VOL_REFLUX
+      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN
+         ! Refluxing on ssh:
+         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)
+      END IF
+#  endif
+      !
 #endif
       !
-   END SUBROUTINE Agrif_Update_Dyn
+   END SUBROUTINE Agrif_Update_ssh
 
 # if defined key_zdftke
@@ -180,9 +192,7 @@
       Agrif_SpecialValueFineGrid = 0.
       ! 
-# if ! defined DECAL_FEEDBACK
-      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)
-# else
-      CALL Agrif_Update_Variable(e3t_id, locupdate=(/1,0/), procname=updatee3t)
-# endif 
+      ! No interface separation here, update vertical grid at T points 
+      ! everywhere over the overlapping regions (one account for refluxing in that case):
+      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) 
       !
       Agrif_UseSpecialValueInUpdate = .FALSE.
@@ -209,6 +219,6 @@
 
       ! Save "old" scale factor (prior update) for subsequent asselin correction
-      ! of prognostic variables (needed to update initial state only)
-      ! -------------------------------------------------------------
+      ! of prognostic variables
+      ! -----------------------
       !
       e3u_a(:,:,:) = e3u_n(:,:,:)
@@ -246,7 +256,8 @@
       ! 2) BEFORE fields:
       !------------------
-      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
-         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
-         & .AND.(.NOT.ln_bt_fw)))) THEN
+!      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
+!         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
+!         & .AND.(.NOT.ln_bt_fw)))) THEN
+      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
          !
          ! Vertical scale factor interpolations
@@ -591,13 +602,15 @@
          END DO
       ELSE
-         IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN
-            IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
-               DO jj=j1,j2
-                  DO ji=i1,i2
-                     sshb(ji,jj) =   sshb(ji,jj) &
-                           & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
-                  END DO
-               END DO
-            ENDIF
+         IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
+            & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
+            & .AND.(.NOT.ln_bt_fw)))) THEN
+! tsplit_new
+!         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN
+            DO jj=j1,j2
+               DO ji=i1,i2
+                  sshb(ji,jj) =   sshb(ji,jj) &
+                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1)
+               END DO
+            END DO
          ENDIF
          !
@@ -618,7 +631,48 @@
 
 
-   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before, nb, ndir )
+   SUBROUTINE updateub2b( tabres, i1, i2, j1, j2, before )
       !!---------------------------------------------
       !!          *** ROUTINE updateub2b ***
+      !!---------------------------------------------
+      INTEGER, INTENT(in) :: i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
+      LOGICAL, INTENT(in) :: before
+      !!
+      INTEGER :: ji, jj
+      REAL(wp) :: zrhoy, za1, zcor
+      !!---------------------------------------------
+      !
+      IF (before) THEN
+         zrhoy = Agrif_Rhoy()
+         DO jj=j1,j2
+            DO ji=i1,i2
+               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj)
+            END DO
+         END DO
+         tabres = zrhoy * tabres
+      ELSE
+         ! 
+         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
+         !
+         za1 = 1._wp / REAL(Agrif_rhot(), wp)
+         DO jj=j1,j2
+            DO ji=i1,i2
+               zcor=tabres(ji,jj) - ub2_b(ji,jj)
+               ! Update time integrated fluxes also in case of multiply nested grids:
+               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 
+               ! Update corrective fluxes:
+! tsplit_new
+!               un_bf(ji,jj)  = un_bf(ji,jj) + zcor
+               ! Update half step back fluxes:
+               ub2_b(ji,jj) = tabres(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE updateub2b
+
+   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir )
+      !!---------------------------------------------
+      !!          *** ROUTINE reflux_sshu ***
       !!---------------------------------------------
       INTEGER, INTENT(in) :: i1, i2, j1, j2
@@ -629,5 +683,5 @@
       LOGICAL :: western_side, eastern_side 
       INTEGER :: ji, jj
-      REAL(wp) :: zrhoy, za1
+      REAL(wp) :: zrhoy, za1, zcor
       !!---------------------------------------------
       !
@@ -644,6 +698,4 @@
          tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2)
          !
-         ! Refluxing here:
-#if defined DECAL_FEEDBACK
          western_side  = (nb == 1).AND.(ndir == 1)
          eastern_side  = (nb == 1).AND.(ndir == 2)
@@ -651,35 +703,66 @@
          IF (western_side) THEN
             DO jj=j1,j2
-               sshn(i1  ,jj) = sshn(i1  ,jj) + rdt * r1_e1e2t(i1  ,jj) &
-                     &         * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))
+               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 
+               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor
             END DO
          ENDIF
          IF (eastern_side) THEN
             DO jj=j1,j2
-               sshn(i2+1,jj) = sshn(i2+1,jj) - rdt * r1_e1e2t(i2+1,jj) &
-                     &         * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
-            END DO
-         ENDIF
+               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj))
+               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor
+            END DO
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE reflux_sshu
+
+
+   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before )
+      !!---------------------------------------------
+      !!          *** ROUTINE updatevb2b ***
+      !!---------------------------------------------
+      INTEGER, INTENT(in) :: i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
+      LOGICAL, INTENT(in) :: before
+      !!
+      INTEGER :: ji, jj
+      REAL(wp) :: zrhox, za1, zcor
+      !!---------------------------------------------
+      !
+      IF (before) THEN
+         zrhox = Agrif_Rhox()
+         DO jj=j1,j2
+            DO ji=i1,i2
+               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj) 
+            END DO
+         END DO
+         tabres = zrhox * tabres
+      ELSE
          ! 
-#endif
+         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
+         !
          za1 = 1._wp / REAL(Agrif_rhot(), wp)
-         !
-         DO jj=j1,j2
-            DO ji=i1,i2
+         DO jj=j1,j2
+            DO ji=i1,i2
+               zcor=tabres(ji,jj) - vb2_b(ji,jj)
                ! Update time integrated fluxes also in case of multiply nested grids:
-               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) & 
-                & + za1 * (tabres(ji,jj) - ub2_b(ji,jj))
+               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 
+               ! Update corrective fluxes:
+! tsplit_new
+!               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor
                ! Update half step back fluxes:
-               ub2_b(ji,jj) = tabres(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      !
-   END SUBROUTINE updateub2b
-
-
-   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before, nb, ndir )
-      !!---------------------------------------------
-      !!          *** ROUTINE updatevb2b ***
+               vb2_b(ji,jj) = tabres(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+   END SUBROUTINE updatevb2b
+
+   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir )
+      !!---------------------------------------------
+      !!          *** ROUTINE reflux_sshv ***
       !!---------------------------------------------
       INTEGER, INTENT(in) :: i1, i2, j1, j2
@@ -690,5 +773,5 @@
       LOGICAL :: southern_side, northern_side 
       INTEGER :: ji, jj
-      REAL(wp) :: zrhox, za1
+      REAL(wp) :: zrhox, za1, zcor
       !!---------------------------------------------
       !
@@ -705,6 +788,4 @@
          tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2)
          !
-         ! Refluxing here:
-#if defined DECAL_FEEDBACK
          southern_side = (nb == 2).AND.(ndir == 1)
          northern_side = (nb == 2).AND.(ndir == 2)
@@ -712,30 +793,20 @@
          IF (southern_side) THEN
             DO ji=i1,i2
-               sshn(ji,j1  ) = sshn(ji,j1  ) + rdt * r1_e1e2t(ji,j1  ) &
-                     &         * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
+               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1))
+               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor
             END DO
          ENDIF
          IF (northern_side) THEN               
             DO ji=i1,i2
-               sshn(ji,j2+1) = sshn(ji,j2+1) - rdt * r1_e1e2t(ji,j2+1) &
-                     &         * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
+               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2))
+               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor
+               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor
             END DO
          ENDIF
          ! 
-#endif
-         za1 = 1._wp / REAL(Agrif_rhot(), wp)
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ! Update time integrated fluxes also in case of multiply nested grids:
-               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) & 
-                & + za1 * (tabres(ji,jj) - vb2_b(ji,jj))
-               ! Update half step back fluxes:
-               vb2_b(ji,jj) = tabres(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      !
-   END SUBROUTINE updatevb2b
-
+      ENDIF
+      !
+   END SUBROUTINE reflux_sshv
 
    SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before )
@@ -844,48 +915,46 @@
 # endif /* key_zdftke */ 
 
-   SUBROUTINE updatee3t( ptab, i1, i2, j1, j2, k1, k2, before )
+   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before )
       !!---------------------------------------------
       !!           *** ROUTINE updatee3t ***
       !!---------------------------------------------
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum
       INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
+      LOGICAL, INTENT(in) :: before
+      !
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab
       INTEGER :: ji,jj,jk
       REAL(wp) :: zcoef
       !!---------------------------------------------
       !
-      IF (before) THEN
-!> jc tmp:
-!         ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2)
-         ptab(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2) / e3t_0(i1:i2,j1:j2,k1:k2) * tmask(i1:i2,j1:j2,k1:k2)
-!< jc tmp:
-      ELSE
-         !
-         ! 1) Updates at BEFORE time step:
-         ! -------------------------------
-         !
-!> jc tmp:
+      IF (.NOT.before) THEN
+         !
+         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk)) 
+         !
+         ! Update e3t from ssh (z* case only)
          DO jk = 1, jpkm1
             DO jj=j1,j2
                DO ji=i1,i2
-                  IF (tmask(ji,jj,jk)==1) THEN
-                     ptab(ji,jj,jk) = ptab(ji,jj,jk) * e3t_0(ji,jj,jk)
-                  ELSE
-                     ptab(ji,jj,jk) = e3t_0(ji,jj,jk)
-                  ENDIF
-               END DO
-            END DO
-         END DO
-!         ptab(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
-!< jc tmp:
-
+                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) &
+                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj)))
+               END DO
+            END DO
+         END DO
+         !
+         ! 1) Updates at BEFORE time step:
+         ! -------------------------------
+         !
          ! Save "old" scale factor (prior update) for subsequent asselin correction
-         ! of prognostic variables (needed to update initial state only)
-         e3t_a(i1:i2,j1:j2,k1:k2) = e3t_n(i1:i2,j1:j2,k1:k2)
-!         hdivn(i1:i2,j1:j2,k1:k2)   = e3t_b(i1:i2,j1:j2,k1:k2)
+         ! of prognostic variables
+         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1)
+
+         ! One should also save e3t_b, but lacking of workspace...
+!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1)
 
          IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) &
             & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    &
             & .AND.(.NOT.ln_bt_fw)))) THEN
+! tsplit_new
+!         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN
 
             DO jk = 1, jpkm1
@@ -923,5 +992,5 @@
          !
          ! Update vertical scale factor at T-points:
-         e3t_n(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2)
+         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1)
          !
          ! Update total depth:
@@ -958,4 +1027,5 @@
          ENDIF
          !
+         DEALLOCATE(ptab)
       ENDIF
       !
Index: /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8802)
+++ /branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8803)
@@ -315,4 +315,5 @@
       IF( Agrif_NbStepint() == 0 ) THEN               ! AGRIF Update 
 !!jc in fact update is useless at last time step, but do it for global diagnostics
+                             CALL Agrif_Update_ssh()      ! Update ssh
          IF(.NOT.ln_linssh)  CALL Agrif_Update_vvl()      ! Update vertical scale factors  
                              CALL Agrif_Update_Tra()      ! Update active tracers
