Index: branches/2015/dev_merge_2015/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
===================================================================
--- branches/2015/dev_merge_2015/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 6059)
+++ branches/2015/dev_merge_2015/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90	(revision 6060)
@@ -22,16 +22,15 @@
    USE oce
    USE dom_oce      
+   USE zdf_oce
    USE agrif_oce
    USE phycst
+   !
    USE in_out_manager
    USE agrif_opa_sponge
    USE lib_mpp
    USE wrk_nemo
-   USE zdf_oce
  
    IMPLICIT NONE
    PRIVATE
-
-   INTEGER :: bdy_tinterp = 0
 
    PUBLIC   Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
@@ -44,12 +43,12 @@
 # endif
 
-#  include "domzgr_substitute.h90"  
+   INTEGER ::   bdy_tinterp = 0
+
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/NST 3.6 , NEMO Consortium (2010)
+   !! NEMO/NST 3.7 , NEMO Consortium (2015)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -60,9 +59,10 @@
       !
       IF( Agrif_Root() )   RETURN
-
-      Agrif_SpecialValue    = 0.e0
+      !
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
-
+      !
       CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
+      !
       Agrif_UseSpecialValue = .FALSE.
       !
@@ -74,50 +74,47 @@
       !!                  ***  ROUTINE Agrif_DYN  ***
       !!----------------------------------------------------------------------  
-      !! 
       INTEGER, INTENT(in) ::   kt
-      !!
-      INTEGER :: ji,jj,jk, j1,j2, i1,i2
-      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb
-      !!----------------------------------------------------------------------  
-
+      !
+      INTEGER ::   ji, jj, jk       ! dummy loop indices
+      INTEGER ::   j1, j2, i1, i2
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zub, zvb
+      !!----------------------------------------------------------------------  
+      !
       IF( Agrif_Root() )   RETURN
-
-      CALL wrk_alloc( jpi, jpj, zub, zvb )
-
-      Agrif_SpecialValue=0.
+      !
+      CALL wrk_alloc( jpi,jpj,   zub, zvb )
+      !
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
-
-      CALL Agrif_Bc_variable(un_interp_id,procname=interpun)
-      CALL Agrif_Bc_variable(vn_interp_id,procname=interpvn)
-
+      !
+      CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
+      CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
+      !
       Agrif_UseSpecialValue = .FALSE.
- 
+      !
       ! prevent smoothing in ghost cells
-      i1=1
-      i2=jpi
-      j1=1
-      j2=jpj
-      IF((nbondj == -1).OR.(nbondj == 2)) j1 = 3
-      IF((nbondj == +1).OR.(nbondj == 2)) j2 = nlcj-2
-      IF((nbondi == -1).OR.(nbondi == 2)) i1 = 3
-      IF((nbondi == +1).OR.(nbondi == 2)) i2 = nlci-2
-
-
-      IF((nbondi == -1).OR.(nbondi == 2)) THEN
-
+      i1 =  1   ;   i2 = jpi
+      j1 =  1   ;   j2 = jpj
+      IF( nbondj == -1 .OR. nbondj == 2 )   j1 = 3
+      IF( nbondj == +1 .OR. nbondj == 2 )   j2 = nlcj-2
+      IF( nbondi == -1 .OR. nbondi == 2 )   i1 = 3
+      IF( nbondi == +1 .OR. nbondi == 2 )   i2 = nlci-2
+
+      IF( nbondi == -1 .OR. nbondi == 2 ) THEN
+         !
          ! Smoothing
          ! ---------
