diff --git a/src/ICE/icectl.F90 b/src/ICE/icectl.F90
index 568b86fbcdd31f58e0eb4ee91a80d8019532cbcd..e66fe77985a8370c07564c29a8c63fb9b154e4c8 100644
--- a/src/ICE/icectl.F90
+++ b/src/ICE/icectl.F90
@@ -47,9 +47,9 @@ MODULE icectl
 
    ! thresold rates for conservation
    !    these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk
-   REAL(wp), PARAMETER ::   zchk_m   = 2.5e-7   ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost
-   REAL(wp), PARAMETER ::   zchk_s   = 2.5e-6   ! g/m2/s  <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
-   REAL(wp), PARAMETER ::   zchk_t   = 7.5e-2   ! W/m2    <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
+   REAL(wp), PARAMETER ::   rchk_m   = 2.5e-7   ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost
+   REAL(wp), PARAMETER ::   rchk_s   = 2.5e-6   ! g/m2/s  <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
+   REAL(wp), PARAMETER ::   rchk_t   = 7.5e-2   ! W/m2    <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
 
    ! for drift outputs
    CHARACTER(LEN=50)   ::   clname="icedrift_diagnostics.ascii"   ! ascii filename
@@ -75,7 +75,7 @@ CONTAINS
       !!
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zchk_m, zchk_s, zchk_t) determine violations
+      !!              The thresholds (rchk_m, rchk_s, rchk_t) determine violations
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10
       !!              and a heat content of 3e5 J/kg (=latent heat of fusion)
       !!-------------------------------------------------------------------
@@ -145,11 +145,11 @@ CONTAINS
 
          IF( lwp ) THEN
             ! check conservation issues
-            IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zchk3(10) ) &
+            IF( ABS(zdiag_mass) > rchk_m * rn_icechk_glo * zchk3(10) ) &
                &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rDt_ice
-            IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zchk3(10) ) &
+            IF( ABS(zdiag_salt) > rchk_s * rn_icechk_glo * zchk3(10) ) &
                &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rDt_ice
-            IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zchk3(10) ) &
+            IF( ABS(zdiag_heat) > rchk_t * rn_icechk_glo * zchk3(10) ) &
                &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rDt_ice
             ! check negative values
             IF( zchk4(1) < 0. )   WRITE(numout,*)   cd_routine,' : violation v_i  < 0        = ',zchk4(1)
