Index: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4853)
+++ branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4854)
@@ -89,4 +89,5 @@
       CALL lim_itd_ini                 ! ice thickness distribution initialization
       !
+      CALL lim_itd_me_init             ! ice thickness distribution initialization
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest
Index: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4853)
+++ branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4854)
@@ -150,6 +150,4 @@
       CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
 
-      IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only)
-
       IF(ln_ctl) THEN
          CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i  : ', tab2d_2=at_i   , clinfo2=' at_i    : ')
@@ -1037,5 +1035,5 @@
             !     / rafting category n1.
             !--------------------------------------------------------------------------
-            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
+            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
             vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por )
             vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por
@@ -1043,5 +1041,5 @@
             vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
             esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
-            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )
+            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
             srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless
 
@@ -1128,5 +1126,5 @@
                jj = indxj(ij)
                ! heat content of ridged ice
-               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 
+               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 
                eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
                e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
Index: branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4853)
+++ branches/2014/dev_4728_CNRS04_coupled_interface/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4854)
@@ -29,4 +29,5 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE limthd_ent
@@ -133,5 +134,5 @@
                   !Energy of melting q(S,T) [J.m-3]
                   zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
-                  e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i )
+                  e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp )
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac
                END DO
@@ -171,7 +172,6 @@
          zgamafr = 0.03
 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
+         DO jj = 2, jpj
+            DO ji = 2, jpi
                IF ( qlead(ji,jj) < 0._wp ) THEN
                   !-------------
@@ -243,4 +243,7 @@
             END DO ! loop on ji ends
          END DO ! loop on jj ends
+      ! 
+      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
+      CALL lbc_lnk( hicol(:,:), 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -255,4 +258,6 @@
       ! This occurs if open water energy budget is negative
       nbpac = 0
+      npac(:) = 0
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -315,7 +320,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_old(:,:) = zv_i_1d(:,:) 
-         za_old(:,:) = za_i_1d(:,:)
-
+         zv_old(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
+         za_old(1:nbpac,:) = za_i_1d(1:nbpac,:)
          !----------------------
          ! Thickness of new ice
@@ -324,5 +328,5 @@
             zh_newice(ji) = hiccrit
          END DO
-         IF( fraz_swi == 1 ) zh_newice(:) = hicol_b(:)
+         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_b(1:nbpac)
 
          !----------------------
@@ -331,5 +335,5 @@
          SELECT CASE ( num_sal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(:) = bulk_sal
+            zs_newice(1:nbpac) = bulk_sal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
             DO ji = 1, nbpac
@@ -339,5 +343,5 @@
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(:) =   2.3
+            zs_newice(1:nbpac) =   2.3
          END SELECT
 
@@ -472,5 +476,4 @@
                za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)               
                zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
-
                ! for remapping
                h_i_old (ji,nlay_i+1) = zv_newfra
@@ -479,7 +482,5 @@
 
             ! --- Ice enthalpy remapping --- !
-            IF( zv_newfra > 0._wp ) THEN
-               CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
-            ENDIF
+            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
 
          ENDDO
@@ -534,5 +535,5 @@
                DO ji = 1, jpi
                   ! heat content in Joules
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac ) 
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 
                END DO
             END DO
