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 8326)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8327)
@@ -83,5 +83,5 @@
       !
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations
+      INTEGER  :: nidx             ! nb of icy pts for vertical thermo calculations
       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
@@ -213,10 +213,10 @@
 
          ! select ice covered grid points
-         nbpb = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
+         nidx = 0 ; idxice(:) = 0
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
                IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
-                  nbpb      = nbpb  + 1
-                  npb(nbpb) = (jj - 1) * jpi + ji
+                  nidx         = nidx  + 1
+                  idxice(nidx) = (jj - 1) * jpi + ji
                ENDIF
             END DO
@@ -234,7 +234,7 @@
          ENDIF
 
-         IF( lk_mpp )         CALL mpp_ini_ice( nbpb , numout )
-
-         IF( nbpb > 0 ) THEN  ! If there is no ice, do nothing.
+         IF( lk_mpp )         CALL mpp_ini_ice( nidx , numout )
+
+         IF( nidx > 0 ) THEN  ! If there is no ice, do nothing.
             !                                                                
             s_i_new   (:) = 0._wp ; dh_s_tot (:) = 0._wp                     ! --- some init --- !
@@ -242,5 +242,5 @@
             dh_snowice(:) = 0._wp ; dh_i_sub (:) = 0._wp
 
-                              CALL lim_thd_1d2d( nbpb, jl, 1 )               ! --- Move to 1D arrays --- !
+                              CALL lim_thd_1d2d( nidx, jl, 1 )               ! --- Move to 1D arrays --- !
             !
             DO jk = 1, nlay_i                                                ! --- Change units from J/m2 to J/m3 --- !
@@ -251,17 +251,17 @@
             ENDDO
             !
-            IF( ln_limdH )    CALL lim_thd_dif( 1, nbpb )                    ! --- Ice/Snow Temperature profile --- !
-            !
-            IF( ln_limdH )    CALL lim_thd_dh( 1, nbpb )                     ! --- Ice/Snow thickness --- !    
-            !
-            IF( ln_limdH )    CALL lim_thd_ent( 1, nbpb, e_i_1d(1:nbpb,:) )  ! --- Ice enthalpy remapping --- !
-            !
-                              CALL lim_thd_sal( 1, nbpb )                    ! --- Ice salinity --- !    
-            !
-                              CALL lim_thd_temp( 1, nbpb )                   ! --- temperature update --- !
+            IF( ln_limdH )    CALL lim_thd_dif( 1, nidx )                    ! --- Ice/Snow Temperature profile --- !
+            !
+            IF( ln_limdH )    CALL lim_thd_dh( 1, nidx )                     ! --- Ice/Snow thickness --- !    
+            !
+            IF( ln_limdH )    CALL lim_thd_ent( 1, nidx, e_i_1d(1:nidx,:) )  ! --- Ice enthalpy remapping --- !
+            !
+                              CALL lim_thd_sal( 1, nidx )                    ! --- Ice salinity --- !    
+            !
+                              CALL lim_thd_temp( 1, nidx )                   ! --- temperature update --- !
             !
             IF( ln_limdH ) THEN
                IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
-                              CALL lim_thd_lam( 1, nbpb )                    ! --- extra lateral melting if monocat --- !
+                              CALL lim_thd_lam( 1, nidx )                    ! --- extra lateral melting if monocat --- !
                END IF
             END IF
@@ -274,5 +274,5 @@
             ENDDO
             !
-                              CALL lim_thd_1d2d( nbpb, jl, 2 )               ! --- Move to 2D arrays --- !
+                              CALL lim_thd_1d2d( nidx, jl, 2 )               ! --- Move to 2D arrays --- !
             !
             IF( lk_mpp )      CALL mpp_comm_free( ncomm_ice ) !RB necessary ??
@@ -282,4 +282,6 @@
 
       IF( ln_limdA)           CALL lim_thd_da                                ! --- lateral melting --- !
+
+      at_i(:,:)    = SUM( a_i(:,:,:), dim=3 )
 
       ! Change thickness to volume
@@ -404,5 +406,5 @@
 
 
-   SUBROUTINE lim_thd_1d2d( nbpb, jl, kn )
+   SUBROUTINE lim_thd_1d2d( nidx, jl, kn )
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE lim_thd_1d2d *** 
@@ -411,5 +413,5 @@
       !!-----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kn       ! 1= from 2D to 1D   ;   2= from 1D to 2D
