Index: NEMO/branches/2018/dev_r10190_ice_thd_zdf/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- NEMO/branches/2018/dev_r10190_ice_thd_zdf/src/ICE/icethd_zdf_bl99.F90	(revision 10191)
+++ NEMO/branches/2018/dev_r10190_ice_thd_zdf/src/ICE/icethd_zdf_bl99.F90	(revision 10192)
@@ -108,23 +108,25 @@
       !
       REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice
-      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
-      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
-      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! 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)   ::   zeta_i      ! Eta factor in the ice
+      REAL(wp), DIMENSION(nlay_i, jpij)     ::   ztiold      ! Old temperature in the ice
+      REAL(wp), DIMENSION(nlay_s, jpij)     ::   ztsold      ! Old temperature in the snow
+      REAL(wp), DIMENSION(nlay_i, jpij)     ::   ztib        ! Temporary temperature in the ice to check the convergence
+      REAL(wp), DIMENSION(nlay_s, jpij)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
+      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   ztcond_i    ! Ice thermal conductivity
+      REAL(wp), DIMENSION(jpij, 0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
+      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   zradab_i    ! Radiation absorbed in the ice
+      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   zkappa_i    ! Kappa factor in the ice
+      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   zeta_i      ! Eta factor 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,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(0:nlay_s, jpij)   ::   zradab_s    ! Radiation absorbed in the snow
+      REAL(wp), DIMENSION(0:nlay_s, jpij)   ::   zkappa_s    ! Kappa factor in the snow
+      REAL(wp), DIMENSION(0:nlay_s, jpij)   ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), DIMENSION(nlay_i+3, jpij)   ::   zindterm    ! 'Ind'ependent term
+      REAL(wp), DIMENSION(nlay_i+3, jpij)   ::   zindtbis    ! Temporary 'ind'ependent term
+      REAL(wp), DIMENSION(nlay_i+3, jpij)   ::   zdiagbis    ! Temporary 'dia'gonal term
+      REAL(wp), DIMENSION(3, nlay_i+3, jpij)::   ztrid       ! Tridiagonal system terms
       REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
+      REAL(wp), DIMENSION(nlay_i, jpij)     ::   tt_i_1d, tsz_i_1d
+      REAL(wp), DIMENSION(nlay_s, jpij)     ::   tt_s_1d
       !
       ! Mono-category
@@ -132,6 +134,8 @@
       REAL(wp) ::   zhe        ! dummy factor
       REAL(wp) ::   zcnd_i     ! mean sea ice thermal conductivity
+      INTEGER  ::   itot
+      LOGICAL, DIMENSION(jpij)  ::   liter
       !!------------------------------------------------------------------     
-
+      itot = npti
       ! --- diag error on heat diffusion - PART 1 --- !
       DO ji = 1, npti
@@ -171,6 +175,20 @@
       ENDIF
       !
-      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
-      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
+      DO jk = 1, nlay_i
+         DO ji = 1, npti
+            tt_i_1d(jk, ji) = t_i_1d(ji, jk)
+            tsz_i_1d(jk, ji) = sz_i_1d(ji, jk)
+         END DO
+      END DO
+
+      DO jk = 1, nlay_s
+         DO ji = 1, npti
+            tt_s_1d(jk, ji) = t_s_1d(ji, jk)
+         END DO
+      END DO
+
+
+      ztsold (:, :) = tt_s_1d(:,:)   ! Old snow temperature
+      ztiold (:, :) = tt_i_1d(:,:)   ! Old ice temperature
 
       !-------------
@@ -184,129 +202,113 @@
             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) )
             !                             ! radiation absorbed by the layer-th snow layer
-            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
+            zradab_s(jk,ji) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
          END DO
       END DO
       !
-      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
+      zradtr_i(1:npti, 0) = zradtr_s(1:npti, nlay_s) * isnow(1:npti) + qtr_ice_top_1d(1:npti) * ( 1._wp - isnow(1:npti) )
       DO jk = 1, nlay_i 
          DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
+            zradtr_i(ji, jk) = zradtr_i(ji, 0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
             !                             ! radiation absorbed by the layer-th ice layer
-            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
+            zradab_i(jk,ji) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
          END DO
       END DO
       !
-      qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
+      qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti, nlay_i)   ! record radiation transmitted below the ice
       !
       iconv    = 0          ! number of iterations
       zdti_max = 1000._wp   ! maximal value of error on all points
       !
-      !                                                          !============================!
-      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
-         !                                                       !============================!
+      liter(:) = .TRUE.                                          !============================!
+!     DO iconv = 1, iconv_max                                    ! Iterative procedure begins !
+      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max)
+         zdti_max = 0._wp
+!        IF(itot<1) EXIT
          iconv = iconv + 1
