Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4672)
@@ -167,5 +167,4 @@
    !                                          !!** ice-dynamic namelist (namicedyn) **
    INTEGER , PUBLIC ::   nevp   = 400          !: number of iterations for subcycling
-   INTEGER , PUBLIC ::   nlay_i = 5            !: number of layers in the ice
 
    !                                          !!** ice-dynamic namelist (namicedyn) **
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4672)
@@ -91,15 +91,8 @@
       REAL(wp) :: zinda, zindb, zareamin 
       REAL(wp) :: zfric_u, zqld, zqfr
-      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx, zqfx
-      REAL(wp)                        :: zhfx_err, ztest
       !
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limthd')
-
-      CALL wrk_alloc( jpij, zdq, zq_ini, zhfx, zqfx )
-   
-      ! init debug
-      zdq(:) = 0._wp ; zq_ini(:) = 0._wp ; zhfx(:) = 0._wp ; zqfx(:) = 0._wp      
 
       ! conservation test
@@ -333,54 +326,18 @@
             ! 4.3) Thermodynamic processes
             !--------------------------------
-            ! --- diag error on heat diffusion - PART 1 --- !
-            DO ji = 1, nbpb
-               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 ) ) 
-            END DO
 
             !---------------------------------!
             ! Ice/Snow Temperature profile    !
             !---------------------------------!
-            CALL lim_thd_dif( 1, nbpb, jl )
-
-            ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
-            CALL lim_thd_enmelt( 1, nbpb )
-
-            DO ji = 1, nbpb
-              ! --- diag error on heat diffusion - PART 2 --- !
-               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) - ftr_ice_1d(ji) - 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 --- !
-               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) )
-            END DO
+            CALL lim_thd_dif( 1, nbpb )
 
             !---------------------------------!
             ! Ice/Snow thicnkess              !
             !---------------------------------!
-            ! --- diag error on heat remapping - PART 1 --- !
-            DO ji = 1, nbpb
-               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 ) ) 
-            END DO
-
-            CALL lim_thd_dh( 1, nbpb, jl )    
+            CALL lim_thd_dh( 1, nbpb )    
 
             ! --- Ice enthalpy remapping --- !
-            CALL lim_thd_ent( 1, nbpb, jl, q_i_b(1:nbpb,:) ) 
-            !                                
-            ! --- diag error on heat remapping - PART 2 --- !
-            DO ji = 1, nbpb
-               zdq(ji)        = - ( zq_ini(ji) + dq_i(ji) + dq_s(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 ) )
-               hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + zdq(ji) * a_i_b(ji) * r1_rdtice
-            END DO
-
+            CALL lim_thd_ent( 1, nbpb, q_i_b(1:nbpb,:) ) 
+                                            
             !---------------------------------!
             ! --- Ice salinity --- !
@@ -528,40 +485,6 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
-      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx, zqfx )
-
       IF( nn_timing == 1 )  CALL timing_stop('limthd')
-   END SUBROUTINE lim_thd
-
- 
-   SUBROUTINE lim_thd_enmelt( kideb, kiut )
-      !!-----------------------------------------------------------------------
-      !!                   ***  ROUTINE lim_thd_enmelt *** 
-      !!                 
-      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
-      !!
-      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
-      !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
-      !!
-      INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts, zindb  ! local scalar 
-      !!-------------------------------------------------------------------
-      !
-      DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = kideb, kiut
-            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt 
-            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
-            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             &
-               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   &
-               &                   - rcp  *                 ( ztmelts-rtt )  ) 
-         END DO
-      END DO
-      DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = kideb, kiut
-            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
-         END DO
-      END DO
-      !
-   END SUBROUTINE lim_thd_enmelt
+   END SUBROUTINE lim_thd 
 
    SUBROUTINE lim_thd_temp( kideb, kiut )
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4672)
@@ -43,5 +43,5 @@
 CONTAINS
 
-   SUBROUTINE lim_thd_dh( kideb, kiut, jl )
+   SUBROUTINE lim_thd_dh( kideb, kiut )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_thd_dh  ***
@@ -68,5 +68,4 @@
       !!------------------------------------------------------------------
       INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