-      INTEGER, INTENT(in) ::   nbpb     ! size of 1D arrays
+      INTEGER, INTENT(in) ::   nidx     ! size of 1D arrays
       INTEGER, INTENT(in) ::   jl       ! ice cat
       !
@@ -421,143 +423,141 @@
       CASE( 1 )            ! from 2D to 1D
          !
-         CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-         !
-         CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nidx, at_i_1d     (1:nidx), at_i            , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, a_i_1d      (1:nidx), a_i(:,:,jl)     , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, ht_i_1d     (1:nidx), ht_i(:,:,jl)    , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, ht_s_1d     (1:nidx), ht_s(:,:,jl)    , jpi, jpj, idxice(1:nidx) )
+         !
+         CALL tab_2d_1d( nidx, t_su_1d     (1:nidx), t_su(:,:,jl)    , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sm_i_1d     (1:nidx), sm_i(:,:,jl)    , jpi, jpj, idxice(1:nidx) )
          DO jk = 1, nlay_s
-            CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, e_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nidx, t_s_1d(1:nidx,jk), t_s(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
+            CALL tab_2d_1d( nidx, e_s_1d(1:nidx,jk), e_s(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, e_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-         END DO
-         !
-         CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, qevap_ice_1d(1:nbpb), qevap_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb), fr2_i0          , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb), qns_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, ftr_ice_1d (1:nbpb), ftr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, evap_ice_1d (1:nbpb), evap_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) ) 
-         CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, qlead_1d   (1:nbpb), qlead           , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, fhld_1d    (1:nbpb), fhld            , jpi, jpj, npb(1:nbpb) )
-         !
-         CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_snw_sum_1d(1:nbpb), wfx_snw_sum  , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_snw_sub_1d(1:nbpb), wfx_snw_sub  , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_ice_sub_1d(1:nbpb), wfx_ice_sub  , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_err_sub_1d(1:nbpb), wfx_err_sub  , jpi, jpj, npb(1:nbpb) )
-         !
-         CALL tab_2d_1d( nbpb, wfx_bog_1d (1:nbpb), wfx_bog         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_bom_1d (1:nbpb), wfx_bom         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) )
-         !
-         CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sfx_bom_1d (1:nbpb), sfx_bom         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sfx_sum_1d (1:nbpb), sfx_sum         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sfx_sub_1d (1:nbpb), sfx_sub         , jpi, jpj,npb(1:nbpb) )
-         !
-         CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_err_1d (1:nbpb), hfx_err         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_res_1d (1:nbpb), hfx_res         , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_err_dif_1d (1:nbpb), hfx_err_dif , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_err_rem_1d (1:nbpb), hfx_err_rem , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, hfx_out_1d (1:nbpb), hfx_out         , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nidx, t_i_1d(1:nidx,jk), t_i(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
+            CALL tab_2d_1d( nidx, e_i_1d(1:nidx,jk), e_i(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
+            CALL tab_2d_1d( nidx, s_i_1d(1:nidx,jk), s_i(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
+         END DO
+         !
+         CALL tab_2d_1d( nidx, qprec_ice_1d(1:nidx), qprec_ice(:,:) , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, qevap_ice_1d(1:nidx), qevap_ice(:,:,jl) , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, qsr_ice_1d (1:nidx), qsr_ice(:,:,jl) , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, fr1_i0_1d  (1:nidx), fr1_i0          , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, fr2_i0_1d  (1:nidx), fr2_i0          , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, qns_ice_1d (1:nidx), qns_ice(:,:,jl) , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, ftr_ice_1d (1:nidx), ftr_ice(:,:,jl) , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, evap_ice_1d (1:nidx), evap_ice(:,:,jl), jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, dqns_ice_1d(1:nidx), dqns_ice(:,:,jl), jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, t_bo_1d     (1:nidx), t_bo            , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sprecip_1d (1:nidx), sprecip         , jpi, jpj, idxice(1:nidx) ) 
+         CALL tab_2d_1d( nidx, fhtur_1d   (1:nidx), fhtur           , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, fhld_1d    (1:nidx), fhld            , jpi, jpj, idxice(1:nidx) )
+         !
+         CALL tab_2d_1d( nidx, wfx_snw_1d (1:nidx), wfx_snw         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_snw_sum_1d(1:nidx), wfx_snw_sum  , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_sub_1d (1:nidx), wfx_sub         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_snw_sub_1d(1:nidx), wfx_snw_sub  , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_ice_sub_1d(1:nidx), wfx_ice_sub  , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_err_sub_1d(1:nidx), wfx_err_sub  , jpi, jpj, idxice(1:nidx) )
+         !
+         CALL tab_2d_1d( nidx, wfx_bog_1d (1:nidx), wfx_bog         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_bom_1d (1:nidx), wfx_bom         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_sum_1d (1:nidx), wfx_sum         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_sni_1d (1:nidx), wfx_sni         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_res_1d (1:nidx), wfx_res         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_spr_1d (1:nidx), wfx_spr         , jpi, jpj, idxice(1:nidx) )
+         !
+         CALL tab_2d_1d( nidx, sfx_bog_1d (1:nidx), sfx_bog         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_bom_1d (1:nidx), sfx_bom         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_sum_1d (1:nidx), sfx_sum         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_sni_1d (1:nidx), sfx_sni         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_bri_1d (1:nidx), sfx_bri         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_res_1d (1:nidx), sfx_res         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_sub_1d (1:nidx), sfx_sub         , jpi, jpj,idxice(1:nidx) )
+         !
+         CALL tab_2d_1d( nidx, hfx_thd_1d (1:nidx), hfx_thd         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_spr_1d (1:nidx), hfx_spr         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_sum_1d (1:nidx), hfx_sum         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_bom_1d (1:nidx), hfx_bom         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_bog_1d (1:nidx), hfx_bog         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_dif_1d (1:nidx), hfx_dif         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_opw_1d (1:nidx), hfx_opw         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_snw_1d (1:nidx), hfx_snw         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_sub_1d (1:nidx), hfx_sub         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_err_1d (1:nidx), hfx_err         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_res_1d (1:nidx), hfx_res         , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_err_dif_1d (1:nidx), hfx_err_dif , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_err_rem_1d (1:nidx), hfx_err_rem , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_out_1d (1:nidx), hfx_out         , jpi, jpj, idxice(1:nidx) )
          !
          ! SIMIP diagnostics
-         CALL tab_2d_1d( nbpb, diag_fc_bo_1d   (1:nbpb), diag_fc_bo  , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, diag_fc_su_1d   (1:nbpb), diag_fc_su  , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nidx, diag_fc_bo_1d   (1:nidx), diag_fc_bo  , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, diag_fc_su_1d   (1:nidx), diag_fc_su  , jpi, jpj, idxice(1:nidx) )
          ! ocean surface fields
-         CALL tab_2d_1d( nbpb, sst_1d(1:nbpb), sst_m, jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, sss_1d(1:nbpb), sss_m, jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nidx, sst_1d(1:nidx), sst_m, jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sss_1d(1:nidx), sss_m, jpi, jpj, idxice(1:nidx) )
          !
       CASE( 2 )            ! from 1D to 2D
          !
-         CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, at_i          , idxice, at_i_1d    (1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, ht_i(:,:,jl)  , idxice, ht_i_1d    (1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, ht_s(:,:,jl)  , idxice, ht_s_1d    (1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, a_i (:,:,jl)  , idxice, a_i_1d     (1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, t_su(:,:,jl)  , idxice, t_su_1d    (1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sm_i(:,:,jl)  , idxice, sm_i_1d    (1:nidx)   , jpi, jpj )
          DO jk = 1, nlay_s
-            CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj)
-            CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, e_s_1d     (1:nbpb,jk), jpi, jpj)
+            CALL tab_1d_2d( nidx, t_s(:,:,jk,jl), idxice, t_s_1d     (1:nidx,jk), jpi, jpj)
+            CALL tab_1d_2d( nidx, e_s(:,:,jk,jl), idxice, e_s_1d     (1:nidx,jk), jpi, jpj)
          END DO
          DO jk = 1, nlay_i
-            CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj)
-            CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, e_i_1d     (1:nbpb,jk), jpi, jpj)
-            CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj)
-         END DO
-         CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj )
-         !
-         CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_snw_sum   , npb, wfx_snw_sum_1d(1:nbpb),jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_snw_sub   , npb, wfx_snw_sub_1d(1:nbpb), jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_ice_sub   , npb, wfx_ice_sub_1d(1:nbpb), jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_err_sub   , npb, wfx_err_sub_1d(1:nbpb), jpi, jpj )
-         !
-         CALL tab_1d_2d( nbpb, wfx_bog       , npb, wfx_bog_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_bom       , npb, wfx_bom_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj )
-         !
-         CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sfx_bom       , npb, sfx_bom_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, sfx_sub       , npb, sfx_sub_1d(1:nbpb)   , jpi, jpj )        
-         !
-         CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_err       , npb, hfx_err_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_err_dif   , npb, hfx_err_dif_1d(1:nbpb), jpi, jpj )
-         CALL tab_1d_2d( nbpb, hfx_out       , npb, hfx_out_1d(1:nbpb)   , jpi, jpj )
-         !
-         CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
-         CALL tab_1d_2d( nbpb, ftr_ice(:,:,jl), npb, ftr_ice_1d(1:nbpb) , jpi, jpj )
+            CALL tab_1d_2d( nidx, t_i(:,:,jk,jl), idxice, t_i_1d     (1:nidx,jk), jpi, jpj)
+            CALL tab_1d_2d( nidx, e_i(:,:,jk,jl), idxice, e_i_1d     (1:nidx,jk), jpi, jpj)
+            CALL tab_1d_2d( nidx, s_i(:,:,jk,jl), idxice, s_i_1d     (1:nidx,jk), jpi, jpj)
+         END DO
+         !
+         CALL tab_1d_2d( nidx, wfx_snw       , idxice, wfx_snw_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_snw_sum   , idxice, wfx_snw_sum_1d(1:nidx),jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_sub       , idxice, wfx_sub_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_snw_sub   , idxice, wfx_snw_sub_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_ice_sub   , idxice, wfx_ice_sub_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_err_sub   , idxice, wfx_err_sub_1d(1:nidx), jpi, jpj )
+         !
+         CALL tab_1d_2d( nidx, wfx_bog       , idxice, wfx_bog_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_bom       , idxice, wfx_bom_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_sum       , idxice, wfx_sum_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_sni       , idxice, wfx_sni_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_res       , idxice, wfx_res_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_spr       , idxice, wfx_spr_1d(1:nidx)   , jpi, jpj )
+         !
+         CALL tab_1d_2d( nidx, sfx_bog       , idxice, sfx_bog_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_bom       , idxice, sfx_bom_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_sum       , idxice, sfx_sum_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_sni       , idxice, sfx_sni_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_res       , idxice, sfx_res_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_bri       , idxice, sfx_bri_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_sub       , idxice, sfx_sub_1d(1:nidx)   , jpi, jpj )        
+         !
+         CALL tab_1d_2d( nidx, hfx_thd       , idxice, hfx_thd_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_spr       , idxice, hfx_spr_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_sum       , idxice, hfx_sum_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_bom       , idxice, hfx_bom_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_bog       , idxice, hfx_bog_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_dif       , idxice, hfx_dif_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_opw       , idxice, hfx_opw_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_snw       , idxice, hfx_snw_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_sub       , idxice, hfx_sub_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_err       , idxice, hfx_err_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_res       , idxice, hfx_res_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_err_rem   , idxice, hfx_err_rem_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_err_dif   , idxice, hfx_err_dif_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_out       , idxice, hfx_out_1d(1:nidx)   , jpi, jpj )
+         !
+         CALL tab_1d_2d( nidx, qns_ice(:,:,jl), idxice, qns_ice_1d(1:nidx) , jpi, jpj)
+         CALL tab_1d_2d( nidx, ftr_ice(:,:,jl), idxice, ftr_ice_1d(1:nidx) , jpi, jpj )
          !
          ! SIMIP diagnostics         
