Index: NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce_interp.F90
===================================================================
--- NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce_interp.F90	(revision 11205)
+++ NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/NST/agrif_oce_interp.F90	(revision 11219)
@@ -37,5 +37,5 @@
    PRIVATE
 
-   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
+   PUBLIC   Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_dyn_ts_flux, Agrif_ssh_ts, Agrif_dta_ts
    PUBLIC   Agrif_tra, Agrif_avm
    PUBLIC   interpun , interpvn
@@ -43,6 +43,4 @@
    PUBLIC   interpunb, interpvnb , interpub2b, interpvb2b
    PUBLIC   interpe3t, interpumsk, interpvmsk
-
-   INTEGER ::   bdy_tinterp = 0
 
 #  include "vectopt_loop_substitute.h90"
@@ -78,5 +76,4 @@
       !
       INTEGER ::   ji, jj, jk       ! dummy loop indices
-      INTEGER ::   j1, j2, i1, i2
       INTEGER ::   ibdy1, jbdy1, ibdy2, jbdy2
       REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb
@@ -93,300 +90,246 @@
       Agrif_UseSpecialValue = .FALSE.
       !
-      ! prevent smoothing in ghost cells
-      i1 =  1   ;   i2 = nlci
-      j1 =  1   ;   j2 = nlcj
-      IF( nbondj == -1 .OR. nbondj == 2 )   j1 = 2 + nbghostcells
-      IF( nbondj == +1 .OR. nbondj == 2 )   j2 = nlcj - nbghostcells - 1
-      IF( nbondi == -1 .OR. nbondi == 2 )   i1 = 2 + nbghostcells 
-      IF( nbondi == +1 .OR. nbondi == 2 )   i2 = nlci - nbghostcells - 1
-
       ! --- West --- !
-      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
-         ibdy1 = 2
-         ibdy2 = 1+nbghostcells 
-         !
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            ua_b(ibdy1:ibdy2,:) = 0._wp
+      ibdy1 = 2
+      ibdy2 = 1+nbghostcells 
+      !
+      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+         DO ji = mi0(ibdy1), mi1(ibdy2)
+            ua_b(ji,:) = 0._wp
+
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) & 
-                      & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
-               END DO
-            END DO
+                  ua_b(ji,jj) = ua_b(ji,jj) + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
+               END DO
+            END DO
+
             DO jj = 1, jpj
-               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
-            END DO
-         ENDIF
-         !
-         IF( .NOT.lk_agrif_clp ) THEN
-            DO jk=1,jpkm1              ! Smooth
-               DO jj=j1,j2
-                  ua(ibdy2,jj,jk) = 0.25_wp*(ua(ibdy2-1,jj,jk)+2._wp*ua(ibdy2,jj,jk)+ua(ibdy2+1,jj,jk))
-               END DO
-            END DO
-         ENDIF
-         !
-         zub(ibdy1:ibdy2,:) = 0._wp    ! Correct transport
+               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      DO ji = mi0(ibdy1), mi1(ibdy2)
+         zub(ji,:) = 0._wp    ! Correct transport
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk)*umask(ibdy1:ibdy2,jj,jk)
+               zub(ji,jj) = zub(ji,jj) & 
+                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk)*umask(ji,jj,jk)
             END DO
          END DO
          DO jj=1,jpj
-            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
+            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
          END DO
             
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) &
-                 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj)) * umask(ibdy1:ibdy2,jj,jk)
-            END DO
-         END DO
+               ua(ji,jj,jk) = ( ua(ji,jj,jk) + ua_b(ji,jj)-zub(ji,jj)) * umask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
             