-      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
       !! 
       INTEGER  ::   ji , jk        ! dummy loop indices
@@ -151,8 +150,4 @@
       zintermelt(:)   = 0._wp
       icount    (:)   = 0
-
-      ! debug
-      dq_i(:) = 0._wp
-      dq_s(:) = 0._wp
 
       ! initialize layer thicknesses and enthalpies
@@ -273,6 +268,4 @@
          zh_s  (ji) = ht_s_b(ji) / REAL( nlay_s )
 
-         ! clem debug: variation of enthalpy (J.m-2)
-         dq_s(ji) = dq_s(ji) + ( zdh_s_pre(ji) + zdh_s_mel(ji) ) * zqprec(ji)  
          ENDIF
       END DO
@@ -297,6 +290,4 @@
             ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )
 
-            ! clem debug: variation of enthalpy (J.m-2)
-            dq_s(ji) = dq_s(ji) + zdeltah(ji,jk) * q_s_b(ji,jk)  
          END DO
       END DO
@@ -325,6 +316,4 @@
             ! new snow thickness
             ht_s_b(ji)     =  MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
-            ! clem debug: variation of enthalpy (J.m-2)
-            dq_s(ji) = dq_s(ji) + zdh_s_sub(ji) * q_s_b(ji,1)  
          END DO
       ENDIF
@@ -397,7 +386,4 @@
             icount(ji)  = icount(ji) + zindh
             zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
-
-            ! clem debug: variation of enthalpy (J.m-2)
-            dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)  
 
             ! update heat content (J.m-2) and layer thickness
@@ -513,7 +499,4 @@
             wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
 
-            ! clem debug: variation of enthalpy (J.m-2)
-            dq_i(ji) = dq_i(ji) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)  
-
             ! update heat content (J.m-2) and layer thickness
             qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
@@ -557,7 +540,4 @@
                   ! Contribution to mass flux
                   wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
-
-                  ! clem debug: variation of enthalpy (J.m-2)
-                  dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)  
 
                   ! update heat content (J.m-2) and layer thickness
@@ -598,7 +578,4 @@
                   ! Contribution to mass flux
                   wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
-
-                  ! clem debug: variation of enthalpy (J.m-2)
-                  dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)  
 
                   ! update heat content (J.m-2) and layer thickness
@@ -664,6 +641,4 @@
 !         ! Contribution to mass flux
 !         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
-!         ! clem debug: variation of enthalpy (J.m-2)
-!         dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)  
 !    
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
@@ -721,8 +696,4 @@
          wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
 
-         ! clem debug: variation of enthalpy (J.m-2)
-         dq_s(ji) = dq_s(ji) - dh_snowice(ji) * q_s_b(ji,1)  
-         dq_i(ji) = dq_i(ji) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw  
-
          ! update heat content (J.m-2) and layer thickness
          qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4672)
@@ -40,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE lim_thd_dif( kideb , kiut , jl )
+   SUBROUTINE lim_thd_dif( kideb , kiut )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_thd_dif  ***
@@ -93,5 +93,4 @@
       !!------------------------------------------------------------------
       INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
-      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
 
       !! * Local variables
@@ -146,4 +145,7 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms
+      ! diag errors on heat
+      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini
+      REAL(wp)                        :: zhfx_err
       !!------------------------------------------------------------------     
       ! 
@@ -155,4 +157,13 @@
       CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  )
       CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
+
+      CALL wrk_alloc( jpij, zdq, zq_ini )
+
+      ! --- diag error on heat diffusion - PART 1 --- !
+      zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
+      DO ji = kideb, kiut
+         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 ) ) 
+      END DO
 
       !------------------------------------------------------------------------------!
@@ -671,6 +682,4 @@
          DO layer  =  1, nlay_s
             DO ji = kideb , kiut
-               ii = MOD( npb(ji) - 1, jpi ) + 1
-               ij = ( npb(ji) - 1 ) / jpi + 1
                t_s_b(ji,layer) = MAX(  MIN( t_s_b(ji,layer), rtt ), 190._wp  )
                zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