-         CALL tab_1d_2d( nbpb, t_si(:,:,jl)   , npb, t_si_1d    (1:nbpb)     , jpi, jpj )
-         CALL tab_1d_2d( nbpb, diag_fc_bo     , npb, diag_fc_bo_1d(1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, diag_fc_su     , npb, diag_fc_su_1d(1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, t_si(:,:,jl)   , idxice, t_si_1d    (1:nidx)     , jpi, jpj )
+         CALL tab_1d_2d( nidx, diag_fc_bo     , idxice, diag_fc_bo_1d(1:nidx)   , jpi, jpj )
+         CALL tab_1d_2d( nidx, diag_fc_su     , idxice, diag_fc_su_1d(1:nidx)   , jpi, jpj )
       END SELECT
       !
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 8326)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90	(revision 8327)
@@ -12,11 +12,14 @@
    !!   lim_thd_da   : sea ice lateral melting
    !!----------------------------------------------------------------------
-   USE par_oce                ! ocean parameters
-   USE phycst                 ! physical constants (ocean directory)
-   USE sbc_oce, ONLY: sst_m   ! Surface boundary condition: ocean fields
-   USE ice                    ! LIM variables
-   USE lib_mpp                ! MPP library
-   USE wrk_nemo               ! work arrays
-   USE lib_fortran            ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   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
+   USE wrk_nemo       ! work arrays
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -97,5 +100,6 @@
       !!              Phil. Trans. R. Soc. A, 373(2052), 20140167.
       !!---------------------------------------------------------------------
-      INTEGER             ::   ji, jj, jl      ! dummy loop indices
+      INTEGER             ::   ji, jj, jk, jl     ! dummy loop indices
+      INTEGER             ::   nidx
       REAL(wp)            ::   zastar, zdfloe, zperi, zwlat, zda
       REAL(wp), PARAMETER ::   zdmax = 300._wp
@@ -104,7 +108,17 @@
       REAL(wp), PARAMETER ::   zm2   = 1.36_wp
       !
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zda_tot
+      REAL(wp), DIMENSION(jpij) ::   zda_tot
       !!---------------------------------------------------------------------
-      CALL wrk_alloc( jpi,jpj, zda_tot )
+
+      ! select ice covered grid points
+      nidx = 0 ; idxice(:) = 0
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            IF ( at_i(ji,jj) > epsi10 ) THEN
+               nidx         = nidx  + 1
+               idxice(nidx) = (jj - 1) * jpi + ji
+            ENDIF
+         END DO
+      END DO
 
       !------------------------------------------------------------!
@@ -112,21 +126,14 @@
       !------------------------------------------------------------!
       zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) )
