Index: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4657)
+++ branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4924)
@@ -25,4 +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
 
    IMPLICIT NONE
@@ -31,5 +32,5 @@
    PUBLIC   lim_thd_dif   ! called by lim_thd
 
-   REAL(wp) ::   epsi10      =  1.e-10_wp    !
+   REAL(wp) ::   epsi10 = 1.e-10_wp    !
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
@@ -39,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE lim_thd_dif( kideb , kiut , jl )
+   SUBROUTINE lim_thd_dif( kideb , kiut )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_thd_dif  ***
@@ -74,10 +75,10 @@
       !!
       !! ** Inputs / Ouputs : (global commons)
-      !!           surface temperature : t_su_b
-      !!           ice/snow temperatures   : t_i_b, t_s_b
-      !!           ice salinities          : s_i_b
+      !!           surface temperature : t_su_1d
+      !!           ice/snow temperatures   : t_i_1d, t_s_1d
+      !!           ice salinities          : s_i_1d
       !!           number of layers in the ice/snow: nlay_i, nlay_s
       !!           profile of the ice/snow layers : z_i, z_s
-      !!           total ice/snow thickness : ht_i_b, ht_s_b
+      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
       !!
       !! ** External : 
@@ -91,7 +92,5 @@
       !!           (04-2007) Energy conservation tested by M. Vancoppenolle
       !!------------------------------------------------------------------
-      INTEGER , INTENT (in) ::   kideb   ! Start point on which the  the computation is applied
-      INTEGER , INTENT (in) ::   kiut    ! End point on which the  the computation is applied
-      INTEGER , INTENT (in) ::   jl      ! Category number
+      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
 
       !! * Local variables
@@ -99,10 +98,10 @@
       INTEGER ::   ii, ij      ! temporary dummy loop index
       INTEGER ::   numeq       ! current reference number of equation
-      INTEGER ::   layer       ! vertical dummy loop index 
+      INTEGER ::   jk       ! vertical dummy loop index 
       INTEGER ::   nconv       ! number of iterations in iterative procedure
       INTEGER ::   minnumeqmin, maxnumeqmax
-      INTEGER, DIMENSION(kiut) ::   numeqmin   ! reference number of top equation
-      INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation
-      INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow
+      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
+      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
+      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow
       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
       REAL(wp) ::   zg1       =  2._wp        !
@@ -111,53 +110,73 @@
       REAL(wp) ::   zraext_s  =  1.e+8_wp     ! extinction coefficient of radiation in the snow
       REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
+      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C 
       REAL(wp) ::   ztmelt_i    ! ice melting temperature
       REAL(wp) ::   zerritmax   ! current maximal error on temperature 
-      REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point
-      REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure )
-      REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration
-      REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness
-      REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness
-      REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function
-      REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function
-      REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature
-      REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4)
-      REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice
-      REAL(wp), DIMENSION(kiut) ::   zihic, zhsu
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat
-      REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow
-      REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow
-      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term
-      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term
-      REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis
-      REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms
+      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point
+      REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure )
+      REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration
+      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
+      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
+      REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface
+      REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function
+      REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function
+      REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature
+      REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4)
+      REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:) ::   zihic, zhsu
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zindterm    ! Independent term
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zindtbis    ! temporary independent term
+      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
       !!------------------------------------------------------------------     
       ! 
+      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow )
+      CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
+      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
+      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
+      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)
+      CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis  )
+      CALL wrk_alloc( jpij, nlay_i+3, 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_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
+            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
+      END DO
+
       !------------------------------------------------------------------------------!
       ! 1) Initialization                                                            !
       !------------------------------------------------------------------------------!
-      !
+      ! clem clean: replace just ztfs by rtt
       DO ji = kideb , kiut
          ! is there snow or not
-         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  )
+         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  )
          ! surface temperature of fusion
-!!gm ???  ztfs(ji) = rtt !!!????
          ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt
          ! layer thickness
-         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
-         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
+         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i )
+         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
       END DO
 
@@ -169,13 +188,13 @@
       z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
 
