Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_eap.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_eap.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_eap.F90	(revision 15020)
@@ -58,14 +58,12 @@
 
    REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) ::   s11r, s12r, s22r, s11s, s12s, s22s
-
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
-#  include "domzgr_substitute.h90"
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice
 
    !! for convergence tests
    INTEGER ::   ncvgid   ! netcdf file id
    INTEGER ::   nvarid   ! netcdf variable id
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   aimsk00
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   eap_res  , aimsk15
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -180,7 +178,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
       !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
 
       REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
@@ -203,55 +201,36 @@
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology'
       !
-      IF( kt == nit000 )  THEN 
-         !
-         ! for diagnostics 
-         ALLOCATE( aimsk00(jpi,jpj) )
-         ! for convergence tests
-         IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) )
-      ENDIF
-      !
+      ! for diagnostics and convergence tests
       DO_2D( 1, 1, 1, 1 )
-         aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
+         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
       END_2D
       IF( nn_rhg_chkcvg > 0 ) THEN
          DO_2D( 1, 1, 1, 1 )
-            aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
          END_2D
       ENDIF
       !
-!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
       !------------------------------------------------------------------------------!
       ! 0) mask at F points for the ice
       !------------------------------------------------------------------------------!
-      ! ocean/land mask
-      DO_2D( 1, 0, 1, 0 )
-         zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
-      END_2D
-      CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1._wp )
-
-      ! Lateral boundary conditions on velocity (modify zfmask)
-      DO_2D( 0, 0, 0, 0 )
-         IF( zfmask(ji,jj) == 0._wp ) THEN
-            zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
-               &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
+      IF( kt == nit000 ) THEN
+         ! ocean/land mask
+         ALLOCATE( fimask(jpi,jpj) )
+         IF( rn_ishlat == 0._wp ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+            END_2D
+         ELSE
+            DO_2D( 0, 0, 0, 0 )
+               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+               ! Lateral boundary conditions on velocity (modify fimask)
+               IF( fimask(ji,jj) == 0._wp ) THEN
+                  fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
+                     &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
+               ENDIF
+            END_2D
          ENDIF
-      END_2D
-      DO jj = 2, jpjm1
-         IF( zfmask(1,jj) == 0._wp ) THEN
-            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
-         ENDIF
-         IF( zfmask(jpi,jj) == 0._wp ) THEN
-            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
-         ENDIF
-      END DO
-      DO ji = 2, jpim1
-         IF( zfmask(ji,1) == 0._wp ) THEN
-            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
-         ENDIF
-         IF( zfmask(ji,jpj) == 0._wp ) THEN
-            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
-         ENDIF
-      END DO
-      CALL lbc_lnk( 'icedyn_rhg_eap', zfmask, 'F', 1.0_wp )
+         CALL lbc_lnk( 'icedyn_rhg_eap', fimask, 'F', 1.0_wp )
+      ENDIF
 
       !------------------------------------------------------------------------------!
@@ -401,5 +380,5 @@
             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
 
          END_2D
@@ -760,5 +739,5 @@
 
          ! convergence test
-         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
          !
          !                                                ! ==================== !
@@ -777,5 +756,5 @@
          zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
             &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-            &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+            &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
 
       END_2D
@@ -830,17 +809,17 @@
             &                            ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
          !
-         CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 )
-         CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 )
-         CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 )
-         CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 )
-         CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 )
-         CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 )
+         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
       ENDIF
 
       ! --- divergence, shear and strength --- !
-      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * aimsk00 )   ! divergence
-      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * aimsk00 )   ! shear
-      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * aimsk00 )   ! delta
-      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * aimsk00 )   ! strength
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
 
       ! --- Stress tensor invariants (SIMIP diags) --- !
@@ -867,6 +846,6 @@
          !
          ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
-         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress
-         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
 
          DEALLOCATE ( zsig_I, zsig_II )
@@ -914,7 +893,7 @@
          CALL lbc_lnk( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp )
 
-         CALL iom_put( 'yield11', zyield11 * aimsk00 )
-         CALL iom_put( 'yield22', zyield22 * aimsk00 )
-         CALL iom_put( 'yield12', zyield12 * aimsk00 )
+         CALL iom_put( 'yield11', zyield11 * zmsk00 )
+         CALL iom_put( 'yield22', zyield22 * zmsk00 )
+         CALL iom_put( 'yield12', zyield12 * zmsk00 )
       ENDIF
 