-      
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            
-            ! Mean floe caliper diameter [m]
-            zdfloe = rn_dmin * ( zastar / ( zastar - at_i(ji,jj) ) )**rn_beta
-            
-            ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2]
-            zperi = at_i(ji,jj) * rpi / ( zcs * zdfloe )
-            
-            ! Melt speed rate [m/s]
-            zwlat = zm1 * ( MAX( 0._wp, sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) )**zm2
-            
-            ! sea ice concentration decrease
-            zda_tot(ji,jj) = - MIN( zwlat * zperi * rdt_ice, at_i(ji,jj) )
-            
-         END DO
+
+      CALL tab_2d_1d( nidx, at_i_1d(1:nidx), at_i , jpi, jpj, idxice(1:nidx) )
+      CALL tab_2d_1d( nidx, t_bo_1d(1:nidx), t_bo , jpi, jpj, idxice(1:nidx) )
+      CALL tab_2d_1d( nidx, sst_1d (1:nidx), sst_m, jpi, jpj, idxice(1:nidx) )
+      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
       
@@ -134,41 +141,49 @@
       ! --- Distribute reduction among ice categories and calculate associated ice-ocean fluxes --- !
       !---------------------------------------------------------------------------------------------!
-      DO jl = jpl, 1, -1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               
-               ! decrease of concentration for the category jl
-               !    1st option: each category contributes to melting in proportion to its concentration
-               rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj) - epsi10 ) )
-               zda     = rswitch * zda_tot(ji,jj) * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi10 )
-               !    2d option: melting of the upper cat first
-               !!zda = MAX( zda_tot(ji,jj), - a_i(ji,jj,jl) )
-               !!zda_tot(ji,jj) = zda_tot(ji,jj) + zda
-               
-               ! Contribution to salt flux
-               sfx_lam(ji,jj) = sfx_lam(ji,jj) - rhoic *  ht_i(ji,jj,jl) * zda * sm_i(ji,jj,jl) * r1_rdtice
-               
-               ! Contribution to heat flux into the ocean [W.m-2], <0  
-!clemX               hfx_thd(ji,jj) = hfx_thd(ji,jj) + zda * r1_rdtice * ( ht_i(ji,jj,jl) * SUM( e_i(ji,jj,:,jl) ) * r1_nlay_i  &
-!                  &                                                + ht_s(ji,jj,jl) *      e_s(ji,jj,1,jl)   * r1_nlay_s )
-               hfx_thd(ji,jj) = hfx_thd(ji,jj) + rswitch * zda_tot(ji,jj) / MAX( at_i(ji,jj), epsi10 )  &
-                  &                                      * r1_rdtice * ( SUM( e_i(ji,jj,:,jl) ) + e_s(ji,jj,1,jl) )
-               
-               ! Contribution to mass flux
-               wfx_lam(ji,jj) =  wfx_lam(ji,jj) - zda * r1_rdtice * ( rhoic * ht_i(ji,jj,jl) + rhosn * ht_s(ji,jj,jl) )
-               
-               ! new concentration
-               a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
-            END DO
-         END DO
+      DO jl = 1, jpl
+         CALL tab_2d_1d( nidx, a_i_1d    (1:nidx), a_i(:,:,jl) , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, ht_i_1d   (1:nidx), ht_i(:,:,jl), jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sm_i_1d   (1:nidx), sm_i(:,:,jl), jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_lam_1d(1:nidx), sfx_lam     , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_thd_1d(1:nidx), hfx_thd     , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_lam_1d(1:nidx), wfx_lam     , jpi, jpj, idxice(1:nidx) )
+         DO jk = 1, nlay_i
+            CALL tab_2d_1d( nidx, e_i_1d(1:nidx,jk), e_i(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
+         END DO
+         DO jk = 1, nlay_s
+            CALL tab_2d_1d( nidx, e_s_1d(1:nidx,jk), e_s(:,:,jk,jl)  , jpi, jpj, idxice(1:nidx) )
+         END DO
+
+         DO ji = 1, nidx
+            ! 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  
+            !clemX               hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda * r1_rdtice * ( ht_i_1d(ji) * SUM( e_i_1d(ji,:) ) * r1_nlay_i  &
+            !                  &                                                     + ht_s_1d(ji) *      e_s_1d(ji,1)   * r1_nlay_s )
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * r1_rdtice * ( SUM( e_i_1d(ji,:) ) + e_s_1d(ji,1) )
+            
+            ! 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) )
+            
+            ! 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 )   ht_i_1d(ji) = 0._wp  
+         END DO
+
+         CALL tab_1d_2d( nidx, a_i (:,:,jl), idxice, a_i_1d     (1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, ht_i(:,:,jl), idxice, ht_i_1d    (1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, sfx_lam     , idxice, sfx_lam_1d(1:nidx) , jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_thd     , idxice, hfx_thd_1d(1:nidx) , jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_lam     , idxice, wfx_lam_1d(1:nidx) , jpi, jpj )
+
       END DO
-      
-      ! total concentration
-      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
-      
-      ! --- ensure that ht_i = 0 where a_i = 0 ---
-      WHERE( a_i == 0._wp ) ht_i = 0._wp
-      !
-      CALL wrk_dealloc( jpi,jpj, zda_tot )
+            
       !
    END SUBROUTINE lim_thd_da
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 8326)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 8327)
@@ -71,10 +71,9 @@
       !!------------------------------------------------------------------------
       INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices
-      INTEGER  ::   nbpac            ! local integers 
-      INTEGER  ::   ii, ij, iter     !   -       -
+      INTEGER  ::   nidx            ! local integers 
+      INTEGER  ::   iter     !   -       -
       REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
-      CHARACTER (len = 15) :: fieldid
 
       REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
@@ -122,18 +121,4 @@
       CALL lim_var_agg(1)
       CALL lim_var_glo2eqv
-      !------------------------------------------------------------------------------|
-      ! 2) Convert units for ice internal energy
-      !------------------------------------------------------------------------------|
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !Energy of melting q(S,T) [J.m-3]
-                  rswitch          = MAX(  0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 )  )   !0 if no ice
-                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp )
-               END DO
-            END DO
-         END DO
-      END DO
 
       !------------------------------------------------------------------------------!
@@ -240,12 +225,10 @@
       !-------------------------------------
       ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
-      nbpac = 0
-      npac(:) = 0
-      !
-      DO jj = 1, jpj
-         DO ji = 1, jpi
+      nidx = 0 ; idxice(:) = 0
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
             IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
-               nbpac = nbpac + 1
-               npac( nbpac ) = (jj - 1) * jpi + ji
+               nidx = nidx + 1
+               idxice( nidx ) = (jj - 1) * jpi + ji
             ENDIF
          END DO
