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 8369)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8370)
@@ -266,4 +266,6 @@
             END IF
             !
+            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
@@ -279,14 +281,4 @@
          !
       END DO
-      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
-
-      ! Change thickness to volume
-      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
-      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
-      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
-      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_thd_da', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      IF( ln_limdA)           CALL lim_thd_da                       ! --- lateral melting --- !
 
       ! Change thickness to volume
@@ -295,5 +287,5 @@
       v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
       smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
-      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_thd_da', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       ! update ice age (in case a_i changed, i.e. becomes 0 or lateral melting in monocat)
@@ -466,4 +458,5 @@
          CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res          )
          CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr          )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam          )
          !
          CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog          )
@@ -474,4 +467,5 @@
          CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res          )
          CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub          )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam          )
          !
          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd          )
@@ -528,4 +522,5 @@
          CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res          )
          CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr          )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam          )
          !
          CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog          )
@@ -536,5 +531,6 @@
          CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res          )
          CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub          )
-
+         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam          )
+         !
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd          )
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr          )
@@ -551,8 +547,7 @@
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem   )
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out          )
-
+         !
          CALL tab_1d_2d( nidx, idxice(1:nidx), qns_ice_1d  (1:nidx), qns_ice(:,:,jl)  )
          CALL tab_1d_2d( nidx, idxice(1:nidx), ftr_ice_1d  (1:nidx), ftr_ice(:,:,jl)  )
-         !
          !
          ! SIMIP diagnostics         
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90	(revision 8369)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90	(revision 8370)
@@ -14,8 +14,6 @@
    USE par_oce        ! ocean parameters
    USE phycst         ! physical constants (ocean directory)
-   USE sbc_oce , ONLY : sst_m
    USE ice            ! LIM variables
    USE thd_ice        ! thermodynamic sea-ice variables
-   USE limtab         ! 1D <==> 2D transformation
    !
    USE lib_mpp        ! MPP library
@@ -100,9 +98,5 @@
       !!              Phil. Trans. R. Soc. A, 373(2052), 20140167.
       !!---------------------------------------------------------------------
-      INTEGER             ::   ji, jj, jk, jl     ! dummy loop indices
-      REAL(wp) ::   ztmelts             ! local scalar
-      REAL(wp) ::   zEi          ! specific enthalpy of sea ice (J/kg)
-      REAL(wp) ::   zEw          ! specific enthalpy of exchanged water (J/kg)
-      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
+      INTEGER  ::   ji     ! dummy loop indices
       REAL(wp)            ::   zastar, zdfloe, zperi, zwlat, zda
       REAL(wp), PARAMETER ::   zdmax = 300._wp
@@ -114,99 +108,45 @@
       !!---------------------------------------------------------------------
 
-      ! select ice covered grid points
-      nidx = 0 ; idxice(:) = 0
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF ( at_i(ji,jj) > epsi10 ) THEN
-               nidx         = nidx  + 1
-               idxice(nidx) = (jj - 1) * jpi + ji
+      !------------------------------------------------------------!
+      !------------------------------------------------------------!
+      zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) )
+      DO ji = 1, nidx   
+         ! --- Calculate reduction of total sea ice concentration --- !
+         zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta         ! Mean floe caliper diameter [m]
+         zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2]
+         zwlat  = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2  ! Melt speed rate [m/s]
+         
+         zda_tot(ji) = MIN( zwlat * zperi * rdt_ice, at_i_1d(ji) )                 ! sea ice concentration decrease (>0)
+      
+         ! --- Distribute reduction among ice categories and calculate associated ice-ocean fluxes --- !
+         IF( a_i_1d(ji) > 0._wp ) THEN
+            ! decrease of concentration for the category jl
+            !    each category contributes to melting in proportion to its concentration
+            zda = MIN( a_i_1d(ji), zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji) )
+            
+            ! Contribution to salt flux
+            sfx_lam_1d(ji) = sfx_lam_1d(ji) + rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice
+            
+            ! Contribution to heat flux into the ocean [W.m-2], (<0)  
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) - zda * r1_rdtice * ( ht_i_1d(ji) * r1_nlay_i * SUM( e_i_1d(ji,1:nlay_i) )  &
+                                                                + ht_s_1d(ji) * r1_nlay_s * SUM( e_s_1d(ji,1:nlay_s) ) ) 
+            
+            ! Contribution to mass flux
+            wfx_lam_1d(ji) =  wfx_lam_1d(ji) + zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) )
+            
+            !! adjust e_i ???
+!!            e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp - zda / a_i_1d(ji) )
+!!            e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp - zda / a_i_1d(ji) )
+            
+            ! new concentration
+            a_i_1d(ji) = a_i_1d(ji) - zda
+            
+            ! ensure that ht_i = 0 where a_i = 0
+            IF( a_i_1d(ji) == 0._wp ) THEN
+               ht_i_1d(ji) = 0._wp
+               ht_s_1d(ji) = 0._wp
             ENDIF