@@ -922,5 +901,5 @@
       IF( iom_use('aniso') ) THEN
          CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp )
-         CALL iom_put( 'aniso' , paniso_11 * aimsk00 )
+         CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
       ENDIF
 
@@ -933,10 +912,10 @@
             &                              zfU, 'U', -1.0_wp,   zfV, 'V', -1.0_wp )
 
-         CALL iom_put( 'dssh_dx' , zspgU * aimsk00 )   ! Sea-surface tilt term in force balance (x)
-         CALL iom_put( 'dssh_dy' , zspgV * aimsk00 )   ! Sea-surface tilt term in force balance (y)
-         CALL iom_put( 'corstrx' , zCorU * aimsk00 )   ! Coriolis force term in force balance (x)
-         CALL iom_put( 'corstry' , zCorV * aimsk00 )   ! Coriolis force term in force balance (y)
-         CALL iom_put( 'intstrx' , zfU   * aimsk00 )   ! Internal force term in force balance (x)
-         CALL iom_put( 'intstry' , zfV   * aimsk00 )   ! Internal force term in force balance (y)
+         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
       ENDIF
 
@@ -949,6 +928,6 @@
          DO_2D( 0, 0, 0, 0 )
             ! 2D ice mass, snow mass, area transport arrays (X, Y)
-            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj)
-            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj)
+            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
+            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
 
             zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
@@ -984,8 +963,8 @@
             IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
                CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
-                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) )
+                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
             ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
                CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
-                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) )
+                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
             ENDIF
          ENDIF
@@ -995,5 +974,5 @@
 
 
-   SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb )
+   SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE rhg_cvg_eap  ***
@@ -1010,4 +989,5 @@
       INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmsk15
       !!
       INTEGER           ::   it, idtime, istatus
@@ -1044,10 +1024,9 @@
          zresm = 0._wp
       ELSE
-         DO_2D( 1, 1, 1, 1 )
-            eap_res(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
-               &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * aimsk15(ji,jj)
+         zresm = 0._wp
+         DO_2D( 0, 0, 0, 0 )
+            zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+               &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
          END_2D
-
-         zresm = MAXVAL( eap_res )
          CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
       ENDIF
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_evp.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_evp.F90	(revision 15020)
@@ -48,12 +48,11 @@
    PUBLIC   rhg_evp_rst       ! called by icedyn_rhg.F90
 
-   !! * Substitutions
-#  include "do_loop_substitute.h90"
-#  include "domzgr_substitute.h90"
-
    !! for convergence tests
    INTEGER ::   ncvgid   ! netcdf file id
    INTEGER ::   nvarid   ! netcdf variable id
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice
+
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -161,7 +160,7 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
       !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
 
       REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
@@ -185,45 +184,35 @@
       !
       ! for diagnostics and convergence tests
-      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
       DO_2D( 1, 1, 1, 1 )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
-         zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
       END_2D
-      !
-      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
+      IF( nn_rhg_chkcvg > 0 ) THEN
+         DO_2D( 1, 1, 1, 1 )
+            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+         END_2D
+      ENDIF
+      !
       !------------------------------------------------------------------------------!
       ! 0) mask at F points for the ice
       !------------------------------------------------------------------------------!
-      ! ocean/land mask
-      DO_2D( 1, 0, 1, 0 )
-         zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
-      END_2D
-      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp)
-
-      ! Lateral boundary conditions on velocity (modify zfmask)
-      DO_2D( 0, 0, 0, 0 )
-         IF( zfmask(ji,jj) == 0._wp ) THEN
-            zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
-               &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
+      IF( kt == nit000 ) THEN
+         ! ocean/land mask
+         ALLOCATE( fimask(jpi,jpj) )
+         IF( rn_ishlat == 0._wp ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+            END_2D
+         ELSE
+            DO_2D( 0, 0, 0, 0 )
+               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+               ! Lateral boundary conditions on velocity (modify fimask)
+               IF( fimask(ji,jj) == 0._wp ) THEN
+                  fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
+                     &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
+               ENDIF
+            END_2D
          ENDIF
-      END_2D
-      DO jj = 2, jpjm1
-         IF( zfmask(1,jj) == 0._wp ) THEN
-            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
-         ENDIF
-         IF( zfmask(jpi,jj) == 0._wp ) THEN
-            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
-        ENDIF
-      END DO
-      DO ji = 2, jpim1
-         IF( zfmask(ji,1) == 0._wp ) THEN
-            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
-         ENDIF
-         IF( zfmask(ji,jpj) == 0._wp ) THEN
-            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
-         ENDIF
-      END DO
-      CALL lbc_lnk( 'icedyn_rhg_evp', zfmask, 'F', 1._wp )
-
+         CALL lbc_lnk( 'icedyn_rhg_evp', fimask, 'F', 1._wp )
+      ENDIF
       !------------------------------------------------------------------------------!
       ! 1) define some variables and initialize arrays
