Index: NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90	(revision 13797)
+++ NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90	(revision 13799)
@@ -13,5 +13,5 @@
    !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
    !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
-   !!                 !  2019     (S. Rynders)  change into eap rheology from
+   !!                 !  2019     (S. Rynders, Y. Aksenov, C. Rousset)  change into eap rheology from
    !!                                           CICE code (Tsamados, Heorton) 
    !!----------------------------------------------------------------------
@@ -450,12 +450,8 @@
 
             ! structure tensor update
-!!$               IF (mod(jter,10) == 0) THEN
                CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11,  zmresult12)
 
-!!$                  paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + zp5kth + zmresult11) * z1dtevpkth ! implicit
-!!$                  paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  + zmresult12) * z1dtevpkth ! implicit
                paniso_11(ji,jj) = (paniso_11(ji,jj)  + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit
                paniso_12(ji,jj) = (paniso_12(ji,jj)                     + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit
-!!$               ENDIF
 
             IF (jter == nn_nevp) THEN
@@ -510,8 +506,5 @@
             
             ! 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) + zstress12tmpF * zalph1  ) * z1_alph1
             zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1
-            !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1
 
          END_2D
@@ -599,5 +592,5 @@
             !
 #if defined key_agrif
-!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+!!          CALL agrif_interp_ice( 'V', jter, nn_nevp )
             CALL agrif_interp_ice( 'V' )
 #endif
@@ -656,5 +649,5 @@
             !
 #if defined key_agrif
-!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+!!          CALL agrif_interp_ice( 'U', jter, nn_nevp )
             CALL agrif_interp_ice( 'U' )
 #endif
@@ -714,5 +707,5 @@
             !
 #if defined key_agrif
-!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+!!          CALL agrif_interp_ice( 'U', jter, nn_nevp )
             CALL agrif_interp_ice( 'U' )
 #endif
@@ -770,5 +763,5 @@
             !
 #if defined key_agrif
-!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+!!          CALL agrif_interp_ice( 'V', jter, nn_nevp )
             CALL agrif_interp_ice( 'V' )
 #endif
@@ -1458,6 +1451,6 @@
                CALL iom_get( numrir, jpdom_auto, 'stress2_i' , stress2_i , cd_type = 'T' )
                CALL iom_get( numrir, jpdom_auto, 'stress12_i', stress12_i, cd_type = 'F' )
-               CALL iom_get( numrir, jpdom_auto, 'aniso_11'  , aniso_11  , cd_type = 'T'   )
-               CALL iom_get( numrir, jpdom_auto, 'aniso_12'  , aniso_12  , cd_type = 'T'   )
+               CALL iom_get( numrir, jpdom_auto, 'aniso_11'  , aniso_11  , cd_type = 'T' )
+               CALL iom_get( numrir, jpdom_auto, 'aniso_12'  , aniso_12  , cd_type = 'T' )
             ELSE                                     ! start rheology from rest
                IF(lwp) WRITE(numout,*)
Index: NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90	(revision 13797)
+++ NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_hgr.F90	(revision 13799)
@@ -28,5 +28,4 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
-#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Index: NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90	(revision 13797)
+++ NEMO/branches/2020/dev_r13787_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/usrdef_sbc.F90	(revision 13799)
@@ -34,4 +34,5 @@
 
    !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -119,11 +120,9 @@
       IF( kt==nit000 .AND. lwp)   WRITE(numout,*)' usrdef_sbc_ice : ICE_RHEO case: analytical stress forcing'
 
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1   
-            ! wind spins up over 6 hours, factor 1000 to balance the units
-            windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
-            windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
-         END DO
-      END DO
+      DO_2D( 0, 0, 0, 0 )
+         ! wind spins up over 6 hours, factor 1000 to balance the units
+         windu(ji,jj) = Umax/sqrt(d*1000)*(d-2*mig(ji)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*min(kt*30./21600,1.)
+         windv(ji,jj) = Umax/sqrt(d*1000)*(d-2*mjg(jj)*res)/((d-2*mig(ji)*res)**2+(d-2*mjg(jj)*res)**2*Rwind**2)**(1/4)*Rwind*min(kt*30./21600,1.)
+      END_2D
       CALL lbc_lnk_multi( 'usrdef_sbc', windu, 'U', -1., windv, 'V', -1. )
 
@@ -134,11 +133,9 @@
       ! ------------------------------------------------------------ !
       ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1   
-            zwndi_t = (  windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
-            zwndj_t = (  windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) )  )
-            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
-         END DO
-      END DO
+      DO_2D( 0, 0, 0, 0 )
+         zwndi_t = (  windu(ji,jj) - r_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
+         zwndj_t = (  windv(ji,jj) - r_vfac * 0.5 * ( v_ice(ji,jj-1) + v_ice(ji,jj) )  )
+         wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
+      END_2D
       CALL lbc_lnk( 'usrdef_sbc', wndm_ice, 'T',  1. )
 
@@ -152,12 +149,10 @@
       ! ------------------------------------------------------------ !
       ! C-grid ice dynamics :   U & V-points (same as ocean)
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1   
-            utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
-               &          * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) )
-            vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
-               &          * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
-         END DO
-      END DO
+      DO_2D( 0, 0, 0, 0 )
+         utau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
+            &          * ( 0.5 * (windu(ji+1,jj) + windu(ji,jj) ) - r_vfac * u_ice(ji,jj) )
+         vtau_ice(ji,jj) = 0.5 * zrhoa * Cd_atm * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
+            &          * ( 0.5 * (windv(ji,jj+1) + windv(ji,jj) ) - r_vfac * v_ice(ji,jj) )
+      END_2D
       CALL lbc_lnk_multi( 'usrdef_sbc', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
       !