@@ -264,5 +247,5 @@
       ENDIF
    
-      IF( ln_limctl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
+      IF( ln_limctl ) WRITE(numout,*) 'lim_thd_lac : nidx = ', nidx
 
       !------------------------------
@@ -271,27 +254,38 @@
       ! If ocean gains heat do nothing. Otherwise compute new ice formation
 
-      IF ( nbpac > 0 ) THEN
-
-         CALL tab_2d_1d( nbpac, zat_i_1d  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
-         DO jl = 1, jpl
-            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
+      IF ( nidx > 0 ) THEN
+
+         CALL tab_2d_1d( nidx, zat_i_1d  (1:nidx)     , at_i         , jpi, jpj, idxice(1:nidx) )
+         DO jl = 1, jpl
+            CALL tab_2d_1d( nidx, za_i_1d  (1:nidx,jl), a_i  (:,:,jl), jpi, jpj, idxice(1:nidx) )
+            CALL tab_2d_1d( nidx, zv_i_1d  (1:nidx,jl), v_i  (:,:,jl), jpi, jpj, idxice(1:nidx) )
+            CALL tab_2d_1d( nidx, zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl), jpi, jpj, idxice(1:nidx) )
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
-            END DO
-         END DO
-
-         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead     , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo      , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw   , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw   , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol     , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel     , jpi, jpj, npac(1:nbpac) )
-
-         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd   , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw   , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, rn_amax_1d(1:nbpac)     , rn_amax_2d, jpi, jpj, npac(1:nbpac) )
-
+               CALL tab_2d_1d( nidx, ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, idxice(1:nidx) )
+            END DO
+         END DO
+
+         CALL tab_2d_1d( nidx, qlead_1d  (1:nidx)     , qlead     , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, t_bo_1d   (1:nidx)     , t_bo      , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sfx_opw_1d(1:nidx)     , sfx_opw   , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, wfx_opw_1d(1:nidx)     , wfx_opw   , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hicol_1d  (1:nidx)     , hicol     , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, zvrel_1d  (1:nidx)     , zvrel     , jpi, jpj, idxice(1:nidx) )
+
+         CALL tab_2d_1d( nidx, hfx_thd_1d(1:nidx)     , hfx_thd   , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, hfx_opw_1d(1:nidx)     , hfx_opw   , jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, rn_amax_1d(1:nidx)     , rn_amax_2d, jpi, jpj, idxice(1:nidx) )
+         CALL tab_2d_1d( nidx, sss_1d    (1:nidx)     , sss_m     , jpi, jpj, idxice(1:nidx) )
+
+         !------------------------------------------------------------------------------|
+         ! 2) Convert units for ice internal energy
+         !------------------------------------------------------------------------------|
+         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
+            END DO
+         END DO
          !------------------------------------------------------------------------------!
          ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
