Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8514)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8518)
@@ -94,11 +94,13 @@
       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
       REAL(wp) ::   ztmelt_i                  ! ice melting temperature
-      REAL(wp) ::   zhsu
+      REAL(wp) ::   z1_hsu
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
       REAL(wp) ::   zdti_bnd = 1.e-4_wp       ! maximal authorized error on temperature 
+      REAL(wp) ::   zcpi                      ! Ice specific heat
+      REAL(wp) ::   zhi                       ! 
+      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
       
       REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
-      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! old surface temperature (before the iterative procedure )
-      REAL(wp), DIMENSION(jpij)     ::   ztsubit     ! surface temperature at previous iteration
+      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij)     ::   zh_i        ! ice layer thickness
       REAL(wp), DIMENSION(jpij)     ::   zh_s        ! snow layer thickness
@@ -106,31 +108,28 @@
       REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   dzf         ! derivative of the surface flux function
+      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
       REAL(wp), DIMENSION(jpij)     ::   zdti        ! current error on temperature
-      REAL(wp), DIMENSION(jpij)     ::   zdifcase    ! case of the equation resolution (1->4)
       REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
-      REAL(wp), DIMENSION(jpij)     ::   zihic
       
+      REAL(wp), DIMENSION(jpij,nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice
+      REAL(wp), DIMENSION(jpij,nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Old temperature in the ice
+      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
+      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztib        ! Old temperature in the ice
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zspeche_i   ! Ice specific heat
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztsb        ! Temporary temperature in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
-      REAL(wp), DIMENSION(jpij)            ::   zdq, zq_ini, zhfx_err ! diag errors on heat
+      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
       
@@ -148,5 +147,4 @@
 
       ! --- diag error on heat diffusion - PART 1 --- !
-      zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
       DO ji = 1, nidx
          zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
@@ -167,17 +165,15 @@
       ! Ice / snow layers
       !--------------------
-
-      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
-      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
-
-      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
-         DO ji = 1 , nidx
-            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_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 = 1 , nidx
-            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
+      z_s(1:nidx,1) = zh_s(1:nidx)
+      z_i(1:nidx,1) = zh_i(1:nidx)
+      DO jk = 2, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
+         DO ji = 1 , nidx
+            z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji)
+         END DO
+      END DO
+
+      DO jk = 2, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
+         DO ji = 1 , nidx
+            z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji)
          END DO
       END DO
@@ -187,34 +183,28 @@
       !------------------------------------------------------------------------------|
       !
-      !-------------------
-      ! Computation of i0
-      !-------------------
-      ! i0 describes the fraction of solar radiation which does not contribute
-      ! to the surface energy budget but rather penetrates inside the ice.
-      ! We assume that no radiation is transmitted through the snow
-      ! If there is no no snow
-      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
-      ! zftrice = io.qsr_ice       is below the surface 
-      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
-      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
-      zhsu = 0.1_wp ! threshold for the computation of i0
+      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
       DO ji = 1 , nidx
-         ! switches
-         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
-         ! hs > 0, isnow = 1
-         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )     
-
-         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
-      END DO
-
-      !-------------------------------------------------------
-      ! Solar radiation absorbed / transmitted at the surface
-      ! Derivative of the non solar flux
-      !-------------------------------------------------------
-      DO ji = 1 , nidx
-         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
-         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
-         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
-         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
+         !-------------------
+         ! Computation of i0
+         !-------------------
+         ! i0 describes the fraction of solar radiation which does not contribute
+         ! to the surface energy budget but rather penetrates inside the ice.
+         ! We assume that no radiation is transmitted through the snow
+         ! If there is no no snow
+         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
+         ! zftrice = io.qsr_ice       is below the surface 
+         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
+         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
+         zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )     
+         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) )
+
+         !-------------------------------------------------------
+         ! Solar radiation absorbed / transmitted at the surface
+         ! Derivative of the non solar flux
+         !-------------------------------------------------------
+         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
+         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
+         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
+         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
       END DO
 
@@ -222,26 +212,19 @@
       ! Transmission - absorption of solar radiation in the ice
       !---------------------------------------------------------
-
-      DO ji = 1, nidx           ! snow initialization
-         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
-      END DO
-
-      DO jk = 1, nlay_s          ! Radiation through snow
+      zradtr_s(1:nidx,0) = zftrice(1:nidx)
+      DO jk = 1, nlay_s
          DO ji = 1, nidx
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) )
             !                             ! radiation absorbed by the layer-th snow layer
             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
          END DO
       END DO
-
-      DO ji = 1, nidx           ! ice initialization
-         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
-      END DO
-
-      DO jk = 1, nlay_i          ! Radiation through ice
+         
+      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
+      DO jk = 1, nlay_i 
          DO ji = 1, nidx
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) )
             !                             ! radiation absorbed by the layer-th ice layer
             zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
@@ -249,7 +232,5 @@
       END DO
 
-      DO ji = 1, nidx           ! Radiation transmitted below the ice
-         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
-      END DO
+      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
 
       !------------------------------------------------------------------------------|
@@ -257,22 +238,9 @@
       !------------------------------------------------------------------------------|
       !
