Index: NEMO/branches/2020/tickets_icb_1900/src/NST/agrif_oce_sponge.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/NST/agrif_oce_sponge.F90	(revision 13277)
+++ NEMO/branches/2020/tickets_icb_1900/src/NST/agrif_oce_sponge.F90	(revision 13899)
@@ -78,8 +78,8 @@
       zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
 
-      Agrif_SpecialValue=0.
+      Agrif_SpecialValue    = 0._wp
       Agrif_UseSpecialValue = ln_spc_dyn
-      use_sign_north = .TRUE.
-      sign_north = -1.
+      use_sign_north        = .TRUE.
+      sign_north            = -1._wp
       !
       tabspongedone_u = .FALSE.
@@ -92,5 +92,5 @@
       !
       Agrif_UseSpecialValue = .FALSE.
-      use_sign_north = .FALSE.
+      use_sign_north        = .FALSE.
 #endif
       !
@@ -109,4 +109,8 @@
       REAL(wp) ::   z1_ispongearea, z1_jspongearea
       REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
+#if defined key_vertical
+      REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu
+      REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv
+#endif
       REAL(wp), DIMENSION(jpjmax)  :: zmskwest,  zmskeast
       REAL(wp), DIMENSION(jpimax)  :: zmsknorth, zmsksouth
@@ -129,50 +133,39 @@
          ! Retrieve masks at open boundaries:
 
-         ! --- West --- !
-         IF( lk_west) THEN
+         IF( lk_west ) THEN                             ! --- West --- !
             ztabramp(:,:) = 0._wp
-            ind1 = 1+nbghostcells
+            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
             DO ji = mi0(ind1), mi1(ind1)                
                ztabramp(ji,:) = ssumask(ji,:)
             END DO
-            !
-            zmskwest(:) = 0._wp
-            zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
+            zmskwest(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1)
+            zmskwest(jpj+1:jpjmax) = 0._wp
          ENDIF
-
-         ! --- East --- !
-         IF( lk_east ) THEN
+         IF( lk_east ) THEN                             ! --- East --- !
             ztabramp(:,:) = 0._wp
-            ind1 = jpiglo - nbghostcells - 1
+            ind1 = jpiglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells
             DO ji = mi0(ind1), mi1(ind1)                 
                ztabramp(ji,:) = ssumask(ji,:)
             END DO
-            !
-            zmskeast(:) = 0._wp
-            zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
+            zmskeast(    1:jpj   ) = MAXVAL(ztabramp(:,:), dim=1)
+            zmskeast(jpj+1:jpjmax) = 0._wp
          ENDIF
-
-         ! --- South --- !
-         IF( lk_south ) THEN
+         IF( lk_south ) THEN                            ! --- South --- !
             ztabramp(:,:) = 0._wp
-            ind1 = 1+nbghostcells
+            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
             DO jj = mj0(ind1), mj1(ind1)                 
                ztabramp(:,jj) = ssvmask(:,jj)
             END DO
-            !
-            zmsksouth(:) = 0._wp
-            zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
+            zmsksouth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2)
+            zmsksouth(jpi+1:jpimax) = 0._wp
          ENDIF
-
-         ! --- North --- !
-         IF( lk_north) THEN
+         IF( lk_north ) THEN                            ! --- North --- !
             ztabramp(:,:) = 0._wp
-            ind1 = jpjglo - nbghostcells - 1
+            ind1 = jpjglo - ( nn_hls + nbghostcells + 1)   ! halo + land + nbghostcells
             DO jj = mj0(ind1), mj1(ind1)                 
                ztabramp(:,jj) = ssvmask(:,jj)
             END DO
-            !
-            zmsknorth(:) = 0._wp
-            zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
+            zmsknorth(    1:jpi   ) = MAXVAL(ztabramp(:,:), dim=2)
+            zmsknorth(jpi+1:jpimax) = 0._wp
          ENDIF
 
@@ -180,6 +173,6 @@
          zmskwest(:)  = 1._wp
          zmskeast(:)  = 1._wp
+         zmsksouth(:) = 1._wp
          zmsknorth(:) = 1._wp
-         zmsksouth(:) = 1._wp
 #if defined key_mpp_mpi
 !         CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax )
@@ -192,8 +185,8 @@
          ! Store it in ztabramp
 
-         ispongearea  = nn_sponge_len * Agrif_irhox()
-         z1_ispongearea = 1._wp / REAL( ispongearea )
-         jspongearea  = nn_sponge_len * Agrif_irhoy()
-         z1_jspongearea = 1._wp / REAL( jspongearea )
+         ispongearea    = nn_sponge_len * Agrif_irhox()
+         z1_ispongearea = 1._wp / REAL( ispongearea, wp )
+         jspongearea    = nn_sponge_len * Agrif_irhoy()
+         z1_jspongearea = 1._wp / REAL( jspongearea, wp )
          
          ztabramp(:,:) = 0._wp
@@ -203,17 +196,15 @@
          IF ( nbcellsy <= 3 ) jspongearea = -1
 
