Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8534)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8559)
@@ -80,5 +80,5 @@
       INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices
       INTEGER  ::   iter     !   -       -
-      REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars
+      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                          ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
@@ -102,5 +102,4 @@
       REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
       REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
-      REAL(wp), DIMENSION(jpij) ::   zat_i_1d    ! total ice fraction    
       REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
       REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
@@ -108,9 +107,6 @@
       REAL(wp), DIMENSION(jpij,jpl) ::   zv_b      ! old volume of ice in category jl
       REAL(wp), DIMENSION(jpij,jpl) ::   za_b      ! old area of ice in category jl
-      REAL(wp), DIMENSION(jpij,jpl) ::   za_i_1d   ! 1-D version of a_i
-      REAL(wp), DIMENSION(jpij,jpl) ::   zv_i_1d   ! 1-D version of v_i
-      REAL(wp), DIMENSION(jpij,jpl) ::   zsmv_i_1d ! 1-D version of smv_i
-
-      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_1d !: 1-D version of e_i
+
+      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i
 
       REAL(wp), DIMENSION(jpi,jpj) ::   zvrel     ! relative ice / frazil velocity
@@ -243,25 +239,24 @@
       IF ( nidx > 0 ) THEN
 
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zat_i_1d  (1:nidx)     , at_i           )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d (1:nidx)      , at_i         )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i  (:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) )
          DO jl = 1, jpl
-            CALL tab_2d_1d( nidx, idxice(1:nidx), za_i_1d  (1:nidx,jl), a_i  (:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), zv_i_1d  (1:nidx,jl), v_i  (:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl)  )
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)   )
-            END DO
-         END DO
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx)     , qlead       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx)     , t_bo        )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx)     , sfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx)     , wfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hicol_1d  (1:nidx)     , hicol       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx)     , zvrel       )
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx)     , hfx_thd     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx)     , hfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx)     , rn_amax_2d  )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx)     , sss_m       )
+               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx) , qlead       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx) , t_bo        )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , hicol       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx) , zvrel       )
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d  )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx) , sss_m       )
 
          !------------------------------------------------------------------------------|
@@ -269,8 +264,10 @@
          !------------------------------------------------------------------------------|
          DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  IF( zv_i_1d(ji,jl) > 0._wp )   ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i )
-               END DO
+            DO jk = 1, nlay_i               
+               WHERE( v_i_2d(1:nidx,jl) > 0._wp )
+                  ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) / v_i_2d(1:nidx,jl) * REAL( nlay_i )
+               ELSEWHERE
+                  ze_i_2d(1:nidx,jk,jl) = 0._wp
+               END WHERE
             END DO
          END DO
@@ -282,11 +279,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_b(1:nidx,:) = zv_i_1d(1:nidx,:) 
-         za_b(1:nidx,:) = za_i_1d(1:nidx,:)
-
-         !----------------------
-         ! Thickness of new ice
-         !----------------------
-         zh_newice(1:nidx) = hicol_1d(1:nidx)
+         zv_b(1:nidx,:) = v_i_2d(1:nidx,:) 
+         za_b(1:nidx,:) = a_i_2d(1:nidx,:)
 
          !----------------------
@@ -309,8 +301,8 @@
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = 1, nidx
-            ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                                         &
-               &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
-               &                       - rcp  *         ( ztmelts - rt0 )  )
+            ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C)
+            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
+               &                       + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+               &                       - rcp  *         ztmelts )
          END DO
 
@@ -318,7 +310,5 @@
          ! Age of new ice
          !----------------
-         DO ji = 1, nidx
-            zo_newice(ji) = 0._wp
-         END DO
+         zo_newice(1:nidx) = 0._wp
 
          !-------------------
@@ -350,9 +340,9 @@
          END DO
          
-         zv_frazb(:) = 0._wp
+         zv_frazb(1:nidx) = 0._wp
          IF( ln_frazil ) THEN
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
             DO ji = 1, nidx
-               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
                zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
                zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
@@ -378,6 +368,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nidx
-            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
+            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -390,15 +380,14 @@
 
          ! find which category to fill
-         zat_i_1d(:) = 0._wp
          DO jl = 1, jpl
             DO ji = 1, nidx
                IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
-                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
-                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
-                  jcat    (ji)    = jl
+                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
+                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
+                  jcat(ji) = jl
                ENDIF
-               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
-            END DO
-         END DO
+            END DO
+         END DO
+         at_i_1d(1:nidx) = SUM( a_i_2d(1:nidx,:), dim=2 )
 
          ! Heat content
@@ -411,8 +400,8 @@
             DO ji = 1, nidx
                jl = jcat(ji)
-               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
-               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
-                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
-                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
+               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
+               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
+                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  &
+                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
             END DO
          END DO
@@ -428,6 +417,6 @@
             DO jk = 1, nlay_i
                DO ji = 1, nidx
-                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
-                  eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
+                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
+                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
                END DO
             END DO
@@ -435,8 +424,8 @@
             ! new volumes including lateral/bottom accretion + residual
             DO ji = 1, nidx
-               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
-               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
-               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)               
-               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
+               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
+               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
+               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
+               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
                ! for remapping
                h_i_old (ji,nlay_i+1) = zv_newfra
@@ -444,5 +433,5 @@
             ENDDO
             ! --- Ice enthalpy remapping --- !
-            CALL ice_thd_ent( ze_i_1d(1:nidx,:,jl) ) 
+            CALL ice_thd_ent( ze_i_2d(1:nidx,:,jl) ) 
          ENDDO
 
@@ -452,6 +441,5 @@
          DO jl = 1, jpl
             DO ji = 1, nidx
-               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
-               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
+               smv_i_2d(ji,jl) = smv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
             END DO
          END DO
@@ -462,7 +450,5 @@
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i 
-               END DO
+               ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) * v_i_2d(1:nidx,jl) * r1_nlay_i 
             END DO
          END DO
@@ -470,16 +456,16 @@
          ! 7) Change 2D vectors to 1D vectors 
          !------------------------------------------------------------------------------!
-         DO jl = 1, jpl
-            CALL tab_1d_2d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl)  )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl)  )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i (:,:,jl)   )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i  (:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) )
+          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)  )
-            END DO
-         END DO
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw  )
+               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw )
          !
       ENDIF ! nidx > 0