-      DO ji = 1, nidx        ! Old surface temperature
-         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), rt0 - ztsu_err )       ! necessary
-         zdti   (ji) =  1000._wp                                 ! initial value of error
-      END DO
-
-      DO jk = 1, nlay_s       ! Old snow temperature
-         DO ji = 1 , nidx
-            ztsb(ji,jk) =  t_s_1d(ji,jk)
-         END DO
-      END DO
-
-      DO jk = 1, nlay_i       ! Old ice temperature
-         DO ji = 1 , nidx
-            ztib(ji,jk) =  t_i_1d(ji,jk)
-         END DO
-      END DO
+      ztsub  (1:nidx)   =  t_su_1d(1:nidx)                          ! temperature at the previous iter
+      t_su_1d(1:nidx)   =  MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
+      zdti   (1:nidx)   =  1000._wp                                 ! initial value of error
+      ztsb   (1:nidx,:) =  t_s_1d(1:nidx,:)                         ! Old snow temperature
+      ztib   (1:nidx,:) =  t_i_1d(1:nidx,:)                         ! Old ice temperature
 
       iconv    =  0           ! number of iterations
@@ -289,12 +257,10 @@
          IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
             DO ji = 1 , nidx
-               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
+               ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1 , nidx
-                  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 * rt0)
-                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
+                  ztcond_i(ji,jk) = MAX( zkimin, 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 * rt0) )
                END DO
             END DO
@@ -302,21 +268,18 @@
          ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
             DO ji = 1 , nidx
-               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
-                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 )  
-               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
+               ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
+                  &                               - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1 , nidx
-                  ztcond_i(ji,jk) = rcdic +                                                                       & 
+                  ztcond_i(ji,jk) = MAX( zkimin, rcdic +                                                           & 
                      &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
                      &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   &
-                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  
-                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
+                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  )
                END DO
             END DO
             DO ji = 1 , nidx
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
-                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )  
-               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
+               ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
+                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) )  
             END DO
          ENDIF
@@ -335,5 +298,5 @@
          SELECT CASE ( nn_monocat )
 
-         CASE (1,3) ! LIM3
+         CASE (1,3)
 
             zepsilon = 0.1_wp
@@ -404,7 +367,7 @@
          DO jk = 1, nlay_i
             DO ji = 1 , nidx
-               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) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
+               ztitemp(ji,jk) = t_i_1d(ji,jk)
+               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
+               zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 )
             END DO
          END DO
@@ -425,5 +388,5 @@
             DO ji = 1 , nidx
                ! 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_1d(ji) - ztsubit(ji) )
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
             END DO
          ENDIF
@@ -507,5 +470,4 @@
                   !  case 1 : no surface melting - snow present                                  |
                   !------------------------------------------------------------------------------|
-                  zdifcase(ji)    =  1.0
                   numeqmin(ji)    =  1
                   numeqmax(ji)    =  nlay_i + nlay_s + 1
@@ -513,7 +475,7 @@
                   !!surface equation
                   ztrid(ji,1,1)  = 0.0
-                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0)
+                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
                   ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
-                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
+                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji)
 
                   !!first layer of snow equation
@@ -529,5 +491,4 @@
                   !------------------------------------------------------------------------------|
                   !
-                  zdifcase(ji)    =  2.0
                   numeqmin(ji)    =  2
                   numeqmax(ji)    =  nlay_i + nlay_s + 1
@@ -552,5 +513,4 @@
                   !------------------------------------------------------------------------------|
                   !
-                  zdifcase(ji)      =  3.0
                   numeqmin(ji)      =  nlay_s + 1
                   numeqmax(ji)      =  nlay_i + nlay_s + 1
@@ -558,7 +518,7 @@
                   !!surface equation	
                   ztrid(ji,numeqmin(ji),1)   =  0.0
-                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1    
+                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
                   ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
+                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji)
 
                   !!first layer of ice equation
@@ -572,5 +532,5 @@
                   IF ( nlay_i == 1 ) THEN
                      ztrid(ji,numeqmin(ji),1)    =  0.0
-                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
+                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
                      ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
                      ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
@@ -589,5 +549,4 @@
                   !------------------------------------------------------------------------------|
                   !
-                  zdifcase(ji)    =  4.0
                   numeqmin(ji)    =  nlay_s + 2
                   numeqmax(ji)    =  nlay_i + nlay_s + 1
@@ -662,5 +621,5 @@
             ! surface temperature
             isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
-            ztsubit(ji) = t_su_1d(ji)
+            ztsub(ji) = t_su_1d(ji)
             IF( t_su_1d(ji) < rt0 ) &
                t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
@@ -676,5 +635,5 @@
          DO ji = 1 , nidx
             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
-            zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
+            zdti   (ji) =  ABS( t_su_1d(ji) - ztsub(ji) )     
          END DO
 
@@ -754,31 +713,20 @@
 
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
+      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
+      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
       DO ji = 1, nidx
-         zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
-            &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
+         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
+            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
+
          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
+            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
          ELSE                          ! case T_su = 0degC
-            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 
+            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
          ENDIF
-         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
+         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
 
          ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
-         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
-      END DO 
-
-      !-----------------------------------------
-      ! Heat flux used to warm/cool ice in W.m-2
-      !-----------------------------------------
-      DO ji = 1, nidx
-         IF( t_su_1d(ji) < rt0 ) 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
-         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
+         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
+
       END DO   
       !