-      DO layer = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
-         DO ji = kideb , kiut
-            z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )
-         END DO
-      END DO
-
-      DO layer = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
-         DO ji = kideb , kiut
-            z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )
+      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
+         DO ji = kideb , kiut
+            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s )
+         END DO
+      END DO
+
+      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
+         DO ji = kideb , kiut
+            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i )
          END DO
       END DO
@@ -194,12 +213,12 @@
       ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
       ! zftrice = io.qsr_ice       is below the surface 
-      ! fstbif  = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
+      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
 
       DO ji = kideb , kiut
          ! switches
-         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  ) 
+         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  ) 
          ! hs > 0, isnow = 1
          zhsu (ji) = hnzst  ! threshold for the computation of i0
-         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )     
+         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) )     
 
          i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
@@ -208,5 +227,5 @@
          !            a function of the cloud cover
          !
-         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0)
+         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0)
          !formula used in Cice
       END DO
@@ -230,10 +249,10 @@
       END DO
 
-      DO layer = 1, nlay_s          ! Radiation through snow
+      DO jk = 1, nlay_s          ! Radiation through snow
          DO ji = kideb, kiut
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
             !                             ! radiation absorbed by the layer-th snow layer
-            zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
+            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
          END DO
       END DO
@@ -243,30 +262,16 @@
       END DO
 
-      DO layer = 1, nlay_i          ! Radiation through ice
+      DO jk = 1, nlay_i          ! Radiation through ice
          DO ji = kideb, kiut
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
             !                             ! radiation absorbed by the layer-th ice layer
-            zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
+            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
          END DO
       END DO
 
       DO ji = kideb, kiut           ! Radiation transmitted below the ice
-         fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif
-      END DO
-
-      ! +++++
-      ! just to check energy conservation
-      DO ji = kideb, kiut
-         ii = MOD( npb(ji) - 1 , jpi ) + 1
-         ij =    ( npb(ji) - 1 ) / jpi + 1
-         fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif
-      END DO
-      ! +++++
-
-      DO layer = 1, nlay_i
-         DO ji = kideb, kiut
-            radab(ji,layer) = zradab_i(ji,layer)
-         END DO
+         !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_1d(ji) / at_i_1d(ji) ! clem modif
+         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
       END DO
 
@@ -277,19 +282,19 @@
       !
       DO ji = kideb, kiut        ! Old surface temperature
-         ztsuold  (ji) =  t_su_b(ji)                              ! temperature at the beg of iter pr.
-         ztsuoldit(ji) =  t_su_b(ji)                              ! temperature at the previous iter
-         t_su_b   (ji) =  MIN( t_su_b(ji), ztfs(ji)-0.00001 )     ! necessary
+         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
+         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
+         t_su_1d   (ji) =  MIN( t_su_1d(ji), ztfs(ji) - ztsu_err )  ! necessary
          zerrit   (ji) =  1000._wp                                ! initial value of error
       END DO
 
-      DO layer = 1, nlay_s       ! Old snow temperature
-         DO ji = kideb , kiut
-            ztsold(ji,layer) =  t_s_b(ji,layer)
-         END DO
-      END DO
-
-      DO layer = 1, nlay_i       ! Old ice temperature
-         DO ji = kideb , kiut
-            ztiold(ji,layer) =  t_i_b(ji,layer)
+      DO jk = 1, nlay_s       ! Old snow temperature
+         DO ji = kideb , kiut
+            ztsb(ji,jk) =  t_s_1d(ji,jk)
+         END DO
+      END DO
+
+      DO jk = 1, nlay_i       ! Old ice temperature
+         DO ji = kideb , kiut
+            ztib(ji,jk) =  t_i_1d(ji,jk)
          END DO
       END DO
@@ -308,12 +313,12 @@
          IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula
             DO ji = kideb , kiut
-               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)
+               ztcond_i(ji,0)        = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt)
                ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin)
             END DO
-            DO layer = 1, nlay_i-1
+            DO jk = 1, nlay_i-1
                DO ji = kideb , kiut
-                  ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  &
-                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
-                  ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
+                  ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
+                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)
+                  ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin)
                END DO
             END DO
@@ -322,19 +327,20 @@
          IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
             DO ji = kideb , kiut