@@ -367,5 +356,5 @@
             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
 
          END_2D
@@ -702,5 +691,5 @@
 
          ! convergence test
-         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice, zmsk15 )
          !
          !                                                ! ==================== !
@@ -717,5 +706,5 @@
          zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
             &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-            &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+            &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
 
       END_2D
@@ -914,10 +903,8 @@
       ENDIF
       !
-      DEALLOCATE( zmsk00, zmsk15 )
-      !
    END SUBROUTINE ice_dyn_rhg_evp
 
 
-   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
+   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb, pmsk15 )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE rhg_cvg  ***
@@ -934,4 +921,5 @@
       INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
       REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pmsk15
       !!
       INTEGER           ::   it, idtime, istatus
@@ -939,5 +927,4 @@
       REAL(wp)          ::   zresm           ! local real
       CHARACTER(len=20) ::   clname
-      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
       !!----------------------------------------------------------------------
 
@@ -969,9 +956,9 @@
          zresm = 0._wp
       ELSE
-         DO_2D( 1, 1, 1, 1 )
-            zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
-               &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
-         END_2D
-         zresm = MAXVAL( zres )
+         zresm = 0._wp
+         DO_2D( 0, 0, 0, 0 )
+            zresm = MAX( zresm, MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+               &                     ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * pmsk15(ji,jj) )
+         END_2D
          CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
       ENDIF
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_vp.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_vp.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/ICE/icedyn_rhg_vp.F90	(revision 15020)
@@ -59,6 +59,8 @@
    INTEGER ::   nvarid_ures_xy, nvarid_vres_xy
 
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
-
+   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   fimask   ! mask at F points for the ice
+   
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -158,5 +160,4 @@
       REAL(wp) ::   zAA3, zw, ztau, zuerr_max, zverr_max
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice
       REAL(wp), DIMENSION(jpi,jpj) ::   za_iU  , za_iV                      ! ice fraction on U/V points
       REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! Acceleration term contribution to RHS
@@ -197,4 +198,5 @@
 !!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00, zmsk15
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! mask for lots of ice (1), little ice (0)
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence (1), no ice (0)
@@ -238,11 +240,12 @@
       
       ! for diagnostics and convergence tests
-      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
+      DO_2D( 1, 1, 1, 1 )
+         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
+      END_2D
+      IF( nn_rhg_chkcvg > 0 ) THEN
+         DO_2D( 1, 1, 1, 1 )
             zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
-         END DO
-      END DO
+         END_2D
+      ENDIF
       
       IF ( lp_zebra_vp ) THEN; nn_zebra_vp = 2
@@ -292,43 +295,28 @@
       ! -- F-mask       (code from EVP)
       !------------------------------
-      ! MartinV: 
-      ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip)
-      ! I am not sure the same recipe applies here
-      
-      ! - ocean/land mask
-      DO jj = 1, jpj - 1
-         DO ji = 1, jpi - 1
-            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
-         END DO
-      END DO
-
-      ! Lateral boundary conditions on velocity (modify zfmask)
-      ! Can be computed once for all, at first time step, for all rheologies
-      DO jj = 2, jpj - 1
-         DO ji = 2, jpi - 1   ! vector opt.
-            IF( zfmask(ji,jj) == 0._wp ) THEN
-               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
-                  &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
-            ENDIF
-         END DO
-      END DO
-      DO jj = 2, jpj - 1
-         IF( zfmask(1,jj) == 0._wp ) THEN
-            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
+      IF( kt == nit000 ) THEN
+         ! MartinV: 
+         ! In EVP routine, fimask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip)
+         ! I am not sure the same recipe applies here
+         
+         ! - ocean/land mask
+         ALLOCATE( fimask(jpi,jpj) )
+         IF( rn_ishlat == 0._wp ) THEN
+            DO_2D( 0, 0, 0, 0 )
+               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+            END_2D
+         ELSE
+            DO_2D( 0, 0, 0, 0 )
+               fimask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+               ! Lateral boundary conditions on velocity (modify fimask)
+               IF( fimask(ji,jj) == 0._wp ) THEN
+                  fimask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), &
+                     &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) )
+               ENDIF
+            END_2D
          ENDIF