-         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport
-            ua_b(2,:)=0._wp
-            DO jk=1,jpkm1
-               DO jj=1,jpj
-                  ua_b(2,jj) = ua_b(2,jj) + fse3u_a(2,jj,jk) * ua(2,jj,jk)
-               END DO
-            END DO
-            DO jj=1,jpj
-               ua_b(2,jj) = ua_b(2,jj) * hur_a(2,jj)            
-            END DO
-         ENDIF
-
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            ua_b(2,:) = 0._wp
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  ua_b(2,jj) = ua_b(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
+               END DO
+            END DO
+            DO jj = 1, jpj
+               ua_b(2,jj) = ua_b(2,jj) * r1_hu_a(2,jj)            
+            END DO
+         ENDIF
+         !
          DO jk=1,jpkm1                 ! Smooth
             DO jj=j1,j2
@@ -126,13 +123,13 @@
             END DO
          END DO
-
-         zub(2,:)=0._wp                ! Correct transport
-         DO jk=1,jpkm1
-            DO jj=1,jpj
-               zub(2,jj) = zub(2,jj) + fse3u_a(2,jj,jk) * ua(2,jj,jk)
+         !
+         zub(2,:) = 0._wp              ! Correct transport
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               zub(2,jj) = zub(2,jj) + e3u_a(2,jj,jk) * ua(2,jj,jk)
             END DO
          END DO
          DO jj=1,jpj
-            zub(2,jj) = zub(2,jj) * hur_a(2,jj)
+            zub(2,jj) = zub(2,jj) * r1_hu_a(2,jj)
          END DO
 
@@ -145,68 +142,136 @@
          ! Set tangential velocities to time splitting estimate
          !-----------------------------------------------------
-         IF ( ln_dynspg_ts) THEN
-            zvb(2,:)=0._wp
+         IF( ln_dynspg_ts ) THEN
+            zvb(2,:) = 0._wp
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  zvb(2,jj) = zvb(2,jj) + e3v_a(2,jj,jk) * va(2,jj,jk)
+               END DO
+            END DO
+            DO jj = 1, jpj
+               zvb(2,jj) = zvb(2,jj) * r1_hv_a(2,jj)
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj)) * vmask(2,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
+         ! Mask domain edges:
+         !-------------------
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               ua(1,jj,jk) = 0._wp
+               va(1,jj,jk) = 0._wp
+            END DO
+         END DO         
+         !
+      ENDIF
+
+      IF( nbondi == 1 .OR. nbondi == 2 ) THEN
+
+         ! Smoothing
+         ! ---------
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            ua_b(nlci-2,:) = 0._wp
             DO jk=1,jpkm1
                DO jj=1,jpj
-                  zvb(2,jj) = zvb(2,jj) + fse3v_a(2,jj,jk) * va(2,jj,jk)
+                  ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
                END DO
             END DO
             DO jj=1,jpj
-               zvb(2,jj) = zvb(2,jj) * hvr_a(2,jj)
-            END DO
-            DO jk=1,jpkm1
-               DO jj=1,jpj
-                  va(2,jj,jk) = (va(2,jj,jk)+va_b(2,jj)-zvb(2,jj))*vmask(2,jj,jk)
-               END DO
-            END DO
-         ENDIF
-
+               ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * r1_hu_a(nlci-2,jj)            
+            END DO
+         ENDIF
+
+         DO jk = 1, jpkm1              ! Smooth
+            DO jj = j1, j2
+               ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk)      &
+                  &             * ( ua(nlci-3,jj,jk) + 2._wp*ua(nlci-2,jj,jk) + ua(nlci-1,jj,jk) )
+            END DO
+         END DO
+
+         zub(nlci-2,:) = 0._wp        ! Correct transport
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               zub(nlci-2,jj) = zub(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
+            END DO
+         END DO
+         DO jj = 1, jpj
+            zub(nlci-2,jj) = zub(nlci-2,jj) * r1_hu_a(nlci-2,jj)
+         END DO
+
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               ua(nlci-2,jj,jk) = ( ua(nlci-2,jj,jk) + ua_b(nlci-2,jj) - zub(nlci-2,jj) ) * umask(nlci-2,jj,jk)
+            END DO
+         END DO
+         !
+         ! Set tangential velocities to time splitting estimate
+         !-----------------------------------------------------
+         IF( ln_dynspg_ts ) THEN
+            zvb(nlci-1,:) = 0._wp
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  zvb(nlci-1,jj) = zvb(nlci-1,jj) + e3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk)
+               END DO
+            END DO
+            DO jj=1,jpj
+               zvb(nlci-1,jj) = zvb(nlci-1,jj) * r1_hv_a(nlci-1,jj)
+            END DO
+            DO jk = 1, jpkm1
+               DO jj = 1, jpj
+                  va(nlci-1,jj,jk) = ( va(nlci-1,jj,jk) + va_b(nlci-1,jj) - zvb(nlci-1,jj) ) * vmask(nlci-1,jj,jk)
+               END DO
+            END DO
+         ENDIF
+         !
          ! Mask domain edges:
          !-------------------
-         DO jk=1,jpkm1
-            DO jj=1,jpj
-               ua(1,jj,jk) = 0._wp
-               va(1,jj,jk) = 0._wp
-            END DO
-         END DO         
-
-      ENDIF
-
-      IF((nbondi == 1).OR.(nbondi == 2)) THEN
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               ua(nlci-1,jj,jk) = 0._wp
+               va(nlci  ,jj,jk) = 0._wp
+            END DO
+         END DO 
+         !
+      ENDIF
+
+      IF( nbondj == -1 .OR. nbondj == 2 ) THEN
 
          ! Smoothing
          ! ---------
-         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport
-            ua_b(nlci-2,:)=0._wp
-            DO jk=1,jpkm1
-               DO jj=1,jpj
-                  ua_b(nlci-2,jj) = ua_b(nlci-2,jj) + fse3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
-               END DO
-            END DO
-            DO jj=1,jpj
-               ua_b(nlci-2,jj) = ua_b(nlci-2,jj) * hur_a(nlci-2,jj)            
-            END DO
-         ENDIF
-
-         DO jk=1,jpkm1                 ! Smooth
-            DO jj=j1,j2
-               ua(nlci-2,jj,jk) = 0.25_wp*(ua(nlci-3,jj,jk)+2._wp*ua(nlci-2,jj,jk)+ua(nlci-1,jj,jk))
-               ua(nlci-2,jj,jk) = ua(nlci-2,jj,jk) * umask(nlci-2,jj,jk)
-            END DO
-         END DO
-
-         zub(nlci-2,:)=0._wp           ! Correct transport
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            va_b(:,2) = 0._wp
+            DO jk = 1, jpkm1
+               DO ji = 1, jpi
+                  va_b(ji,2) = va_b(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk)
+               END DO
+            END DO
+            DO ji=1,jpi
+               va_b(ji,2) = va_b(ji,2) * r1_hv_a(ji,2)            
+            END DO
+         ENDIF
+         !
+         DO jk = 1, jpkm1              ! Smooth
+            DO ji = i1, i2
+               va(ji,2,jk) = 0.25_wp * vmask(ji,2,jk)    &
+                  &        * ( va(ji,1,jk) + 2._wp*va(ji,2,jk) + va(ji,3,jk) )
+            END DO
+         END DO
+         !
+         zvb(:,2) = 0._wp              ! Correct transport
          DO jk=1,jpkm1
-            DO jj=1,jpj
-               zub(nlci-2,jj) = zub(nlci-2,jj) + fse3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk)
-            END DO
-         END DO
-         DO jj=1,jpj
-            zub(nlci-2,jj) = zub(nlci-2,jj) * hur_a(nlci-2,jj)
-         END DO
-
-         DO jk=1,jpkm1
-            DO jj=1,jpj
-               ua(nlci-2,jj,jk) = (ua(nlci-2,jj,jk)+ua_b(nlci-2,jj)-zub(nlci-2,jj))*umask(nlci-2,jj,jk)
+            DO ji=1,jpi
+               zvb(ji,2) = zvb(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk)
+            END DO
+         END DO
+         DO ji = 1, jpi
+            zvb(ji,2) = zvb(ji,2) * r1_hv_a(ji,2)
+         END DO
+         DO jk = 1, jpkm1
+            DO ji = 1, jpi
+               va(ji,2,jk) = ( va(ji,2,jk) + va_b(ji,2) - zvb(ji,2) ) * vmask(ji,2,jk)
             END DO
          END DO