-         !
-         ztib(1:npti,:) = t_i_1d(1:npti,:)
-         ztsb(1:npti,:) = t_s_1d(1:npti,:)
-         !
-         !--------------------------------
-         ! 3) Sea ice thermal conductivity
-         !--------------------------------
-         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
-            !
-            DO ji = 1, npti
-               ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
-            END DO
-            DO jk = 1, nlay_i-1
-               DO ji = 1, npti
-                  ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  &
-                     &                       MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+         DO ji = 1, npti
+          IF(liter(ji)) THEN
+!            zdti_max = 0._wp
+            !
+            ztib(:, ji) = tt_i_1d(:, ji)
+            ztsb(:, ji) = tt_s_1d(:, ji)
+            !
+            !--------------------------------
+            ! 3) Sea ice thermal conductivity
+            !--------------------------------
+            IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
+            !
+               ztcond_i(0, ji)      = rcnd_i + zbeta * tsz_i_1d(1, ji)      / MIN( -epsi10, tt_i_1d(1, ji) - rt0 )
+               ztcond_i(nlay_i, ji) = rcnd_i + zbeta * tsz_i_1d(nlay_i, ji) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
+               DO jk = 1, nlay_i-1
+                    ztcond_i(jk, ji) = rcnd_i + zbeta * 0.5_wp * ( tsz_i_1d(jk, ji) + tsz_i_1d(jk+1, ji) ) /  &
+                        &                       MIN( -epsi10, 0.5_wp * (tt_i_1d(jk, ji) + tt_i_1d(jk+1, ji)) - rt0 )
+                  END DO
+               !
+            ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
+               !
+                  ztcond_i(0, ji)      = rcnd_i + 0.09_wp  *  tsz_i_1d(1, ji)      / MIN( -epsi10, tt_i_1d(1, ji) - rt0 )  &
+                     &                         - 0.011_wp * ( tt_i_1d(1, ji) - rt0 )
+                  ztcond_i(nlay_i, ji) = rcnd_i + 0.09_wp  *  tsz_i_1d(nlay_i, ji) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
+                     &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 )
+               DO jk = 1, nlay_i-1
+                  ztcond_i(jk,ji) = rcnd_i + 0.09_wp  *   0.5_wp * ( tsz_i_1d(jk,ji) + tsz_i_1d(jk+1,ji) ) /        &
+                     &                      MIN( -epsi10, 0.5_wp * ( tt_i_1d (jk,ji) + tt_i_1d (jk+1,ji) ) - rt0 )  &
+                     &                     - 0.011_wp * ( 0.5_wp * ( tt_i_1d (jk,ji) + tt_i_1d (jk+1,ji) ) - rt0 )
                END DO
-            END DO
-            !
-         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
-            !
-            DO ji = 1, npti
-               ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
-                  &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  &
-                  &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 )
-            END DO
-            DO jk = 1, nlay_i-1
-               DO ji = 1, npti
-                  ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        &
-                     &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  &
-                     &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )
-               END DO
-            END DO
-            !
-         ENDIF
-         ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )        
-         !
-         !--- G(he) : enhancement of thermal conductivity in mono-category case
-         ! Computation of effective thermal conductivity G(h)
-         ! Used in mono-category case only to simulate an ITD implicitly
-         ! Fichefet and Morales Maqueda, JGR 1997
-         zghe(1:npti) = 1._wp
-         !
-         SELECT CASE ( nn_virtual_itd )
-         !
-         CASE ( 1 , 2 )
+            !
+            ENDIF
+            ztcond_i(:, ji) = MAX( zkimin, ztcond_i(:, ji) )        
+            !
+            !--- G(he) : enhancement of thermal conductivity in mono-category case
+            ! Computation of effective thermal conductivity G(h)
+            ! Used in mono-category case only to simulate an ITD implicitly
+            ! Fichefet and Morales Maqueda, JGR 1997
+            zghe(ji) = 1._wp
+            !
+            SELECT CASE ( nn_virtual_itd )
+            !
+            CASE ( 1 , 2 )
             !
             zepsilon = 0.1_wp
-            DO ji = 1, npti
-               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
+               zcnd_i = SUM( ztcond_i(:,ji) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity
                zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe)
                IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  &
                   &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he)
+            !
+            END SELECT
+            !
+            !-----------------
+            ! 4) kappa factors
+            !-----------------
+            !--- Snow
+            DO jk = 0, nlay_s-1
+               zkappa_s(jk,ji) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
             END DO