-         IF( zfmask(jpi,jj) == 0._wp ) THEN
-            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpi - 1,jj,1), umask(jpi,jj-1,1) ) )
-        ENDIF
-      END DO
-      DO ji = 2, jpi - 1
-         IF( zfmask(ji,1) == 0._wp ) THEN
-            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
-         ENDIF
-         IF( zfmask(ji,jpj) == 0._wp ) THEN
-            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpj - 1,1) ) )
-         ENDIF
-      END DO
-      
-      CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp )
+         
+         CALL lbc_lnk( 'icedyn_rhg_vp', fimask, 'F', 1._wp )
+      ENDIF
       
       !----------------------------------------------------------------------------------------------------------
@@ -455,5 +443,5 @@
                zds(ji,jj) = ( ( zu_c(ji,jj+1) * r1_e1u(ji,jj+1) - zu_c(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                   &         + ( zv_c(ji+1,jj) * r1_e2v(ji+1,jj) - zv_c(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+                  &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
 
             END DO
@@ -521,5 +509,5 @@
                
                ! Temporary zef factor at F-point
-               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj)
+               zef(ji,jj)      = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * fimask(ji,jj)
 
             END DO
@@ -611,5 +599,5 @@
          
          ! --- Stress contributions at f-points         
-         ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing
+         ! MV NOTE: I applied fimask on zds, by mimetism on EVP, but without deep understanding of what I was doing
          ! My guess is that this is the way to enforce boundary conditions on strain rate tensor
 
@@ -620,8 +608,8 @@
                
                ! sig12 contribution to RHS of U equation at F-points 
-               zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)
+               zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * fimask(ji,jj)
                
                ! sig12 contribution to RHS of V equation at F-points
-               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)
+               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * fimask(ji,jj)
 
             END DO
@@ -1181,5 +1169,5 @@
             zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+               &         ) * r1_e1e2f(ji,jj) * fimask(ji,jj)
 
          END DO
@@ -1501,6 +1489,4 @@
       ENDIF
 
-      DEALLOCATE( zmsk00, zmsk15 )
-
    END SUBROUTINE ice_dyn_rhg_vp
    
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/BDY/bdyini.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/BDY/bdyini.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/BDY/bdyini.F90	(revision 15020)
@@ -44,4 +44,7 @@
    INTEGER, DIMENSION(jp_nseg) ::   jpjnob, jpindt, jpinft, npckgn   !
    INTEGER, DIMENSION(jp_nseg) ::   jpjsob, jpisdt, jpisft, npckgs   !
+   
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -140,4 +143,5 @@
       !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
       !!----------------------------------------------------------------------      
+      INTEGER  ::   ji, jj                                 ! dummy loop indices
       INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices
       INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers
@@ -630,10 +634,8 @@
       ! For the flagu/flagv calculation below we require a version of fmask without
       ! the land boundary condition (shlat) included:
-      DO ij = 1, jpjm1
-         DO ii = 1, jpim1
-            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
-               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
-         END DO
-      END DO
+      DO_2D( 0, 0, 0, 0 )
+         zfmask(ji,jj) =  ztmask(ji,jj  ) * ztmask(ji+1,jj  )   &
+            &           * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
+      END_2D
       CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
 
@@ -646,10 +648,8 @@
 
       ! Derive mask on U and V grid from mask on T grid
-      DO ij = 1, jpjm1
-         DO ii = 1, jpim1
-            bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij  )
-            bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1)  
-         END DO
-      END DO
+      DO_2D( 0, 0, 0, 0 )
+            bdyumask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji+1,jj  )
+            bdyvmask(ji,jj) = bdytmask(ji,jj) * bdytmask(ji  ,jj+1)  
+      END_2D
       CALL lbc_lnk( 'bdyini', bdyumask, 'U', 1.0_wp , bdyvmask, 'V', 1.0_wp )   ! Lateral boundary cond. 
 
@@ -687,10 +687,8 @@
 
       ! Recompute zfmask