-         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-            zvb(ibdy1:ibdy2,:) = 0._wp
+      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+         DO ji = mi0(ibdy1), mi1(ibdy2)
+            zvb(ji,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) & 
-                     & + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * vmask(ibdy1:ibdy2,jj,jk)
+                  zvb(ji,jj) = zvb(ji,jj) + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
+               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
             END DO
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 
-                    & + va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj))*vmask(ibdy1:ibdy2,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         DO jk = 1, jpkm1              ! Mask domain edges
-            DO jj = 1, jpj
-               ua(1,jj,jk) = 0._wp
-               va(1,jj,jk) = 0._wp
-            END DO
-         END DO 
+                  va(ji,jj,jk) = ( va(ji,jj,jk) + va_b(ji,jj)-zvb(ji,jj))*vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 
       ! --- East --- !
-      IF( nbondi ==  1 .OR. nbondi == 2 ) THEN
-         ibdy1 = nlci-1-nbghostcells
-         ibdy2 = nlci-2 
-         !
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            ua_b(ibdy1:ibdy2,:) = 0._wp
+      ibdy1 = jpiglo-1-nbghostcells
+      ibdy2 = jpiglo-2 
+      !
+      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+         DO ji = mi0(ibdy1), mi1(ibdy2)
+            ua_b(ji,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) & 
-                      & + e3u_a(ibdy1:ibdy2,jj,jk) * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
+                  ua_b(ji,jj) = ua_b(ji,jj) & 
+                      & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               ua_b(ibdy1:ibdy2,jj) = ua_b(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
-            END DO
-         ENDIF
-         !
-         IF( .NOT.lk_agrif_clp ) THEN
-            DO jk=1,jpkm1              ! Smooth
-               DO jj=j1,j2
-                  ua(ibdy1,jj,jk) = 0.25_wp*(ua(ibdy1-1,jj,jk)+2._wp*ua(ibdy1,jj,jk)+ua(ibdy1+1,jj,jk))
-               END DO
-            END DO
-         ENDIF
-         !
-         zub(ibdy1:ibdy2,:) = 0._wp    ! Correct transport
+               ua_b(ji,jj) = ua_b(ji,jj) * r1_hu_a(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      DO ji = mi0(ibdy1), mi1(ibdy2)
+         zub(ji,:) = 0._wp    ! Correct transport
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) & 
-                  & + e3u_a(ibdy1:ibdy2,jj,jk)  * ua(ibdy1:ibdy2,jj,jk) * umask(ibdy1:ibdy2,jj,jk)
+               zub(ji,jj) = zub(ji,jj) & 
+                  & + e3u_a(ji,jj,jk)  * ua(ji,jj,jk) * umask(ji,jj,jk)
             END DO
          END DO
          DO jj=1,jpj
-            zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
+            zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
          END DO
             
          DO jk = 1, jpkm1
             DO jj = 1, jpj
-               ua(ibdy1:ibdy2,jj,jk) = ( ua(ibdy1:ibdy2,jj,jk) & 
-                 & + ua_b(ibdy1:ibdy2,jj)-zub(ibdy1:ibdy2,jj))*umask(ibdy1:ibdy2,jj,jk)
-            END DO
-         END DO
+               ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
+                 & + ua_b(ji,jj)-zub(ji,jj))*umask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
             
-         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-            ibdy1 = ibdy1 + 1
-            ibdy2 = ibdy2 + 1 
-            zvb(ibdy1:ibdy2,:) = 0._wp
+      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+         ibdy1 = jpiglo-nbghostcells
+         ibdy2 = jpiglo-1 
+         DO ji = mi0(ibdy1), mi1(ibdy2)
+            zvb(ji,:) = 0._wp
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) &
-                     & + e3v_a(ibdy1:ibdy2,jj,jk) * va(ibdy1:ibdy2,jj,jk) * vmask(ibdy1:ibdy2,jj,jk)
+                  zvb(ji,jj) = zvb(ji,jj) &
+                     & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
             DO jj = 1, jpj
-               zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
+               zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
             END DO
             DO jk = 1, jpkm1
                DO jj = 1, jpj
-                  va(ibdy1:ibdy2,jj,jk) = ( va(ibdy1:ibdy2,jj,jk) & 
-                      & + va_b(ibdy1:ibdy2,jj)-zvb(ibdy1:ibdy2,jj)) * vmask(ibdy1:ibdy2,jj,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         DO jk = 1, jpkm1              ! Mask domain edges
-            DO jj = 1, jpj
-               ua(nlci-1,jj,jk) = 0._wp
-               va(nlci  ,jj,jk) = 0._wp
-            END DO
-         END DO 
+                  va(ji,jj,jk) = ( va(ji,jj,jk) & 
+                      & + va_b(ji,jj)-zvb(ji,jj)) * vmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 
       ! --- South --- !
-      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
-         jbdy1 = 2
-         jbdy2 = 1+nbghostcells 
-         !
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            va_b(:,jbdy1:jbdy2) = 0._wp
+      jbdy1 = 2
+      jbdy2 = 1+nbghostcells 
+      !
+      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+         DO jj = mj0(jbdy1), mj1(jbdy2)
+            va_b(:,jj) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) & 
-                      & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+                  va_b(ji,jj) = va_b(ji,jj) & 
+                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
             DO ji=1,jpi
-               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
-            END DO
-         ENDIF
-         !
-         IF ( .NOT.lk_agrif_clp ) THEN
-            DO jk = 1, jpkm1           ! Smooth
-               DO ji = i1, i2
-                  va(ji,jbdy2,jk) = 0.25_wp*(va(ji,jbdy2-1,jk)+2._wp*va(ji,jbdy2,jk)+va(ji,jbdy2+1,jk))
-               END DO
-            END DO
-         ENDIF
-         !
-         zvb(:,jbdy1:jbdy2) = 0._wp    ! Correct transport
+               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)     
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jj = mj0(jbdy1), mj1(jbdy2)
+         zvb(:,jj) = 0._wp    ! Correct transport
          DO jk=1,jpkm1
             DO ji=1,jpi
-               zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+               zvb(ji,jj) = zvb(ji,jj) & 
+                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
             END DO
          END DO
          DO ji = 1, jpi
-            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
+            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
          END DO
 
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 
-                 & + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
-            END DO
-         END DO
+               va(ji,jj,jk) = ( va(ji,jj,jk) & 
+                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
             
-         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-            zub(:,jbdy1:jbdy2) = 0._wp
+      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+         DO jj = mj0(jbdy1), mj1(jbdy2)
+            zub(:,jj) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)
+                  zub(ji,jj) = zub(ji,jj) & 
+                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
             DO ji = 1, jpi
-               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
+               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
             END DO
                
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 
-                    & + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         DO jk = 1, jpkm1              ! Mask domain edges
-            DO ji = 1, jpi
-               ua(ji,1,jk) = 0._wp
-               va(ji,1,jk) = 0._wp
-            END DO
-         END DO 
+                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
+                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
 
       ! --- North --- !
-      IF( nbondj ==  1 .OR. nbondj == 2 ) THEN
-         jbdy1 = nlcj-1-nbghostcells
-         jbdy2 = nlcj-2 
-         !
-         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
-            va_b(:,jbdy1:jbdy2) = 0._wp
+      jbdy1 = nlcj-1-nbghostcells
+      jbdy2 = nlcj-2 
+      !
+      IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+         DO jj = mj0(jbdy1), mj1(jbdy2)
+            va_b(:,jj) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) & 
-                      & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+                  va_b(ji,jj) = va_b(ji,jj) & 
+                      & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
                END DO
             END DO
             DO ji=1,jpi