-            !
-         END SELECT
-         !
-         !-----------------
-         ! 4) kappa factors
-         !-----------------
-         !--- Snow
-         DO jk = 0, nlay_s-1
-            DO ji = 1, npti
-               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
+            ! Snow-ice interface
+            zfac = 0.5_wp * ( ztcond_i(0, ji) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
+            IF( zfac > epsi10 ) THEN
+               zkappa_s(nlay_s, ji) = zghe(ji) * rn_cnd_s * ztcond_i(0, ji) / zfac
+            ELSE
+               zkappa_s(nlay_s, ji) = 0._wp
+            ENDIF
+
+            !--- Ice
+            DO jk = 0, nlay_i
+               zkappa_i(jk,ji) = zghe(ji) * ztcond_i(jk,ji) * z1_h_i(ji)
             END DO
-         END DO
-         DO ji = 1, npti   ! Snow-ice interface
-            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
-            IF( zfac > epsi10 ) THEN
-               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
-            ELSE
-               zkappa_s(ji,nlay_s) = 0._wp
-            ENDIF
-         END DO
-
-         !--- Ice
-         DO jk = 0, nlay_i
-            DO ji = 1, npti
-               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
+            zkappa_i(0, ji) = zkappa_s(nlay_s, ji) * isnow(ji) + zkappa_i(0, ji) * ( 1._wp - isnow(ji) )
+            !
+            !--------------------------------------
+            ! 5) Sea ice specific heat, eta factors
+            !--------------------------------------
+            DO jk = 1, nlay_i
+               zcpi = rcpi + zgamma * tsz_i_1d(jk,ji) / MAX( ( tt_i_1d(jk,ji) - rt0 ) * ( ztiold(jk,ji) - rt0 ), epsi10 )
+               zeta_i(jk,ji) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )
             END DO
-         END DO
-         DO ji = 1, npti   ! Snow-ice interface
-            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
-         END DO
-         !
-         !--------------------------------------
-         ! 5) Sea ice specific heat, eta factors
-         !--------------------------------------
-         DO jk = 1, nlay_i
-            DO ji = 1, npti
-               zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
+
+            DO jk = 1, nlay_s
+               zeta_s(jk,ji) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
             END DO
-         END DO
-
-         DO jk = 1, nlay_s
-            DO ji = 1, npti
-               zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji)
-            END DO
-         END DO
          !
          !----------------------------------------!
@@ -325,11 +327,7 @@
             !----------------------------
             ! update of the non solar flux according to the update in T_su
-            DO ji = 1, npti
                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
-            END DO
-
-            DO ji = 1, npti
+
                zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar
-            END DO
             !
             !----------------------------
@@ -342,30 +340,25 @@
 
             ! ice interior terms (top equation has the same form as the others)
-            ztrid   (1:npti,:,:) = 0._wp
-            zindterm(1:npti,:)   = 0._wp
-            zindtbis(1:npti,:)   = 0._wp
-            zdiagbis(1:npti,:)   = 0._wp
+            ztrid   (:, :, ji) = 0._wp
+            zindterm(:, ji)   = 0._wp
+            zindtbis(:, ji)   = 0._wp
+            zdiagbis(:, ji)   = 0._wp
 
             DO jm = nlay_s + 2, nlay_s + nlay_i 
-               DO ji = 1, npti
                   jk = jm - nlay_s - 1
-                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
-                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
-                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-               END DO
+                  ztrid   (1, jm, ji) =       - zeta_i(jk,ji) *   zkappa_i(jk-1, ji)
+                  ztrid   (2, jm, ji) = 1._wp + zeta_i(jk,ji) * ( zkappa_i(jk-1, ji) + zkappa_i(jk, ji) )
+                  ztrid   (3, jm, ji) =       - zeta_i(jk,ji) *                       zkappa_i(jk, ji)
+                  zindterm(jm,ji)   = ztiold(jk,ji) + zeta_i(jk,ji) * zradab_i(jk, ji)
             END DO
 
             jm =  nlay_s + nlay_i + 1
-            DO ji = 1, npti
                ! ice bottom term
-               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
-               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
-               ztrid   (ji,jm,3) = 0._wp
-               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
-            END DO
-
-            DO ji = 1, npti
+               ztrid   (1, jm, ji) =       - zeta_i(nlay_i, ji) *   zkappa_i(nlay_i-1, ji)   
+               ztrid   (2, jm, ji) = 1._wp + zeta_i(nlay_i, ji) * ( zkappa_i(nlay_i-1, ji) + zkappa_i(nlay_i, ji) * zg1 )
+               ztrid   (3, jm, ji) = 0._wp
+               zindterm(jm,ji)   = ztiold(nlay_i, ji) + zeta_i(nlay_i, ji) *  &
+                  &              ( zradab_i(nlay_i, ji) + zkappa_i(nlay_i, ji) * zg1 * t_bo_1d(ji) ) 
+
                !                               !---------------------!
                IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
@@ -374,14 +367,14 @@
                   DO jm = 3, nlay_s + 1
                      jk = jm - 1
-                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
-                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
-                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+                     ztrid   (1, jm, ji) =       - zeta_s(jk,ji) *   zkappa_s(jk-1,ji)
+                     ztrid   (2, jm, ji) = 1._wp + zeta_s(jk,ji) * ( zkappa_s(jk-1,ji) + zkappa_s(jk,ji) )
+                     ztrid   (3, jm, ji) =       - zeta_s(jk,ji) *                       zkappa_s(jk,ji)
+                     zindterm(jm,ji)   = ztsold(jk,ji) + zeta_s(jk,ji) * zradab_s(jk,ji)
                   END DO
                   
                   ! case of only one layer in the ice (ice equation is altered)
                   IF( nlay_i == 1 ) THEN
