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 14914)
+++ 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