@@ -214,17 +279,18 @@
          ! Set tangential velocities to time splitting estimate
          !-----------------------------------------------------
-         IF ( ln_dynspg_ts) THEN
-            zvb(nlci-1,:)=0._wp
-            DO jk=1,jpkm1
-               DO jj=1,jpj
-                  zvb(nlci-1,jj) = zvb(nlci-1,jj) + fse3v_a(nlci-1,jj,jk) * va(nlci-1,jj,jk)
-               END DO
-            END DO
-            DO jj=1,jpj
-               zvb(nlci-1,jj) = zvb(nlci-1,jj) * hvr_a(nlci-1,jj)
-            END DO
-            DO jk=1,jpkm1
-               DO jj=1,jpj
-                  va(nlci-1,jj,jk) = (va(nlci-1,jj,jk)+va_b(nlci-1,jj)-zvb(nlci-1,jj))*vmask(nlci-1,jj,jk)
+         IF( ln_dynspg_ts ) THEN
+            zub(:,2) = 0._wp
+            DO jk = 1, jpkm1
+               DO ji = 1, jpi
+                  zub(ji,2) = zub(ji,2) + e3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk)
+               END DO
+            END DO
+            DO ji = 1, jpi
+               zub(ji,2) = zub(ji,2) * r1_hu_a(ji,2)
+            END DO
+
+            DO jk = 1, jpkm1
+               DO ji = 1, jpi
+                  ua(ji,2,jk) = ( ua(ji,2,jk) + ua_b(ji,2) - zub(ji,2) ) * umask(ji,2,jk)
                END DO
             END DO
@@ -233,8 +299,8 @@
          ! Mask domain edges:
          !-------------------
-         DO jk=1,jpkm1
-            DO jj=1,jpj
-               ua(nlci-1,jj,jk) = 0._wp
-               va(nlci  ,jj,jk) = 0._wp
+         DO jk = 1, jpkm1
+            DO ji = 1, jpi
+               ua(ji,1,jk) = 0._wp
+               va(ji,1,jk) = 0._wp
             END DO
          END DO 
@@ -242,144 +308,77 @@
       ENDIF
 
-      IF((nbondj == -1).OR.(nbondj == 2)) THEN
-
+      IF( nbondj == 1 .OR. nbondj == 2 ) THEN
+         !
          ! Smoothing
          ! ---------
-         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport
-            va_b(:,2)=0._wp
-            DO jk=1,jpkm1
-               DO ji=1,jpi
-                  va_b(ji,2) = va_b(ji,2) + fse3v_a(ji,2,jk) * va(ji,2,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               va_b(ji,2) = va_b(ji,2) * hvr_a(ji,2)            
-            END DO
-         ENDIF
-
-         DO jk=1,jpkm1                 ! Smooth
-            DO ji=i1,i2
-               va(ji,2,jk)=0.25_wp*(va(ji,1,jk)+2._wp*va(ji,2,jk)+va(ji,3,jk))
-               va(ji,2,jk)=va(ji,2,jk)*vmask(ji,2,jk)
-            END DO
-         END DO
-
-         zvb(:,2)=0._wp                ! Correct transport
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               zvb(ji,2) = zvb(ji,2) + fse3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk)
-            END DO
-         END DO
-         DO ji=1,jpi
-            zvb(ji,2) = zvb(ji,2) * hvr_a(ji,2)
-         END DO
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               va(ji,2,jk) = (va(ji,2,jk)+va_b(ji,2)-zvb(ji,2))*vmask(ji,2,jk)
-            END DO
-         END DO
-
+         IF( .NOT.ln_dynspg_ts ) THEN  ! Store transport
+            va_b(:,nlcj-2) = 0._wp
+            DO jk = 1, jpkm1
+               DO ji = 1, jpi
+                  va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk)
+               END DO
+            END DO
+            DO ji = 1, jpi
+               va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)            
+            END DO
+         ENDIF
+         !
+         DO jk = 1, jpkm1              ! Smooth
+            DO ji = i1, i2
+               va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk)   &
+                  &             * ( va(ji,nlcj-3,jk) + 2._wp * va(ji,nlcj-2,jk) + va(ji,nlcj-1,jk) )
+            END DO
+         END DO
+         !
+         zvb(:,nlcj-2) = 0._wp         ! Correct transport
+         DO jk = 1, jpkm1
+            DO ji = 1, jpi
+               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
+            END DO
+         END DO
+         DO ji = 1, jpi
+            zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * r1_hv_a(ji,nlcj-2)
+         END DO
+         DO jk = 1, jpkm1
+            DO ji = 1, jpi
+               va(ji,nlcj-2,jk) = ( va(ji,nlcj-2,jk) + va_b(ji,nlcj-2) - zvb(ji,nlcj-2) ) * vmask(ji,nlcj-2,jk)
+            END DO
+         END DO
+         !
          ! Set tangential velocities to time splitting estimate
          !-----------------------------------------------------