-         ! --- West --- !
-         IF(lk_west) THEN
-            ind1 = 1+nbghostcells
-            ind2 = 1+nbghostcells + ispongearea 
+         IF( lk_west ) THEN                             ! --- West --- !
+            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
+            ind2 = nn_hls + 1 + nbghostcells + ispongearea 
             DO ji = mi0(ind1), mi1(ind2)   
                DO jj = 1, jpj               
-                  ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)
-               END DO
-            END DO         
-
+                  ztabramp(ji,jj) =                       REAL(ind2 - mig(ji), wp) * z1_ispongearea   * zmskwest(jj)
+               END DO
+            END DO
             ! ghost cells:
             ind1 = 1
-            ind2 = nbghostcells + 1
+            ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
             DO ji = mi0(ind1), mi1(ind2)   
                DO jj = 1, jpj               
@@ -222,40 +213,32 @@
             END DO
          ENDIF
-
-         ! --- East --- !
-         IF(lk_east) THEN
-            ind1 = jpiglo - nbghostcells - ispongearea
-            ind2 = jpiglo - nbghostcells
+         IF( lk_east ) THEN                             ! --- East --- !
+            ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea
+            ind2 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
             DO ji = mi0(ind1), mi1(ind2)
-
                DO jj = 1, jpj
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)
-               ENDDO
-            END DO
-
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea ) * zmskeast(jj)
+               END DO
+            END DO
             ! ghost cells:
-            ind1 = jpiglo - nbghostcells
+            ind1 = jpiglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
             ind2 = jpiglo
             DO ji = mi0(ind1), mi1(ind2)
-
                DO jj = 1, jpj
                   ztabramp(ji,jj) = zmskeast(jj)
-               ENDDO
-            END DO
-         ENDIF
-
-         ! --- South --- !
-         IF( lk_south ) THEN 
-            ind1 = 1+nbghostcells
-            ind2 = 1+nbghostcells + jspongearea
+               END DO
+            END DO
+         ENDIF      
+         IF( lk_south ) THEN                            ! --- South --- !
+            ind1 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
+            ind2 = nn_hls + 1 + nbghostcells + jspongearea 
             DO jj = mj0(ind1), mj1(ind2) 
                DO ji = 1, jpi
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)
-               END DO
-            END DO
-
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea ) * zmsksouth(ji)
+               END DO
+            END DO
             ! ghost cells:
             ind1 = 1
-            ind2 = nbghostcells + 1
+            ind2 = nn_hls + 1 + nbghostcells               ! halo + land + nbghostcells
             DO jj = mj0(ind1), mj1(ind2) 
                DO ji = 1, jpi
@@ -264,17 +247,14 @@
             END DO
          ENDIF
-
-         ! --- North --- !
-         IF( lk_north ) THEN  
-            ind1 = jpjglo - nbghostcells - jspongearea
-            ind2 = jpjglo - nbghostcells
+         IF( lk_north ) THEN                            ! --- North --- !
+            ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea
+            ind2 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
             DO jj = mj0(ind1), mj1(ind2)
                DO ji = 1, jpi
-                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)
-               END DO
-            END DO
-
+                  ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea ) * zmsknorth(ji)
+               END DO
+            END DO
             ! ghost cells:
-            ind1 = jpjglo - nbghostcells
+            ind1 = jpjglo - ( nn_hls + nbghostcells )      ! halo + land + nbghostcells - 1
             ind2 = jpjglo
             DO jj = mj0(ind1), mj1(ind2)
@@ -284,5 +264,5 @@
             END DO
          ENDIF
-      
+         !
       ENDIF
 
@@ -291,12 +271,8 @@
          fspu(:,:) = 0._wp
          fspv(:,:) = 0._wp
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj  ) ) * ssumask(ji,jj)
             fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) ) * ssvmask(ji,jj)
          END_2D
-         CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1.0_wp )   ! Lateral boundary conditions
-         CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1.0_wp )
-
-         spongedoneT = .TRUE.
       ENDIF
 
@@ -305,5 +281,5 @@
          fspt(:,:) = 0._wp
          fspf(:,:) = 0._wp
-         DO_2D_00_00
+         DO_2D( 0, 0, 0, 0 )
             fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj)
             fspf(ji,jj) = 0.25_wp * ( ztabramp(ji  ,jj  ) + ztabramp(ji  ,jj+1)   &
@@ -311,13 +287,23 @@
                                   &  * ssvmask(ji,jj) * ssvmask(ji,jj+1)
          END_2D
-         CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1.0_wp )   ! Lateral boundary conditions
-         CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1.0_wp )
-         
+      ENDIF
+      
+      IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
+         CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp )
+         spongedoneT = .TRUE.
          spongedoneU = .TRUE.
       ENDIF
+      IF( .NOT. spongedoneT ) THEN
+         CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp )
+         spongedoneT = .TRUE.
+      ENDIF
+      IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
+         CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1._wp, fspf, 'F', 1._wp )
+         spongedoneU = .TRUE.
+      ENDIF
 
 #if defined key_vertical
       ! Remove vertical interpolation where not needed:
