Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8378)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8379)
@@ -83,5 +83,5 @@
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
       REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
       REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
@@ -218,15 +218,4 @@
          END DO
 
-         ! debug point to follow
-         jiindex_1d = 0
-         IF( ln_limctl ) THEN
-            DO ji = mi0(iiceprt), mi1(iiceprt)
-               DO jj = mj0(jiceprt), mj1(jiceprt)
-                  jiindex_1d = (jj - 1) * jpi + ji
-                  WRITE(numout,*) ' lim_thd : Category no : ', jl 
-               END DO
-            END DO
-         ENDIF
-
          IF( lk_mpp )         CALL mpp_ini_ice( nidx , numout )
 
@@ -234,11 +223,5 @@
             !                                                                
                               CALL lim_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- !
-            !
-            DO jk = 1, nlay_i                                       ! --- Change units from J/m2 to J/m3 --- !
-               WHERE( ht_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (ht_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i
-            ENDDO
-            DO jk = 1, nlay_s
-               WHERE( ht_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (ht_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s
-            ENDDO
+            !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
             !
             s_i_new   (1:nidx) = 0._wp ; dh_s_tot (1:nidx) = 0._wp  ! --- some init --- !  (important to have them here) 
@@ -264,17 +247,6 @@
             IF( ln_limdA )    CALL lim_thd_da                       ! --- lateral melting --- !
             !
-            DO jk = 1, nlay_i                                       ! --- Change units from J/m3 to J/m2 --- !
-               e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * ht_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i
-            ENDDO
-            DO jk = 1, nlay_s
-               e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * ht_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s
-            ENDDO
-            !
-            ! Change thickness to volume
-            v_i_1d(1:nidx)   = ht_i_1d(1:nidx) * a_i_1d(1:nidx)
-            v_s_1d(1:nidx)   = ht_s_1d(1:nidx) * a_i_1d(1:nidx)
-            smv_i_1d(1:nidx) = sm_i_1d(1:nidx) * v_i_1d(1:nidx)
-            !
-                              CALL lim_thd_1d2d( jl, 2 )            ! --- Move to 2D arrays --- !
+                              CALL lim_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
+            !                                                       ! --- & Move to 2D arrays --- !
             !
             IF( lk_mpp )      CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
@@ -282,41 +254,18 @@
          !
       END DO
-      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
-
       ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting)
-      DO jl  = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i_b(ji,jj,jl) - epsi10 ) )
-               oa_i(ji,jj,jl) = rswitch * oa_i(ji,jj,jl) * a_i(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 )
-            END DO
-         END DO
-      END DO
+      oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
 
       IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      CALL lim_var_zapsmall
-
+      !
+                         CALL lim_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
+      !                                                  !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
+      !                   
+      IF( jpl > 1 )      CALL lim_itd_th_rem( kt )       ! --- Transport ice between thickness categories --- !
+      !
+      IF( ln_limdO )     CALL lim_thd_lac                ! --- frazil ice growing in leads --- !
+      !
       IF( ln_limctl )    CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print
-      !
-      !------------------------------------------------!
-      !  Transport ice between thickness categories
-      !------------------------------------------------!
-      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      IF( jpl > 1 )      CALL lim_itd_th_rem( kt )
-
-      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      !------------------------------------------------!
-      !  Add frazil ice growing in leads
-      !------------------------------------------------!
-      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      IF( ln_limdO )     CALL lim_thd_lac
-      
-      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      IF( ln_ctl )       CALL lim_prt3D( 'limthd' )  ! Control print
+      IF( ln_ctl )       CALL lim_prt3D( 'limthd' )      ! Control print
       !
       IF( nn_timing == 1 )  CALL timing_stop('limthd')
@@ -480,7 +429,28 @@
          CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m )
          CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m )
+
+         ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
+         DO jk = 1, nlay_i
+            WHERE( ht_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (ht_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i
+         ENDDO
+         DO jk = 1, nlay_s
+            WHERE( ht_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (ht_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s
+         ENDDO
          !
       CASE( 2 )            ! from 1D to 2D
          !
+         ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
+         DO jk = 1, nlay_i
+            e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * ht_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i
+         ENDDO
+         DO jk = 1, nlay_s
+            e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * ht_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s
+         ENDDO
+         !
+         ! Change thickness to volume
+         v_i_1d(1:nidx)   = ht_i_1d(1:nidx) * a_i_1d(1:nidx)
+         v_s_1d(1:nidx)   = ht_s_1d(1:nidx) * a_i_1d(1:nidx)
+         smv_i_1d(1:nidx) = sm_i_1d(1:nidx) * v_i_1d(1:nidx)
+         
          CALL tab_1d_2d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i             )
          CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl)      )