@@ -164,9 +164,9 @@ CONTAINS
             IF( zchk3(7)>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &
                &                  WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zchk3(7)
             ! check if advection scheme is conservative
-            IF( ABS(zchk3(8)) > zchk_m * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
+            IF( ABS(zchk3(8)) > rchk_m * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
                &                  WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zchk3(8) * rDt_ice
-            IF( ABS(zchk3(9)) > zchk_t * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
+            IF( ABS(zchk3(9)) > rchk_t * rn_icechk_glo * zchk3(10) .AND. cd_routine == 'icedyn_adv' ) &
                &                  WRITE(numout,*)   cd_routine,' : violation adv scheme [J]  = ',zchk3(9) * rDt_ice
          ENDIF
          !
@@ -182,7 +182,7 @@ CONTAINS
       !!
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zchk_m, zchk_s, zchk_t) determine the violations
+      !!              The thresholds (rchk_m, rchk_s, rchk_t) determine the violations
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10
       !!              and a heat content of 3e5 J/kg (=latent heat of fusion)
       !!-------------------------------------------------------------------
@@ -204,11 +204,11 @@ CONTAINS
       zchk(1:4)   = glob_sum_vec( 'icectl', ztmp(:,:,1:4) )
       
       IF( lwp ) THEN
-         IF( ABS(zchk(1)) > zchk_m * rn_icechk_glo * zchk(4) ) &
+         IF( ABS(zchk(1)) > rchk_m * rn_icechk_glo * zchk(4) ) &
             &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zchk(1) * rDt_ice
-         IF( ABS(zchk(2)) > zchk_s * rn_icechk_glo * zchk(4) ) &
+         IF( ABS(zchk(2)) > rchk_s * rn_icechk_glo * zchk(4) ) &
             &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zchk(2) * rDt_ice
-         IF( ABS(zchk(3)) > zchk_t * rn_icechk_glo * zchk(4) ) &
+         IF( ABS(zchk(3)) > rchk_t * rn_icechk_glo * zchk(4) ) &
             &                   WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zchk(3) * rDt_ice
       ENDIF
       !
@@ -259,20 +259,20 @@ CONTAINS
             &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
             &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           &
             &         - pdiag_fv
-         IF( MAXVAL( ABS(zdiag_mass) ) > zchk_m * rn_icechk_cel )   ll_stop_m = .TRUE.
+         IF( MAXVAL( ABS(zdiag_mass) ) > rchk_m * rn_icechk_cel )   ll_stop_m = .TRUE.
          !
          ! -- salt diag -- !
          zdiag_salt =   ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_Dt_ice                                                  &
             &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) &
             &         - pdiag_fs
-         IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel )   ll_stop_s = .TRUE.
+         IF( MAXVAL( ABS(zdiag_salt) ) > rchk_s * rn_icechk_cel )   ll_stop_s = .TRUE.
          !
          ! -- heat diag -- !
          zdiag_heat =   ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_Dt_ice &
             &         + (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                                &
             &            - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr )                                        &
             &         - pdiag_ft
-         IF( MAXVAL( ABS(zdiag_heat) ) > zchk_t * rn_icechk_cel )   ll_stop_t = .TRUE.
+         IF( MAXVAL( ABS(zdiag_heat) ) > rchk_t * rn_icechk_cel )   ll_stop_t = .TRUE.
          !
          ! -- other diags -- !
          ! a_i < 0
diff --git a/src/ICE/icedia.F90 b/src/ICE/icedia.F90
index 290d77f6935484a257de6336c057159e301709db..154d89de1e4311b194a31b1e645aa65cf523a6d9 100644
--- a/src/ICE/icedia.F90
+++ b/src/ICE/icedia.F90
@@ -33,7 +33,7 @@ MODULE icedia
    PUBLIC   ice_dia        ! called by icestp.F90
    PUBLIC   ice_dia_init   ! called in icestp.F90
 
-   REAL(wp), SAVE ::   z1_e1e2  ! inverse of the ocean area
+   REAL(wp), SAVE ::   r1_area  ! inverse of the ocean area
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini                    ! initial volume, salt and heat contents
    REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends
 
@@ -76,7 +76,7 @@ CONTAINS
       ENDIF
 
       IF( kt == nit000 ) THEN
-         z1_e1e2 = 1._wp / glob_sum( 'icedia', e1e2t(:,:) )
+         r1_area = 1._wp / glob_sum( 'icedia', e1e2t(:,:) )
       ENDIF
 
       ztmp(:,:,:) = 0._wp ! should be better coded
@@ -152,8 +152,8 @@ CONTAINS
       CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal  forcing                      (psu*km3 equivalent ocean water)
       CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)
       CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)
-      CALL iom_put( 'ibgfrchfxtop' , frc_temtop * z1_e1e2 * 1.e-20 * kt*rn_Dt ) ! heat on top of ice/snw/ocean      (W/m2)
-      CALL iom_put( 'ibgfrchfxbot' , frc_tembot * z1_e1e2 * 1.e-20 * kt*rn_Dt ) ! heat on top of ocean(below ice)   (W/m2)
+      CALL iom_put( 'ibgfrchfxtop' , frc_temtop * r1_area * 1.e-20 * kt*rn_Dt ) ! heat on top of ice/snw/ocean      (W/m2)
+      CALL iom_put( 'ibgfrchfxbot' , frc_tembot * r1_area * 1.e-20 * kt*rn_Dt ) ! heat on top of ocean(below ice)   (W/m2)
 
       CALL iom_put( 'ibgvol_tot'  , zbg(6)  )
       CALL iom_put( 'sbgvol_tot'  , zbg(7)  )