-         IF ( ln_dynspg_ts ) THEN
-            zub(:,2)=0._wp
-            DO jk=1,jpkm1
-               DO ji=1,jpi
-                  zub(ji,2) = zub(ji,2) + fse3u_a(ji,2,jk) * ua(ji,2,jk) * umask(ji,2,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               zub(ji,2) = zub(ji,2) * hur_a(ji,2)
-            END DO
-
-            DO jk=1,jpkm1
-               DO ji=1,jpi
-                  ua(ji,2,jk) = (ua(ji,2,jk)+ua_b(ji,2)-zub(ji,2))*umask(ji,2,jk)
-               END DO
-            END DO
-         ENDIF
-
+         IF( ln_dynspg_ts ) THEN
+            zub(:,nlcj-1) = 0._wp
+            DO jk = 1, jpkm1
+               DO ji = 1, jpi
+                  zub(ji,nlcj-1) = zub(ji,nlcj-1) + e3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk)
+               END DO
+            END DO
+            DO ji = 1, jpi
+               zub(ji,nlcj-1) = zub(ji,nlcj-1) * r1_hu_a(ji,nlcj-1)
+            END DO
+            !
+            DO jk = 1, jpkm1
+               DO ji = 1, jpi
+                  ua(ji,nlcj-1,jk) = ( ua(ji,nlcj-1,jk) + ua_b(ji,nlcj-1) - zub(ji,nlcj-1) ) * umask(ji,nlcj-1,jk)
+               END DO
+            END DO
+         ENDIF
+         !
          ! Mask domain edges:
          !-------------------
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               ua(ji,1,jk) = 0._wp
-               va(ji,1,jk) = 0._wp
-            END DO
-         END DO 
-
-      ENDIF
-
-      IF((nbondj == 1).OR.(nbondj == 2)) THEN
-         ! Smoothing
-         ! ---------
-         IF ( .NOT.ln_dynspg_ts ) THEN ! Store transport
-            va_b(:,nlcj-2)=0._wp
-            DO jk=1,jpkm1
-               DO ji=1,jpi
-                  va_b(ji,nlcj-2) = va_b(ji,nlcj-2) + fse3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               va_b(ji,nlcj-2) = va_b(ji,nlcj-2) * hvr_a(ji,nlcj-2)            
-            END DO
-         ENDIF
-
-         DO jk=1,jpkm1                 ! Smooth
-            DO ji=i1,i2
-               va(ji,nlcj-2,jk)=0.25_wp*(va(ji,nlcj-3,jk)+2._wp*va(ji,nlcj-2,jk)+va(ji,nlcj-1,jk))
-               va(ji,nlcj-2,jk)=va(ji,nlcj-2,jk)*vmask(ji,nlcj-2,jk)
-            END DO
-         END DO
-
-         zvb(:,nlcj-2)=0._wp           ! Correct transport
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + fse3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk)
-            END DO
-         END DO
-         DO ji=1,jpi
-            zvb(ji,nlcj-2) = zvb(ji,nlcj-2) * hvr_a(ji,nlcj-2)
-         END DO
-         DO jk=1,jpkm1
-            DO ji=1,jpi
-               va(ji,nlcj-2,jk) = (va(ji,nlcj-2,jk)+va_b(ji,nlcj-2)-zvb(ji,nlcj-2))*vmask(ji,nlcj-2,jk)
-            END DO
-         END DO
-
-         ! Set tangential velocities to time splitting estimate
-         !-----------------------------------------------------
-         IF ( ln_dynspg_ts ) THEN
-            zub(:,nlcj-1)=0._wp
-            DO jk=1,jpkm1
-               DO ji=1,jpi
-                  zub(ji,nlcj-1) = zub(ji,nlcj-1) + fse3u_a(ji,nlcj-1,jk) * ua(ji,nlcj-1,jk) * umask(ji,nlcj-1,jk)
-               END DO
-            END DO
-            DO ji=1,jpi
-               zub(ji,nlcj-1) = zub(ji,nlcj-1) * hur_a(ji,nlcj-1)
-            END DO
-
-            DO jk=1,jpkm1
-               DO ji=1,jpi
-                  ua(ji,nlcj-1,jk) = (ua(ji,nlcj-1,jk)+ua_b(ji,nlcj-1)-zub(ji,nlcj-1))*umask(ji,nlcj-1,jk)
-               END DO
-            END DO
-         ENDIF
-
-         ! Mask domain edges:
-         !-------------------
-         DO jk=1,jpkm1
-            DO ji=1,jpi
+         DO jk = 1, jpkm1
+            DO ji = 1, jpi
                ua(ji,nlcj  ,jk) = 0._wp
                va(ji,nlcj-1,jk) = 0._wp
             END DO
          END DO 
-
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zub, zvb )
+         !
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi,jpj,   zub, zvb )
       !
    END SUBROUTINE Agrif_dyn
+
 
    SUBROUTINE Agrif_dyn_ts( jn )
@@ -392,7 +391,7 @@
       INTEGER :: ji, jj
       !!----------------------------------------------------------------------  
-
+      !
       IF( Agrif_Root() )   RETURN
-
+      !
       IF((nbondi == -1).OR.(nbondi == 2)) THEN
          DO jj=1,jpj
@@ -405,5 +404,5 @@
          END DO
       ENDIF
-
+      !
       IF((nbondi == 1).OR.(nbondi == 2)) THEN
          DO jj=1,jpj
@@ -416,5 +415,5 @@
          END DO
       ENDIF
-
+      !
       IF((nbondj == -1).OR.(nbondj == 2)) THEN
          DO ji=1,jpi
@@ -427,5 +426,5 @@
          END DO
       ENDIF
-
+      !
       IF((nbondj == 1).OR.(nbondj == 2)) THEN
          DO ji=1,jpi
@@ -441,4 +440,5 @@
    END SUBROUTINE Agrif_dyn_ts
 
+
    SUBROUTINE Agrif_dta_ts( kt )
       !!----------------------------------------------------------------------
@@ -452,47 +452,46 @@
       REAL(wp) :: zrhot, zt
       !!----------------------------------------------------------------------  
-
+      !
       IF( Agrif_Root() )   RETURN
-
-      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in
-      ! the forward case only
-
+      !
+      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
+      !
       zrhot = Agrif_rhot()
-
+      !
       ! "Central" time index for interpolation:
-      IF (ln_bt_fw) THEN
-         zt = REAL(Agrif_NbStepint()+0.5_wp,wp) / zrhot
+      IF( ln_bt_fw ) THEN
+         zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot
       ELSE
-         zt = REAL(Agrif_NbStepint(),wp) / zrhot
-      ENDIF
-
+         zt = REAL( Agrif_NbStepint()       , wp ) / zrhot
+      ENDIF
+      !
       ! Linear interpolation of sea level
-      Agrif_SpecialValue    = 0.e0
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = .TRUE.
-      CALL Agrif_Bc_variable(sshn_id,calledweight=zt, procname=interpsshn )
+      CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn )
       Agrif_UseSpecialValue = .FALSE.
