Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4901)
@@ -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_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4901)
@@ -29,5 +29,4 @@
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE wrk_nemo         ! work arrays
-   USE cpl_oasis3, ONLY : lk_cpl
 
    IMPLICIT NONE
@@ -113,8 +112,4 @@
 
       CALL lim_istate_init     !  reading the initials parameters of the ice
-
-# if defined key_coupled
-      albege(:,:)   = 0.8 * tms(:,:)
-# endif
 
       ! surface temperature
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4901)
@@ -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_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4901)
@@ -32,6 +32,5 @@
    USE sbc_oce          ! Surface boundary condition: ocean fields
    USE sbccpl
-   USE cpl_oasis3, ONLY : lk_cpl
-   USE oce       , ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
+   USE oce       , ONLY : fraqsr_1lev, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass
    USE albedo           ! albedo parameters
    USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
@@ -98,19 +97,21 @@
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
-      !!              - alb_ice : sea-ice alberdo (lk_cpl=T)
+      !!              - alb_ice : sea-ice albedo (lk_cpl=T)
       !!
       !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
+      !!              These refs are now obsolete since everything has been revised
+      !!              The ref should be Rousset et al., 2015?
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
-      !
-      INTEGER  ::   ji, jj, jl, jk           ! dummy loop indices
-      REAL(wp) ::   zinda, zemp      ! local scalars
-      REAL(wp) ::   zf_mass         ! Heat flux associated with mass exchange ice->ocean (W.m-2)
-      REAL(wp) ::   zfcm1           ! New solar flux received by the ocean
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
+      INTEGER, INTENT(in) ::   kt                                   ! number of iteration
+      !
+      INTEGER  ::   ji, jj, jl, jk                                  ! dummy loop indices
+      !
+      REAL(wp) ::   zinda, zemp                                     !  local scalars
+      REAL(wp) ::   zf_mass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2)
+      REAL(wp) ::   zfcm1                                           ! New solar flux received by the ocean
+      !
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace
       !!---------------------------------------------------------------------
-      
-      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp )
 
       ! make calls for heat fluxes before it is modified
@@ -134,14 +135,12 @@
             ! Solar heat flux reaching the ocean = zfcm1 (W.m-2) 
             !---------------------------------------------------
-            IF( lk_cpl ) THEN ! be carfeful: not been tested yet
-               ! original line
+            IF( lk_cpl ) THEN 
+               !!! LIM2 version zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
                zfcm1 = qsr_tot(ji,jj)
-               !!!zfcm1 = qsr_tot(ji,jj) + ftr_ice(ji,jj) * ( 1._wp - pfrld(ji,jj) ) / ( 1._wp - zinda + zinda * iatte(ji,jj) )
                DO jl = 1, jpl
-                  zfcm1 = zfcm1 - ( qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)
+                  zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)
                END DO
             ELSE
-               !!!zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + &
-               !!!     &    ( 1._wp - pfrld(ji,jj) ) * ftr_ice(ji,jj) / ( 1._wp - zinda + zinda * iatte(ji,jj) )
+               !!! LIM2 version zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
                zfcm1   = pfrld(ji,jj) * qsr(ji,jj)
                DO jl = 1, jpl
@@ -215,9 +214,16 @@
 
       !------------------------------------------------!
-      !    Computation of snow/ice and ocean albedo    !
+      !    Snow/ice albedo (only if sent to coupler)   !
       !------------------------------------------------!
       IF( lk_cpl ) THEN          ! coupled case
-         CALL albedo_ice( t_su, ht_i, ht_s, zalbp, zalb )                  ! snow/ice albedo
-         alb_ice(:,:,:) =  0.5_wp * zalbp(:,:,:) + 0.5_wp * zalb (:,:,:)   ! Ice albedo (mean clear and overcast skys)
+
+            CALL wrk_alloc( jpi, jpj, jpl, zalb_cs, zalb_os )
+
+            CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos
+
+            alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
+
+            CALL wrk_dealloc( jpi, jpj, jpl, zalb_cs, zalb_os )
+
       ENDIF
 
@@ -229,7 +235,5 @@
          CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
       ENDIF
-      !
-      IF( lk_cpl )   CALL wrk_dealloc( jpi, jpj, jpl, zalb, zalbp )
-      ! 
+
    END SUBROUTINE lim_sbc_flx
 
@@ -344,6 +348,5 @@
       ! clem modif
       IF( .NOT. ln_rstart ) THEN
-         iatte(:,:) = 1._wp
-         oatte(:,:) = 1._wp
+         fraqsr_1lev(:,:) = 1._wp
       ENDIF
       !
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4901)
@@ -22,5 +22,5 @@
    USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain variables