diff --git a/src/ICE/icedyn_adv_umx.F90 b/src/ICE/icedyn_adv_umx.F90
index f0a4aab4f92b0173bc2e211e3371dc1586f66036..e5113ac90de67aa1d107f9c1aea3275f2c71f9a2 100644
--- a/src/ICE/icedyn_adv_umx.F90
+++ b/src/ICE/icedyn_adv_umx.F90
@@ -44,8 +44,8 @@ MODULE icedyn_adv_umx
    !                                                 then interpolate T at u/v points using the upstream scheme
    LOGICAL            ::   ll_prelim = .FALSE.  ! prelimiter from: Zalesak(1979) eq. 14 => not well defined in 2D
    !
-   REAL(wp)           ::   z1_6   = 1._wp /   6._wp   ! =1/6
-   REAL(wp)           ::   z1_120 = 1._wp / 120._wp   ! =1/120
+   REAL(wp)           ::   r1_6   = 1._wp /   6._wp   ! =1/6
+   REAL(wp)           ::   r1_120 = 1._wp / 120._wp   ! =1/120
    !
    INTEGER, ALLOCATABLE, DIMENSION(:,:,:) ::   imsk_small, jmsk_small
    !
@@ -936,7 +936,7 @@ CONTAINS
 !!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
                   &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
-                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
+                  &        + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
                   &                                               - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
             END_2D
          END DO
@@ -950,7 +950,7 @@ CONTAINS
 !!rachid          zdx2 = e1u(ji,jj) * e1t(ji,jj)
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (         (                      pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
                   &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
-                  &        + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
+                  &        + r1_6 * zdx2 * ( zcu*zcu - 1._wp ) *    (                      ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
                   &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) )
             END_2D
          END DO
@@ -967,9 +967,9 @@ CONTAINS
                zdx4 = zdx2 * zdx2
                pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * (        (                       pt  (ji+1,jj,jl) + pt  (ji,jj,jl)     &
                   &                                                            - zcu   * ( pt  (ji+1,jj,jl) - pt  (ji,jj,jl) ) ) &
-                  &        + z1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
+                  &        + r1_6   * zdx2 * ( zcu*zcu - 1._wp ) * (                       ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl)     &
                   &                                                   - 0.5_wp * zcu   * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) &
-                  &        + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     &
+                  &        + r1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl)     &
                   &                                               - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) )
             END_2D
          END DO
@@ -1072,7 +1072,7 @@ CONTAINS
 !!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
                   &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
-                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
+                  &        + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
                   &                                               - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
             END_2D
          END DO
@@ -1085,7 +1085,7 @@ CONTAINS
 !!rachid          zdy2 = e2v(ji,jj) * e2t(ji,jj)
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (      (                         pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
                   &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
-                  &        + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
+                  &        + r1_6 * zdy2 * ( zcv*zcv - 1._wp ) * (                         ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
                   &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) )
             END_2D
          END DO
@@ -1102,9 +1102,9 @@ CONTAINS
                zdy4 = zdy2 * zdy2
                pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * (                              ( pt  (ji,jj+1,jl) + pt  (ji,jj,jl)     &
                   &                                                            - zcv   * ( pt  (ji,jj+1,jl) - pt  (ji,jj,jl) ) ) &
-                  &        + z1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
+                  &        + r1_6   * zdy2 * ( zcv*zcv - 1._wp ) * (                       ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl)     &
                   &                                                   - 0.5_wp * zcv   * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) &
-                  &        + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
+                  &        + r1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl)     &
                   &                                               - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) )
             END_2D
          END DO