-               va_b(ji,jbdy1:jbdy2) = va_b(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
-            END DO
-         ENDIF
-         !
-         IF ( .NOT.lk_agrif_clp ) THEN
-            DO jk = 1, jpkm1           ! Smooth
-               DO ji = i1, i2
-                  va(ji,jbdy1,jk) = 0.25_wp*(va(ji,jbdy1-1,jk)+2._wp*va(ji,jbdy1,jk)+va(ji,jbdy1+1,jk))
-               END DO
-            END DO
-         ENDIF
-         !
-         zvb(:,jbdy1:jbdy2) = 0._wp    ! Correct transport
+               va_b(ji,jj) = va_b(ji,jj) * r1_hv_a(ji,jj)
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jj = mj0(jbdy1), mj1(jbdy2)
+         zvb(:,jj) = 0._wp    ! Correct transport
          DO jk=1,jpkm1
             DO ji=1,jpi
-               zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) & 
-                  & + e3v_a(ji,jbdy1:jbdy2,jk) * va(ji,jbdy1:jbdy2,jk) * vmask(ji,jbdy1:jbdy2,jk)
+               zvb(ji,jj) = zvb(ji,jj) & 
+                  & + e3v_a(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
             END DO
          END DO
          DO ji = 1, jpi
-            zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
+            zvb(ji,jj) = zvb(ji,jj) * r1_hv_a(ji,jj)
          END DO
 
          DO jk = 1, jpkm1
             DO ji = 1, jpi
-               va(ji,jbdy1:jbdy2,jk) = ( va(ji,jbdy1:jbdy2,jk) & 
-                 & + va_b(ji,jbdy1:jbdy2) - zvb(ji,jbdy1:jbdy2) ) * vmask(ji,jbdy1:jbdy2,jk)
-            END DO
-         END DO
+               va(ji,jj,jk) = ( va(ji,jj,jk) & 
+                 & + va_b(ji,jj) - zvb(ji,jj) ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
             
-         IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
-            jbdy1 = jbdy1 + 1
-            jbdy2 = jbdy2 + 1 
-            zub(:,jbdy1:jbdy2) = 0._wp
+      IF( ln_dynspg_ts ) THEN       ! Set tangential velocities to time splitting estimate
+         jbdy1 = jbdy1 + 1
+         jbdy2 = jbdy2 + 1 
+         DO jj = mj0(jbdy1), mj1(jbdy2)
+            zub(:,jj) = 0._wp
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) & 
-                     & + e3u_a(ji,jbdy1:jbdy2,jk) * ua(ji,jbdy1:jbdy2,jk) * umask(ji,jbdy1:jbdy2,jk)
+                  zub(ji,jj) = zub(ji,jj) & 
+                     & + e3u_a(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
                END DO
             END DO
             DO ji = 1, jpi
-               zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
+               zub(ji,jj) = zub(ji,jj) * r1_hu_a(ji,jj)
             END DO
                
             DO jk = 1, jpkm1
                DO ji = 1, jpi
-                  ua(ji,jbdy1:jbdy2,jk) = ( ua(ji,jbdy1:jbdy2,jk) & 
-                    & + ua_b(ji,jbdy1:jbdy2) - zub(ji,jbdy1:jbdy2) ) * umask(ji,jbdy1:jbdy2,jk)
-               END DO
-            END DO
-         ENDIF
-         !
-         DO jk = 1, jpkm1              ! Mask domain edges
-            DO ji = 1, jpi
-               ua(ji,nlcj  ,jk) = 0._wp
-               va(ji,nlcj-1,jk) = 0._wp
-            END DO
-         END DO 
+                  ua(ji,jj,jk) = ( ua(ji,jj,jk) & 
+                    & + ua_b(ji,jj) - zub(ji,jj) ) * umask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
       ENDIF
       !
@@ -412,17 +355,7 @@
          DO jj=1,jpj
             va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
-            ! Specified fluxes:
             ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
          END DO
       END DO
-      ! Characteristics method (only at boundary point):
-!      istart = 2
-!      iend   = 2
-!      DO ji = mi0(istart), mi1(iend)
-!         DO jj=1,jpj
-!            ua_e(ji,jj) = 0.5_wp * ( ubdy(ji,jj) * hur_e(ji,jj) + ua_e(ji+1,jj) &
-!                        &            - sqrt(grav * hur_e(ji,jj)) * (sshn_e(ji+1,jj) - hbdy(ji,jj)) )
-!         END DO
-!      END DO
       !
       !--- East ---!
@@ -441,13 +374,4 @@
          END DO
       END DO
-      ! Characteristics method (only at boundary point):
-!      istart = jpiglo-2
-!      iend   = jpiglo-2
-!      DO ji = mi0(istart), mi1(iend)
-!         DO jj=1,jpj
-!            ua_e(ji,jj) = 0.5_wp * ( ubdy(ji,jj) * hur_e(ji,jj) + ua_e(ji-1,jj) &
-!                        &            + sqrt(grav * hur_e(ji,jj)) * (sshn_e(ji,jj) - hbdy(ji+1,jj)) )
-!         END DO
-!      END DO
       !
       !--- South ---!
@@ -457,17 +381,7 @@
          DO ji=1,jpi
             ua_e(ji,jj) = ubdy(ji,jj) * hur_e(ji,jj)
-            ! Specified fluxes:
             va_e(ji,jj) = vbdy(ji,jj) * hvr_e(ji,jj)
          END DO
       END DO
-      ! Characteristics method (only at boundary point):
-!      jstart = 2
-!      jend   = 2
-!      DO jj = mj0(jstart), mj1(jend)
-!         DO ji=1,jpi
-!            va_e(ji,jj) = 0.5_wp * ( vbdy(ji,jj) * hvr_e(ji,jj) + va_e(ji,jj+1) &
-!                        &            - sqrt(grav * hvr_e(ji,jj)) * (sshn_e(ji,jj+1) - hbdy(ji,jj)) )
-!         END DO
-!      END DO
       !
       !--- North ---!
@@ -486,16 +400,73 @@
          END DO
       END DO
-      ! Characteristics method (only at boundary point):
-!      jstart = jpjglo-2
-!      jend   = jpjglo-2
-!      DO jj = mj0(jstart), mj1(jend)
-!         DO ji=1,jpi
-!            va_e(ji,jj) = 0.5_wp * ( vbdy(ji,jj) * hvr_e(ji,jj) + va_e(ji,jj-1) &
-!                        &            + sqrt(grav * hvr_e(ji,jj)) * (sshn_e(ji,jj) - hbdy(ji,jj+1)) )
-!         END DO
-!      END DO
       !
    END SUBROUTINE Agrif_dyn_ts
 
+   SUBROUTINE Agrif_dyn_ts_flux( jn, zu, zv )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE Agrif_dyn_ts_flux  ***
+      !!----------------------------------------------------------------------  
+      INTEGER, INTENT(in) ::   jn
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zu, zv
+      !!
+      INTEGER :: ji, jj
+      INTEGER :: istart, iend, jstart, jend
+      !!----------------------------------------------------------------------  
+      !
+      IF( Agrif_Root() )   RETURN
+      !
+      !--- West ---!
+      istart = 2
+      iend   = nbghostcells+1
+      DO ji = mi0(istart), mi1(iend)
+         DO jj=1,jpj
+            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+         END DO
+      END DO
+      !
+      !--- East ---!
+      istart = jpiglo-nbghostcells
+      iend   = jpiglo-1
+      DO ji = mi0(istart), mi1(iend)
+         DO jj=1,jpj
+            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+         END DO
+      END DO
+      istart = jpiglo-nbghostcells-1
+      iend   = jpiglo-2
+      DO ji = mi0(istart), mi1(iend)
+         DO jj=1,jpj
+            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+         END DO
+      END DO
+      !
+      !--- South ---!
+      jstart = 2
+      jend   = nbghostcells+1
+      DO jj = mj0(jstart), mj1(jend)
+         DO ji=1,jpi
+            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+         END DO
+      END DO
+      !
+      !--- North ---!
+      jstart = jpjglo-nbghostcells
+      jend   = jpjglo-1
+      DO jj = mj0(jstart), mj1(jend)
+         DO ji=1,jpi
+            zu(ji,jj) = ubdy(ji,jj) * e2u(ji,jj)
+         END DO
+      END DO
+      jstart = jpjglo-nbghostcells-1
+      jend   = jpjglo-2
+      DO jj = mj0(jstart), mj1(jend)
+         DO ji=1,jpi
+            zv(ji,jj) = vbdy(ji,jj) * e1v(ji,jj)
+         END DO
+      END DO
+      !
+   END SUBROUTINE Agrif_dyn_ts_flux
 
    SUBROUTINE Agrif_dta_ts( kt )
@@ -517,19 +488,23 @@
       !
       ! Interpolate barotropic fluxes
-      Agrif_SpecialValue=0._wp
+      Agrif_SpecialValue = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
+      !
+      ! Set bdy time interpolation stage to 0 (latter incremented locally do deal with corners)
+      utint_stage(:,:) = 0
+      vtint_stage(:,:) = 0
       !
       IF( ll_int_cons ) THEN  ! Conservative interpolation
          ! order matters here !!!!!!
          CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
-         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )
-         bdy_tinterp = 1
+         CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b ) 
+         !
          CALL Agrif_Bc_variable( unb_id        , calledweight=1._wp, procname=interpunb  ) ! After
          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  )         
       ELSE ! Linear interpolation