-                     ztrid   (ji,nlay_s+2,3) = 0._wp
-                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
+                     ztrid   (3, nlay_s+2, ji) = 0._wp
+                     zindterm(nlay_s+2, ji)   = zindterm(nlay_s+2, ji) + zeta_i(1, ji) * zkappa_i(1, ji) * t_bo_1d(ji) 
                   ENDIF
                   
@@ -392,14 +385,14 @@
                      
                      ! surface equation
-                     ztrid   (ji,1,1) = 0._wp
-                     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)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
+                     ztrid   (1, 1, ji) = 0._wp
+                     ztrid   (2, 1, ji) = zdqns_ice_b(ji) - zg1s * zkappa_s(0, ji)
+                     ztrid   (3, 1, ji) =                   zg1s * zkappa_s(0, ji)
+                     zindterm(1, ji)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
                      
                      ! first layer of snow equation
-                     ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s
-                     ztrid   (ji,2,2) = 1._wp + 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)
+                     ztrid   (1, 2, ji) =       - zeta_s(1, ji) *                    zkappa_s(0, ji) * zg1s
+                     ztrid   (2, 2, ji) = 1._wp + zeta_s(1, ji) * ( zkappa_s(1, ji) + zkappa_s(0, ji) * zg1s )
+                     ztrid   (3, 2, ji) =       - zeta_s(1, ji) *   zkappa_s(1, ji)
+                     zindterm(2, ji)   = ztsold(1, ji) + zeta_s(1, ji) * zradab_s(1, ji)
                      
                   ELSE                            !--  case 2 : surface is melting
@@ -409,8 +402,8 @@
                      
                      ! first layer of snow equation
-                     ztrid   (ji,2,1) = 0._wp
-                     ztrid   (ji,2,2) = 1._wp + 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) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
+                     ztrid   (1, 2, ji) = 0._wp
+                     ztrid   (2, 2, ji) = 1._wp + zeta_s(1, ji) * ( zkappa_s(1, ji) + zkappa_s(0, ji) * zg1s )
+                     ztrid   (3, 2, ji) =       - zeta_s(1, ji) *   zkappa_s(1, ji) 
+                     zindterm(2, ji)   = ztsold(1, ji) + zeta_s(1, ji) * ( zradab_s(1, ji) + zkappa_s(0, ji) * zg1s * t_su_1d(ji) ) 
                   ENDIF
                   !                            !---------------------!
@@ -424,24 +417,24 @@
                      
                      ! surface equation	
-                     ztrid   (ji,jm_min(ji),1) = 0._wp
-                     ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1    
-                     ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1
-                     zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
+                     ztrid   (1, jm_min(ji), ji) = 0._wp
+                     ztrid   (2, jm_min(ji), ji) = zdqns_ice_b(ji) - zkappa_i(0, ji) * zg1    
+                     ztrid   (3, jm_min(ji), ji) =                   zkappa_i(0, ji) * zg1
+                     zindterm(jm_min(ji), ji)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
                      
                      ! first layer of ice equation
-                     ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1
-                     ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
-                     ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)  
-                     zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
+                     ztrid   (1, jm_min(ji)+1, ji) =       - zeta_i(1, ji) *                    zkappa_i(0, ji) * zg1
+                     ztrid   (2, jm_min(ji)+1, ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(1, ji) + zkappa_i(0, ji) * zg1 )
+                     ztrid   (3, jm_min(ji)+1, ji) =       - zeta_i(1, ji) *   zkappa_i(1, ji)  
+                     zindterm(jm_min(ji)+1, ji)   = ztiold(1, ji) + zeta_i(1, ji) * zradab_i(1, ji)  
                      
                      ! case of only one layer in the ice (surface & ice equations are altered)
                      IF( nlay_i == 1 ) THEN
-                        ztrid   (ji,jm_min(ji),1)   = 0._wp
-                        ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp
-                        ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp
-                        ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp
-                        ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
-                        ztrid   (ji,jm_min(ji)+1,3) = 0._wp
-                        zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji))
+                        ztrid   (1, jm_min(ji), ji)   = 0._wp
+                        ztrid   (2, jm_min(ji), ji)   = zdqns_ice_b(ji)      -   zkappa_i(0, ji) * 2._wp
+                        ztrid   (3, jm_min(ji), ji)   =                          zkappa_i(0, ji) * 2._wp
+                        ztrid   (1, jm_min(ji)+1, ji) =       - zeta_i(1, ji) *   zkappa_i(0, ji) * 2._wp
+                        ztrid   (2, jm_min(ji)+1, ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(0, ji) * 2._wp + zkappa_i(1, ji) )
+                        ztrid   (3, jm_min(ji)+1, ji) = 0._wp
+                        zindterm(jm_min(ji)+1, ji)   = ztiold(1, ji) + zeta_i(1, ji) * (zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji))
                      ENDIF
                      