@@ -722,9 +731,26 @@
          IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC
             hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
-         ELSE                                    ! case T_su = 0degC
+         ELSE                         ! case T_su = 0degC
             hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
          ENDIF
       END DO
 
+      ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
+      CALL lim_thd_enmelt( kideb, kiut )
+
+      ! --- diag error 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 --- !
+         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) )
+      END DO
+   
       !
       CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
@@ -735,6 +761,38 @@
       CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
       CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
+      CALL wrk_dealloc( jpij, zdq, zq_ini )
 
    END SUBROUTINE lim_thd_dif
+
+   SUBROUTINE lim_thd_enmelt( kideb, kiut )
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE lim_thd_enmelt *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
+      !
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts, zindb  ! local scalar 
+      !!-------------------------------------------------------------------
+      !
+      DO jk = 1, nlay_i             ! Sea ice energy of melting
+         DO ji = kideb, kiut
+            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt 
+            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
+            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             &
+               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   &
+               &                   - rcp  *                 ( ztmelts-rtt )  ) 
+         END DO
+      END DO
+      DO jk = 1, nlay_s             ! Snow energy of melting
+         DO ji = kideb, kiut
+            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
+         END DO
+      END DO
+      !
+   END SUBROUTINE lim_thd_enmelt
 
 #else
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4672)
@@ -36,5 +36,5 @@
    PRIVATE
 
-   PUBLIC   lim_thd_ent         ! called by lim_thd
+   PUBLIC   lim_thd_ent         ! called by limthd and limthd_lac
 
    REAL(wp) :: epsi20 = 1.e-20   ! constant values
@@ -48,5 +48,5 @@
 CONTAINS
  
-   SUBROUTINE lim_thd_ent( kideb, kiut, jl, qnew )
+   SUBROUTINE lim_thd_ent( kideb, kiut, qnew )
       !!-------------------------------------------------------------------
       !!               ***   ROUTINE lim_thd_ent  ***
@@ -74,12 +74,10 @@
       !!-------------------------------------------------------------------
       INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
-      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number
 
-      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (remapped)
+      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (J.m-3, remapped)
 
-      INTEGER  :: ji,ii,ij   !  dummy loop indices
+      INTEGER  :: ji         !  dummy loop indices
       INTEGER  :: jk0, jk1   !  old/new layer indices
-      REAL(wp) :: ztmelts    ! temperature of melting
-      REAL(wp) :: zswitch, zaaa, zbbb, zccc, zdiscrim ! converting enthalpy to temperature
+      REAL(wp) :: zswitch
       !
       REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces
@@ -139,8 +137,17 @@
       DO jk1 = 1, nlay_i
          DO ji = kideb, kiut
-            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) ) 
+            zswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) ) 
             qnew(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 )
          ENDDO
       ENDDO
+
+      ! --- diag error on heat remapping --- !
+      ! comment: if input h_i_old and qh_i_old are already multiplied by a_i (as in limthd_lac), 
+      ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
+      DO ji = kideb, kiut
+         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_b(ji) * r1_rdtice *  &
+            &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) ) 
+      END DO
+      
       !
       CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 )
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4672)
@@ -329,5 +329,4 @@
          ! Salinity of new ice 
          !----------------------
-
          SELECT CASE ( num_sal )
          CASE ( 1 )                    ! Sice = constant 
@@ -343,5 +342,4 @@
          END SELECT
 
-
          !-------------------------
          ! Heat content of new ice
@@ -354,4 +352,5 @@
                &                       - rcp  *         ( ztmelts - rtt )  )
          END DO ! ji
+
          !----------------
          ! Age of new ice
@@ -395,5 +394,4 @@
          END DO
 
-
          !-----------------
          ! Area of new ice
@@ -407,7 +405,7 @@
          !------------------------------------------------------------------------------!
 
-         !-------------------------------------------
-         ! Compute excessive new ice area and volume
-         !-------------------------------------------
+         !------------------------
+         ! 6.1) lateral ice growth
+         !------------------------
          ! 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
