Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8325)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8342)
@@ -37,5 +37,5 @@
 CONTAINS
 
-   SUBROUTINE lim_thd_dif( kideb , kiut )
+   SUBROUTINE lim_thd_dif
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_thd_dif  ***
@@ -67,7 +67,4 @@
       !!           of temperature
       !!
-      !! ** Arguments :
-      !!           kideb , kiut : Starting and ending points on which the 
-      !!                         the computation is applied
       !!
       !! ** Inputs / Ouputs : (global commons)
@@ -89,6 +86,4 @@
       !!           (04-2007) Energy conservation tested by M. Vancoppenolle
       !!------------------------------------------------------------------
-      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied
-
       !! * Local variables
       INTEGER ::   ji             ! spatial loop index
@@ -180,5 +175,5 @@
       ! --- diag error on heat diffusion - PART 1 --- !
       zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
-      DO ji = kideb, kiut
+      DO ji = 1, nidx
          zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
             &           SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
@@ -188,5 +183,5 @@
       ! 1) Initialization                                                            !
       !------------------------------------------------------------------------------!
-      DO ji = kideb , kiut
+      DO ji = 1 , nidx
          isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
          ! layer thickness
@@ -203,5 +198,5 @@
 
       DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
          END DO
@@ -209,5 +204,5 @@
 
       DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
          END DO
@@ -230,5 +225,5 @@
       ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
       zhsu = 0.1_wp ! threshold for the computation of i0
-      DO ji = kideb , kiut
+      DO ji = 1 , nidx
          ! switches
          isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
@@ -243,5 +238,5 @@
       ! Derivative of the non solar flux
       !-------------------------------------------------------
-      DO ji = kideb , kiut
+      DO ji = 1 , nidx
          zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
          zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
@@ -254,10 +249,10 @@
       !---------------------------------------------------------
 
-      DO ji = kideb, kiut           ! snow initialization
+      DO ji = 1, nidx           ! snow initialization
          zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
       END DO
 
       DO jk = 1, nlay_s          ! Radiation through snow
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             !                             ! radiation transmitted below the layer-th snow layer
             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
@@ -267,10 +262,10 @@
       END DO
 
-      DO ji = kideb, kiut           ! ice initialization
+      DO ji = 1, nidx           ! ice initialization
          zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
       END DO
 
       DO jk = 1, nlay_i          ! Radiation through ice
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             !                             ! radiation transmitted below the layer-th ice layer
             zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
@@ -280,5 +275,5 @@
       END DO
 
-      DO ji = kideb, kiut           ! Radiation transmitted below the ice
+      DO ji = 1, nidx           ! Radiation transmitted below the ice
          ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
       END DO
@@ -288,5 +283,5 @@
       !------------------------------------------------------------------------------|
       !
-      DO ji = kideb, kiut        ! Old surface temperature
+      DO ji = 1, nidx        ! Old surface temperature
          ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
          ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
@@ -296,5 +291,5 @@
 
       DO jk = 1, nlay_s       ! Old snow temperature
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             ztsb(ji,jk) =  t_s_1d(ji,jk)
          END DO
@@ -302,5 +297,5 @@
 
       DO jk = 1, nlay_i       ! Old ice temperature
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             ztib(ji,jk) =  t_i_1d(ji,jk)
          END DO
@@ -319,10 +314,10 @@
          !
          IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
             END DO
             DO jk = 1, nlay_i-1
-               DO ji = kideb , kiut
+               DO ji = 1 , nidx
                   ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
                      MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0)
@@ -333,5 +328,5 @@
 
          IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
                   &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 )  
@@ -339,5 +334,5 @@
             END DO
             DO jk = 1, nlay_i-1
-               DO ji = kideb , kiut
+               DO ji = 1 , nidx
                   ztcond_i(ji,jk) = rcdic +                                                                       & 
                      &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
@@ -347,5 +342,5 @@
                END DO
             END DO
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
                   &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )  
@@ -372,5 +367,5 @@
             zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp
 
-            DO ji = kideb, kiut
+            DO ji = 1, nidx
    
                ! Mean sea ice thermal conductivity
@@ -400,5 +395,5 @@
          !
          !--- Snow
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
             zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac
@@ -407,5 +402,5 @@
 
          DO jk = 1, nlay_s-1
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) )
             END DO
@@ -414,5 +409,5 @@
          !--- Ice
          DO jk = 1, nlay_i-1
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) )
             END DO
@@ -420,5 +415,5 @@
 
          !--- Snow-ice interface
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             zfac                  = 1./ MAX( epsi10 , zh_i(ji) )
             zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