-         bdy_tinterp = 0
+         !
          ubdy(:,:) = 0._wp   ;   vbdy(:,:) = 0._wp 
          CALL Agrif_Bc_variable( unb_id, procname=interpunb )
@@ -671,5 +646,5 @@
    
 
-   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
+   SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!----------------------------------------------------------------------
       !!                  *** ROUTINE interptsn ***
@@ -678,10 +653,7 @@
       INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
       LOGICAL                                     , INTENT(in   ) ::   before
-      INTEGER                                     , INTENT(in   ) ::   nb , ndir
-      !
-      INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices
-      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
-      REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+      !
+      INTEGER  ::   ji, jj, jk, jn  ! dummy loop indices
+      INTEGER  ::   N_in, N_out
       ! vertical interpolation:
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
@@ -689,5 +661,4 @@
       REAL(wp), DIMENSION(k1:k2) :: h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
-      REAL(wp) :: h_diff
 
       IF( before ) THEN         
@@ -713,16 +684,7 @@
       ELSE 
 
-         western_side  = (nb == 1).AND.(ndir == 1)   ;   eastern_side  = (nb == 1).AND.(ndir == 2)
-         southern_side = (nb == 2).AND.(ndir == 1)   ;   northern_side = (nb == 2).AND.(ndir == 2)
-
 # if defined key_vertical              
          DO jj=j1,j2
             DO ji=i1,i2