-               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   &
-                  &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )  
+               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt )   &
+                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rtt )  
                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
             END DO
-            DO layer = 1, nlay_i-1
+            DO jk = 1, nlay_i-1
                DO ji = kideb , kiut
-                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
-                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   &
-                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )  
-                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
+                  ztcond_i(ji,jk) = rcdic +                                                                     & 
+                     &                 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                          &
+                     &                 / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)   &
+                     &               - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt )  
+                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
                END DO
             END DO
             DO ji = kideb , kiut
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   &
-                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )  
+               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt)   &
+                  &                        - 0.011_wp * ( t_bo_1d(ji) - rtt )  
                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
             END DO
@@ -352,15 +358,15 @@
          END DO
 
-         DO layer = 1, nlay_s-1
-            DO ji = kideb , kiut
-               zkappa_s(ji,layer)  = 2.0 * rcdsn / &
+         DO jk = 1, nlay_s-1
+            DO ji = kideb , kiut
+               zkappa_s(ji,jk)  = 2.0 * rcdsn / &
                   MAX(epsi10,2.0*zh_s(ji))
             END DO
          END DO
 
-         DO layer = 1, nlay_i-1
+         DO jk = 1, nlay_i-1
             DO ji = kideb , kiut
                !-- Ice kappa factors
-               zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ &
+               zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ &
                   MAX(epsi10,2.0*zh_i(ji)) 
             END DO
@@ -381,18 +387,18 @@
          !------------------------------------------------------------------------------|
          !
-         DO layer = 1, nlay_i
-            DO ji = kideb , kiut
-               ztitemp(ji,layer)   = t_i_b(ji,layer)
-               zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
-                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
-               zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &
+         DO jk = 1, nlay_i
+            DO ji = kideb , kiut
+               ztitemp(ji,jk)   = t_i_1d(ji,jk)
+               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ &
+                  MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10)
+               zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), &
                   epsi10)
             END DO
          END DO
 
-         DO layer = 1, nlay_s
-            DO ji = kideb , kiut
-               ztstemp(ji,layer) = t_s_b(ji,layer)
-               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
+         DO jk = 1, nlay_s
+            DO ji = kideb , kiut
+               ztstemp(ji,jk) = t_s_1d(ji,jk)
+               zeta_s(ji,jk)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
             END DO
          END DO
@@ -403,14 +409,11 @@
          !
          DO ji = kideb , kiut
-
             ! update of the non solar flux according to the update in T_su
-            qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & 
-               ( t_su_b(ji) - ztsuoldit(ji) )
+            qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
 
             ! update incoming flux
             zf(ji)    =   zfsw(ji)              & ! net absorbed solar radiation
-               + qnsr_ice_1d(ji)           ! non solar total flux 
+               + qns_ice_1d(ji)                  ! non solar total flux 
             ! (LWup, LWdw, SH, LH)
-
          END DO
 
@@ -427,5 +430,5 @@
          !!ice interior terms (top equation has the same form as the others)
 
-         DO numeq=1,jkmax+2
+         DO numeq=1,nlay_i+3
             DO ji = kideb , kiut
                ztrid(ji,numeq,1) = 0.
@@ -440,11 +443,11 @@
          DO numeq = nlay_s + 2, nlay_s + nlay_i 
             DO ji = kideb , kiut
-               layer              = numeq - nlay_s - 1
-               ztrid(ji,numeq,1)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
-               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
-                  zkappa_i(ji,layer))
-               ztrid(ji,numeq,3)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer)
-               zindterm(ji,numeq) =  ztiold(ji,layer) + zeta_i(ji,layer)* &
-                  zradab_i(ji,layer)
+               jk              = numeq - nlay_s - 1
+               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk-1)
+               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + &
+                  zkappa_i(ji,jk))
+               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk)
+               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* &
+                  zradab_i(ji,jk)
             END DO
          ENDDO
@@ -457,12 +460,12 @@
                +  zkappa_i(ji,nlay_i-1) )
             ztrid(ji,numeq,3)  =  0.0
-            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
+            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &
                ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
-               *  t_bo_b(ji) ) 
+               *  t_bo_1d(ji) ) 
          ENDDO
 
 
          DO ji = kideb , kiut