-
+      !
       ! Interpolate barotropic fluxes
       Agrif_SpecialValue=0.
       Agrif_UseSpecialValue = ln_spc_dyn
-
-      IF (ll_int_cons) THEN ! Conservative interpolation
+      !
+      IF( ll_int_cons ) THEN  ! Conservative interpolation
          ! orders 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)
+         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(unb_id ,calledweight=1._wp, procname=interpunb) ! After
-         CALL Agrif_Bc_variable(vnb_id ,calledweight=1._wp, procname=interpvnb)
+         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)         
+         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_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 
-         ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 
-         ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 
-         ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 
-         CALL Agrif_Bc_variable(unb_id,calledweight=zt, procname=interpunb)
-         CALL Agrif_Bc_variable(vnb_id,calledweight=zt, procname=interpvnb)
+         ubdy_w(:) = 0._wp   ;   vbdy_w(:) = 0._wp 
+         ubdy_e(:) = 0._wp   ;   vbdy_e(:) = 0._wp 
+         ubdy_n(:) = 0._wp   ;   vbdy_n(:) = 0._wp 
+         ubdy_s(:) = 0._wp   ;   vbdy_s(:) = 0._wp
+         CALL Agrif_Bc_variable( unb_id, calledweight=zt, procname=interpunb )
+         CALL Agrif_Bc_variable( vnb_id, calledweight=zt, procname=interpvnb )
       ENDIF
       Agrif_UseSpecialValue = .FALSE.
@@ -500,4 +499,5 @@
    END SUBROUTINE Agrif_dta_ts
 
+
    SUBROUTINE Agrif_ssh( kt )
       !!----------------------------------------------------------------------
@@ -507,28 +507,29 @@
       !!
       !!----------------------------------------------------------------------  
-
+      !
       IF( Agrif_Root() )   RETURN
-
+      !
       IF((nbondi == -1).OR.(nbondi == 2)) THEN
          ssha(2,:)=ssha(3,:)
          sshn(2,:)=sshn(3,:)
       ENDIF
-
+      !
       IF((nbondi == 1).OR.(nbondi == 2)) THEN
          ssha(nlci-1,:)=ssha(nlci-2,:)
          sshn(nlci-1,:)=sshn(nlci-2,:)
       ENDIF
-
+      !
       IF((nbondj == -1).OR.(nbondj == 2)) THEN
          ssha(:,2)=ssha(:,3)
          sshn(:,2)=sshn(:,3)
       ENDIF
-
+      !
       IF((nbondj == 1).OR.(nbondj == 2)) THEN
          ssha(:,nlcj-1)=ssha(:,nlcj-2)
          sshn(:,nlcj-1)=sshn(:,nlcj-2)
       ENDIF
-
+      !
    END SUBROUTINE Agrif_ssh
+
 
    SUBROUTINE Agrif_ssh_ts( jn )
@@ -540,32 +541,33 @@
       INTEGER :: ji,jj
       !!----------------------------------------------------------------------  
-
+      !
       IF((nbondi == -1).OR.(nbondi == 2)) THEN
-         DO jj=1,jpj
+         DO jj = 1, jpj
             ssha_e(2,jj) = hbdy_w(jj)
          END DO
       ENDIF
-
+      !
       IF((nbondi == 1).OR.(nbondi == 2)) THEN
-         DO jj=1,jpj
+         DO jj = 1, jpj
             ssha_e(nlci-1,jj) = hbdy_e(jj)
          END DO
       ENDIF
-
+      !
       IF((nbondj == -1).OR.(nbondj == 2)) THEN
-         DO ji=1,jpi
+         DO ji = 1, jpi
             ssha_e(ji,2) = hbdy_s(ji)
          END DO
       ENDIF
-
+      !
       IF((nbondj == 1).OR.(nbondj == 2)) THEN
-         DO ji=1,jpi
+         DO ji = 1, jpi
             ssha_e(ji,nlcj-1) = hbdy_n(ji)
          END DO
       ENDIF
-
+      !
    END SUBROUTINE Agrif_ssh_ts
 
 # if defined key_zdftke
+
    SUBROUTINE Agrif_tke
       !!----------------------------------------------------------------------
@@ -573,33 +575,36 @@
       !!----------------------------------------------------------------------  
       REAL(wp) ::   zalpha
+      !!----------------------------------------------------------------------  
       !
       zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp )
       IF( zalpha > 1. )   zalpha = 1.
-      
+      !
       Agrif_SpecialValue    = 0.e0
       Agrif_UseSpecialValue = .TRUE.
-      
+      !
       CALL Agrif_Bc_variable(avm_id ,calledweight=zalpha, procname=interpavm)       
-              
+      !
       Agrif_UseSpecialValue = .FALSE.
       !
    END SUBROUTINE Agrif_tke
+   
 # endif
 
-   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, nb, ndir )
+      !!----------------------------------------------------------------------
       !!   *** ROUTINE interptsn ***
-      !!---------------------------------------------
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
+      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   ! dummy loop indices
-      INTEGER :: imin, imax, jmin, jmax
+      INTEGER  ::   imin, imax, jmin, jmax
       REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3
       REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
-
+      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
+      !!----------------------------------------------------------------------
+      !
       IF (before) THEN         
          ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
@@ -634,14 +639,14 @@
          IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2        
          !