@@ -452,16 +445,16 @@
                      
                      ! first layer of ice equation
-                     ztrid   (ji,jm_min(ji),1) = 0._wp
-                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
-                     ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1)
-                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * (zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji)) 
+                     ztrid   (1, jm_min(ji), ji) = 0._wp
+                     ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(1, ji) + zkappa_i(0, ji) * zg1 )  
+                     ztrid   (3, jm_min(ji), ji) =       - zeta_i(1, ji) *   zkappa_i(1, ji)
+                     zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) * (zradab_i(1, ji) + zkappa_i(0, ji) * zg1 * t_su_1d(ji)) 
                      
                      ! case of only one layer in the ice (surface & ice equations are altered)
                      IF( nlay_i == 1 ) THEN
-                        ztrid   (ji,jm_min(ji),1) = 0._wp
-                        ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) )
-                        ztrid   (ji,jm_min(ji),3) = 0._wp
-                        zindterm(ji,jm_min(ji))   = ztiold(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._wp
+                        ztrid   (1, jm_min(ji), ji) = 0._wp
+                        ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(0, ji) * 2._wp + zkappa_i(1, ji) )
+                        ztrid   (3, jm_min(ji), ji) = 0._wp
+                        zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) * ( zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji) ) &
+                           &                      + t_su_1d(ji) * zeta_i(1, ji) * zkappa_i(0, ji) * 2._wp
                      ENDIF
                      
@@ -469,8 +462,7 @@
                ENDIF
                !
-               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
-               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
+               zindtbis(jm_min(ji), ji) = zindterm(jm_min(ji), ji)
+               zdiagbis(jm_min(ji), ji) = ztrid   (2, jm_min(ji), ji)
                !
-            END DO
             !
             !------------------------------
@@ -479,43 +471,32 @@
             ! Solve the tridiagonal system with Gauss elimination method.
             ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
-            jm_maxt = 0
-            jm_mint = nlay_i+5
-            DO ji = 1, npti
-               jm_mint = MIN(jm_min(ji),jm_mint)
-               jm_maxt = MAX(jm_max(ji),jm_maxt)
-            END DO
-
-            DO jk = jm_mint+1, jm_maxt
-               DO ji = 1, npti
-                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
-                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
-                  zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1)
+
+
+
+
+
+               DO jm = jm_min(ji)+1, jm_max(ji)
+                  zdiagbis(jm,ji) = ztrid   (2, jm, ji) - ztrid(1, jm, ji) * ztrid   (3, jm-1, ji) / zdiagbis(jm-1,ji)
+                  zindtbis(jm,ji) = zindterm(jm,ji  ) - ztrid(1, jm, ji) * zindtbis(jm-1,ji  ) / zdiagbis(jm-1,ji)
                END DO
-            END DO
 
             ! ice temperatures
-            DO ji = 1, npti
-               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
-            END DO
-
-            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-               DO ji = 1, npti
+               tt_i_1d(nlay_i, ji) = zindtbis(jm_max(ji), ji) / zdiagbis(jm_max(ji), ji)
+
+               DO jm = nlay_i + nlay_s, nlay_s + 2, -1
                   jk = jm - nlay_s - 1
-                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+                  tt_i_1d(jk, ji) = ( zindtbis(jm,ji) - ztrid(3, jm,ji) * tt_i_1d(jk+1, ji) ) / zdiagbis(jm,ji)
                END DO
-            END DO
-
-            DO ji = 1, npti
+
                ! snow temperatures      
                IF( h_s_1d(ji) > 0._wp ) THEN
-                  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)
+                  tt_s_1d(nlay_s, ji) = ( zindtbis(nlay_s+1, ji) - ztrid(3, nlay_s+1, ji) * tt_i_1d(1, ji) ) / zdiagbis(nlay_s+1, ji)
                ENDIF
                ! surface temperature
                ztsub(ji) = t_su_1d(ji)
                IF( t_su_1d(ji) < rt0 ) THEN
-                  t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
-                     &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
+                  t_su_1d(ji) = ( zindtbis(jm_min(ji), ji) - ztrid(3,jm_min(ji),ji) *  &
+                     &          ( isnow(ji) * tt_s_1d(1, ji) + ( 1._wp - isnow(ji) ) * tt_i_1d(1, ji) ) ) / zdiagbis(jm_min(ji), ji)
                ENDIF
-            END DO
             !
             !--------------------------------------------------------------
@@ -524,28 +505,18 @@
             ! check that nowhere it has started to melt
             ! zdti_max is a measure of error, it has to be under zdti_bnd