-            IF ( ht_s_b(ji).gt.0.0 ) THEN
+            IF ( ht_s_1d(ji).gt.0.0 ) THEN
                !
                !------------------------------------------------------------------------------|
@@ -472,11 +475,11 @@
                !!snow interior terms (bottom equation has the same form as the others)
                DO numeq = 3, nlay_s + 1
-                  layer =  numeq - 1
-                  ztrid(ji,numeq,1)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
-                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
-                     zkappa_s(ji,layer) )
-                  ztrid(ji,numeq,3)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer)
-                  zindterm(ji,numeq)  =  ztsold(ji,layer) + zeta_s(ji,layer)* &
-                     zradab_s(ji,layer)
+                  jk =  numeq - 1
+                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk-1)
+                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + &
+                     zkappa_s(ji,jk) )
+                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
+                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* &
+                     zradab_s(ji,jk)
                END DO
 
@@ -485,8 +488,8 @@
                   ztrid(ji,nlay_s+2,3)    =  0.0
                   zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
-                     t_bo_b(ji) 
+                     t_bo_1d(ji) 
                ENDIF
 
-               IF ( t_su_b(ji) .LT. rtt ) THEN
+               IF ( t_su_1d(ji) .LT. rtt ) THEN
 
                   !------------------------------------------------------------------------------|
@@ -501,5 +504,5 @@
                   ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
                   ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
-                  zindterm(ji,1) = dzf(ji)*t_su_b(ji)   - zf(ji)
+                  zindterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji)
 
                   !!first layer of snow equation
@@ -507,5 +510,5 @@
                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)
                   ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
+                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
 
                ELSE 
@@ -524,7 +527,7 @@
                      zkappa_s(ji,0) * zg1s )
                   ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *            &
+                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            &
                      ( zradab_s(ji,1) +                         &
-                     zkappa_s(ji,0) * zg1s * t_su_b(ji) ) 
+                     zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
                ENDIF
             ELSE
@@ -534,5 +537,5 @@
                !------------------------------------------------------------------------------|
                !
-               IF (t_su_b(ji) .LT. rtt) THEN
+               IF (t_su_1d(ji) .LT. rtt) THEN
                   !
                   !------------------------------------------------------------------------------|
@@ -548,5 +551,5 @@
                   ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1    
                   ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_b(ji) - zf(ji)
+                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
 
                   !!first layer of ice equation
@@ -555,5 +558,5 @@
                      + zkappa_i(ji,0) * zg1 )
                   ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)  
-                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)  
+                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)  
 
                   !!case of only one layer in the ice (surface & ice equations are altered)
@@ -568,6 +571,6 @@
                      ztrid(ji,numeqmin(ji)+1,3)  =  0.0
 
-                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1)* &
-                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
+                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* &
+                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) )
                   ENDIF
 
@@ -588,6 +591,6 @@
                      zg1)  
                   ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
-                     zkappa_i(ji,0) * zg1 * t_su_b(ji) ) 
+                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
+                     zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
 
                   !!case of only one layer in the ice (surface & ice equations are altered)
@@ -597,7 +600,7 @@
                         zkappa_i(ji,1))
                      ztrid(ji,numeqmin(ji),3)  =  0.0
-                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1)* &
-                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
-                        + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
+                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* &
+                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) &
+                        + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
                   ENDIF
 
@@ -618,5 +621,5 @@
 
          maxnumeqmax = 0
-         minnumeqmin = jkmax+4
+         minnumeqmin = nlay_i+5
 
          DO ji = kideb , kiut
@@ -627,7 +630,7 @@
          END DO
 
-         DO layer = minnumeqmin+1, maxnumeqmax
-            DO ji = kideb , kiut
-               numeq               =  min(max(numeqmin(ji)+1,layer),numeqmax(ji))
+         DO jk = minnumeqmin+1, maxnumeqmax
+            DO ji = kideb , kiut
+               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
                   ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
@@ -639,12 +642,12 @@
          DO ji = kideb , kiut
             ! ice temperatures