-   USE oce     , ONLY :  iatte, oatte
+   USE oce     , ONLY : fraqsr_1lev
    USE ice            ! LIM: sea-ice variables
    USE par_ice        ! LIM: sea-ice parameters
@@ -43,5 +43,4 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE timing         ! Timing
-   USE cpl_oasis3, ONLY : lk_cpl
    USE limcons        ! conservation tests
 
@@ -68,5 +67,5 @@
       !!                ***  ROUTINE lim_thd  ***       
       !!  
-      !! ** Purpose : This routine manages the ice thermodynamic.
+      !! ** Purpose : This routine manages ice thermodynamics
       !!         
       !! ** Action : - Initialisation of some variables
@@ -74,11 +73,12 @@
       !!               at the ice base, snow acc.,heat budget of the leads)
       !!             - selection of the icy points and put them in an array
-      !!             - call lim_vert_ther for vert ice thermodynamic
-      !!             - back to the geographic grid
-      !!             - selection of points for lateral accretion
-      !!             - call lim_lat_acc  for the ice accretion
+      !!             - call lim_thd_dif  for vertical heat diffusion
+      !!             - call lim_thd_dh   for vertical ice growth and melt
+      !!             - call lim_thd_ent  for enthalpy remapping
+      !!             - call lim_thd_sal  for ice desalination
+      !!             - call lim_thd_temp to  retrieve temperature from ice enthalpy
       !!             - back to the geographic grid
       !!     
-      !! ** References : H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
+      !! ** References : 
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt    ! number of iteration
@@ -93,5 +93,9 @@
       !
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      !
+      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns
       !!-------------------------------------------------------------------
+      CALL wrk_alloc( jpi, jpj, zqsr, zqns )
+
       IF( nn_timing == 1 )  CALL timing_start('limthd')
 
@@ -137,4 +141,24 @@
       !-----------------------------------------------------------------------------!
 
+      !--- Ocean solar and non solar fluxes to be used in zqld
+      IF ( .NOT. lk_cpl ) THEN   ! --- forced case, fluxes to the lead are the same as over the ocean
+         !
+         zqsr(:,:) = qsr(:,:)      ; zqns(:,:) = qns(:,:)
+         !
+      ELSE                       ! --- coupled case, fluxes to the lead are total - intercepted
+         !
+         zqsr(:,:) = qsr_tot(:,:)  ; zqns(:,:) = qns_tot(:,:)
+         !
+         DO jl = 1, jpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * old_a_i(ji,jj,jl)
+                  zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * old_a_i(ji,jj,jl)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+
 !CDIR NOVERRCHK
       DO jj = 1, jpj
@@ -149,11 +173,13 @@
             !           !  temperature and turbulent mixing (McPhee, 1992)
             !
+
             ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
-            zqld =  tms(ji,jj) * rdt_ice *                                       &
-               &  ( pfrld(ji,jj)         * ( qsr(ji,jj) * oatte(ji,jj)           &   ! solar heat + clem modif
-               &                           + qns(ji,jj) )                        &   ! non solar heat
-               ! latent heat of precip (note that precip is included in qns but not in qns_ice)
-               &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  &
+            zqld =  tms(ji,jj) * rdt_ice *                                          &
+               &      ( pfrld(ji,jj) * ( zqsr(ji,jj) * fraqsr_1lev(ji,jj) + zqns(ji,jj) ) &
+               &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) *     & ! heat content of precip
+               &      ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  &
                &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) )
+               ! REMARK valid at least in forced mode from clem
+               ! precip is included in qns but not in qns_ice
 
             !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
@@ -185,5 +211,5 @@
             hfx_in(ji,jj) = hfx_in(ji,jj)                                                                                         & 
                ! heat flux above the ocean
-               &    +             pfrld(ji,jj)   * ( qns(ji,jj) + qsr(ji,jj) )                                                    &
+               &    +             pfrld(ji,jj)   * ( zqns(ji,jj) + zqsr(ji,jj) )                                                  &
                ! latent heat of precip (note that precip is included in qns but not in qns_ice)
                &    +   ( 1._wp - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  &
@@ -306,7 +332,4 @@
             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, iatte_1d   (1:nbpb), iatte           , jpi, jpj, npb(1:nbpb) ) 
-            CALL tab_2d_1d( nbpb, oatte_1d   (1:nbpb), oatte           , jpi, jpj, npb(1:nbpb) ) 
 
             CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
@@ -482,8 +505,13 @@
       ENDIF
       !
+      !
+      CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
+
+      !
       ! conservation test
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
       IF( nn_timing == 1 )  CALL timing_stop('limthd')