-            zdti_max = 0._wp
-            DO ji = 1, npti
+
                t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
                zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
-            END DO
-
-            DO jk = 1, nlay_s
-               DO ji = 1, npti
-                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
-                  zdti_max      = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+
+               DO jk = 1, nlay_s
+                  tt_s_1d(jk,ji) = MAX( MIN( tt_s_1d(jk,ji), rt0 ), rt0 - 100._wp )
+                  zdti_max      = MAX( zdti_max, ABS( tt_s_1d(jk,ji) - ztsb(jk,ji) ) )
                END DO
-            END DO
-
-            DO jk = 1, nlay_i
-               DO ji = 1, npti
-                  ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
-                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
-                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+
+               DO jk = 1, nlay_i
+                  ztmelts       = -rTmlt * tsz_i_1d(jk, ji) + rt0 
+                  tt_i_1d(jk, ji) =  MAX( MIN( tt_i_1d(jk, ji), ztmelts ), rt0 - 100._wp )
+                  zdti_max      =  MAX( zdti_max, ABS( tt_i_1d(jk, ji) - ztib(jk,ji) ) )
                END DO
-            END DO
-
-            ! Compute spatial maximum over all errors
-            ! note that this could be optimized substantially by iterating only the non-converging points
-            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
          !
          !----------------------------------------!
@@ -568,30 +539,25 @@
 
             ! ice interior terms (top equation has the same form as the others)
-            ztrid   (1:npti,:,:) = 0._wp
-            zindterm(1:npti,:)   = 0._wp
-            zindtbis(1:npti,:)   = 0._wp
-            zdiagbis(1:npti,:)   = 0._wp
-
-            DO jm = nlay_s + 2, nlay_s + nlay_i 
-               DO ji = 1, npti
+               ztrid   (:, :, ji) = 0._wp
+               zindterm(:, ji)   = 0._wp
+               zindtbis(:, ji)   = 0._wp
+               zdiagbis(:, ji)   = 0._wp
+
+               DO jm = nlay_s + 2, nlay_s + nlay_i 
                   jk = jm - nlay_s - 1
-                  ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1)
-                  ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-                  ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk)
-                  zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-               END DO
-            ENDDO
-
-            jm =  nlay_s + nlay_i + 1
-            DO ji = 1, npti
+                  ztrid   (1, jm, ji) =       - zeta_i(jk,ji) *   zkappa_i(jk-1, ji)
+                  ztrid   (2, jm, ji) = 1._wp + zeta_i(jk,ji) * ( zkappa_i(jk-1, ji) + zkappa_i(jk,ji) )
+                  ztrid   (3, jm, ji) =       - zeta_i(jk,ji) *                       zkappa_i(jk,ji)
+                  zindterm(jm,ji)   = ztiold(jk,ji) + zeta_i(jk,ji) * zradab_i(jk,ji)
+               ENDDO
+
+               jm =  nlay_s + nlay_i + 1
                ! ice bottom term
-               ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)   
-               ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 )
-               ztrid   (ji,jm,3) = 0._wp
-               zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-                  &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 
-            ENDDO
-
-            DO ji = 1, npti
+               ztrid   (1, jm, ji) =       - zeta_i(nlay_i, ji) *   zkappa_i(nlay_i-1, ji)   
+               ztrid   (2, jm, ji) = 1._wp + zeta_i(nlay_i, ji) * ( zkappa_i(nlay_i-1, ji) + zkappa_i(nlay_i, ji) * zg1 )
+               ztrid   (3, jm, ji) = 0._wp
+               zindterm(jm,ji)   = ztiold(nlay_i, ji) + zeta_i(nlay_i, ji) *  &
+                  &              ( zradab_i(nlay_i, ji) + zkappa_i(nlay_i, ji) * zg1 * t_bo_1d(ji) ) 
+
                !                               !---------------------!
                IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells !
@@ -600,14 +566,14 @@
                   DO jm = 3, nlay_s + 1
                      jk = jm - 1
-                     ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1)
-                     ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                     ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk)
-                     zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+                     ztrid   (1, jm, ji) =       - zeta_s(jk,ji) *   zkappa_s(jk-1,ji)
+                     ztrid   (2, jm, ji) = 1._wp + zeta_s(jk,ji) * ( zkappa_s(jk-1,ji) + zkappa_s(jk,ji) )
+                     ztrid   (3, jm, ji) =       - zeta_s(jk,ji) *                       zkappa_s(jk,ji)
+                     zindterm(jm,ji)   = ztsold(jk,ji) + zeta_s(jk,ji) * zradab_s(jk,ji)
                   END DO
                   
                   ! case of only one layer in the ice (ice equation is altered)
                   IF ( nlay_i == 1 ) THEN
