Changeset 10192


Ignore:
Timestamp:
2018-10-15T12:22:07+02:00 (20 months ago)
Author:
andmirek
Message:

Ticket #2139. Changes to the solver reproducing reference solution

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2018/dev_r10190_ice_thd_zdf/src/ICE/icethd_zdf_bl99.F90

    r10069 r10192  
    108108      ! 
    109109      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    110       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    111       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
    112       REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
    113       REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
    114       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity 
    115       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
    116       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice 
    117       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice 
    118       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice 
     110      REAL(wp), DIMENSION(nlay_i, jpij)     ::   ztiold      ! Old temperature in the ice 
     111      REAL(wp), DIMENSION(nlay_s, jpij)     ::   ztsold      ! Old temperature in the snow 
     112      REAL(wp), DIMENSION(nlay_i, jpij)     ::   ztib        ! Temporary temperature in the ice to check the convergence 
     113      REAL(wp), DIMENSION(nlay_s, jpij)     ::   ztsb        ! Temporary temperature in the snow to check the convergence 
     114      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   ztcond_i    ! Ice thermal conductivity 
     115      REAL(wp), DIMENSION(jpij, 0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice 
     116      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   zradab_i    ! Radiation absorbed in the ice 
     117      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   zkappa_i    ! Kappa factor in the ice 
     118      REAL(wp), DIMENSION(0:nlay_i, jpij)   ::   zeta_i      ! Eta factor in the ice 
    119119      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow 
    120       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow 
    121       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow 
    122       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow 
    123       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term 
    124       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term 
    125       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term 
    126       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms 
     120      REAL(wp), DIMENSION(0:nlay_s, jpij)   ::   zradab_s    ! Radiation absorbed in the snow 
     121      REAL(wp), DIMENSION(0:nlay_s, jpij)   ::   zkappa_s    ! Kappa factor in the snow 
     122      REAL(wp), DIMENSION(0:nlay_s, jpij)   ::   zeta_s      ! Eta factor in the snow 
     123      REAL(wp), DIMENSION(nlay_i+3, jpij)   ::   zindterm    ! 'Ind'ependent term 
     124      REAL(wp), DIMENSION(nlay_i+3, jpij)   ::   zindtbis    ! Temporary 'ind'ependent term 
     125      REAL(wp), DIMENSION(nlay_i+3, jpij)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     126      REAL(wp), DIMENSION(3, nlay_i+3, jpij)::   ztrid       ! Tridiagonal system terms 
    127127      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    128128      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
     129      REAL(wp), DIMENSION(nlay_i, jpij)     ::   tt_i_1d, tsz_i_1d 
     130      REAL(wp), DIMENSION(nlay_s, jpij)     ::   tt_s_1d 
    129131      ! 
    130132      ! Mono-category 
     
    132134      REAL(wp) ::   zhe        ! dummy factor 
    133135      REAL(wp) ::   zcnd_i     ! mean sea ice thermal conductivity 
     136      INTEGER  ::   itot 
     137      LOGICAL, DIMENSION(jpij)  ::   liter 
    134138      !!------------------------------------------------------------------      
    135  
     139      itot = npti 
    136140      ! --- diag error on heat diffusion - PART 1 --- ! 
    137141      DO ji = 1, npti 
     
    171175      ENDIF 
    172176      ! 
    173       ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature 
    174       ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature 
     177      DO jk = 1, nlay_i 
     178         DO ji = 1, npti 
     179            tt_i_1d(jk, ji) = t_i_1d(ji, jk) 
     180            tsz_i_1d(jk, ji) = sz_i_1d(ji, jk) 
     181         END DO 
     182      END DO 
     183 
     184      DO jk = 1, nlay_s 
     185         DO ji = 1, npti 
     186            tt_s_1d(jk, ji) = t_s_1d(ji, jk) 
     187         END DO 
     188      END DO 
     189 
     190 
     191      ztsold (:, :) = tt_s_1d(:,:)   ! Old snow temperature 
     192      ztiold (:, :) = tt_i_1d(:,:)   ! Old ice temperature 
    175193 
    176194      !------------- 
     
    184202            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * h_s_1d(ji) * r1_nlay_s * REAL(jk) ) 
    185203            !                             ! radiation absorbed by the layer-th snow layer 
    186             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
     204            zradab_s(jk,ji) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 
    187205         END DO 
    188206      END DO 
    189207      ! 
    190       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) ) 
     208      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) ) 
    191209      DO jk = 1, nlay_i  
    192210         DO ji = 1, npti 
    193211            !                             ! radiation transmitted below the layer-th ice layer 
    194             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 
     212            zradtr_i(ji, jk) = zradtr_i(ji, 0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) ) 
    195213            !                             ! radiation absorbed by the layer-th ice layer 
    196             zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
     214            zradab_i(jk,ji) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 
    197215         END DO 
    198216      END DO 
    199217      ! 
    200       qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice 
     218      qtr_ice_bot_1d(1:npti) = zradtr_i(1:npti, nlay_i)   ! record radiation transmitted below the ice 
    201219      ! 
    202220      iconv    = 0          ! number of iterations 
    203221      zdti_max = 1000._wp   ! maximal value of error on all points 
    204222      ! 
    205       !                                                          !============================! 
    206       DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins ! 
    207          !                                                       !============================! 
     223      liter(:) = .TRUE.                                          !============================! 
     224!     DO iconv = 1, iconv_max                                    ! Iterative procedure begins ! 
     225      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max) 
     226         zdti_max = 0._wp 
     227!        IF(itot<1) EXIT 
    208228         iconv = iconv + 1 
    209          ! 
    210          ztib(1:npti,:) = t_i_1d(1:npti,:) 
    211          ztsb(1:npti,:) = t_s_1d(1:npti,:) 
    212          ! 
    213          !-------------------------------- 
    214          ! 3) Sea ice thermal conductivity 
    215          !-------------------------------- 
    216          IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T 
    217             ! 
    218             DO ji = 1, npti 
    219                ztcond_i(ji,0)      = rcnd_i + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 
    220                ztcond_i(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
    221             END DO 
    222             DO jk = 1, nlay_i-1 
    223                DO ji = 1, npti 
    224                   ztcond_i(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /  & 
    225                      &                       MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 
     229         DO ji = 1, npti 
     230          IF(liter(ji)) THEN 
     231!            zdti_max = 0._wp 
     232            ! 
     233            ztib(:, ji) = tt_i_1d(:, ji) 
     234            ztsb(:, ji) = tt_s_1d(:, ji) 
     235            ! 
     236            !-------------------------------- 
     237            ! 3) Sea ice thermal conductivity 
     238            !-------------------------------- 
     239            IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T 
     240            ! 
     241               ztcond_i(0, ji)      = rcnd_i + zbeta * tsz_i_1d(1, ji)      / MIN( -epsi10, tt_i_1d(1, ji) - rt0 ) 
     242               ztcond_i(nlay_i, ji) = rcnd_i + zbeta * tsz_i_1d(nlay_i, ji) / MIN( -epsi10, t_bo_1d(ji)  - rt0 ) 
     243               DO jk = 1, nlay_i-1 
     244                    ztcond_i(jk, ji) = rcnd_i + zbeta * 0.5_wp * ( tsz_i_1d(jk, ji) + tsz_i_1d(jk+1, ji) ) /  & 
     245                        &                       MIN( -epsi10, 0.5_wp * (tt_i_1d(jk, ji) + tt_i_1d(jk+1, ji)) - rt0 ) 
     246                  END DO 
     247               ! 
     248            ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 
     249               ! 
     250                  ztcond_i(0, ji)      = rcnd_i + 0.09_wp  *  tsz_i_1d(1, ji)      / MIN( -epsi10, tt_i_1d(1, ji) - rt0 )  & 
     251                     &                         - 0.011_wp * ( tt_i_1d(1, ji) - rt0 ) 
     252                  ztcond_i(nlay_i, ji) = rcnd_i + 0.09_wp  *  tsz_i_1d(nlay_i, ji) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
     253                     &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
     254               DO jk = 1, nlay_i-1 
     255                  ztcond_i(jk,ji) = rcnd_i + 0.09_wp  *   0.5_wp * ( tsz_i_1d(jk,ji) + tsz_i_1d(jk+1,ji) ) /        & 
     256                     &                      MIN( -epsi10, 0.5_wp * ( tt_i_1d (jk,ji) + tt_i_1d (jk+1,ji) ) - rt0 )  & 
     257                     &                     - 0.011_wp * ( 0.5_wp * ( tt_i_1d (jk,ji) + tt_i_1d (jk+1,ji) ) - rt0 ) 
    226258               END DO 
    227             END DO 
    228             ! 
    229          ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 
    230             ! 
    231             DO ji = 1, npti 
    232                ztcond_i(ji,0)      = rcnd_i + 0.09_wp  *  sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  & 
    233                   &                         - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 
    234                ztcond_i(ji,nlay_i) = rcnd_i + 0.09_wp  *  sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )  & 
    235                   &                         - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 
    236             END DO 
    237             DO jk = 1, nlay_i-1 
    238                DO ji = 1, npti 
    239                   ztcond_i(ji,jk) = rcnd_i + 0.09_wp  *   0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) /        & 
    240                      &                      MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )  & 
    241                      &                     - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 
    242                END DO 
    243             END DO 
    244             ! 
    245          ENDIF 
    246          ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )         
    247          ! 
    248          !--- G(he) : enhancement of thermal conductivity in mono-category case 
    249          ! Computation of effective thermal conductivity G(h) 
    250          ! Used in mono-category case only to simulate an ITD implicitly 
    251          ! Fichefet and Morales Maqueda, JGR 1997 
    252          zghe(1:npti) = 1._wp 
    253          ! 
    254          SELECT CASE ( nn_virtual_itd ) 
    255          ! 
    256          CASE ( 1 , 2 ) 
     259            ! 
     260            ENDIF 
     261            ztcond_i(:, ji) = MAX( zkimin, ztcond_i(:, ji) )         
     262            ! 
     263            !--- G(he) : enhancement of thermal conductivity in mono-category case 
     264            ! Computation of effective thermal conductivity G(h) 
     265            ! Used in mono-category case only to simulate an ITD implicitly 
     266            ! Fichefet and Morales Maqueda, JGR 1997 
     267            zghe(ji) = 1._wp 
     268            ! 
     269            SELECT CASE ( nn_virtual_itd ) 
     270            ! 
     271            CASE ( 1 , 2 ) 
    257272            ! 
    258273            zepsilon = 0.1_wp 
    259             DO ji = 1, npti 
    260                zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity 
     274               zcnd_i = SUM( ztcond_i(:,ji) ) / REAL( nlay_i+1, wp )                                ! Mean sea ice thermal conductivity 
    261275               zhe = ( rn_cnd_s * h_i_1d(ji) + zcnd_i * h_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )        ! Effective thickness he (zhe) 
    262276               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) )  & 
    263277                  &   zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )   ! G(he) 
     278            ! 
     279            END SELECT 
     280            ! 
     281            !----------------- 
     282            ! 4) kappa factors 
     283            !----------------- 
     284            !--- Snow 
     285            DO jk = 0, nlay_s-1 
     286               zkappa_s(jk,ji) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
    264287            END DO 
    265             ! 
    266          END SELECT 
    267          ! 
    268          !----------------- 
    269          ! 4) kappa factors 
    270          !----------------- 
    271          !--- Snow 
    272          DO jk = 0, nlay_s-1 
    273             DO ji = 1, npti 
    274                zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 
     288            ! Snow-ice interface 
     289            zfac = 0.5_wp * ( ztcond_i(0, ji) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
     290            IF( zfac > epsi10 ) THEN 
     291               zkappa_s(nlay_s, ji) = zghe(ji) * rn_cnd_s * ztcond_i(0, ji) / zfac 
     292            ELSE 
     293               zkappa_s(nlay_s, ji) = 0._wp 
     294            ENDIF 
     295 
     296            !--- Ice 
     297            DO jk = 0, nlay_i 
     298               zkappa_i(jk,ji) = zghe(ji) * ztcond_i(jk,ji) * z1_h_i(ji) 
    275299            END DO 
    276          END DO 
    277          DO ji = 1, npti   ! Snow-ice interface 
    278             zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 
    279             IF( zfac > epsi10 ) THEN 
    280                zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 
    281             ELSE 
    282                zkappa_s(ji,nlay_s) = 0._wp 
    283             ENDIF 
    284          END DO 
    285  
    286          !--- Ice 
    287          DO jk = 0, nlay_i 
    288             DO ji = 1, npti 
    289                zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 
     300            zkappa_i(0, ji) = zkappa_s(nlay_s, ji) * isnow(ji) + zkappa_i(0, ji) * ( 1._wp - isnow(ji) ) 
     301            ! 
     302            !-------------------------------------- 
     303            ! 5) Sea ice specific heat, eta factors 
     304            !-------------------------------------- 
     305            DO jk = 1, nlay_i 
     306               zcpi = rcpi + zgamma * tsz_i_1d(jk,ji) / MAX( ( tt_i_1d(jk,ji) - rt0 ) * ( ztiold(jk,ji) - rt0 ), epsi10 ) 
     307               zeta_i(jk,ji) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi ) 
    290308            END DO 
    291          END DO 
    292          DO ji = 1, npti   ! Snow-ice interface 
    293             zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 
    294          END DO 
    295          ! 
    296          !-------------------------------------- 
    297          ! 5) Sea ice specific heat, eta factors 
    298          !-------------------------------------- 
    299          DO jk = 1, nlay_i 
    300             DO ji = 1, npti 
    301                zcpi = rcpi + zgamma * sz_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 ) 
    302                zeta_i(ji,jk) = rdt_ice * r1_rhoi * z1_h_i(ji) / MAX( epsi10, zcpi )  
     309 
     310            DO jk = 1, nlay_s 
     311               zeta_s(jk,ji) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 
    303312            END DO 
    304          END DO 
    305  
    306          DO jk = 1, nlay_s 
    307             DO ji = 1, npti 
    308                zeta_s(ji,jk) = rdt_ice * r1_rhos * r1_rcpi * z1_h_s(ji) 
    309             END DO 
    310          END DO 
    311313         ! 
    312314         !----------------------------------------! 
     
    325327            !---------------------------- 
    326328            ! update of the non solar flux according to the update in T_su 
    327             DO ji = 1, npti 
    328329               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    329             END DO 
    330  
    331             DO ji = 1, npti 
     330 
    332331               zfnet(ji) = qsr_ice_1d(ji) - qtr_ice_top_1d(ji) + qns_ice_1d(ji) ! net heat flux = net - transmitted solar + non solar 
    333             END DO 
    334332            ! 
    335333            !---------------------------- 
     
    342340 
    343341            ! ice interior terms (top equation has the same form as the others) 
    344             ztrid   (1:npti,:,:) = 0._wp 
    345             zindterm(1:npti,:)   = 0._wp 
    346             zindtbis(1:npti,:)   = 0._wp 
    347             zdiagbis(1:npti,:)   = 0._wp 
     342            ztrid   (:, :, ji) = 0._wp 
     343            zindterm(:, ji)   = 0._wp 
     344            zindtbis(:, ji)   = 0._wp 
     345            zdiagbis(:, ji)   = 0._wp 
    348346 
    349347            DO jm = nlay_s + 2, nlay_s + nlay_i  
    350                DO ji = 1, npti 
    351348                  jk = jm - nlay_s - 1 
    352                   ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1) 
    353                   ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
    354                   ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk) 
    355                   zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    356                END DO 
     349                  ztrid   (1, jm, ji) =       - zeta_i(jk,ji) *   zkappa_i(jk-1, ji) 
     350                  ztrid   (2, jm, ji) = 1._wp + zeta_i(jk,ji) * ( zkappa_i(jk-1, ji) + zkappa_i(jk, ji) ) 
     351                  ztrid   (3, jm, ji) =       - zeta_i(jk,ji) *                       zkappa_i(jk, ji) 
     352                  zindterm(jm,ji)   = ztiold(jk,ji) + zeta_i(jk,ji) * zradab_i(jk, ji) 
    357353            END DO 
    358354 
    359355            jm =  nlay_s + nlay_i + 1 
    360             DO ji = 1, npti 
    361356               ! ice bottom term 
    362                ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)    
    363                ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 
    364                ztrid   (ji,jm,3) = 0._wp 
    365                zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    366                   &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )  
    367             END DO 
    368  
    369             DO ji = 1, npti 
     357               ztrid   (1, jm, ji) =       - zeta_i(nlay_i, ji) *   zkappa_i(nlay_i-1, ji)    
     358               ztrid   (2, jm, ji) = 1._wp + zeta_i(nlay_i, ji) * ( zkappa_i(nlay_i-1, ji) + zkappa_i(nlay_i, ji) * zg1 ) 
     359               ztrid   (3, jm, ji) = 0._wp 
     360               zindterm(jm,ji)   = ztiold(nlay_i, ji) + zeta_i(nlay_i, ji) *  & 
     361                  &              ( zradab_i(nlay_i, ji) + zkappa_i(nlay_i, ji) * zg1 * t_bo_1d(ji) )  
     362 
    370363               !                               !---------------------! 
    371364               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells ! 
     
    374367                  DO jm = 3, nlay_s + 1 
    375368                     jk = jm - 1 
    376                      ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1) 
    377                      ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    378                      ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk) 
    379                      zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     369                     ztrid   (1, jm, ji) =       - zeta_s(jk,ji) *   zkappa_s(jk-1,ji) 
     370                     ztrid   (2, jm, ji) = 1._wp + zeta_s(jk,ji) * ( zkappa_s(jk-1,ji) + zkappa_s(jk,ji) ) 
     371                     ztrid   (3, jm, ji) =       - zeta_s(jk,ji) *                       zkappa_s(jk,ji) 
     372                     zindterm(jm,ji)   = ztsold(jk,ji) + zeta_s(jk,ji) * zradab_s(jk,ji) 
    380373                  END DO 
    381374                   
    382375                  ! case of only one layer in the ice (ice equation is altered) 
    383376                  IF( nlay_i == 1 ) THEN 
    384                      ztrid   (ji,nlay_s+2,3) = 0._wp 
    385                      zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)  
     377                     ztrid   (3, nlay_s+2, ji) = 0._wp 
     378                     zindterm(nlay_s+2, ji)   = zindterm(nlay_s+2, ji) + zeta_i(1, ji) * zkappa_i(1, ji) * t_bo_1d(ji)  
    386379                  ENDIF 
    387380                   
     
    392385                      
    393386                     ! surface equation 
    394                      ztrid   (ji,1,1) = 0._wp 
    395                      ztrid   (ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    396                      ztrid   (ji,1,3) =                   zg1s * zkappa_s(ji,0) 
    397                      zindterm(ji,1)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
     387                     ztrid   (1, 1, ji) = 0._wp 
     388                     ztrid   (2, 1, ji) = zdqns_ice_b(ji) - zg1s * zkappa_s(0, ji) 
     389                     ztrid   (3, 1, ji) =                   zg1s * zkappa_s(0, ji) 
     390                     zindterm(1, ji)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
    398391                      
    399392                     ! first layer of snow equation 
    400                      ztrid   (ji,2,1) =       - zeta_s(ji,1) *                    zkappa_s(ji,0) * zg1s 
    401                      ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    402                      ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1) 
    403                      zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
     393                     ztrid   (1, 2, ji) =       - zeta_s(1, ji) *                    zkappa_s(0, ji) * zg1s 
     394                     ztrid   (2, 2, ji) = 1._wp + zeta_s(1, ji) * ( zkappa_s(1, ji) + zkappa_s(0, ji) * zg1s ) 
     395                     ztrid   (3, 2, ji) =       - zeta_s(1, ji) *   zkappa_s(1, ji) 
     396                     zindterm(2, ji)   = ztsold(1, ji) + zeta_s(1, ji) * zradab_s(1, ji) 
    404397                      
    405398                  ELSE                            !--  case 2 : surface is melting 
     
    409402                      
    410403                     ! first layer of snow equation 
    411                      ztrid   (ji,2,1) = 0._wp 
    412                      ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    413                      ztrid   (ji,2,3) =       - zeta_s(ji,1) *   zkappa_s(ji,1)  
    414                      zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     404                     ztrid   (1, 2, ji) = 0._wp 
     405                     ztrid   (2, 2, ji) = 1._wp + zeta_s(1, ji) * ( zkappa_s(1, ji) + zkappa_s(0, ji) * zg1s ) 
     406                     ztrid   (3, 2, ji) =       - zeta_s(1, ji) *   zkappa_s(1, ji)  
     407                     zindterm(2, ji)   = ztsold(1, ji) + zeta_s(1, ji) * ( zradab_s(1, ji) + zkappa_s(0, ji) * zg1s * t_su_1d(ji) )  
    415408                  ENDIF 
    416409                  !                            !---------------------! 
     
    424417                      
    425418                     ! surface equation    
    426                      ztrid   (ji,jm_min(ji),1) = 0._wp 
    427                      ztrid   (ji,jm_min(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * zg1     
    428                      ztrid   (ji,jm_min(ji),3) =                   zkappa_i(ji,0) * zg1 
    429                      zindterm(ji,jm_min(ji))   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
     419                     ztrid   (1, jm_min(ji), ji) = 0._wp 
     420                     ztrid   (2, jm_min(ji), ji) = zdqns_ice_b(ji) - zkappa_i(0, ji) * zg1     
     421                     ztrid   (3, jm_min(ji), ji) =                   zkappa_i(0, ji) * zg1 
     422                     zindterm(jm_min(ji), ji)   = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
    430423                      
    431424                     ! first layer of ice equation 
    432                      ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *                    zkappa_i(ji,0) * zg1 
    433                      ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
    434                      ztrid   (ji,jm_min(ji)+1,3) =       - zeta_i(ji,1) *   zkappa_i(ji,1 
    435                      zindterm(ji,jm_min(ji)+1)   = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1 
     425                     ztrid   (1, jm_min(ji)+1, ji) =       - zeta_i(1, ji) *                    zkappa_i(0, ji) * zg1 
     426                     ztrid   (2, jm_min(ji)+1, ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(1, ji) + zkappa_i(0, ji) * zg1 ) 
     427                     ztrid   (3, jm_min(ji)+1, ji) =       - zeta_i(1, ji) *   zkappa_i(1, ji 
     428                     zindterm(jm_min(ji)+1, ji)   = ztiold(1, ji) + zeta_i(1, ji) * zradab_i(1, ji 
    436429                      
    437430                     ! case of only one layer in the ice (surface & ice equations are altered) 
    438431                     IF( nlay_i == 1 ) THEN 
    439                         ztrid   (ji,jm_min(ji),1)   = 0._wp 
    440                         ztrid   (ji,jm_min(ji),2)   = zdqns_ice_b(ji)      -   zkappa_i(ji,0) * 2._wp 
    441                         ztrid   (ji,jm_min(ji),3)   =                          zkappa_i(ji,0) * 2._wp 
    442                         ztrid   (ji,jm_min(ji)+1,1) =       - zeta_i(ji,1) *   zkappa_i(ji,0) * 2._wp 
    443                         ztrid   (ji,jm_min(ji)+1,2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) ) 
    444                         ztrid   (ji,jm_min(ji)+1,3) = 0._wp 
    445                         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)) 
     432                        ztrid   (1, jm_min(ji), ji)   = 0._wp 
     433                        ztrid   (2, jm_min(ji), ji)   = zdqns_ice_b(ji)      -   zkappa_i(0, ji) * 2._wp 
     434                        ztrid   (3, jm_min(ji), ji)   =                          zkappa_i(0, ji) * 2._wp 
     435                        ztrid   (1, jm_min(ji)+1, ji) =       - zeta_i(1, ji) *   zkappa_i(0, ji) * 2._wp 
     436                        ztrid   (2, jm_min(ji)+1, ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(0, ji) * 2._wp + zkappa_i(1, ji) ) 
     437                        ztrid   (3, jm_min(ji)+1, ji) = 0._wp 
     438                        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)) 
    446439                     ENDIF 
    447440                      
     
    452445                      
    453446                     ! first layer of ice equation 
    454                      ztrid   (ji,jm_min(ji),1) = 0._wp 
    455                      ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    456                      ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) *   zkappa_i(ji,1) 
    457                      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))  
     447                     ztrid   (1, jm_min(ji), ji) = 0._wp 
     448                     ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(1, ji) + zkappa_i(0, ji) * zg1 )   
     449                     ztrid   (3, jm_min(ji), ji) =       - zeta_i(1, ji) *   zkappa_i(1, ji) 
     450                     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))  
    458451                      
    459452                     ! case of only one layer in the ice (surface & ice equations are altered) 
    460453                     IF( nlay_i == 1 ) THEN 
    461                         ztrid   (ji,jm_min(ji),1) = 0._wp 
    462                         ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2._wp + zkappa_i(ji,1) ) 
    463                         ztrid   (ji,jm_min(ji),3) = 0._wp 
    464                         zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 
    465                            &                      + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2._wp 
     454                        ztrid   (1, jm_min(ji), ji) = 0._wp 
     455                        ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * ( zkappa_i(0, ji) * 2._wp + zkappa_i(1, ji) ) 
     456                        ztrid   (3, jm_min(ji), ji) = 0._wp 
     457                        zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) * ( zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji) ) & 
     458                           &                      + t_su_1d(ji) * zeta_i(1, ji) * zkappa_i(0, ji) * 2._wp 
    466459                     ENDIF 
    467460                      
     
    469462               ENDIF 
    470463               ! 
    471                zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
    472                zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2) 
     464               zindtbis(jm_min(ji), ji) = zindterm(jm_min(ji), ji) 
     465               zdiagbis(jm_min(ji), ji) = ztrid   (2, jm_min(ji), ji) 
    473466               ! 
    474             END DO 
    475467            ! 
    476468            !------------------------------ 
     
    479471            ! Solve the tridiagonal system with Gauss elimination method. 
    480472            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    481             jm_maxt = 0 
    482             jm_mint = nlay_i+5 
    483             DO ji = 1, npti 
    484                jm_mint = MIN(jm_min(ji),jm_mint) 
    485                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    486             END DO 
    487  
    488             DO jk = jm_mint+1, jm_maxt 
    489                DO ji = 1, npti 
    490                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
    491                   zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    492                   zindtbis(ji,jm) = zindterm(ji,jm  ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1  ) / zdiagbis(ji,jm-1) 
     473 
     474 
     475 
     476 
     477 
     478               DO jm = jm_min(ji)+1, jm_max(ji) 
     479                  zdiagbis(jm,ji) = ztrid   (2, jm, ji) - ztrid(1, jm, ji) * ztrid   (3, jm-1, ji) / zdiagbis(jm-1,ji) 
     480                  zindtbis(jm,ji) = zindterm(jm,ji  ) - ztrid(1, jm, ji) * zindtbis(jm-1,ji  ) / zdiagbis(jm-1,ji) 
    493481               END DO 
    494             END DO 
    495482 
    496483            ! ice temperatures 
    497             DO ji = 1, npti 
    498                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    499             END DO 
    500  
    501             DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    502                DO ji = 1, npti 
     484               tt_i_1d(nlay_i, ji) = zindtbis(jm_max(ji), ji) / zdiagbis(jm_max(ji), ji) 
     485 
     486               DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    503487                  jk = jm - nlay_s - 1 
    504                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     488                  tt_i_1d(jk, ji) = ( zindtbis(jm,ji) - ztrid(3, jm,ji) * tt_i_1d(jk+1, ji) ) / zdiagbis(jm,ji) 
    505489               END DO 
    506             END DO 
    507  
    508             DO ji = 1, npti 
     490 
    509491               ! snow temperatures       
    510492               IF( h_s_1d(ji) > 0._wp ) THEN 
    511                   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) 
     493                  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) 
    512494               ENDIF 
    513495               ! surface temperature 
    514496               ztsub(ji) = t_su_1d(ji) 
    515497               IF( t_su_1d(ji) < rt0 ) THEN 
    516                   t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    517                      &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     498                  t_su_1d(ji) = ( zindtbis(jm_min(ji), ji) - ztrid(3,jm_min(ji),ji) *  & 
     499                     &          ( isnow(ji) * tt_s_1d(1, ji) + ( 1._wp - isnow(ji) ) * tt_i_1d(1, ji) ) ) / zdiagbis(jm_min(ji), ji) 
    518500               ENDIF 
    519             END DO 
    520501            ! 
    521502            !-------------------------------------------------------------- 
     
    524505            ! check that nowhere it has started to melt 
    525506            ! zdti_max is a measure of error, it has to be under zdti_bnd 
    526             zdti_max = 0._wp 
    527             DO ji = 1, npti 
     507 
    528508               t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    529509               zdti_max    = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    530             END DO 
    531  
    532             DO jk = 1, nlay_s 
    533                DO ji = 1, npti 
    534                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    535                   zdti_max      = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     510 
     511               DO jk = 1, nlay_s 
     512                  tt_s_1d(jk,ji) = MAX( MIN( tt_s_1d(jk,ji), rt0 ), rt0 - 100._wp ) 
     513                  zdti_max      = MAX( zdti_max, ABS( tt_s_1d(jk,ji) - ztsb(jk,ji) ) ) 
    536514               END DO 
    537             END DO 
    538  
    539             DO jk = 1, nlay_i 
    540                DO ji = 1, npti 
    541                   ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    542                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    543                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     515 
     516               DO jk = 1, nlay_i 
     517                  ztmelts       = -rTmlt * tsz_i_1d(jk, ji) + rt0  
     518                  tt_i_1d(jk, ji) =  MAX( MIN( tt_i_1d(jk, ji), ztmelts ), rt0 - 100._wp ) 
     519                  zdti_max      =  MAX( zdti_max, ABS( tt_i_1d(jk, ji) - ztib(jk,ji) ) ) 
    544520               END DO 
    545             END DO 
    546  
    547             ! Compute spatial maximum over all errors 
    548             ! note that this could be optimized substantially by iterating only the non-converging points 
    549             IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    550521         ! 
    551522         !----------------------------------------! 
     
    568539 
    569540            ! ice interior terms (top equation has the same form as the others) 
    570             ztrid   (1:npti,:,:) = 0._wp 
    571             zindterm(1:npti,:)   = 0._wp 
    572             zindtbis(1:npti,:)   = 0._wp 
    573             zdiagbis(1:npti,:)   = 0._wp 
    574  
    575             DO jm = nlay_s + 2, nlay_s + nlay_i  
    576                DO ji = 1, npti 
     541               ztrid   (:, :, ji) = 0._wp 
     542               zindterm(:, ji)   = 0._wp 
     543               zindtbis(:, ji)   = 0._wp 
     544               zdiagbis(:, ji)   = 0._wp 
     545 
     546               DO jm = nlay_s + 2, nlay_s + nlay_i  
    577547                  jk = jm - nlay_s - 1 
    578                   ztrid   (ji,jm,1) =       - zeta_i(ji,jk) *   zkappa_i(ji,jk-1) 
    579                   ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
    580                   ztrid   (ji,jm,3) =       - zeta_i(ji,jk) *                       zkappa_i(ji,jk) 
    581                   zindterm(ji,jm)   = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    582                END DO 
    583             ENDDO 
    584  
    585             jm =  nlay_s + nlay_i + 1 
    586             DO ji = 1, npti 
     548                  ztrid   (1, jm, ji) =       - zeta_i(jk,ji) *   zkappa_i(jk-1, ji) 
     549                  ztrid   (2, jm, ji) = 1._wp + zeta_i(jk,ji) * ( zkappa_i(jk-1, ji) + zkappa_i(jk,ji) ) 
     550                  ztrid   (3, jm, ji) =       - zeta_i(jk,ji) *                       zkappa_i(jk,ji) 
     551                  zindterm(jm,ji)   = ztiold(jk,ji) + zeta_i(jk,ji) * zradab_i(jk,ji) 
     552               ENDDO 
     553 
     554               jm =  nlay_s + nlay_i + 1 
    587555               ! ice bottom term 
    588                ztrid   (ji,jm,1) =       - zeta_i(ji,nlay_i) *   zkappa_i(ji,nlay_i-1)    
    589                ztrid   (ji,jm,2) = 1._wp + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i-1) + zkappa_i(ji,nlay_i) * zg1 ) 
    590                ztrid   (ji,jm,3) = 0._wp 
    591                zindterm(ji,jm)   = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    592                   &              ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) )  
    593             ENDDO 
    594  
    595             DO ji = 1, npti 
     556               ztrid   (1, jm, ji) =       - zeta_i(nlay_i, ji) *   zkappa_i(nlay_i-1, ji)    
     557               ztrid   (2, jm, ji) = 1._wp + zeta_i(nlay_i, ji) * ( zkappa_i(nlay_i-1, ji) + zkappa_i(nlay_i, ji) * zg1 ) 
     558               ztrid   (3, jm, ji) = 0._wp 
     559               zindterm(jm,ji)   = ztiold(nlay_i, ji) + zeta_i(nlay_i, ji) *  & 
     560                  &              ( zradab_i(nlay_i, ji) + zkappa_i(nlay_i, ji) * zg1 * t_bo_1d(ji) )  
     561 
    596562               !                               !---------------------! 
    597563               IF( h_s_1d(ji) > 0._wp ) THEN   !  snow-covered cells ! 
     
    600566                  DO jm = 3, nlay_s + 1 
    601567                     jk = jm - 1 
    602                      ztrid   (ji,jm,1) =       - zeta_s(ji,jk) *   zkappa_s(ji,jk-1) 
    603                      ztrid   (ji,jm,2) = 1._wp + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    604                      ztrid   (ji,jm,3) =       - zeta_s(ji,jk) *                       zkappa_s(ji,jk) 
    605                      zindterm(ji,jm)   = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     568                     ztrid   (1, jm, ji) =       - zeta_s(jk,ji) *   zkappa_s(jk-1,ji) 
     569                     ztrid   (2, jm, ji) = 1._wp + zeta_s(jk,ji) * ( zkappa_s(jk-1,ji) + zkappa_s(jk,ji) ) 
     570                     ztrid   (3, jm, ji) =       - zeta_s(jk,ji) *                       zkappa_s(jk,ji) 
     571                     zindterm(jm,ji)   = ztsold(jk,ji) + zeta_s(jk,ji) * zradab_s(jk,ji) 
    606572                  END DO 
    607573                   
    608574                  ! case of only one layer in the ice (ice equation is altered) 
    609575                  IF ( nlay_i == 1 ) THEN 
    610                      ztrid   (ji,nlay_s+2,3) = 0._wp 
    611                      zindterm(ji,nlay_s+2)   = zindterm(ji,nlay_s+2) + zeta_i(ji,1) * zkappa_i(ji,1) * t_bo_1d(ji)  
     576                     ztrid   (3, nlay_s+2, ji) = 0._wp 
     577                     zindterm(nlay_s+2, ji)   = zindterm(nlay_s+2, ji) + zeta_i(1, ji) * zkappa_i(1, ji) * t_bo_1d(ji)  
    612578                  ENDIF 
    613579                   
     
    616582                   
    617583                  ! first layer of snow equation 
    618                   ztrid   (ji,2,1) = 0._wp 
    619                   ztrid   (ji,2,2) = 1._wp + zeta_s(ji,1) * zkappa_s(ji,1) 
    620                   ztrid   (ji,2,3) =       - zeta_s(ji,1) * zkappa_s(ji,1)  
    621                   zindterm(ji,2)   = ztsold(ji,1) + zeta_s(ji,1) * ( zradab_s(ji,1) + qcn_ice_1d(ji) )  
     584                  ztrid   (1, 2, ji) = 0._wp 
     585                  ztrid   (2, 2, ji) = 1._wp + zeta_s(1, ji) * zkappa_s(1, ji) 
     586                  ztrid   (3, 2, ji) =       - zeta_s(1, ji) * zkappa_s(1, ji)  
     587                  zindterm(2, ji)   = ztsold(1, ji) + zeta_s(1, ji) * ( zradab_s(1, ji) + qcn_ice_1d(ji) )  
    622588                   
    623589                  !                            !---------------------! 
     
    628594                   
    629595                  ! first layer of ice equation 
    630                   ztrid   (ji,jm_min(ji),1) = 0._wp 
    631                   ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1 
    632                   ztrid   (ji,jm_min(ji),3) =       - zeta_i(ji,1) * zkappa_i(ji,1) 
    633                   zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 
     596                  ztrid   (1, jm_min(ji), ji) = 0._wp 
     597                  ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * zkappa_i(1, ji 
     598                  ztrid   (3, jm_min(ji), ji) =       - zeta_i(1, ji) * zkappa_i(1, ji) 
     599                  zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) * ( zradab_i(1, ji) + qcn_ice_1d(ji) ) 
    634600                   
    635601                  ! case of only one layer in the ice (surface & ice equations are altered) 
    636602                  IF( nlay_i == 1 ) THEN 
    637                      ztrid   (ji,jm_min(ji),1) = 0._wp 
    638                      ztrid   (ji,jm_min(ji),2) = 1._wp + zeta_i(ji,1) * zkappa_i(ji,1) 
    639                      ztrid   (ji,jm_min(ji),3) = 0._wp 
    640                      zindterm(ji,jm_min(ji))   = ztiold(ji,1) + zeta_i(ji,1) *  & 
    641                         &                                     ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) + qcn_ice_1d(ji) ) 
     603                     ztrid   (1, jm_min(ji), ji) = 0._wp 
     604                     ztrid   (2, jm_min(ji), ji) = 1._wp + zeta_i(1, ji) * zkappa_i(1, ji) 
     605                     ztrid   (3, jm_min(ji), ji) = 0._wp 
     606                     zindterm(jm_min(ji), ji)   = ztiold(1, ji) + zeta_i(1, ji) *  & 
     607                        &                                     ( zradab_i(1, ji) + zkappa_i(1, ji) * t_bo_1d(ji) + qcn_ice_1d(ji) ) 
    642608                  ENDIF 
    643609                   
    644610               ENDIF 
    645611               ! 
    646                zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
    647                zdiagbis(ji,jm_min(ji)) = ztrid   (ji,jm_min(ji),2) 
     612               zindtbis(jm_min(ji), ji) = zindterm(jm_min(ji), ji) 
     613               zdiagbis(jm_min(ji), ji) = ztrid   (2, jm_min(ji), ji) 
    648614               ! 
    649             END DO 
    650615            ! 
    651616            !------------------------------ 
     
    654619            ! Solve the tridiagonal system with Gauss elimination method. 
    655620            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    656             jm_maxt = 0 
    657             jm_mint = nlay_i+5 
    658             DO ji = 1, npti 
    659                jm_mint = MIN(jm_min(ji),jm_mint) 
    660                jm_maxt = MAX(jm_max(ji),jm_maxt) 
    661             END DO 
    662              
    663             DO jk = jm_mint+1, jm_maxt 
    664                DO ji = 1, npti 
    665                   jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 
    666                   zdiagbis(ji,jm) = ztrid   (ji,jm,2) - ztrid(ji,jm,1) * ztrid   (ji,jm-1,3) / zdiagbis(ji,jm-1) 
    667                   zindtbis(ji,jm) = zindterm(ji,jm)   - ztrid(ji,jm,1) * zindtbis(ji,jm-1)   / zdiagbis(ji,jm-1) 
    668                END DO 
     621 
     622 
     623 
     624 
     625 
     626            DO jm = jm_min(ji)+1, jm_max(ji) 
     627                  zdiagbis(jm,ji) = ztrid   (2, jm, ji) - ztrid(1, jm, ji) * ztrid   (3, jm-1, ji) / zdiagbis(jm-1,ji) 
     628                  zindtbis(jm,ji) = zindterm(jm,ji)   - ztrid(1, jm,ji) * zindtbis(jm-1,ji)   / zdiagbis(jm-1,ji) 
    669629            END DO 
    670630             
    671631            ! ice temperatures 
    672            DO ji = 1, npti 
    673                t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    674             END DO 
    675  
    676             DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    677                DO ji = 1, npti 
     632               tt_i_1d(nlay_i, ji) = zindtbis(jm_max(ji), ji) / zdiagbis(jm_max(ji), ji) 
     633 
     634               DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    678635                  jk = jm - nlay_s - 1 
    679                   t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     636                  tt_i_1d(jk, ji) = ( zindtbis(jm,ji) - ztrid(3, jm, ji) * tt_i_1d(jk+1, ji) ) / zdiagbis(jm,ji) 
    680637               END DO 
    681             END DO 
    682638             
    683639            ! snow temperatures       
    684             DO ji = 1, npti 
    685640               IF( h_s_1d(ji) > 0._wp ) THEN 
    686                   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) 
     641                  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) 
    687642               ENDIF 
    688             END DO 
    689643            ! 
    690644            !-------------------------------------------------------------- 
     
    693647            ! check that nowhere it has started to melt 
    694648            ! zdti_max is a measure of error, it has to be under zdti_bnd 
    695             zdti_max = 0._wp 
    696  
    697             DO jk = 1, nlay_s 
    698                DO ji = 1, npti 
    699                   t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    700                   zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     649 
     650 
     651               DO jk = 1, nlay_s 
     652                     tt_s_1d(jk,ji) = MAX( MIN( tt_s_1d(jk,ji), rt0 ), rt0 - 100._wp ) 
     653                     zdti_max = MAX( zdti_max, ABS( tt_s_1d(jk,ji) - ztsb(jk,ji) ) ) 
    701654               END DO 
    702             END DO 
    703655             
    704             DO jk = 1, nlay_i 
    705                DO ji = 1, npti 
    706                   ztmelts       = -rTmlt * sz_i_1d(ji,jk) + rt0  
    707                   t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 
    708                   zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     656               DO jk = 1, nlay_i 
     657                     ztmelts       = -rTmlt * tsz_i_1d(jk,ji) + rt0  
     658                     tt_i_1d(jk,ji) =  MAX( MIN( tt_i_1d(jk,ji), ztmelts ), rt0 - 100._wp ) 
     659                     zdti_max      =  MAX( zdti_max, ABS( tt_i_1d(jk,ji) - ztib(jk,ji) ) ) 
    709660               END DO 
    710             END DO 
    711  
    712             ! Compute spatial maximum over all errors 
    713             ! note that this could be optimized substantially by iterating only the non-converging points 
    714             IF( lk_mpp )   CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    715  
    716          ENDIF ! k_jules 
    717           
     661 
     662        ENDIF ! k_jules 
     663 
     664!       IF(zdti_max.le.zdti_bnd) THEN 
     665!         liter(ji) = .FALSE. 
     666!         itot = itot - 1 
     667!       ENDIF 
     668 
     669        ENDIF ! liter 
     670        END DO ! ji 
     671        IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    718672      END DO  ! End of the do while iterative procedure 
    719        
    720       IF( ln_icectl .AND. lwp ) THEN 
    721          WRITE(numout,*) ' zdti_max : ', zdti_max 
    722          WRITE(numout,*) ' iconv    : ', iconv 
    723       ENDIF 
     673!need to move this part where all processors are available       
     674!     IF( ln_icectl .AND. lwp ) THEN 
     675!        WRITE(numout,*) ' zdti_max : ', zdti_max 
     676!        WRITE(numout,*) ' iconv    : ', iconv 
     677!     ENDIF 
    724678       
    725679      ! 
     
    729683      ! 
    730684      ! --- calculate conduction fluxes (positive downward) 
     685      DO jk = 1, nlay_s 
     686         t_s_1d(1:npti, jk) = tt_s_1d(jk, 1:npti) 
     687      ENDDO 
     688      DO jk = 1, nlay_i 
     689         t_i_1d(1:npti, jk) = tt_i_1d(jk, 1:npti) 
     690      ENDDO 
    731691 
    732692      DO ji = 1, npti 
    733693         !                                ! surface ice conduction flux 
    734          qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(ji,0)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
    735             &                  - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
     694         qcn_ice_top_1d(ji) =  -           isnow(ji)   * zkappa_s(0, ji)      * zg1s * ( t_s_1d(ji,1) - t_su_1d(ji) )  & 
     695            &                  - ( 1._wp - isnow(ji) ) * zkappa_i(0, ji)      * zg1  * ( t_i_1d(ji,1) - t_su_1d(ji) ) 
    736696         !                                ! bottom ice conduction flux 
    737          qcn_ice_bot_1d(ji) =                          - zkappa_i(ji,nlay_i) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
     697         qcn_ice_bot_1d(ji) =                          - zkappa_i(nlay_i, ji) * zg1  * ( t_bo_1d(ji ) - t_i_1d (ji,nlay_i) ) 
    738698      END DO 
    739699       
     
    770730                
    771731               IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    772                   zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  & 
     732                  zhfx_err = ( qns_ice_1d(ji)     + qsr_ice_1d(ji)     - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji)  & 
    773733                     &       + zdq * r1_rdtice ) * a_i_1d(ji) 
    774734               ELSE                          ! case T_su = 0degC 
    775                   zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  & 
     735                  zhfx_err = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji)  & 
    776736                     &       + zdq * r1_rdtice ) * a_i_1d(ji) 
    777737               ENDIF 
     
    779739            ELSEIF( k_jules == np_jules_ACTIVE ) THEN 
    780740             
    781                zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji,nlay_i) - qcn_ice_bot_1d(ji)  & 
     741               zhfx_err    = ( qcn_ice_top_1d(ji) + qtr_ice_top_1d(ji) - zradtr_i(ji, nlay_i) - qcn_ice_bot_1d(ji)  & 
    782742                  &          + zdq * r1_rdtice ) * a_i_1d(ji) 
    783743             
     
    800760      DO ji = 1, npti 
    801761         IF( h_s_1d(ji) > 0.1_wp ) THEN  
    802             cnd_ice_1d(ji) = 2._wp * zkappa_s(ji,0) 
     762            cnd_ice_1d(ji) = 2._wp * zkappa_s(0, ji) 
    803763         ELSE 
    804764            IF( h_i_1d(ji) > 0.1_wp ) THEN 
    805                cnd_ice_1d(ji) = 2._wp * zkappa_i(ji,0) 
     765               cnd_ice_1d(ji) = 2._wp * zkappa_i(0, ji) 
    806766            ELSE 
    807                cnd_ice_1d(ji) = 2._wp * ztcond_i(ji,0) * 10._wp 
     767               cnd_ice_1d(ji) = 2._wp * ztcond_i(0, ji) * 10._wp 
    808768            ENDIF 
    809769         ENDIF 
     
    813773      IF( k_jules == np_jules_EMULE ) THEN 
    814774         ! Restore temperatures to their initial values 
    815          t_s_1d    (1:npti,:) = ztsold        (1:npti,:) 
    816          t_i_1d    (1:npti,:) = ztiold        (1:npti,:) 
     775         DO jk = 1, nlay_s 
     776            t_s_1d    (1:npti,jk) = ztsold        (jk, 1:npti) 
     777         ENDDO 
     778         DO jk = 1, nlay_i 
     779            t_i_1d    (1:npti,jk) = ztiold        (jk, 1:npti) 
     780         ENDDO 
    817781         qcn_ice_1d(1:npti)   = qcn_ice_top_1d(1:npti) 
    818782      ENDIF 
     
    822786      DO ji = 1, npti          
    823787         !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    824          zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
     788         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(1, ji) * zh_s(ji) 
    825789         IF( h_s_1d(ji) >= zhs_min ) THEN 
    826790            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
    827                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 
     791               &            ztcond_i(1, ji) * zh_s(ji) * t_i_1d(ji,1) ) / MAX( epsi10, zfac ) 
    828792         ELSE 
    829793            t_si_1d(ji) = t_su_1d(ji) 
Note: See TracChangeset for help on using the changeset viewer.