diff --git a/src/ICE/icedyn_rdgrft.F90 b/src/ICE/icedyn_rdgrft.F90
index 2a8437e06da2a9a3b3c09806d64f8c3087bbd8b8..306fc5e5c31cc80c0dc63b6d7068eeea184ad1c1 100644
--- a/src/ICE/icedyn_rdgrft.F90
+++ b/src/ICE/icedyn_rdgrft.F90
@@ -788,12 +788,17 @@ CONTAINS
       !!----------------------------------------------------------------------
       INTEGER             ::   ji, jj, jl  ! dummy loop indices
       REAL(wp)            ::   z1_3        ! local scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk, zworka           ! temporary array used here
       !!clem
       LOGICAL                   ::   ln_str_R75
       REAL(wp)                  ::   zhi, zcp, rn_pe_rdg
       REAL(wp), DIMENSION(jpij) ::   zstrength, zaksum   ! strength in 1D      
       !!----------------------------------------------------------------------
+      ! prepare the mask
+      at_i(:,:) = SUM( a_i, dim=3 )
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         zmsk(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10  ) ) ! 1 if ice    , 0 if no ice
+      END_2D
       !
       SELECT CASE( nice_str )          !--- Set which ice strength is chosen
 
@@ -805,7 +810,6 @@ CONTAINS
          strength(:,:) = 0._wp
          !
          ! Identify grid cells with ice
-         at_i(:,:) = SUM( a_i, dim=3 )
          npti = 0   ;   nptidx(:) = 0
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
             IF ( at_i(ji,jj) > epsi10 ) THEN
@@ -859,10 +863,10 @@ CONTAINS
          ENDIF
          !
       CASE ( np_strh79 )           !== Hibler(1979)'s method ==!
-         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
+         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - at_i(:,:) ) ) * zmsk(:,:)
          !
       CASE ( np_strcst )           !== Constant strength ==!
-         strength(:,:) = rn_str
+         strength(:,:) = rn_str * zmsk(:,:)
          !
       END SELECT
       !
@@ -879,7 +883,7 @@ CONTAINS
          END_2D
 
          DO_2D( 0, 0, 0, 0 )
-            strength(ji,jj) = zworka(ji,jj)
+            strength(ji,jj) = zworka(ji,jj) * zmsk(ji,jj)
          END_2D
          CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp )
          !
diff --git a/src/ICE/icedyn_rhg_evp.F90 b/src/ICE/icedyn_rhg_evp.F90
index 21c9f53a9cb5809ee3c9581b7d4519b7cbf76c53..2efe027c4f8e0cc917a521ec1f35a20e6f226c70 100644
--- a/src/ICE/icedyn_rhg_evp.F90
+++ b/src/ICE/icedyn_rhg_evp.F90
@@ -158,7 +158,7 @@ CONTAINS
       REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
       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) ::   zmsk, zmsk00, zmsk15
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence
 
@@ -189,6 +189,7 @@ CONTAINS
       ! for diagnostics and convergence tests
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
+         zmsk  (ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10  ) ) ! 1 if ice    , 0 if no ice
       END_2D
       IF( nn_rhg_chkcvg > 0 ) THEN
          DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
@@ -383,10 +384,10 @@ CONTAINS
             zdt2 = zdt * zdt
 
             ! delta at T points
-            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
+            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj)        ! zmsk is for reducing cpu
 
             ! P/delta at T points
-            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) * zmsk(ji,jj) ! zmsk is for reducing cpu
 
          END_2D
          CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp, zp_delt, 'T', 1.0_wp )
@@ -421,8 +422,10 @@ CONTAINS
             ENDIF
 
             ! stress at T points (zkt/=0 if landfast)
-            zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
-            zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+            zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) &
+               &         * z1_alph1 * zmsk(ji,jj) ! zmsk is for reducing cpu
+            zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) &
+               &         * z1_alph2 * zmsk(ji,jj) ! zmsk is for reducing cpu
 
          END_2D
 
@@ -449,7 +452,8 @@ CONTAINS
             zp_delf = 0.25_wp * ( (zp_delt(ji,jj) + zp_delt(ji+1,jj)) + (zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1)) )
 
             ! stress at F points (zkt/=0 if landfast)
-            zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
+            zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) &
+               &         * z1_alph2 * zmsk(ji,jj) ! zmsk is for reducing cpu
 
          END_2D
 