@@ -301,11 +295,11 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
-         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
+         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:nbpac) = hicol_1d(1:nbpac)
+         zh_newice(1:nidx) = hicol_1d(1:nidx)
 
          !----------------------
@@ -314,13 +308,11 @@
          SELECT CASE ( nn_icesal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(1:nbpac) = rn_icesal
+            zs_newice(1:nidx) = rn_icesal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
-            DO ji = 1, nbpac
-               ii =   MOD( npac(ji) - 1 , jpi ) + 1
-               ij =      ( npac(ji) - 1 ) / jpi + 1
-               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij)  )
+            DO ji = 1, nidx
+               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(1:nbpac) =   2.3
+            zs_newice(1:nidx) =   2.3
          END SELECT
 
@@ -329,5 +321,5 @@
          !-------------------------
          ! We assume that new ice is formed at the seawater freezing point
-         DO ji = 1, nbpac
+         DO ji = 1, nidx
             ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K)
             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                                         &
@@ -339,5 +331,5 @@
          ! Age of new ice
          !----------------
-         DO ji = 1, nbpac
+         DO ji = 1, nidx
             zo_newice(ji) = 0._wp
          END DO
@@ -346,5 +338,5 @@
          ! Volume of new ice
          !-------------------
-         DO ji = 1, nbpac
+         DO ji = 1, nidx
 
             zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
@@ -374,5 +366,5 @@
          IF( ln_frazil ) THEN
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            DO ji = 1, nbpac
+            DO ji = 1, nidx
                rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
                zfrazb        = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb
@@ -385,5 +377,5 @@
          ! Area of new ice
          !-----------------
-         DO ji = 1, nbpac
+         DO ji = 1, nidx
             za_newice(ji) = zv_newice(ji) / zh_newice(ji)
          END DO
@@ -398,5 +390,5 @@
          ! If lateral ice growth gives an ice concentration gt 1, then
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
-         DO ji = 1, nbpac
+         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) )
@@ -413,5 +405,5 @@
          zat_i_1d(:) = 0._wp
          DO jl = 1, jpl
-            DO ji = 1, nbpac
+            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)
@@ -424,5 +416,5 @@
 
          ! Heat content
-         DO ji = 1, nbpac
+         DO ji = 1, nidx
             jl = jcat(ji)                                                    ! categroy in which new ice is put
             zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
@@ -430,5 +422,5 @@
 
          DO jk = 1, nlay_i
-            DO ji = 1, nbpac
+            DO ji = 1, nidx
                jl = jcat(ji)
                rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
@@ -445,8 +437,8 @@
 
             ! for remapping