-                     ztrid   (ji,nlay_s+2,3) = 0._wp
-                     zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji) 
+                     ztrid   (3, nlay_s+2, ji) = 0._wp
+                     zindterm(nlay_s+2, ji)   = zindterm(nlay_s+2, ji) + zeta_i(1, ji) * zkappa_i(1, ji) * t_bo_1d(ji) 
                   ENDIF
                   
@@ -616,8 +582,8 @@
                   
                   ! first layer of snow equation
-                  ztrid   (ji,2,1) = 0._wp
-                  ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1)
-                  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) + qcn_ice_1d(ji) ) 
+                  ztrid   (1, 2, ji) = 0._wp
+                  ztrid   (2, 2, ji) = 1._wp + zeta_s(1, ji) * zkappa_s(1, ji)
+                  ztrid   (3, 2, ji) =       - zeta_s(1, ji) * zkappa_s(1, ji) 
+                  zindterm(2, ji)   = ztsold(1, ji) + zeta_s(1, ji) * ( zradab_s(1, ji) + qcn_ice_1d(ji) ) 
                   
                   !                            !---------------------!
@@ -628,24 +594,23 @@
                   
                   ! first layer of ice equation
-                  ztrid   (ji,jm_min(ji),1) = 0._wp
-                  ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)  
-                  ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) )
+                  ztrid   (1, jm_min(ji), ji) = 0._wp
+                  ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * zkappa_i(1, ji)  
+                  ztrid   (3, jm_min(ji), ji) =       - zeta_i(1, ji) * zkappa_i(1, ji)
+                  zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) * ( zradab_i(1, ji) + qcn_ice_1d(ji) )
                   
                   ! case of only one layer in the ice (surface & ice equations are altered)
                   IF( nlay_i == 1 ) THEN
-                     ztrid   (ji,jm_min(ji),1) = 0._wp
-                     ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1)
-                     ztrid   (ji,jm_min(ji),3) = 0._wp
-                     zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  &
-                        &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) )
+                     ztrid   (1, jm_min(ji), ji) = 0._wp
+                     ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * zkappa_i(1, ji)
+                     ztrid   (3, jm_min(ji), ji) = 0._wp
+                     zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) *  &
+                        &                                     ( zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji) + qcn_ice_1d(ji) )
                   ENDIF
                   
                ENDIF
                !
-               zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
-               zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2)
+               zindtbis(jm_min(ji), ji) = zindterm(jm_min(ji), ji)
+               zdiagbis(jm_min(ji), ji) = ztrid   (2, jm_min(ji), ji)
                !
-            END DO
             !
             !------------------------------
@@ -654,37 +619,26 @@
             ! Solve the tridiagonal system with Gauss elimination method.
             ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
-            jm_maxt = 0
-            jm_mint = nlay_i+5
-            DO ji = 1, npti
-               jm_mint = MIN(jm_min(ji),jm_mint)
-               jm_maxt = MAX(jm_max(ji),jm_maxt)
-            END DO
-            
-            DO jk = jm_mint+1, jm_maxt
-               DO ji = 1, npti
-                  jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji))
-                  zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1)
-                  zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1)
-               END DO
+
+
+
+
+
+            DO jm = jm_min(ji)+1, jm_max(ji)
+                  zdiagbis(jm,ji) = ztrid   (2, jm, ji) - ztrid(1, jm, ji) * ztrid   (3, jm-1, ji) / zdiagbis(jm-1,ji)
+                  zindtbis(jm,ji) = zindterm(jm,ji)   - ztrid(1, jm,ji) * zindtbis(jm-1,ji)   / zdiagbis(jm-1,ji)
             END DO
             
             ! ice temperatures
-           DO ji = 1, npti
-               t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
-            END DO
-
-            DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-               DO ji = 1, npti
+               tt_i_1d(nlay_i, ji) = zindtbis(jm_max(ji), ji) / zdiagbis(jm_max(ji), ji)
+
+               DO jm = nlay_i + nlay_s, nlay_s + 2, -1
                   jk = jm - nlay_s - 1
-                  t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
+                  tt_i_1d(jk, ji) = ( zindtbis(jm,ji) - ztrid(3, jm, ji) * tt_i_1d(jk+1, ji) ) / zdiagbis(jm,ji)
                END DO
-            END DO
             
             ! snow temperatures      
-            DO ji = 1, npti
                IF( h_s_1d(ji) > 0._wp ) THEN
-                  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)
+                  tt_s_1d(nlay_s, ji) = ( zindtbis(nlay_s+1, ji) - ztrid(3, nlay_s+1, ji) * tt_i_1d(1, ji) ) / zdiagbis(nlay_s+1, ji)
                ENDIF
-            END DO
             !
             !--------------------------------------------------------------
@@ -693,33 +647,33 @@
             ! check that nowhere it has started to melt
             ! zdti_max is a measure of error, it has to be under zdti_bnd