@@ -746,15 +750,15 @@ CONTAINS
             &   ) * 0.25_wp * r1_e1e2t(ji,jj)
 
          ! shear at T points
-         pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+         pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) * zmsk(ji,jj)
 
          ! divergence at T points
          pdivu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
             &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-            &             ) * r1_e1e2t(ji,jj)
+            &             ) * r1_e1e2t(ji,jj) * zmsk(ji,jj)
 
          ! delta at T points
-         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta
+         zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) * zmsk(ji,jj) ! delta
          rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
          pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
 
@@ -851,8 +855,8 @@ CONTAINS
             zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
          END_2D
          !
-         CALL iom_put( 'sig1_pnorm' , zsig1_p )
-         CALL iom_put( 'sig2_pnorm' , zsig2_p )
+         CALL iom_put( 'sig1_pnorm' , zsig1_p * zmsk00 )
+         CALL iom_put( 'sig2_pnorm' , zsig2_p * zmsk00 )
 
          DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
 
diff --git a/src/ICE/iceupdate.F90 b/src/ICE/iceupdate.F90
index 9484a039044977c5590f17cc351db80aa003dd58..a973312f04364781465403618f175cf52c6170ed 100644
--- a/src/ICE/iceupdate.F90
+++ b/src/ICE/iceupdate.F90
@@ -339,7 +339,7 @@ CONTAINS
       REAL(wp) ::   zat_v, zvtau_ice, zv_t, zrhoco  !   -      -
       REAL(wp) ::   zflagi                          !   -      -
       !!---------------------------------------------------------------------
-      IF( ln_timing )   CALL timing_start('ice_update')
+      IF( ln_timing )   CALL timing_start('iceupdate')
 
       IF( kt == nit000 .AND. lwp ) THEN
          WRITE(numout,*)
@@ -391,7 +391,7 @@ CONTAINS
       END_2D
       CALL lbc_lnk( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp )   ! lateral boundary condition
       !
-      IF( ln_timing )   CALL timing_stop('ice_update')
+      IF( ln_timing )   CALL timing_stop('iceupdate')
       !
    END SUBROUTINE ice_update_tau
 
diff --git a/tests/ICE_AGRIF/EXPREF/1_namelist_cfg b/tests/ICE_AGRIF/EXPREF/1_namelist_cfg
index ffbec2ea8d93641b7948c63504218313540ba16c..d21e97691895853f49473251884ed819eb9c6880 100644
--- a/tests/ICE_AGRIF/EXPREF/1_namelist_cfg
+++ b/tests/ICE_AGRIF/EXPREF/1_namelist_cfg
@@ -104,6 +104,8 @@
 !-----------------------------------------------------------------------
 &namagrif      !  AGRIF zoom                                            ("key_agrif")
 !-----------------------------------------------------------------------
+   ln_init_chfrpar = .true.  !  initialize child grids from parent
+   ln_chk_bathy    = .false. !  =T  check the parent bathymetry
 /
 !!======================================================================
 !!                ***  Top/Bottom boundary condition  ***             !!
diff --git a/tests/ICE_AGRIF/EXPREF/namelist_cfg b/tests/ICE_AGRIF/EXPREF/namelist_cfg
index 4c24668745bbd9ae0d486d8c664885155bb573c9..4e54a36bc4f169be6b026b7d16ca7ba55f92446e 100644
--- a/tests/ICE_AGRIF/EXPREF/namelist_cfg
+++ b/tests/ICE_AGRIF/EXPREF/namelist_cfg
@@ -104,6 +104,8 @@
 !-----------------------------------------------------------------------
 &namagrif      !  AGRIF zoom                                            ("key_agrif")
 !-----------------------------------------------------------------------
+   ln_init_chfrpar = .true.  !  initialize child grids from parent
+   ln_chk_bathy    = .false. !  =T  check the parent bathymetry
 /
 !!======================================================================
 !!                ***  Top/Bottom boundary condition  ***             !!