-      DO ij = 1, jpjm1
-         DO ii = 1, jpim1
-            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
-               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
-         END DO
-      END DO
+      DO_2D( 0, 0, 0, 0 )
+         zfmask(ji,jj) =  ztmask(ji,jj  ) * ztmask(ji+1,jj  )   &
+            &           * ztmask(ji,jj+1) * ztmask(ji+1,jj+1)
+      END_2D
       CALL lbc_lnk( 'bdyini', zfmask, 'F', 1.0_wp )
 
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domain.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domain.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domain.F90	(revision 15020)
@@ -148,7 +148,7 @@
       END DO
       !
-      DO jk = 1, jpkm1
-         hf_0(1:jpim1,:) = hf_0(1:jpim1,:) + e3f_0(1:jpim1,:,jk)*vmask(1:jpim1,:,jk)*vmask(2:jpi,:,jk)
-      END DO
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         hf_0(ji,jj) = hf_0(ji,jj) + e3f_0(ji,jj,jk)*vmask(ji,jj,jk)*vmask(ji+1,jj,jk)
+      END_3D
       CALL lbc_lnk('domain', hf_0, 'F', 1._wp)
       !
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/dommsk.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/dommsk.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/dommsk.F90	(revision 15020)
@@ -182,5 +182,5 @@
       ssfmask(:,:) = MAXVAL( fmask(:,:,:), DIM=3 )
       IF( lk_SWE ) THEN      ! Shallow Water Eq. case : redefine ssfmask
-         DO_2D( 0,0, 0,0  )
+         DO_2D( 0, 0, 0, 0 )
             ssfmask(ji,jj) = MAX(  ssmask(ji,jj+1), ssmask(ji+1,jj+1),  & 
                &                   ssmask(ji,jj  ), ssmask(ji+1,jj  )   )
@@ -202,31 +202,12 @@
       ! Lateral boundary conditions on velocity (modify fmask)
       ! ---------------------------------------  
-      IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition
+      IF( rn_shlat /= 0._wp ) THEN      ! Not free-slip lateral boundary condition
          !
-         DO jk = 1, jpk
-            DO_2D( 0, 0, 0, 0 )
-               IF( fmask(ji,jj,jk) == 0._wp ) THEN
-                  fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk), &
-                     &                                           vmask(ji,jj,jk), vmask(ji+1,jj,jk) ) )
-               ENDIF
-            END_2D
-            DO jj = 2, jpjm1
-               IF( fmask(1,jj,jk) == 0._wp ) THEN
-                  fmask(1  ,jj,jk) = rn_shlat * MIN( 1._wp , MAX( vmask(2,jj,jk), umask(1,jj+1,jk), umask(1,jj,jk) ) )
-               ENDIF
-               IF( fmask(jpi,jj,jk) == 0._wp ) THEN
-                  fmask(jpi,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(jpi,jj+1,jk), vmask(jpim1,jj,jk), umask(jpi,jj-1,jk) ) )
-               ENDIF
-            END DO         
-            DO ji = 2, jpim1
-               IF( fmask(ji,1,jk) == 0._wp ) THEN
-                  fmask(ji, 1 ,jk) = rn_shlat * MIN( 1._wp , MAX( vmask(ji+1,1,jk), umask(ji,2,jk), vmask(ji,1,jk) ) )
-               ENDIF
-               IF( fmask(ji,jpj,jk) == 0._wp ) THEN
-                  fmask(ji,jpj,jk) = rn_shlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,jk), vmask(ji-1,jpj,jk), umask(ji,jpjm1,jk) ) )
-               ENDIF
-            END DO
-         END DO
-         !
+         DO_3D( 0, 0, 0, 0, 1, jpk )
+            IF( fmask(ji,jj,jk) == 0._wp ) THEN
+               fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk), &
+                  &                                           vmask(ji,jj,jk), vmask(ji+1,jj,jk) ) )
+            ENDIF
+         END_3D
          CALL lbc_lnk( 'dommsk', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask
          !
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domqco.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domqco.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domqco.F90	(revision 15020)
@@ -184,5 +184,5 @@
          !                                ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
 
-      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )                               ! start from 1 since lbc_lnk('F') doesn't update the 1st row/line
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ! round brackets added to fix the order of floating point operations
          ! needed to ensure halo 1 - halo 2 compatibility
@@ -197,9 +197,9 @@
 !!st         ELSE                                      !- Flux Form   (simple averaging)
 #else