-         END DO
+         ENDIF
       END DO
-
-      IF( nidx > 0 ) THEN
-         !------------------------------------------------------------!
-         ! --- Calculate reduction of total sea ice concentration --- !
-         !------------------------------------------------------------!
-         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i  )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d(1:nidx), t_bo  )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d (1:nidx), sst_m )
-         
-         zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) )
-         DO ji = 1, nidx   
-            zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta         ! Mean floe caliper diameter [m]
-            zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2]
-            zwlat  = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2  ! Melt speed rate [m/s]
-            
-            zda_tot(ji) = - MIN( zwlat * zperi * rdt_ice, at_i_1d(ji) )               ! sea ice concentration decrease
-         END DO
-         
-         !---------------------------------------------------------------------------------------------!
-         ! --- Distribute reduction among ice categories and calculate associated ice-ocean fluxes --- !
-         !---------------------------------------------------------------------------------------------!
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_lam_1d(1:nidx), sfx_lam )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_lam_1d(1:nidx), wfx_lam )
-         
-         DO jl = 1, jpl
-            
-            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d(1:nidx), ht_i(:,:,jl) )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), ht_s_1d(1:nidx), ht_s(:,:,jl) )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,jl) )
-            DO jk = 1, nlay_i
-               CALL tab_2d_1d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,jl) )
-            END DO
-            DO jk = 1, nlay_s
-               CALL tab_2d_1d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,jl) )
-            END DO
-            
-            DO ji = 1, nidx
-               IF( a_i_1d(ji) > 0._wp ) THEN
-                  ! decrease of concentration for the category jl
-                  !    each category contributes to melting in proportion to its concentration
-                  zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji)
-                  
-                  ! Contribution to salt flux
-                  sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic *  ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice
-                  
-                  ! Contribution to heat flux into the ocean [W.m-2], (<0)  
-                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice 
-                  
-                  ! Contribution to mass flux
-                  wfx_lam_1d(ji) =  wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) )
-                  
-                  !! adjust e_i ???
-                  e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
-                  e_s_1d(ji,1)        = e_s_1d(ji,1)        * ( 1._wp + zda_tot(ji) / at_i_1d(ji) )
-                  
-                  ! new concentration
-                  a_i_1d(ji) = a_i_1d(ji) + zda
-                  
-                  ! ensure that ht_i = 0 where a_i = 0
-                  IF( a_i_1d(ji) == 0._wp ) THEN
-                     ht_i_1d(ji) = 0._wp
-                     ht_s_1d(ji) = 0._wp
-                  ENDIF
-               ENDIF
-            END DO
-            
-            CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d    (1:nidx), a_i (:,:,jl) )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d   (1:nidx), ht_i(:,:,jl) )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), ht_s_1d   (1:nidx), ht_s(:,:,jl) )
-            DO jk = 1, nlay_s
-               CALL tab_1d_2d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,jl) )
-            END DO
-            DO jk = 1, nlay_i
-               CALL tab_1d_2d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,jl) )
-            END DO
-            
-         END DO
-         
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_lam_1d(1:nidx), sfx_lam )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_lam_1d(1:nidx), wfx_lam )
-         
-      ENDIF
       !
    END SUBROUTINE lim_thd_da
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8369)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8370)
@@ -20,4 +20,5 @@
    USE ice            ! LIM variables
    USE thd_ice        ! LIM thermodynamics
+   !
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