@@ -422,9 +420,7 @@
                zdv_res(ji) = 0._wp
             ENDIF
-         END DO ! ji
-
-         !------------------------------------------------
-         ! Laterally redistribute new ice volume and area
-         !------------------------------------------------
+         END DO
+
+         ! find which category to fill
          zat_i_1d(:) = 0._wp
          DO jl = 1, jpl
@@ -433,5 +429,5 @@
                   za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
                   zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
-                  jcat  (ji)    = jl
+                  jcat    (ji)    = jl
                ENDIF
                zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
@@ -439,9 +435,7 @@
          END DO
 
-         !----------------------------------
-         ! Heat content - lateral accretion
-         !----------------------------------
-         DO ji = 1, nbpac
-            jl = jcat(ji)                                                  ! categroy in which new ice is put
+         ! Heat content
+         DO ji = 1, nbpac
+            jl = jcat(ji)                                                    ! categroy in which new ice is put
             zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice
          END DO
@@ -457,11 +451,12 @@
          END DO
 
-         !-----------------------------------------------
-         ! Add excessive volume of new ice at the bottom
-         !-----------------------------------------------
+         !------------------------------------------------
+         ! 6.2) bottom ice growth + ice enthalpy remapping
+         !------------------------------------------------
          DO jl = 1, jpl
+
+            ! for remapping
             h_i_old (1:nbpac,0:nlay_i+1) = 0._wp
             qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp
-
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
@@ -471,17 +466,20 @@
             END DO
 
+            ! new volumes including lateral/bottom accretion + residual
             DO ji = 1, nbpac
                zinda          = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
                zv_newfra      = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
                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
                qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
-
-               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
             ENDDO
 
             ! --- Ice enthalpy remapping --- !
-            CALL lim_thd_ent( 1, nbpac, jl, ze_i_1d(1:nbpac,:,jl) ) 
+            IF( zv_newfra > 0._wp ) THEN
+               CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
+            ENDIF
 
          ENDDO
@@ -500,15 +498,13 @@
          ! Update salinity
          !-----------------
-         !clem IF(  num_sal == 2  ) THEN
-            DO jl = 1, jpl
-               DO ji = 1, nbpac
-                  zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl)
-                  zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
-               END DO
-            END DO   
-         !clem ENDIF
+         DO jl = 1, jpl
+            DO ji = 1, nbpac
+               zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl)
+               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
+            END DO
+         END DO
 
          !------------------------------------------------------------------------------!
-         ! 8) Change 2D vectors to 1D vectors 
+         ! 7) Change 2D vectors to 1D vectors 
          !------------------------------------------------------------------------------!
          DO jl = 1, jpl
@@ -531,5 +527,5 @@
 
       !------------------------------------------------------------------------------!
-      ! 9) Change units for e_i
+      ! 8) Change units for e_i
       !------------------------------------------------------------------------------!    
       DO jl = 1, jpl
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 4672)
@@ -13,4 +13,5 @@
    !                                             !!! ice thermodynamics
    INTEGER, PUBLIC, PARAMETER ::   jkmax    = 6   !: maximum number of ice layers
+   INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 5   !: number of ice layers
    INTEGER, PUBLIC, PARAMETER ::   nlay_s   = 1   !: number of snow layers
 
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4659)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4672)
@@ -124,7 +124,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_s_b   !:    Snow enthalpy per unit volume
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dq_i   !: variation of ice enthalpy (debug)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dq_s   !: variation of snw enthalpy (debug)
-
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qh_i_old  !: ice heat content (q*h, J.m-2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old   !: ice thickness layer (m)
@@ -171,5 +168,5 @@
          &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &    
          &      dh_snowice(jpij) , sm_i_b   (jpij) , s_i_new  (jpij) ,    &
-         &      dq_i      (jpij) , dq_s     (jpij),  t_s_b(jpij,nlay_s),  &
+         &      t_s_b(jpij,nlay_s),                                       &
          &      t_i_b(jpij,jkmax), s_i_b(jpij,jkmax)                ,     &            
          &      q_i_b(jpij,jkmax), q_s_b(jpij,jkmax)                ,     &