+
    END SUBROUTINE lim_thd 
 
@@ -552,4 +580,6 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
       IF(lwm) WRITE ( numoni, namicethd )
+
+      IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
       !
       IF(lwp) THEN                          ! control print
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4901)
@@ -26,5 +26,4 @@
    USE wrk_nemo       ! work arrays
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE cpl_oasis3, ONLY : lk_cpl
    
    IMPLICIT NONE
@@ -166,9 +165,9 @@
       !
       DO ji = kideb, kiut
-         zinda         = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
-         ztmelts       = zinda * rtt + ( 1._wp - zinda ) * rtt
-
-         zfdum     = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
-         zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
+         zinda      = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
+         ztmelts    = zinda * rtt + ( 1._wp - zinda ) * rtt
+
+         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
 
          zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4901)
@@ -25,5 +25,5 @@
    USE wrk_nemo       ! work arrays
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE cpl_oasis3, ONLY : lk_cpl
+   USE sbc_oce, ONLY : lk_cpl
 
    IMPLICIT NONE
@@ -146,6 +146,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms
       ! diag errors on heat
-      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini
-      REAL(wp)                        :: zhfx_err
+      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err
       !!------------------------------------------------------------------     
       ! 
@@ -158,5 +157,5 @@
       CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
 
-      CALL wrk_alloc( jpij, zdq, zq_ini )
+      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
 
       ! --- diag error on heat diffusion - PART 1 --- !
@@ -272,5 +271,4 @@
 
       DO ji = kideb, kiut           ! Radiation transmitted below the ice
-         !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif
          ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
       END DO
@@ -407,11 +405,16 @@
          !------------------------------------------------------------------------------|
          !
-         DO ji = kideb , kiut
-            ! update of the non solar flux according to the update in T_su
-            qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
-
+         IF( .NOT. lk_cpl ) THEN   !--- forced atmosphere case
+            DO ji = kideb , kiut
+               ! update of the non solar flux according to the update in T_su
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) )
+            END DO
+         ENDIF
+
+         ! Update incoming flux
+         DO ji = kideb , kiut
             ! update incoming flux
             zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation
-               + qns_ice_1d(ji)                  ! non solar total flux 
+               + qns_ice_1d(ji)                   ! non solar total flux 
             ! (LWup, LWdw, SH, LH)
          END DO
@@ -737,14 +740,30 @@
       CALL lim_thd_enmelt( kideb, kiut )
 
-      ! --- diag error on heat diffusion - PART 2 --- !
+      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
       DO ji = kideb, kiut
          zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  &
             &                              SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
-         zhfx_err    = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
-         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji)
-         ! --- correction of qns_ice and surface conduction flux --- !
-         qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 
-         fc_su     (ji) = fc_su     (ji) - zhfx_err 
-         ! --- Heat flux at the ice surface in W.m-2 --- !
+         zhfx_err(ji)   = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
+         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_b(ji)
+      END DO 
+
+      ! diagnose external surface (forced case) or bottom (forced case) from heat conservation
+      IF( .NOT. lk_cpl ) THEN   ! --- forced case: qns_ice and fc_su are diagnosed
+         !
+         DO ji = kideb, kiut
+            qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji)
+            fc_su     (ji) = fc_su(ji)      - zhfx_err(ji)
+         END DO
+         !
+      ELSE                      ! --- coupled case: ocean turbulent heat flux is diagnosed
+         !
+         DO ji = kideb, kiut
+            fhtur_1d  (ji) = fhtur_1d(ji)   - zhfx_err(ji)
+         END DO
+         !
+      ENDIF
+
+      ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2)
+      DO ji = kideb, kiut
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
          hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
@@ -759,5 +778,5 @@
       CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
       CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
-      CALL wrk_dealloc( jpij, zdq, zq_ini )
+      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
 
    END SUBROUTINE lim_thd_dif
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4901)
@@ -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
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4900)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4901)
@@ -115,7 +115,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new     !: Salinity of new ice at the bottom
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   iatte_1d   !: clem attenuation coef of the input solar flux (unitless)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oatte_1d   !: clem attenuation coef of the input solar flux (unitless)
-
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_b   !: corresponding to the 2D var  t_s
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_b   !: corresponding to the 2D var  t_i
@@ -149,5 +146,5 @@
          &      qsr_ice_1d (jpij) ,     &
          &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) ,     &
-         &      t_bo_b   (jpij) , iatte_1d  (jpij) , oatte_1d (jpij) ,     &
+         &      t_bo_b   (jpij) ,                                          &
          &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , &
