Index: NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/iceistate.F90	(revision 11960)
+++ NEMO/branches/2019/dev_r11943_MERGE_2019/src/ICE/iceistate.F90	(revision 12340)
@@ -61,4 +61,6 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
    !   
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -269,12 +271,10 @@
          ! select ice covered grid points
          npti = 0 ; nptidx(:) = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
-                  npti         = npti  + 1
-                  nptidx(npti) = (jj - 1) * jpi + ji
-               ENDIF
-            END DO
-         END DO
+         DO_2D_11_11
+            IF ( zht_i_ini(ji,jj) > 0._wp ) THEN
+               npti         = npti  + 1
+               nptidx(npti) = (jj - 1) * jpi + ji
+            ENDIF
+         END_2D
 
          ! move to 1D arrays: (jpi,jpj) -> (jpi*jpj)
@@ -321,38 +321,28 @@
          CALL ice_var_salprof ! for sz_i
          DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
-                  v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
-                  sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
-               END DO
-            END DO
+            DO_2D_11_11
+               v_i (ji,jj,jl) = h_i(ji,jj,jl) * a_i(ji,jj,jl)
+               v_s (ji,jj,jl) = h_s(ji,jj,jl) * a_i(ji,jj,jl)
+               sv_i(ji,jj,jl) = MIN( MAX( rn_simin , s_i(ji,jj,jl) ) , rn_simax ) * v_i(ji,jj,jl)
+            END_2D
          END DO
          !
          DO jl = 1, jpl
-            DO jk = 1, nlay_s
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
-                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
-                        &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
-                  END DO
-               END DO
-            END DO
+            DO_3D_11_11( 1, nlay_s )
+               t_s(ji,jj,jk,jl) = zts_3d(ji,jj,jl)
+               e_s(ji,jj,jk,jl) = zswitch(ji,jj) * v_s(ji,jj,jl) * r1_nlay_s * &
+                  &               rhos * ( rcpi * ( rt0 - t_s(ji,jj,jk,jl) ) + rLfus )
+            END_3D
          END DO
          !
          DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
-                     ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
-                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
-                        &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
-                        &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
-                        &                       - rcp   * ( ztmelts - rt0 ) )
-                  END DO
-               END DO
-            END DO
+            DO_3D_11_11( 1, nlay_i )
+               t_i (ji,jj,jk,jl) = zti_3d(ji,jj,jl) 
+               ztmelts          = - rTmlt * sz_i(ji,jj,jk,jl) + rt0 ! melting temperature in K
+               e_i(ji,jj,jk,jl) = zswitch(ji,jj) * v_i(ji,jj,jl) * r1_nlay_i * &
+                  &               rhoi * (  rcpi  * ( ztmelts - t_i(ji,jj,jk,jl) ) + &
+                  &                         rLfus * ( 1._wp - (ztmelts-rt0) / MIN( (t_i(ji,jj,jk,jl)-rt0), -epsi20 ) ) &
+                  &                       - rcp   * ( ztmelts - rt0 ) )
+            END_3D
          END DO
 