-            zdti_max = 0._wp
-
-            DO jk = 1, nlay_s
-               DO ji = 1, npti
-                  t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
-                  zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+
+
+               DO jk = 1, nlay_s
+                     tt_s_1d(jk,ji) = MAX( MIN( tt_s_1d(jk,ji), rt0 ), rt0 - 100._wp )
+                     zdti_max = MAX( zdti_max, ABS( tt_s_1d(jk,ji) - ztsb(jk,ji) ) )
                END DO
-            END DO
             
-            DO jk = 1, nlay_i
-               DO ji = 1, npti
-                  ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0 
-                  t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp )
-                  zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+               DO jk = 1, nlay_i
+                     ztmelts       = -rTmlt * tsz_i_1d(jk,ji) + rt0 
+                     tt_i_1d(jk,ji) =  MAX( MIN( tt_i_1d(jk,ji), ztmelts ), rt0 - 100._wp )
+                     zdti_max      =  MAX( zdti_max, ABS( tt_i_1d(jk,ji) - ztib(jk,ji) ) )
                END DO
-            END DO
-
-            ! Compute spatial maximum over all errors
-            ! note that this could be optimized substantially by iterating only the non-converging points
-            IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice )
-
-         ENDIF ! k_jules
-         
+
+        ENDIF ! k_jules
+
+!       IF(zdti_max.le.zdti_bnd) THEN
+!         liter(ji) = .FALSE.
+!         itot = itot - 1
+!       ENDIF
+
+        ENDIF ! liter
+        END DO ! ji
+        IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
       END DO  ! End of the do while iterative procedure
-      
-      IF( ln_icectl .AND. lwp ) THEN
-         WRITE(numout,*) ' zdti_max : ', zdti_max
-         WRITE(numout,*) ' iconv    : ', iconv
-      ENDIF
+!need to move this part where all processors are available      
+!     IF( ln_icectl .AND. lwp ) THEN
+!        WRITE(numout,*) ' zdti_max : ', zdti_max
+!        WRITE(numout,*) ' iconv    : ', iconv
+!     ENDIF
       
       !
@@ -729,11 +683,17 @@
       !
       ! --- calculate conduction fluxes (positive downward)
+      DO jk = 1, nlay_s
+         t_s_1d(1:npti, jk) = tt_s_1d(jk, 1:npti)
+      ENDDO
+      DO jk = 1, nlay_i
+         t_i_1d(1:npti, jk) = tt_i_1d(jk, 1:npti)
+      ENDDO
 
       DO ji = 1, npti
          !                                ! surface ice conduction flux
-         qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
-            &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
+         qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(0, ji)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  &
+            &                  - ( 1._wp - isnow(ji) ) * zkappa_i(0, ji)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) )
          !                                ! bottom ice conduction flux
-         qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
+         qcn_ice_bot_1d(ji) =                          - zkappa_i(nlay_i, ji) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) )
       END DO
       
@@ -770,8 +730,8 @@
                
                IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
+                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji)  &
                      &       + zdq * r1_rdtice ) * a_i_1d(ji)
                ELSE                          ! case T_su = 0degC
-                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
+                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji)  &
                      &       + zdq * r1_rdtice ) * a_i_1d(ji)
                ENDIF
@@ -779,5 +739,5 @@
             ELSEIF( k_jules == np_jules_ACTIVE ) THEN
             
-               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  &
+               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji)  &
                   &          + zdq * r1_rdtice ) * a_i_1d(ji)
             
@@ -800,10 +760,10 @@
       DO ji = 1, npti
          IF( h_s_1d(ji) > 0.1_wp ) THEN 
-            cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0)
+            cnd_ice_1d(ji) = 2._wp * zkappa_s(0, ji)
          ELSE
             IF( h_i_1d(ji) > 0.1_wp ) THEN
-               cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0)
+               cnd_ice_1d(ji) = 2._wp * zkappa_i(0, ji)
             ELSE
-               cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp
+               cnd_ice_1d(ji) = 2._wp * ztcond_i(0, ji) * 10._wp
             ENDIF
          ENDIF
@@ -813,6 +773,10 @@
       IF( k_jules == np_jules_EMULE ) THEN
          ! Restore temperatures to their initial values
-         t_s_1d    (1:npti,:) = ztsold        (1:npti,:)
-         t_i_1d    (1:npti,:) = ztiold        (1:npti,:)
+         DO jk = 1, nlay_s
+            t_s_1d    (1:npti,jk) = ztsold        (jk, 1:npti)
+         ENDDO
+         DO jk = 1, nlay_i
+            t_i_1d    (1:npti,jk) = ztiold        (jk, 1:npti)
+         ENDDO
          qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti)
       ENDIF
@@ -822,8 +786,8 @@
       DO ji = 1, npti         
          !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
+         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(1, ji) * zh_s(ji)
          IF( h_s_1d(ji) >= zhs_min ) THEN
             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
-               &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
+               &            ztcond_i(1, ji) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac )
          ELSE
             t_si_1d(ji) = t_su_1d(ji)
