Index: NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- NEMO/trunk/src/ICE/icevar.F90	(revision 11732)
+++ NEMO/trunk/src/ICE/icevar.F90	(revision 12377)
@@ -82,4 +82,6 @@
    END INTERFACE
 
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -115,13 +117,18 @@
       !
       ato_i(:,:) = 1._wp - at_i(:,:)         ! open water fraction  
-
+      !
+      !!GS: tm_su always needed by ABL over sea-ice
+      ALLOCATE( z1_at_i(jpi,jpj) )
+      WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
+      ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
+      END WHERE
+      tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
+      WHERE( at_i(:,:)<=epsi20 ) tm_su(:,:) = rt0
+      !
       ! The following fields are calculated for diagnostics and outputs only
       ! ==> Do not use them for other purposes
       IF( kn > 1 ) THEN
          !
-         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
-         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
-         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
-         END WHERE
+         ALLOCATE( z1_vt_i(jpi,jpj) , z1_vt_s(jpi,jpj) )
          WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
          ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
@@ -136,5 +143,4 @@
          !         
          !                          ! mean temperature (K), salinity and age
-         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:)
          om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(:,:)
@@ -154,5 +160,4 @@
          !                           ! put rt0 where there is no ice
          WHERE( at_i(:,:)<=epsi20 )
-            tm_su(:,:) = rt0
             tm_si(:,:) = rt0
             tm_i (:,:) = rt0
@@ -165,7 +170,9 @@
          END WHERE         
          !
-         DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s )
+         DEALLOCATE( z1_vt_i , z1_vt_s )
          !
       ENDIF
+      !
+      DEALLOCATE( z1_at_i )
       !
    END SUBROUTINE ice_var_agg
@@ -236,22 +243,18 @@
       zlay_i   = REAL( nlay_i , wp )    ! number of layers
       DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
-                     !
-                     ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3]
-                     ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C]
-                     ! Conversion q(S,T) -> T (second order equation)
-                     zbbb             = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
-                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
-                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
-                     !
-                  ELSE                                   !--- no ice
-                     t_i(ji,jj,jk,jl) = rt0
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, nlay_i )
+            IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
+               !
+               ze_i             =   e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i             ! Energy of melting e(S,T) [J.m-3]
+               ztmelts          = - sz_i(ji,jj,jk,jl) * rTmlt                                 ! Ice layer melt temperature [C]
+               ! Conversion q(S,T) -> T (second order equation)
+               zbbb             = ( rcp - rcpi ) * ztmelts + ze_i * r1_rhoi - rLfus
+               zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * rcpi * rLfus * ztmelts , 0._wp) )
+               t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
+               !
+            ELSE                                   !--- no ice
+               t_i(ji,jj,jk,jl) = rt0
+            ENDIF
+         END_3D
       END DO
 
@@ -344,25 +347,19 @@
          z1_dS = 1._wp / ( zsi1 - zsi0 )
          DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
-                  !                             ! force a constant profile when SSS too low (Baltic Sea)
-                  IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp  
-               END DO
-            END DO
+            DO_2D_11_11
+               zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - s_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
+               !                             ! force a constant profile when SSS too low (Baltic Sea)
+               IF( 2._wp * s_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp  
+            END_2D
          END DO
          !
          ! Computation of the profile
          DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     !                          ! linear profile with 0 surface value
-                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
-                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
-                     sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
-                  END DO
-               END DO
-            END DO
+            DO_3D_11_11( 1, nlay_i )
+               !                          ! linear profile with 0 surface value
+               zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * h_i(ji,jj,jl) * r1_nlay_i
+               zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * s_i(ji,jj,jl)     ! weighting the profile
+               sz_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
+            END_3D
          END DO
          !
@@ -489,51 +486,41 @@
          ! Zap ice energy and use ocean heat to melt ice
          !-----------------------------------------------------------------
-         DO jk = 1, nlay_i
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  ! update exchanges with ocean
-                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
-                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
-               END DO
-            END DO
-         END DO
-         !
-         DO jk = 1, nlay_s
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  ! update exchanges with ocean
-                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
-                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
-                  t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, nlay_i )
+            ! update exchanges with ocean
+            hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+            e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
+            t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
+         END_3D
+         !
+         DO_3D_11_11( 1, nlay_s )
+            ! update exchanges with ocean
+            hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+            e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * zswitch(ji,jj)
+            t_s(ji,jj,jk,jl) = t_s(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
+         END_3D
          !
          !-----------------------------------------------------------------
          ! zap ice and snow volume, add water and salt to ocean
          !-----------------------------------------------------------------
-         DO jj = 1 , jpj
-            DO ji = 1 , jpi
-               ! update exchanges with ocean
-               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice
-               !
-               a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
-               v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
-               v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
-               t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
-               oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
-               sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
-               !
-               h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
-               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
-               !
-               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
-               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
-               !
-            END DO
-         END DO
+         DO_2D_11_11
+            ! update exchanges with ocean
+            sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl)   * rhoi * r1_rdtice
+            wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl)   * rhoi * r1_rdtice
+            wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhos * r1_rdtice
+            !
+            a_i  (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj)
+            v_i  (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj)
+            v_s  (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj)
+            t_su (ji,jj,jl) = t_su(ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
+            oa_i (ji,jj,jl) = oa_i(ji,jj,jl) * zswitch(ji,jj)
+            sv_i (ji,jj,jl) = sv_i(ji,jj,jl) * zswitch(ji,jj)
+            !
+            h_i (ji,jj,jl) = h_i (ji,jj,jl) * zswitch(ji,jj)
+            h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
+            !
+            a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
+            v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+            !
+         END_2D
          !
       END DO 
@@ -587,45 +574,35 @@
          ! zap ice energy and send it to the ocean
          !----------------------------------------
-         DO jk = 1, nlay_i
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
-                     pe_i(ji,jj,jk,jl) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END DO
-         !
-         DO jk = 1, nlay_s
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
-                     hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
-                     pe_s(ji,jj,jk,jl) = 0._wp
-                  ENDIF
-               END DO
-            END DO
-         END DO
+         DO_3D_11_11( 1, nlay_i )
+            IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+               hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0
+               pe_i(ji,jj,jk,jl) = 0._wp
+            ENDIF
+         END_3D
+         !
+         DO_3D_11_11( 1, nlay_s )
+            IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+               hfx_res(ji,jj)   = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0
+               pe_s(ji,jj,jk,jl) = 0._wp
+            ENDIF
+         END_3D
          !
          !-----------------------------------------------------
          ! zap ice and snow volume, add water and salt to ocean
          !-----------------------------------------------------
-         DO jj = 1 , jpj
-            DO ji = 1 , jpi
-               IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
-                  pv_i   (ji,jj,jl) = 0._wp
-               ENDIF
-               IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
-                  wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
-                  pv_s   (ji,jj,jl) = 0._wp
-               ENDIF
-               IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
-                  sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
-                  psv_i  (ji,jj,jl) = 0._wp
-               ENDIF
-            END DO
-         END DO
+         DO_2D_11_11
+            IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+               wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt
+               pv_i   (ji,jj,jl) = 0._wp
+            ENDIF
+            IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN
+               wfx_res(ji,jj)    = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt
+               pv_s   (ji,jj,jl) = 0._wp
+            ENDIF
+            IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN
+               sfx_res(ji,jj)    = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt
+               psv_i  (ji,jj,jl) = 0._wp
+            ENDIF
+         END_2D
          !
       END DO 