-            h_i_old (1:nbpac,0:nlay_i+1) = 0._wp
-            eh_i_old(1:nbpac,0:nlay_i+1) = 0._wp
+            h_i_old (1:nidx,0:nlay_i+1) = 0._wp
+            eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
             DO jk = 1, nlay_i
-               DO ji = 1, nbpac
+               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)
@@ -455,5 +447,5 @@
 
             ! new volumes including lateral/bottom accretion + residual
-            DO ji = 1, nbpac
+            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 )
@@ -465,5 +457,5 @@
             ENDDO
             ! --- Ice enthalpy remapping --- !
-            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
+            CALL lim_thd_ent( 1, nidx, ze_i_1d(1:nidx,:,jl) ) 
          ENDDO
 
@@ -472,5 +464,5 @@
          !-----------------
          DO jl = 1, jpl
-            DO ji = 1, nbpac
+            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)
@@ -479,36 +471,30 @@
 
          !------------------------------------------------------------------------------!
+         ! 8) Change units for e_i
+         !------------------------------------------------------------------------------!    
+         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
+            END DO
+         END DO
+         !------------------------------------------------------------------------------!
          ! 7) Change 2D vectors to 1D vectors 
          !------------------------------------------------------------------------------!
          DO jl = 1, jpl
-            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )
+            CALL tab_1d_2d( nidx, a_i (:,:,jl), idxice(1:nidx), za_i_1d (1:nidx,jl), jpi, jpj )
+            CALL tab_1d_2d( nidx, v_i (:,:,jl), idxice(1:nidx), zv_i_1d (1:nidx,jl), jpi, jpj )
+            CALL tab_1d_2d( nidx, smv_i (:,:,jl), idxice(1:nidx), zsmv_i_1d(1:nidx,jl) , jpi, jpj )
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj )
-            END DO
-         END DO
-         CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj )
-         CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
-
-         CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
-         CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj )
+               CALL tab_1d_2d( nidx, e_i(:,:,jk,jl), idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), jpi, jpj )
+            END DO
+         END DO
+         CALL tab_1d_2d( nidx, sfx_opw, idxice(1:nidx), sfx_opw_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, wfx_opw, idxice(1:nidx), wfx_opw_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_thd, idxice(1:nidx), hfx_thd_1d(1:nidx), jpi, jpj )
+         CALL tab_1d_2d( nidx, hfx_opw, idxice(1:nidx), hfx_opw_1d(1:nidx), jpi, jpj )
          !
-      ENDIF ! nbpac > 0
-
-      !------------------------------------------------------------------------------!
-      ! 8) Change units for e_i
-      !------------------------------------------------------------------------------!    
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  ! heat content in J/m2
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
-               END DO
-            END DO
-         END DO
-      END DO
-
+      ENDIF ! nidx > 0
       !
       CALL wrk_dealloc( jpij, jcat )   ! integer
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8326)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8327)
@@ -476,6 +476,6 @@
          DO jk = 1, nlay_i
             DO ji = kideb, kiut
-               ii =  MOD( npb(ji) - 1 , jpi ) + 1
-               ij =     ( npb(ji) - 1 ) / jpi + 1
+               ii =  MOD( idxice(ji) - 1 , jpi ) + 1
+               ij =     ( idxice(ji) - 1 ) / jpi + 1
                ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
                zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8326)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8327)
@@ -26,7 +26,5 @@
    !: are the variables corresponding to 2d vectors
 
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npb    !: address vector for 1d vertical thermo computations
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   nplm   !: address vector for mono-category lateral melting
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: address vector for new ice formation
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   idxice !: selected points for ice thermo
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlead_1d     
@@ -64,4 +62,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_sub_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_err_sub_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_lam_1d 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bog_1d    
@@ -80,6 +79,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_opw_1d   
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d  
-
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sub_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_lam_1d
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
@@ -146,5 +145,5 @@
 
       ii = 1
-      ALLOCATE( npb      (jpij) , nplm      (jpij) , npac       (jpij) ,   &
+      ALLOCATE( idxice   (jpij) ,   &
          &      qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   &
          &      fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij)  ,   &
@@ -163,9 +162,9 @@
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
          &      wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,              &
-         &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         &
+         &      wfx_lam_1d(jpij)  , dqns_ice_1d(jpij) , evap_ice_1d (jpij),                     &
          &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &  
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
-         &      hicol_1d   (jpij)                     , STAT=ierr(ii) )
+         &      sfx_lam_1d (jpij) , hicol_1d   (jpij)                     , STAT=ierr(ii) )
       !
       ii = ii + 1