-               iref = ji
-               jref = jj
-               if(western_side) iref=MAX(2,ji)
-               if(eastern_side) iref=MIN(nlci-1,ji)
-               if(southern_side) jref=MAX(2,jj)
-               if(northern_side) jref=MIN(nlcj-1,jj)
                N_in = 0
                DO jk=k1,k2 !k2 = jpk of parent grid
@@ -734,7 +696,7 @@
                N_out = 0
                DO jk=1,jpk ! jpk of child grid
-                  IF (tmask(iref,jref,jk) == 0) EXIT 
+                  IF (tmask(ji,jj,jk) == 0) EXIT 
                   N_out = N_out + 1
-                  h_out(jk) = e3t_n(iref,jref,jk)
+                  h_out(jk) = e3t_n(ji,jj,jk)
                ENDDO
                IF (N_in > 0) THEN
@@ -753,132 +715,4 @@
          END DO
 
-         IF ( .NOT.lk_agrif_clp ) THEN 
-            !
-            imin = i1 ; imax = i2
-            jmin = j1 ; jmax = j2
-            ! 
-            ! Remove CORNERS
-            IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells
-            IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1
-            IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells
-            IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1      
-            !
-            IF( eastern_side ) THEN
-               zrho = Agrif_Rhox()
-               z1 = ( zrho - 1._wp ) * 0.5_wp                    
-               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
-               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
-               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
-               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
-               !
-               ibdy = nlci-nbghostcells
-               DO jn = 1, jpts
-                  tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
-                  DO jk = 1, jpkm1
-                     DO jj = jmin,jmax
-                        IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
-                           tsa(ibdy,jj,jk,jn) = tsa(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)
-                        ELSE
-                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy+1,jj,jk,jn)+z3*tsa(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)
-                           IF( un(ibdy-1,jj,jk) > 0._wp ) THEN
-                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy-1,jj,jk,jn)+z5*tsa(ibdy+1,jj,jk,jn) & 
-                                                 + z7*tsa(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tsa(ibdy+1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
-               END DO
-            ENDIF
-            ! 
-            IF( northern_side ) THEN
-               zrho = Agrif_Rhoy()
-               z1 = ( zrho - 1._wp ) * 0.5_wp                    
-               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
-               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
-               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
-               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
-               !
-               jbdy = nlcj-nbghostcells         
-               DO jn = 1, jpts
-                  tsa(imin:imax,jbdy+1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
-                  DO jk = 1, jpkm1
-                     DO ji = imin,imax
-                        IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
-                           tsa(ji,jbdy,jk,jn) = tsa(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy+1,jk,jn)+z3*tsa(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)        
-                           IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN
-                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy-1,jk,jn)+z5*tsa(ji,jbdy+1,jk,jn)  &
-                                                 + z7*tsa(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tsa(imin:imax,jbdy+1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
-               END DO
-            ENDIF
-            !
-            IF( western_side ) THEN
-               zrho = Agrif_Rhox()
-               z1 = ( zrho - 1._wp ) * 0.5_wp                    
-               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
-               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
-               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
-               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
-               !    
-               ibdy = 1+nbghostcells       
-               DO jn = 1, jpts
-                  tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
-                  DO jk = 1, jpkm1
-                     DO jj = jmin,jmax
-                        IF( umask(ibdy,jj,jk) == 0._wp ) THEN
-                           tsa(ibdy,jj,jk,jn) = tsa(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)
-                        ELSE
-                           tsa(ibdy,jj,jk,jn)=(z4*tsa(ibdy-1,jj,jk,jn)+z3*tsa(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)        
-                           IF( un(ibdy,jj,jk) < 0._wp ) THEN
-                              tsa(ibdy,jj,jk,jn)=( z6*tsa(ibdy+1,jj,jk,jn)+z5*tsa(ibdy-1,jj,jk,jn) &
-                                                 + z7*tsa(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tsa(ibdy-1,jmin:jmax,1:jpkm1,jn) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
-               END DO
-            ENDIF
-            !
-            IF( southern_side ) THEN
-               zrho = Agrif_Rhoy()
-               z1 = ( zrho - 1._wp ) * 0.5_wp                    
-               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
-               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
-               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
-               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
-               !  
-               jbdy=1+nbghostcells        
-               DO jn = 1, jpts
-                  tsa(imin:imax,jbdy-1,1:jpkm1,jn) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
-                  DO jk = 1, jpkm1      
-                     DO ji = imin,imax
-                        IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
-                           tsa(ji,jbdy,jk,jn)=tsa(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)
-                        ELSE
-                           tsa(ji,jbdy,jk,jn)=(z4*tsa(ji,jbdy-1,jk,jn)+z3*tsa(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)
-                           IF( vn(ji,jbdy,jk) < 0._wp ) THEN
-                              tsa(ji,jbdy,jk,jn)=( z6*tsa(ji,jbdy+1,jk,jn)+z5*tsa(ji,jbdy-1,jk,jn) & 
-                                                 + z7*tsa(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk)
-                           ENDIF
-                        ENDIF
-                     END DO
-                  END DO
-                  ! Restore ghost points:
-                  tsa(imin:imax,jbdy-1,1:jpkm1,jn) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
-               END DO
-            ENDIF
-            !
-         ENDIF
       ENDIF
       !
@@ -1092,21 +926,30 @@
          zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
          zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot      
-         ! Polynomial interpolation coefficients:
-         IF( bdy_tinterp == 1 ) THEN
-            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-               &               - zt0**2._wp * (       zt0 - 1._wp)        )
-         ELSEIF( bdy_tinterp == 2 ) THEN
-            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-               &               - zt0        * (       zt0 - 1._wp)**2._wp )
-         ELSE
-            ztcoeff = 1
-         ENDIF
-         !   
-         ubdy(i1:i2,j1:j2) = ubdy(i1:i2,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)
-         !            
-         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
-            ubdy(i1:i2,j1:j2) = ubdy(i1:i2,j1:j2) / (zrhoy*e2u(i1:i2,j1:j2)) * umask(i1:i2,j1:j2,1)
-         ENDIF
-      ENDIF
+         ! 
+         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
+                  utint_stage(ji,jj) = utint_stage(ji,jj) + 1
+               ENDIF
+            END DO
+         END DO
+      END IF
       ! 
    END SUBROUTINE interpunb
@@ -1132,19 +975,30 @@
          ! Time indexes bounds for integration
          zt0 = REAL(Agrif_NbStepint()  , wp) / zrhot
-         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot      
-         IF( bdy_tinterp == 1 ) THEN
-            ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-               &               - zt0**2._wp * (       zt0 - 1._wp)        )
-         ELSEIF( bdy_tinterp == 2 ) THEN
-            ztcoeff = zrhot * (  zt1        * (       zt1 - 1._wp)**2._wp &
-               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
-         ELSE
-            ztcoeff = 1
-         ENDIF
-         vbdy(i1:i2,j1:j2) = vbdy(i1:i2,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)
-         !            
-         IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
-            vbdy(i1:i2,j1:j2) = vbdy(i1:i2,j1:j2) / (zrhox*e1v(i1:i2,j1:j2)) * vmask(i1:i2,j1:j2,1)
-         ENDIF            
+         zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot 
+         !     
+         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
+                  vtint_stage(ji,jj) = vtint_stage(ji,jj) + 1
+               ENDIF
+            END DO
+         END DO          
       ENDIF
       !
@@ -1179,4 +1033,7 @@
          !
          ubdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2) 
+         !
+         ! Update interpolation stage:
+         utint_stage(i1:i2,j1:j2) = 1
       ENDIF
       ! 
@@ -1193,5 +1050,5 @@
       !
       INTEGER ::   ji,jj
-      REAL(wp) ::   zrhot, zt0, zt1,zat
+      REAL(wp) ::   zrhot, zt0, zt1, zat
       !!----------------------------------------------------------------------  
       !
@@ -1212,4 +1069,7 @@
          !
          vbdy(i1:i2,j1:j2) = zat * ptab(i1:i2,j1:j2)
+         !
+         ! update interpolation stage:
+         vtint_stage(i1:i2,j1:j2) = 1
       ENDIF
       !      