-            t_i_b(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
+            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
          END DO
 
          DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
             DO ji = kideb , kiut
-               layer    =  numeq - nlay_s - 1
-               t_i_b(ji,layer)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
-                  t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
+               jk    =  numeq - nlay_s - 1
+               t_i_1d(ji,jk)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
+                  t_i_1d(ji,jk+1))/zdiagbis(ji,numeq)
             END DO
          END DO
@@ -652,15 +655,15 @@
          DO ji = kideb , kiut
             ! snow temperatures      
-            IF (ht_s_b(ji).GT.0._wp) &
-               t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
-               *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
-               *        MAX(0.0,SIGN(1.0,ht_s_b(ji))) 
+            IF (ht_s_1d(ji).GT.0._wp) &
+               t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
+               *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) &
+               *        MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 
 
             ! surface temperature
-            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  )
-            ztsuoldit(ji) = t_su_b(ji)
-            IF( t_su_b(ji) < ztfs(ji) ) &
-               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1)   &
-               &          + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
+            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) )  )  )
+            ztsubit(ji) = t_su_1d(ji)
+            IF( t_su_1d(ji) < ztfs(ji) ) &
+               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   &
+               &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
          END DO
          !
@@ -672,22 +675,20 @@
          ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
          DO ji = kideb , kiut
-            t_su_b(ji) =  MAX(  MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp  )
-            zerrit(ji) =  ABS( t_su_b(ji) - ztsuoldit(ji) )     
-         END DO
-
-         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)))
-            END DO
-         END DO
-
-         DO layer  =  1, nlay_i
-            DO ji = kideb , kiut
-               ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt 
-               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
-               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
+            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp  )
+            zerrit(ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
+         END DO
+
+         DO jk  =  1, nlay_s
+            DO ji = kideb , kiut
+               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rtt ), 190._wp  )
+               zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk)))
+            END DO
+         END DO
+
+         DO jk  =  1, nlay_i
+            DO ji = kideb , kiut
+               ztmelt_i        = -tmut * s_i_1d(ji,jk) + rtt 
+               t_i_1d(ji,jk) =  MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp)
+               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk)))
             END DO
          END DO
@@ -713,45 +714,87 @@
       !-------------------------------------------------------------------------!
       DO ji = kideb, kiut
-#if ! defined key_coupled
          ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 
-         qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )
-#endif
+         IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) )
          !                                ! surface ice conduction flux
-         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  )
-         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   &
-            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji))
+         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )  )
+         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
+            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
          !                                ! bottom ice conduction flux
-         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
-      END DO
-
-      !-------------------------!
-      ! Heat conservation       !
-      !-------------------------!
-      IF( con_i .AND. jiindex_1d > 0 ) THEN
+         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
+      END DO
+
+      !-----------------------------------------
+      ! Heat flux used to warm/cool ice in W.m-2
+      !-----------------------------------------
+      DO ji = kideb, kiut
+         IF( t_su_1d(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_1d(ji)
+         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_1d(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_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
+            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(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_1d(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_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
+      END DO
+   
+      !
+      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
+      CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
+      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
+      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   &
+         &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
+      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
+      CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis )
+      CALL wrk_dealloc( jpij, nlay_i+3, 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
-            ! Upper snow value
-            fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 
-            ! Bott. snow value
-            fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 
-         END DO
-         DO ji = kideb, kiut         ! Upper ice layer
-            fc_i(ji,0) = - REAL( isnow(ji) ) * &  ! interface flux if there is snow
-               ( zkappa_i(ji,0)  * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) &
-               - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * & 
-               zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not
-         END DO
-         DO layer = 1, nlay_i - 1         ! Internal ice layers
-            DO ji = kideb, kiut
-               fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) )
-               ii = MOD( npb(ji) - 1, jpi ) + 1
-               ij = ( npb(ji) - 1 ) / jpi + 1
-            END DO
-         END DO
-         DO ji = kideb, kiut         ! Bottom ice layers
-            fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
-         END DO
-      ENDIF
+            ztmelts      = - tmut  * s_i_1d(ji,jk) + rtt 
+            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) )
+            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                             &
+               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_1d(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_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus )
+         END DO
+      END DO
       !
-   END SUBROUTINE lim_thd_dif
+   END SUBROUTINE lim_thd_enmelt
 
 #else