-      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ! round brackets added to fix the order of floating point operations
          ! needed to ensure halo 1 - halo 2 compatibility
          pr3f(ji,jj) = 0.25_wp * ( ( pssh(ji,jj  ) + pssh(ji+1,jj  ) ) &
-            &                     + ( pssh(ji,jj+1) + pssh(ji+1,jj+1)  &
+            &                    + ( pssh(ji,jj+1) + pssh(ji+1,jj+1)   &
             &                       )                                  & ! bracket for halo 1 - halo 2 compatibility
             &                    ) * r1_hf_0(ji,jj)
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domvvl.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domvvl.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domvvl.F90	(revision 15020)
@@ -713,5 +713,5 @@
          !
       CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
-         DO_3D( 1, 0, 1, 0, 1, jpk )
+         DO_3D( 0, 0, 0, 0, 1, jpk )
             pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
                &                       *    r1_e1e2f(ji,jj)                                                  &
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domzgr.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domzgr.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/DOM/domzgr.F90	(revision 15020)
@@ -340,5 +340,5 @@
       !                                    ! N.B.  top     k-index of W-level = mikt
       !                                    !       bottom  k-index of W-level = mbkt+1
-      DO_2D( 1, 0, 1, 0 )
+      DO_2D( 0, 0, 0, 0 )
          miku(ji,jj) = MAX(  mikt(ji+1,jj  ) , mikt(ji,jj)  )
          mikv(ji,jj) = MAX(  mikt(ji  ,jj+1) , mikt(ji,jj)  )
@@ -349,10 +349,33 @@
       END_2D
       ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
-      zk(:,:) = REAL( miku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1.0_wp )   ;   miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
-      zk(:,:) = REAL( mikv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1.0_wp )   ;   mikv(:,:) = MAX( NINT( zk(:,:) ), 1 )
-      zk(:,:) = REAL( mikf(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'F', 1.0_wp )   ;   mikf(:,:) = MAX( NINT( zk(:,:) ), 1 )
-      !
-      zk(:,:) = REAL( mbku(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'U', 1.0_wp )   ;   mbku(:,:) = MAX( NINT( zk(:,:) ), 1 )
-      zk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk( 'domzgr', zk, 'V', 1.0_wp )   ;   mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 )
+      DO_2D( 0, 0, 0, 0 )
+         zk(ji,jj) = REAL( miku(ji,jj), wp )
+      END_2D
+      CALL lbc_lnk( 'domzgr', zk, 'U', 1.0_wp )
+      miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
+
+      DO_2D( 0, 0, 0, 0 )
+         zk(ji,jj) = REAL( mikv(ji,jj), wp )
+      END_2D
+      CALL lbc_lnk( 'domzgr', zk, 'V', 1.0_wp )
+      mikv(:,:) = MAX( NINT( zk(:,:) ), 1 )
+      
+      DO_2D( 0, 0, 0, 0 )
+         zk(ji,jj) = REAL( mikf(ji,jj), wp )
+      END_2D
+      CALL lbc_lnk( 'domzgr', zk, 'F', 1.0_wp )
+      mikf(:,:) = MAX( NINT( zk(:,:) ), 1 )
+      !
+      DO_2D( 0, 0, 0, 0 )
+         zk(ji,jj) = REAL( mbku(ji,jj), wp )
+      END_2D
+      CALL lbc_lnk( 'domzgr', zk, 'U', 1.0_wp )
+      mbku(:,:) = MAX( NINT( zk(:,:) ), 1 )
+      
+      DO_2D( 0, 0, 0, 0 )
+         zk(ji,jj) = REAL( mbkv(ji,jj), wp )
+      END_2D
+      CALL lbc_lnk( 'domzgr', zk, 'V', 1.0_wp )
+      mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 )
       !
    END SUBROUTINE zgr_top_bot
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/LDF/ldfc1d_c2d.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/LDF/ldfc1d_c2d.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/LDF/ldfc1d_c2d.F90	(revision 15020)
@@ -80,5 +80,5 @@
             pah1(:,:,jk) = pahs1(:,:) * (  zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) )  )
          END DO
-         DO_3DS( 1, 0, 1, 0, jpkm1, 1, -1 )  ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
+         DO_3DS( 0, 0, 0, 0, jpkm1, 1, -1 )  ! pah2 at F-point (zdep2 is an approximation in zps-coord.)
             zdep2 = (  gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk)   &
                &     + gdept_0(ji,jj  ,jk) + gdept_0(ji+1,jj  ,jk)  ) * r1_4