-      DO_2D_00_00
+      DO_2D( 0, 0, 0, 0 )
          IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. &
          &   (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0
@@ -334,10 +320,11 @@
       END_2D
       !
-      ztabramp(:,:) = REAL( mbkt_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1.0_wp )
-      mbkt_parent(:,:) = NINT( ztabramp(:,:) )
-      ztabramp(:,:) = REAL( mbku_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1.0_wp )
-      mbku_parent(:,:) = NINT( ztabramp(:,:) )
-      ztabramp(:,:) = REAL( mbkv_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1.0_wp )
-      mbkv_parent(:,:) = NINT( ztabramp(:,:) )
+      ztabramp (:,:) = REAL( mbkt_parent(:,:), wp )
+      ztabrampu(:,:) = REAL( mbku_parent(:,:), wp )
+      ztabrampv(:,:) = REAL( mbkv_parent(:,:), wp )
+      CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1._wp, ztabrampu, 'U', 1._wp, ztabrampv, 'V', 1._wp )
+      mbkt_parent(:,:) = NINT( ztabramp (:,:) )
+      mbku_parent(:,:) = NINT( ztabrampu(:,:) )
+      mbkv_parent(:,:) = NINT( ztabrampv(:,:) )
 #endif
       !
@@ -346,4 +333,5 @@
    END SUBROUTINE Agrif_Sponge
 
+   
    SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
       !!----------------------------------------------------------------------
@@ -433,5 +421,5 @@
                   N_out = N_out + 1
                   h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
-               ENDDO
+               END DO
 
                ! Account for small differences in free-surface
@@ -444,6 +432,6 @@
                   CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
                ENDIF
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 # endif
 
@@ -456,7 +444,7 @@
                   tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
 # endif
-               ENDDO
-            ENDDO
-         ENDDO
+               END DO
+            END DO
+         END DO
 
          DO jn = 1, jpts            
@@ -513,4 +501,5 @@
    END SUBROUTINE interptsn_sponge
 
+   
    SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
       !!---------------------------------------------
@@ -521,6 +510,6 @@
       LOGICAL, INTENT(in) :: before
 
-      INTEGER :: ji,jj,jk,jmax
-
+      INTEGER  :: ji,jj,jk,jmax
+      INTEGER  :: ind1
       ! sponge parameters 
       REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
@@ -586,5 +575,5 @@
                   zhtot = zhtot + h_in(jk)
                   tabin(jk) = tabres(ji,jj,jk,m1)
-               ENDDO
+               END DO
                !         
                N_out = 0
@@ -593,5 +582,5 @@
                   N_out = N_out + 1
                   h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
-               ENDDO
+               END DO
 
                ! Account for small differences in free-surface
@@ -605,6 +594,6 @@
                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
                ENDIF 
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 
          ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
@@ -659,6 +648,8 @@
 
          jmax = j2-1
-        ! IF (lk_north) jmax = MIN(jmax,nlcj-nbghostcells-2)   ! North
-         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2)   ! North
+         ind1 = jpjglo - ( nn_hls + nbghostcells + 2 )   ! North
+         DO jj = mj0(ind1), mj1(ind1)                 
+            jmax = MIN(jmax,jj)
+         END DO
 
          DO jj = j1+1, jmax
@@ -688,5 +679,6 @@
    END SUBROUTINE interpun_sponge
 
-   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
+   
+   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
       !!---------------------------------------------
       !!   *** ROUTINE interpvn_sponge ***
@@ -695,7 +687,7 @@
       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
       LOGICAL, INTENT(in) :: before
-      INTEGER, INTENT(in) :: nb , ndir
       !
       INTEGER  ::   ji, jj, jk, imax
+      INTEGER  :: ind1
       REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, zhtot
       REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
@@ -759,5 +751,5 @@
                   zhtot = zhtot + h_in(jk)
                   tabin(jk) = tabres(ji,jj,jk,m1)
-               ENDDO
+               END DO
                !          
                N_out = 0
@@ -766,5 +758,5 @@
                   N_out = N_out + 1
                   h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
-               ENDDO
+               END DO
 
                ! Account for small differences in free-surface
@@ -778,6 +770,6 @@
                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
                ENDIF
-            ENDDO
-         ENDDO
+            END DO
+         END DO
 
          vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)  
@@ -812,7 +804,9 @@
 
          imax = i2 - 1
-      !   IF(lk_east) imax = MIN(imax,nlci-nbghostcells-2)   ! East
-         IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-nbghostcells-2)   ! East
-
+         ind1 = jpiglo - ( nn_hls + nbghostcells + 2 )   ! East
+         DO ji = mi0(ind1), mi1(ind1)                
+            imax = MIN(imax,ji)
+         END DO
+         
          DO jj = j1+1, j2
             DO ji = i1+1, imax   ! vector opt.