@@ -435,5 +430,5 @@
          !
          DO jk = 1, nlay_i
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztitemp(ji,jk)   = t_i_1d(ji,jk)
                zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
@@ -443,5 +438,5 @@
 
          DO jk = 1, nlay_s
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztstemp(ji,jk) = t_s_1d(ji,jk)
                zeta_s(ji,jk)  = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 )
@@ -455,5 +450,5 @@
          !
          IF ( ln_dqnsice ) THEN 
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ! update of the non solar flux according to the update in T_su
                qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
@@ -462,5 +457,5 @@
 
          ! Update incoming flux
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             ! update incoming flux
             zf(ji)    =          zfsw(ji)              & ! net absorbed solar radiation
@@ -481,5 +476,5 @@
 
          DO numeq=1,nlay_i+3
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztrid(ji,numeq,1) = 0.
                ztrid(ji,numeq,2) = 0.
@@ -492,5 +487,5 @@
 
          DO numeq = nlay_s + 2, nlay_s + nlay_i 
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                jk                 = numeq - nlay_s - 1
                ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
@@ -502,5 +497,5 @@
 
          numeq =  nlay_s + nlay_i + 1
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             !!ice bottom term
             ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
@@ -512,5 +507,5 @@
 
 
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             IF ( ht_s_1d(ji) > 0.0 ) THEN
                !
@@ -659,5 +654,5 @@
          minnumeqmin = nlay_i+5
 
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
             zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
@@ -667,5 +662,5 @@
 
          DO jk = minnumeqmin+1, maxnumeqmax
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
@@ -674,5 +669,5 @@
          END DO
 
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             ! ice temperatures
             t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
@@ -680,5 +675,5 @@
 
          DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                jk    =  numeq - nlay_s - 1
                t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
@@ -686,5 +681,5 @@
          END DO
 
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             ! snow temperatures      
             IF (ht_s_1d(ji) > 0._wp) &
@@ -706,5 +701,5 @@
          ! check that nowhere it has started to melt
          ! zdti(ji) is a measure of error, it has to be under zdti_bnd
-         DO ji = kideb , kiut
+         DO ji = 1 , nidx
             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
             zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
@@ -712,5 +707,5 @@
 
          DO jk  =  1, nlay_s
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  )
                zdti  (ji)    = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
@@ -719,5 +714,5 @@
 
          DO jk  =  1, nlay_i
-            DO ji = kideb , kiut
+            DO ji = 1 , nidx
                ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
                t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
@@ -729,5 +724,5 @@
          ! note that this could be optimized substantially by iterating only the non-converging points
          zdti_max = 0._wp
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             zdti_max = MAX( zdti_max, zdti(ji) )   
          END DO
@@ -738,5 +733,5 @@
       ! MV SIMIP 2016
       !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-      DO ji = kideb, kiut
+      DO ji = 1, nidx
          zfac        = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
          t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + &
@@ -755,5 +750,5 @@
       !   12) Fluxes at the interfaces                                          !
       !-------------------------------------------------------------------------!
-      DO ji = kideb, kiut
+      DO ji = 1, nidx
          !                                ! surface ice conduction flux
          isnow(ji)       = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )
@@ -771,9 +766,9 @@
 
       ! --- computes sea ice energy of melting compulsory for limthd_dh --- !
-      CALL lim_thd_enmelt( kideb, kiut )
+      CALL lim_thd_enmelt
 
       ! --- diagnose the change in non-solar flux due to surface temperature change --- !
       IF ( ln_dqnsice ) THEN
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
          END DO
@@ -781,5 +776,5 @@
 
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
-      DO ji = kideb, kiut
+      DO ji = 1, nidx
          zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
             &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
@@ -798,5 +793,5 @@
       ! Heat flux used to warm/cool ice in W.m-2
       !-----------------------------------------
-      DO ji = kideb, kiut
+      DO ji = 1, nidx
          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
             hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
@@ -821,5 +816,5 @@
    END SUBROUTINE lim_thd_dif
 
-   SUBROUTINE lim_thd_enmelt( kideb, kiut )
+   SUBROUTINE lim_thd_enmelt
       !!-----------------------------------------------------------------------
       !!                   ***  ROUTINE lim_thd_enmelt *** 
@@ -829,6 +824,4 @@
       !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
-      !
       INTEGER  ::   ji, jk   ! dummy loop indices
       REAL(wp) ::   ztmelts  ! local scalar 
@@ -836,5 +829,5 @@
       !
       DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             ztmelts      = - tmut  * s_i_1d(ji,jk) + rt0
             t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point
@@ -846,5 +839,5 @@
       END DO
       DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = kideb, kiut
+         DO ji = 1, nidx
             e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
          END DO