-         IF( eastern_side) THEN
+         IF( eastern_side ) THEN
             DO jn = 1, jpts
                tsa(nlci,j1:j2,k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
-                     IF( umask(nlci-2,jj,jk) == 0.e0 ) THEN
+                     IF( umask(nlci-2,jj,jk) == 0._wp ) THEN
                         tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
                      ELSE
                         tsa(nlci-1,jj,jk,jn)=(zalpha4*tsa(nlci,jj,jk,jn)+zalpha3*tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
-                        IF( un(nlci-2,jj,jk) > 0.e0 ) THEN
+                        IF( un(nlci-2,jj,jk) > 0._wp ) THEN
                            tsa(nlci-1,jj,jk,jn)=( zalpha6*tsa(nlci-2,jj,jk,jn)+zalpha5*tsa(nlci,jj,jk,jn) & 
                                  + zalpha7*tsa(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk)
@@ -651,5 +656,5 @@
                END DO
                tsa(nlci,j1:j2,k1:k2,jn) = 0._wp
-            ENDDO
+            END DO
          ENDIF
          ! 
@@ -659,9 +664,9 @@
                DO jk = 1, jpkm1
                   DO ji = imin,imax
-                     IF( vmask(ji,nlcj-2,jk) == 0.e0 ) THEN
+                     IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN
                         tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
                      ELSE
                         tsa(ji,nlcj-1,jk,jn)=(zalpha4*tsa(ji,nlcj,jk,jn)+zalpha3*tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)        
-                        IF (vn(ji,nlcj-2,jk) > 0.e0 ) THEN
+                        IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN
                            tsa(ji,nlcj-1,jk,jn)=( zalpha6*tsa(ji,nlcj-2,jk,jn)+zalpha5*tsa(ji,nlcj,jk,jn)  &
                                  + zalpha7*tsa(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk)
@@ -671,17 +676,17 @@
                END DO
                tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp
-            ENDDO
-         ENDIF
-         !
-         IF( western_side) THEN            
+            END DO
+         ENDIF
+         !
+         IF( western_side ) THEN            
             DO jn = 1, jpts
                tsa(1,j1:j2,k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)
                DO jk = 1, jpkm1
                   DO jj = jmin,jmax
-                     IF( umask(2,jj,jk) == 0.e0 ) THEN
+                     IF( umask(2,jj,jk) == 0._wp ) THEN
                         tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
                      ELSE
                         tsa(2,jj,jk,jn)=(zalpha4*tsa(1,jj,jk,jn)+zalpha3*tsa(3,jj,jk,jn))*tmask(2,jj,jk)        
-                        IF( un(2,jj,jk) < 0.e0 ) THEN
+                        IF( un(2,jj,jk) < 0._wp ) THEN
                            tsa(2,jj,jk,jn)=(zalpha6*tsa(3,jj,jk,jn)+zalpha5*tsa(1,jj,jk,jn)+zalpha7*tsa(4,jj,jk,jn))*tmask(2,jj,jk)
                         ENDIF
@@ -696,11 +701,11 @@
             DO jn = 1, jpts
                tsa(i1:i2,1,k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)
-               DO jk=1,jpk      
+               DO jk = 1, jpk      
                   DO ji=imin,imax
-                     IF( vmask(ji,2,jk) == 0.e0 ) THEN
+                     IF( vmask(ji,2,jk) == 0._wp ) THEN
                         tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
                      ELSE
                         tsa(ji,2,jk,jn)=(zalpha4*tsa(ji,1,jk,jn)+zalpha3*tsa(ji,3,jk,jn))*tmask(ji,2,jk)
-                        IF( vn(ji,2,jk) < 0.e0 ) THEN
+                        IF( vn(ji,2,jk) < 0._wp ) THEN
                            tsa(ji,2,jk,jn)=(zalpha6*tsa(ji,3,jk,jn)+zalpha5*tsa(ji,1,jk,jn)+zalpha7*tsa(ji,4,jk,jn))*tmask(ji,2,jk)
                         ENDIF
@@ -709,5 +714,5 @@
                END DO
                tsa(i1:i2,1,k1:k2,jn) = 0._wp
-            ENDDO
+            END DO
          ENDIF
          !
@@ -735,12 +740,14 @@
    END SUBROUTINE interptsn
 
-   SUBROUTINE interpsshn(ptab,i1,i2,j1,j2,before,nb,ndir)
+
+   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpsshn  ***
       !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) :: i1,i2,j1,j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      INTEGER                         , INTENT(in   ) ::   nb , ndir
+      !
       LOGICAL :: western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
@@ -761,32 +768,26 @@
    END SUBROUTINE interpsshn
 
-   SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2, before)
-      !!---------------------------------------------
+
+   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before )
+      !!----------------------------------------------------------------------
       !!   *** ROUTINE interpun ***
-      !!---------------------------------------------    
-      !!
-      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
-      !!
-      INTEGER :: ji,jj,jk
-      REAL(wp) :: zrhoy	
-      !!---------------------------------------------    
-      !
-      IF (before) THEN 
-         DO jk=1,jpk
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk)
-                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u(ji,jj,jk)
-               END DO
-            END DO
+      !!----------------------------------------------------------------------
+      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
+      !
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zrhoy	
+      !!----------------------------------------------------------------------
+      !
+      IF( before ) THEN 
+         DO jk = k1, jpk
+            ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)
          END DO
       ELSE
          zrhoy = Agrif_Rhoy()
-         DO jk=1,jpkm1
+         DO jk = 1, jpkm1
             DO jj=j1,j2
-               ua(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj)))
-               ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u(i1:i2,jj,jk)
+               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) )
             END DO
          END DO
@@ -795,59 +796,47 @@
    END SUBROUTINE interpun
 
-   SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before)
-      !!---------------------------------------------
+
+   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before )
+      !!----------------------------------------------------------------------
       !!   *** ROUTINE interpvn ***
-      !!---------------------------------------------    
-      !
-      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
-      !
-      INTEGER :: ji,jj,jk
-      REAL(wp) :: zrhox	
-      !!---------------------------------------------    
+      !!----------------------------------------------------------------------
+      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
+      !
+      INTEGER  ::   ji, jj, jk
+      REAL(wp) ::   zrhox	
+      !!----------------------------------------------------------------------
       !      
-      IF (before) THEN          
-         !interpv entre 1 et k2 et interpv2d en jpkp1
-         DO jk=k1,jpk
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk)
-                  ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v(ji,jj,jk)
-               END DO
-            END DO
+      IF( before ) THEN       !interpv entre 1 et k2 et interpv2d en jpkp1
+         DO jk = k1, jpk
+            ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk)
          END DO
       ELSE          
          zrhox= Agrif_Rhox()