@@ -88,5 +88,5 @@
          !
       CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
-         DO_3DS( 1, 0, 1, 0, jpkm1, 1, -1 )
+         DO_3DS( 0, 0, 0, 0, jpkm1, 1, -1 )
             zdep1 = (  gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk)  ) * 0.5_wp
             zdep2 = (  gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk)  ) * 0.5_wp
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/LDF/ldfdyn.F90
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/LDF/ldfdyn.F90	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/src/OCE/LDF/ldfdyn.F90	(revision 15020)
@@ -385,10 +385,7 @@
                   zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
                   zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
+                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk,Kbb) * uu(ji  ,jj+1,jk,Kbb) + vv(ji+1,jj  ,jk,Kbb) * vv(ji+1,jj  ,jk,Kbb)
                   zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
                   ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk)      ! 288= 12*12 * 2
-               END_2D
-               DO_2D( 1, 0, 1, 0 )
-                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, Kbb) * uu(ji  ,jj+1,jk, Kbb) + vv(ji+1,jj  ,jk, Kbb) * vv(ji+1,jj  ,jk, Kbb)
-                  zu2pv2_ij    = uu(ji  ,jj  ,jk, Kbb) * uu(ji  ,jj  ,jk, Kbb) + vv(ji  ,jj  ,jk, Kbb) * vv(ji  ,jj  ,jk, Kbb)
                   zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
                   ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk)      ! 288= 12*12 * 2
@@ -400,10 +397,7 @@
                   zu2pv2_ij    = uu(ji  ,jj  ,jk,Kbb) * uu(ji  ,jj  ,jk,Kbb) + vv(ji  ,jj  ,jk,Kbb) * vv(ji  ,jj  ,jk,Kbb)
                   zu2pv2_ij_m1 = uu(ji-1,jj  ,jk,Kbb) * uu(ji-1,jj  ,jk,Kbb) + vv(ji  ,jj-1,jk,Kbb) * vv(ji  ,jj-1,jk,Kbb)
+                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk,Kbb) * uu(ji  ,jj+1,jk,Kbb) + vv(ji+1,jj  ,jk,Kbb) * vv(ji+1,jj  ,jk,Kbb)
                   zemax = MAX( e1t(ji,jj) , e2t(ji,jj) )
                   ahmt(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax  ) * zemax * tmask(ji,jj,jk)
-               END_2D
-               DO_2D( 1, 0, 1, 0 )
-                  zu2pv2_ij_p1 = uu(ji  ,jj+1,jk, Kbb) * uu(ji  ,jj+1,jk, Kbb) + vv(ji+1,jj  ,jk, Kbb) * vv(ji+1,jj  ,jk, Kbb)
-                  zu2pv2_ij    = uu(ji  ,jj  ,jk, Kbb) * uu(ji  ,jj  ,jk, Kbb) + vv(ji  ,jj  ,jk, Kbb) * vv(ji  ,jj  ,jk, Kbb)
                   zemax = MAX( e1f(ji,jj) , e2f(ji,jj) )
                   ahmf(ji,jj,jk) = SQRT(  SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax  ) * zemax * fmask(ji,jj,jk)
@@ -487,6 +481,4 @@
                DO_2D( 0, 0, 0, 0 )
                   ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) )
-               END_2D
-               DO_2D( 1, 0, 1, 0 )
                   ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) )
                END_2D
Index: NEMO/branches/2021/ticket2680_C1D_PAPA/tests/C1D_ASICS/EXP00/namelist_cfg
===================================================================
--- NEMO/branches/2021/ticket2680_C1D_PAPA/tests/C1D_ASICS/EXP00/namelist_cfg	(revision 15018)
+++ NEMO/branches/2021/ticket2680_C1D_PAPA/tests/C1D_ASICS/EXP00/namelist_cfg	(revision 15020)
@@ -49,7 +49,8 @@
 &namdom        !   time and space domain
 !-----------------------------------------------------------------------
-   ln_linssh   = .true.   !  =T  linear free surface  ==>>  model level are fixed in time
+   ln_linssh   = .true.    !  =T  linear free surface  ==>>  model level are fixed in time
    !
-   rn_Dt      =  360.     !  time step for the dynamics and tracer
+   rn_Dt       =  360.     !  time step for the dynamics and tracer
+   ln_c1d      = .true.    !  Single column domain (1x1pt)              (T => fill namc1d)
 /
 !-----------------------------------------------------------------------