-         DO jk=1,jpkm1
-            DO jj=j1,j2
-               va(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj)))
-               va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v(i1:i2,jj,jk)
-            END DO
+         DO jk = 1, jpkm1
+            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) )
          END DO
       ENDIF
       !        
    END SUBROUTINE interpvn
-
-   SUBROUTINE interpunb(ptab,i1,i2,j1,j2,before,nb,ndir)
+   
+
+   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpunb  ***
       !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) :: i1,i2,j1,j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !!
-      INTEGER :: ji,jj
-      REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
-      !!----------------------------------------------------------------------  
-      !
-      IF (before) THEN 
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ptab(ji,jj) = un_b(ji,jj) * e2u(ji,jj) * hu(ji,jj) 
-            END DO
-         END DO
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      INTEGER                         , INTENT(in   ) ::   nb , ndir
+      !
+      INTEGER  ::   ji, jj
+      REAL(wp) ::   zrhoy, zrhot, zt0, zt1, ztcoeff
+      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
+      !!----------------------------------------------------------------------  
+      !
+      IF( before ) THEN 
+         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -863,8 +852,8 @@
          IF( bdy_tinterp == 1 ) THEN
             ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-                  &      - zt0**2._wp * (       zt0 - 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 ) 
+               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
 
          ELSE
@@ -887,18 +876,14 @@
          IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
             IF(western_side) THEN
-               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
-                     &                                  * umask(i1,j1:j2,1)
+               ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
             ENDIF
             IF(eastern_side) THEN
-               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2))   &
-                     &                                  * umask(i1,j1:j2,1)
+               ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoy*e2u(i1,j1:j2)) * umask(i1,j1:j2,1)
             ENDIF
             IF(southern_side) THEN
-               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
-                     &                                  * umask(i1:i2,j1,1)
+               ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
             ENDIF
             IF(northern_side) THEN
-               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1))   &
-                     &                                  * umask(i1:i2,j1,1)
+               ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoy*e2u(i1:i2,j1)) * umask(i1:i2,j1,1)
             ENDIF
          ENDIF
@@ -907,24 +892,21 @@
    END SUBROUTINE interpunb
 
-   SUBROUTINE interpvnb(ptab,i1,i2,j1,j2,before,nb,ndir)
+
+   SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpvnb  ***
       !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) :: i1,i2,j1,j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !!
-      INTEGER :: ji,jj
-      REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff   
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      INTEGER                         , INTENT(in   ) ::   nb , ndir
+      !
+      INTEGER  ::   ji,jj
+      REAL(wp) ::   zrhox, zrhot, zt0, zt1, ztcoeff   
+      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
       ! 
-      IF (before) THEN 
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ptab(ji,jj) = vn_b(ji,jj) * e1v(ji,jj) * hv(ji,jj) 
-            END DO
-         END DO
+      IF( before ) THEN 
+         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -939,9 +921,8 @@
          IF( bdy_tinterp == 1 ) THEN
             ztcoeff = zrhot * (  zt1**2._wp * (       zt1 - 1._wp)        &
-                  &      - zt0**2._wp * (       zt0 - 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 ) 
-
+               &               - zt0        * (       zt0 - 1._wp)**2._wp ) 
          ELSE
             ztcoeff = 1
@@ -983,23 +964,20 @@
    END SUBROUTINE interpvnb
 
-   SUBROUTINE interpub2b(ptab,i1,i2,j1,j2,before,nb,ndir)
+
+   SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpub2b  ***
       !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) :: i1,i2,j1,j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !!
-      INTEGER :: ji,jj
-      REAL(wp) :: zrhot, zt0, zt1,zat
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      INTEGER                         , INTENT(in   ) ::   nb , ndir
+      !
+      INTEGER  ::   ji,jj
+      REAL(wp) ::   zrhot, zt0, zt1,zat
+      LOGICAL  ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
       IF( before ) THEN
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ptab(ji,jj) = ub2_b(ji,jj) * e2u(ji,jj)
-            END DO
-         END DO
+         ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1012,6 +990,6 @@
          zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
          ! Polynomial interpolation coefficients:
-         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
-               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
+         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
+            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
          ! 
          IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
@@ -1022,25 +1000,22 @@
       ! 
    END SUBROUTINE interpub2b
-
-   SUBROUTINE interpvb2b(ptab,i1,i2,j1,j2,before,nb,ndir)
+   
+
+   SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpvb2b  ***
       !!----------------------------------------------------------------------  
-      INTEGER, INTENT(in) :: i1,i2,j1,j2
-      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
-      LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !!
-      INTEGER :: ji,jj
-      REAL(wp) :: zrhot, zt0, zt1,zat
-      LOGICAL :: western_side, eastern_side,northern_side,southern_side
+      INTEGER                         , INTENT(in   ) ::   i1, i2, j1, j2
+      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) ::   ptab
+      LOGICAL                         , INTENT(in   ) ::   before
+      INTEGER                         , INTENT(in   ) ::   nb , ndir
+      !
+      INTEGER ::   ji,jj
+      REAL(wp) ::   zrhot, zt0, zt1,zat
+      LOGICAL ::   western_side, eastern_side,northern_side,southern_side
       !!----------------------------------------------------------------------  
       !
       IF( before ) THEN
-         DO jj=j1,j2
-            DO ji=i1,i2
-               ptab(ji,jj) = vb2_b(ji,jj) * e1v(ji,jj)
-            END DO
-         END DO
+         ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
       ELSE      
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1053,24 +1028,24 @@
          zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
          ! Polynomial interpolation coefficients:
-         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)        &
-               &      - zt0**2._wp * (-2._wp*zt0 + 3._wp)        ) 
-         !
-         IF(western_side ) vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(eastern_side ) vbdy_e(j1:j2) = zat * ptab(i1,j1:j2)  
-         IF(southern_side) vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
-         IF(northern_side) vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
+         zat = zrhot * (  zt1**2._wp * (-2._wp*zt1 + 3._wp)    &
+            &           - zt0**2._wp * (-2._wp*zt0 + 3._wp)    ) 
+         !
+         IF(western_side )   vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)  
+         IF(eastern_side )   vbdy_e(j1:j2) = zat * ptab(i1,j1:j2)  
+         IF(southern_side)   vbdy_s(i1:i2) = zat * ptab(i1:i2,j1) 
+         IF(northern_side)   vbdy_n(i1:i2) = zat * ptab(i1:i2,j1) 
       ENDIF
       !      
    END SUBROUTINE interpvb2b
 
-   SUBROUTINE interpe3t(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
+
+   SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpe3t  ***
       !!----------------------------------------------------------------------  
-      ! 
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
+      INTEGER                              , INTENT(in   ) :: i1, i2, j1, j2, k1, k2
       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL :: before
-      INTEGER, INTENT(in) :: nb , ndir
+      LOGICAL                              , INTENT(in   ) :: before
+      INTEGER                              , INTENT(in   ) :: nb , ndir
       !
       INTEGER :: ji, jj, jk
@@ -1079,12 +1054,6 @@
       !!----------------------------------------------------------------------  
       !    
-      IF (before) THEN
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  ptab(ji,jj,jk) = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+      IF( before ) THEN
+         ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
       ELSE
          western_side  = (nb == 1).AND.(ndir == 1)
@@ -1093,14 +1062,14 @@
          northern_side = (nb == 2).AND.(ndir == 2)
 
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
+         DO jk = k1, k2
+            DO jj = j1, j2
+               DO ji = i1, i2
                   ! Get velocity mask at boundary edge points:
-                  IF (western_side)  ztmpmsk = umask(ji    ,jj    ,1)
-                  IF (eastern_side)  ztmpmsk = umask(nlci-2,jj    ,1)
-                  IF (northern_side) ztmpmsk = vmask(ji    ,nlcj-2,1)
-                  IF (southern_side) ztmpmsk = vmask(ji    ,2     ,1)
-
-                  IF (ABS(ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk))*ztmpmsk > 1.D-2) THEN
+                  IF( western_side )   ztmpmsk = umask(ji    ,jj    ,1)
+                  IF( eastern_side )   ztmpmsk = umask(nlci-2,jj    ,1)
+                  IF( northern_side)   ztmpmsk = vmask(ji    ,nlcj-2,1)
+                  IF( southern_side)   ztmpmsk = vmask(ji    ,2     ,1)
+                  !
+                  IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
                      IF (western_side) THEN
                         WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
@@ -1112,5 +1081,5 @@
                         WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
                      ENDIF
-                     WRITE(numout,*) '      ptab(ji,jj,jk), fse3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
+                     WRITE(numout,*) '      ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
                      kindic_agr = kindic_agr + 1
                   ENDIF
@@ -1118,38 +1087,31 @@
             END DO
          END DO
-
+         !
       ENDIF
       ! 
    END SUBROUTINE interpe3t
 
-   SUBROUTINE interpumsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
+
+   SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpumsk  ***
       !!----------------------------------------------------------------------  
-      ! 
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !
-      INTEGER :: ji, jj, jk
-      LOGICAL :: western_side, eastern_side   
+      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
+      INTEGER                              , INTENT(in   ) ::   nb , ndir
+      !
+      INTEGER ::   ji, jj, jk
+      LOGICAL ::   western_side, eastern_side   
       !!----------------------------------------------------------------------  
       !    
-      IF (before) THEN
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  ptab(ji,jj,jk) = umask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+      IF( before ) THEN
+         ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
       ELSE
-
-         western_side  = (nb == 1).AND.(ndir == 1)
-         eastern_side  = (nb == 1).AND.(ndir == 2)
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
+         western_side = (nb == 1).AND.(ndir == 1)
+         eastern_side = (nb == 1).AND.(ndir == 2)
+         DO jk = k1, k2
+            DO jj = j1, j2
+               DO ji = i1, i2
                    ! Velocity mask at boundary edge points:
                   IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
@@ -1167,38 +1129,31 @@
             END DO
          END DO
-
+         !
       ENDIF
       ! 
    END SUBROUTINE interpumsk
 
-   SUBROUTINE interpvmsk(ptab,i1,i2,j1,j2,k1,k2,before,nb,ndir)
+
+   SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interpvmsk  ***
       !!----------------------------------------------------------------------  
-      ! 
-      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
-      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
-      LOGICAL :: before
-      INTEGER, INTENT(in) :: nb , ndir
-      !
-      INTEGER :: ji, jj, jk
-      LOGICAL :: northern_side, southern_side     
+      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
+      INTEGER                              , INTENT(in   ) :: nb , ndir
+      !
+      INTEGER ::   ji, jj, jk
+      LOGICAL ::   northern_side, southern_side     
       !!----------------------------------------------------------------------  
       !    
-      IF (before) THEN
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
-                  ptab(ji,jj,jk) = vmask(ji,jj,jk)
-               END DO
-            END DO
-         END DO
+      IF( before ) THEN
+         ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
       ELSE
-
          southern_side = (nb == 2).AND.(ndir == 1)
          northern_side = (nb == 2).AND.(ndir == 2)
-         DO jk=k1,k2
-            DO jj=j1,j2
-               DO ji=i1,i2
+         DO jk = k1, k2
+            DO jj = j1, j2
+               DO ji = i1, i2
                    ! Velocity mask at boundary edge points:
                   IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
@@ -1216,5 +1171,5 @@
             END DO
          END DO
-
+         !
       ENDIF
       ! 
@@ -1223,14 +1178,14 @@
 # if defined key_zdftke
 
-   SUBROUTINE interpavm(ptab,i1,i2,j1,j2,k1,k2,before)
+   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE interavm  ***
       !!----------------------------------------------------------------------  
-      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
+      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
       !!----------------------------------------------------------------------  
       !      
-      IF( before) THEN
+      IF( before ) THEN
          ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
       ELSE
