Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8565)
@@ -68,6 +68,4 @@
       rn_pstar      =   2.0e+04       !     ice strength thickness parameter [N/m2]
       rn_crhg       =   20.0          !     ice strength conc. parameter (-)
-   ln_str_R75       = .false.         !  ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2)    
-      rn_perdg      =   17.0          !     ridging work divided by pot. energy change in ridging
                    ! -- ice_rdgrft -- !
    rn_csrdg         =   0.5           !  fraction of shearing energy contributing to ridging
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8565)
@@ -27,8 +27,6 @@
    !: are the variables corresponding to 2d vectors
 
-!! please, DOCTOR naming convention   starting by i means LOCAL integer
-!!         ===>>>  rename idxice  as nidice or nidthd, or nidx_thd or midx ... ?
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   idxice !: selected points for ice thermo
-   INTEGER , PUBLIC                                  ::   nidx   !  number of selected points
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   nptidx !: selected points for ice thermo
+   INTEGER , PUBLIC                                  ::   npti   !  number of selected points
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlead_1d     
@@ -175,5 +173,5 @@
 
       ii = 1
-      ALLOCATE( idxice   (jpij) ,   &
+      ALLOCATE( nptidx   (jpij) ,   &
          &      qlead_1d (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   &
          &      fr1_i0_1d(jpij) , fr2_i0_1d (jpij) , qns_ice_1d(jpij)  ,   &
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8565)
@@ -29,9 +29,9 @@
    REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
 
-   REAL(wp) , PARAMETER ::   rc1    = 0.05    ! snow thickness (only for nn_ice_alb=0)
-   REAL(wp) , PARAMETER ::   rc2    = 0.10    !  "        "
-   REAL(wp) , PARAMETER ::   rcloud = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
-   REAL(wp) , PARAMETER ::   r1_c1 = 1. / rc1
-   REAL(wp) , PARAMETER ::   r1_c2 = 1. / rc2
+   REAL(wp) , PARAMETER ::   ppc1    = 0.05    ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   ppc2    = 0.10    !  "        "
+   REAL(wp) , PARAMETER ::   ppcloud = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
+   REAL(wp) , PARAMETER ::   pp1_c1 = 1. / ppc1
+   REAL(wp) , PARAMETER ::   pp1_c2 = 1. / ppc2
    !
    ! ** albedo namelist (namalb)
@@ -154,14 +154,14 @@
                DO ji = 1, jpi
                   !                    ! Freezing snow
-                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > rc2
-                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - rc1 ) ) )
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > ppc2
+                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - ppc1 ) ) )
                   zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
-                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1  )   &
+                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * pp1_c1  )   &
                      &    +           zswitch   * rn_alb_sdry  
                      !
                   !                    ! Melting snow
-                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > rc2
-                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - rc2 ) )
-                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 )   &
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > ppc2
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - ppc2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * pp1_c2 )   &
                      &      +         zswitch   *   rn_alb_smlt 
                      !
@@ -178,5 +178,5 @@
          END DO
          !
-         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + ppcloud       ! Oberhuber correction for overcast sky
          !
          !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8565)
@@ -71,5 +71,5 @@
       REAL(wp) ::   zvtrp, zetrp
       REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill
-      REAL(wp), PARAMETER ::   zconv = 1.e-9 ! convert W to GW and kg to Mt
+      REAL(wp), PARAMETER ::   ppconv = 1.e-9 ! convert W to GW and kg to Mt
       !!-------------------------------------------------------------------
       !
@@ -79,22 +79,22 @@
             &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
             &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    &
-            &                  ) * e1e2t(:,:) ) * zconv
+            &                  ) * e1e2t(:,:) ) * ppconv
          !
          !                          ! salt flux
          pdiag_fs = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
             &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
-            &                  ) *  e1e2t(:,:) ) * zconv 
+            &                  ) *  e1e2t(:,:) ) * ppconv 
          !
          !                          ! heat flux
          pdiag_ft = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
             &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                  ) *  e1e2t(:,:) ) * zconv
-
-         pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )
-
-         pdiag_s = glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t * zconv )
+            &                  ) *  e1e2t(:,:) ) * ppconv
+
+         pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * ppconv )
+
+         pdiag_s = glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t * ppconv )
 
          pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
-            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv
+            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * ppconv
 
       ELSEIF( icount == 1 ) THEN
@@ -104,30 +104,30 @@
             &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv
+            &              ) * e1e2t(:,:) ) * ppconv - pdiag_fv
 
          ! salt flux
          zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
             &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_fs
+            &              ) * e1e2t(:,:) ) * ppconv - pdiag_fs
 
          ! heat flux
          zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
             &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_ft
+            &              ) * e1e2t(:,:) ) * ppconv - pdiag_ft
  
          ! outputs
-         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  &
+         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * ppconv  &
             &     - pdiag_v ) * r1_rdtice - zfv ) * rday
 
-         zs = ( ( glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  &
+         zs = ( ( glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t ) * ppconv  &
             &     - pdiag_s ) * r1_rdtice + zfs ) * rday
 
          zt = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   &
-            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   &
+            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * ppconv   &
             &   - pdiag_t ) * r1_rdtice + zft
 
          ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
-         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t  ) * zconv * rday 
-         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t  ) * zconv
+         zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t  ) * ppconv * rday 
+         zetrp = glob_sum( ( diag_trp_ei         + diag_trp_es         ) * e1e2t  ) * ppconv
 
          zvmin = glob_min( v_i )
@@ -136,5 +136,5 @@
 
          ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
-         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2
+         zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * ppconv ! in 1.e9 m2
          zv_sill = zarea * 2.5e-5
          zs_sill = zarea * 25.e-5
@@ -176,20 +176,20 @@
       REAL(wp)                        :: zhfx, zsfx, zvfx
       REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill
-      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
+      REAL(wp), PARAMETER             :: ppconv = 1.e-9 ! convert W to GW and kg to Mt
       !!-------------------------------------------------------------------
 
       ! water flux
-      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday
+      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * ppconv * rday
 
       ! salt flux
-      zsfx  = glob_sum( ( sfx + diag_sice ) * e1e2t ) * zconv * rday
+      zsfx  = glob_sum( ( sfx + diag_sice ) * e1e2t ) * ppconv * rday
 
       ! heat flux
       zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es   &
       !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
-         &              ) * e1e2t ) * zconv
+         &              ) * e1e2t ) * ppconv
 
       ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
-      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2
+      zarea   = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * ppconv ! in 1.e9 m2
       zv_sill = zarea * 2.5e-5
       zs_sill = zarea * 25.e-5
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8565)
@@ -61,6 +61,4 @@
    REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
    REAL(wp) ::   rn_crhg          ! determines changes in ice strength
-   LOGICAL  ::   ln_str_R75       ! ice strength parameterization (Rothrock75)
-   REAL(wp) ::   rn_perdg         ! ridging work divided by pot. energy change in ridging
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
@@ -130,12 +128,12 @@
       !!
       INTEGER  ::   ji, jj, jk, jl             ! dummy loop index
-      INTEGER  ::   niter, iterate_ridging     ! local integer 
-      INTEGER  ::   nidx2                      ! local integer
+      INTEGER  ::   iter, iterate_ridging      ! local integer 
+      INTEGER  ::   ipti                       ! local integer
       REAL(wp) ::   zfac                       ! local scalar
-      INTEGER , DIMENSION(jpij) ::   idxice2       ! compute ridge/raft or not
+      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
       REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       !
-      INTEGER, PARAMETER ::   nitermax = 20    
+      INTEGER, PARAMETER ::   jp_itermax = 20    
       !!-------------------------------------------------------------------
       ! controls
@@ -154,25 +152,25 @@
       ! 0) Identify grid cells with ice
       !--------------------------------
-      nidx = 0   ;   idxice(:) = 0
+      npti = 0   ;   nptidx(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN
-               nidx           = nidx + 1
-               idxice( nidx ) = (jj - 1) * jpi + ji
+               npti           = npti + 1
+               nptidx( npti ) = (jj - 1) * jpi + ji
             ENDIF
          END DO
       END DO
 
-      IF( nidx > 0 ) THEN
+      IF( npti > 0 ) THEN
         
          ! just needed here
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) )
          ! needed here and in the iteration loop
-         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i(:,:,:) )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i(:,:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx)      , ato_i(:,:) )
-
-         DO ji = 1, nidx
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i(:,:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i(:,:) )
+
+         DO ji = 1, npti
             !-----------------------------------------------------------------------------!
             ! 2) Dynamical inputs (closing rate, divu_adv, opning)
@@ -219,19 +217,19 @@
          !------------------------------------------------
          CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
-         nidx2 = 0 ; idxice2(:) = 0
-         DO ji = 1, nidx
+         ipti = 0 ; iptidx(:) = 0
+         DO ji = 1, npti
             IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
-               nidx2 = nidx2 + 1
-               idxice2    (nidx2)   = idxice     (ji)
-               a_i_2d     (nidx2,:) = a_i_2d     (ji,:)  ! adjust to new indices
-               v_i_2d     (nidx2,:) = v_i_2d     (ji,:)  ! adjust to new indices
-               ato_i_1d   (nidx2)   = ato_i_1d   (ji)    ! adjust to new indices
-               closing_net(nidx2)   = closing_net(ji)    ! adjust to new indices
-               zdivu_adv  (nidx2)   = zdivu_adv  (ji)    ! adjust to new indices
-               opning     (nidx2)   = opning     (ji)    ! adjust to new indices
+               ipti = ipti + 1
+               iptidx     (ipti)   = nptidx     (ji)
+               a_i_2d     (ipti,:) = a_i_2d     (ji,:)  ! adjust to new indices
+               v_i_2d     (ipti,:) = v_i_2d     (ji,:)  ! adjust to new indices
+               ato_i_1d   (ipti)   = ato_i_1d   (ji)    ! adjust to new indices
+               closing_net(ipti)   = closing_net(ji)    ! adjust to new indices
+               zdivu_adv  (ipti)   = zdivu_adv  (ji)    ! adjust to new indices
+               opning     (ipti)   = opning     (ji)    ! adjust to new indices
             ENDIF
          END DO
-         idxice(:) = idxice2(:)
-         nidx      = nidx2
+         nptidx(:) = iptidx(:)
+         npti      = ipti
 
       ENDIF
@@ -240,5 +238,5 @@
       ! Start ridging
       !--------------
-      IF( nidx > 0 ) THEN
+      IF( npti > 0 ) THEN
          
          !-----------
@@ -246,32 +244,32 @@
          !-----------
          ! fields used but not modified
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) )
          ! the following fields are modified in this routine
-         !!CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )
-         !!CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) )
-         !!CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i(:,:,:) )
+         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
+         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) )
+         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
          IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )
-            CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )
+            CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+            CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
          ENDIF
          DO jl = 1, jpl
             DO jk = 1, nlay_s
-               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )
+               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
             END DO
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
-            END DO
-         END DO
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_dyn_1d    (1:nidx), sfx_dyn    (:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d    (1:nidx), sfx_bri    (:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_dyn_1d    (1:nidx), wfx_dyn    (:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dyn_1d    (1:nidx), hfx_dyn    (:,:) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )
+               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
          IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
          ENDIF
          
@@ -279,7 +277,7 @@
          ! 3) Ridging iteration
          !-----------------------------------------------------------------------------!
-         niter           = 1
+         iter            = 1
          iterate_ridging = 1      
-         DO WHILE( iterate_ridging > 0 .AND. niter < nitermax )
+         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )
 
             CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )
@@ -294,5 +292,5 @@
             ! rates were reduced above), ridge again with new rates.
             iterate_ridging = 0
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) )
                IF( ABS( zfac ) < epsi10 ) THEN
@@ -308,6 +306,6 @@
             END DO
             !
-            niter = niter + 1
-            IF( niter  >  nitermax )    CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )
+            iter = iter + 1
+            IF( iter  >  jp_itermax )    CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )
             !
          END DO
@@ -316,32 +314,32 @@
          ! 1D <==> 2D
          !-----------
-         CALL tab_1d_2d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i(:,:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
          IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )
-            CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )
+            CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+            CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
          ENDIF
          DO jl = 1, jpl
             DO jk = 1, nlay_s
-               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )
+               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) )
             END DO
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
-            END DO
-         END DO
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_dyn_1d    (1:nidx), sfx_dyn    (:,:) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d    (1:nidx), sfx_bri    (:,:) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_dyn_1d    (1:nidx), wfx_dyn    (:,:) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dyn_1d    (1:nidx), hfx_dyn    (:,:) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )
+               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
          IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
          ENDIF
 
-      ENDIF ! nidx > 0
+      ENDIF ! npti > 0
    
       CALL ice_var_agg( 1 ) 
@@ -378,6 +376,6 @@
 
       !                       ! Ice thickness needed for rafting
-      WHERE( pa_i(1:nidx,:) > 0._wp )   ;   zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:)
-      ELSEWHERE                         ;   zhi(1:nidx,:) = 0._wp
+      WHERE( pa_i(1:npti,:) > 0._wp )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+      ELSEWHERE                         ;   zhi(1:npti,:) = 0._wp
       END WHERE
 
@@ -398,8 +396,8 @@
       ! Compute total area of ice plus open water.
       ! This is in general not equal to one because of divergence during transport
-      zasum(1:nidx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 )
-      !
-      WHERE( zasum(1:nidx) > 0._wp )   ;   z1_asum(1:nidx) = 1._wp / zasum(1:nidx)
-      ELSEWHERE                        ;   z1_asum(1:nidx) = 0._wp
+      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
+      !
+      WHERE( zasum(1:npti) > 0._wp )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
+      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp
       END WHERE
       ! Compute cumulative thickness distribution function
@@ -407,13 +405,13 @@
       ! where zGsum(n) is the fractional area in categories 0 to n.
       ! initial value (in h = 0) equals open water area
-      zGsum(1:nidx,-1) = 0._wp
-      zGsum(1:nidx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx)
+      zGsum(1:npti,-1) = 0._wp
+      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti)
       DO jl = 1, jpl
-         zGsum(1:nidx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx)  ! sum(1:jl) is ok (and not jpl)
+         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is ok (and not jpl)
       END DO
       !
       IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
          DO jl = 0, jpl    
-            DO ji = 1, nidx
+            DO ji = 1, npti
                IF    ( zGsum(ji,jl)   < rn_gstar ) THEN
                   apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * &
@@ -432,10 +430,10 @@
          zfac = 1._wp / ( 1._wp - EXP(-z1_astar) )
          DO jl = -1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac
             END DO
          END DO
          DO jl = 0, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl)
             END DO
@@ -447,5 +445,5 @@
       IF( ln_rafting .AND. ln_ridging ) THEN             !- ridging & rafting
          DO jl = 1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl)
                araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl)
@@ -454,5 +452,5 @@
       ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
          DO jl = 1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                aridge(ji,jl) = apartf(ji,jl)
                araft (ji,jl) = 0._wp
@@ -461,5 +459,5 @@
       ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
          DO jl = 1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                aridge(ji,jl) = 0._wp
                araft (ji,jl) = apartf(ji,jl)
@@ -468,5 +466,5 @@
       ELSE                                               !- no ridging & no rafting
          DO jl = 1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                aridge(ji,jl) = 0._wp
                araft (ji,jl) = 0._wp         
@@ -502,8 +500,8 @@
 
       zfac = 1._wp / hi_hrft
-      zaksum(1:nidx) = apartf(1:nidx,0)
+      zaksum(1:npti) = apartf(1:npti,0)
       ! Transfer function
       DO jl = 1, jpl !all categories have a specific transfer function
-         DO ji = 1, nidx
+         DO ji = 1, npti
             IF ( apartf(ji,jl) > 0._wp ) THEN
                zhmean         = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min )
@@ -530,9 +528,9 @@
       !  closing rate to a gross closing rate.  
       ! NOTE: 0 < aksum <= 1
-      WHERE( zaksum(1:nidx) > 0._wp )   ;   closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx)
-      ELSEWHERE                         ;   closing_gross(1:nidx) = 0._wp
+      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti)
+      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp
       END WHERE
       
-      DO ji = 1, nidx  
+      DO ji = 1, npti  
          ! correction to closing rate and opening if closing rate is excessive
          !---------------------------------------------------------------------
@@ -552,5 +550,5 @@
       ! would be removed.  Reduce the opening rate proportionately.
       DO jl = 1, jpl
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice
 !!            IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN
@@ -594,5 +592,5 @@
       ! 1) Compute change in open water area due to closing and opening.
       !-------------------------------------------------------------------------------
-      DO ji = 1, nidx
+      DO ji = 1, npti
          ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice )
       END DO
@@ -603,7 +601,7 @@
       DO jl1 = 1, jpl
 
-         CALL tab_2d_1d( nidx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,jl1) )
-
-         DO ji = 1, nidx
+         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) )
+
+         DO ji = 1, npti
 
             !------------------------------------------------
@@ -716,5 +714,5 @@
          ! special loop for e_i because of layers jk
          DO jk = 1, nlay_i
-            DO ji = 1, nidx
+            DO ji = 1, npti
                IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
                   ! Compute ridging /rafting fractions
@@ -736,5 +734,5 @@
          DO jl2  = 1, jpl 
             !
-            DO ji = 1, nidx
+            DO ji = 1, npti
 
                IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN
@@ -783,5 +781,5 @@
 
             DO jk = 1, nlay_i
-               DO ji = 1, nidx
+               DO ji = 1, npti
                   IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp )   &
                      &   ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji)                  
@@ -794,6 +792,6 @@
       !
       ! In case ridging/rafting lead to very small negative values (sometimes it happens)
-      WHERE( a_i_2d(1:nidx,:) < 0._wp )   a_i_2d(1:nidx,:) = 0._wp
-      WHERE( v_i_2d(1:nidx,:) < 0._wp )   v_i_2d(1:nidx,:) = 0._wp
+      WHERE( a_i_2d(1:npti,:) < 0._wp )   a_i_2d(1:npti,:) = 0._wp
+      WHERE( v_i_2d(1:npti,:) < 0._wp )   v_i_2d(1:npti,:) = 0._wp
       !
    END SUBROUTINE rdgrft_shift
@@ -820,13 +818,13 @@
       !!----------------------------------------------------------------------
       !                              !--------------------------------------------------!
-      IF( ln_str_R75 ) THEN          ! Ice strength => Rothrock (1975) method           !
-      !                              !--------------------------------------------------!
-      !                              !--------------------------------------------------!
-      ELSEIF( ln_str_H79 ) THEN      ! Ice strength => Hibler (1979) method             !
+      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             !
       !                              !--------------------------------------------------!
          strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
-         !
          ismooth = 1
-         !
+         !                           !--------------------------------------------------!
+      ELSE                           ! Zero strength                                    !
+         !                           !--------------------------------------------------!
+         strength(:,:) = 0._wp
+         ismooth = 0
       ENDIF
       !                              !--------------------------------------------------!
@@ -896,5 +894,4 @@
       !!
       NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
-         &                    ln_str_R75, rn_perdg,          &
          &                    rn_csrdg  ,                    &
          &                    ln_partf_lin, rn_gstar,        &
@@ -921,6 +918,4 @@
          WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
          WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
-         WRITE(numout,*) '      ice strength parameterization Rothrock (1975)            ln_str_R75   = ', ln_str_R75 
-         WRITE(numout,*) '            Ratio of ridging work to PotEner change in ridging rn_perdg     = ', rn_perdg 
          WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_csrdg     = ', rn_csrdg 
          WRITE(numout,*) '      linear ridging participation function                    ln_partf_lin = ', ln_partf_lin
@@ -940,8 +935,4 @@
       ENDIF
       !
-      IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN
-         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one formulation for ice strength (ln_str_H79 or ln_str_R75)' )
-      ENDIF
-      !
       IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN
          CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8565)
@@ -60,5 +60,5 @@
       !
       INTEGER  ::   ji, jj, jl, jcat     ! dummy loop index
-      INTEGER  ::   nidx2                ! local integer
+      INTEGER  ::   ipti                 ! local integer
       REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars
       REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      -
@@ -66,5 +66,5 @@
       REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
 
-      INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not
+      INTEGER , DIMENSION(jpij)       ::   iptidx          ! compute remapping or not
       INTEGER , DIMENSION(jpij,jpl-1) ::   jdonor          ! donor category index
       REAL(wp), DIMENSION(jpij,jpl)   ::   zdhice          ! ice thickness increment
@@ -83,10 +83,10 @@
       !  1) Identify grid cells with ice
       !-----------------------------------------------------------------------------------------------
-      nidx = 0   ;   idxice(:) = 0
+      npti = 0   ;   nptidx(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF ( at_i(ji,jj) > epsi10 ) THEN
-               nidx         = nidx + 1
-               idxice( nidx ) = (jj - 1) * jpi + ji
+               npti         = npti + 1
+               nptidx( npti ) = (jj - 1) * jpi + ji
             ENDIF
          END DO
@@ -96,17 +96,17 @@
       !  2) Compute new category boundaries
       !-----------------------------------------------------------------------------------------------
-      IF( nidx > 0 ) THEN
+      IF( npti > 0 ) THEN
 
          zdhice(:,:) = 0._wp
          zhbnew(:,:) = 0._wp
 
-         CALL tab_3d_2d( nidx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i   )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), h_ib_2d(1:nidx,1:jpl), h_i_b )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i    )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b  )
+         CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i   )
+         CALL tab_3d_2d( npti, nptidx(1:npti), h_ib_2d(1:npti,1:jpl), h_i_b )
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i    )
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_ib_2d (1:npti,1:jpl), a_i_b  )
 
          DO jl = 1, jpl
             ! Compute thickness change in each ice category
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zdhice(ji,jl) = h_i_2d(ji,jl) - h_ib_2d(ji,jl)
             END DO
@@ -116,5 +116,5 @@
          DO jl = 1, jpl - 1
             !
-            DO ji = 1, nidx
+            DO ji = 1, npti
                !
                ! --- New boundary: Hn* = Hn + Fn*dt --- !
@@ -136,10 +136,10 @@
                !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
                !          in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
-               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   idxice(ji) = 0
-               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   idxice(ji) = 0
+               IF( a_i_2d(ji,jl  ) > epsi10 .AND. h_i_2d(ji,jl  ) > ( zhbnew(ji,jl) - epsi10 ) )   nptidx(ji) = 0
+               IF( a_i_2d(ji,jl+1) > epsi10 .AND. h_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) )   nptidx(ji) = 0
                
                ! 2) Hn-1 < Hn* < Hn+1  
-               IF( zhbnew(ji,jl) < hi_max(jl-1) )   idxice(ji) = 0
-               IF( zhbnew(ji,jl) > hi_max(jl+1) )   idxice(ji) = 0
+               IF( zhbnew(ji,jl) < hi_max(jl-1) )   nptidx(ji) = 0
+               IF( zhbnew(ji,jl) > hi_max(jl+1) )   nptidx(ji) = 0
                
             END DO
@@ -147,5 +147,5 @@
          !
          ! --- New boundaries for category jpl --- !
-         DO ji = 1, nidx
+         DO ji = 1, npti
             IF( a_i_2d(ji,jpl) > epsi10 ) THEN
                zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * h_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) )
@@ -158,6 +158,6 @@
             !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 
             !    in ice_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice)
-            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   idxice(ji) = 0
-            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   idxice(ji) = 0
+            IF( h_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) )   nptidx(ji) = 0
+            IF( h_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) )   nptidx(ji) = 0
          END DO
          !
@@ -165,14 +165,14 @@
          !  3) Identify cells where remapping
          !-----------------------------------------------------------------------------------------------
-         nidx2 = 0 ; idxice2(:) = 0
-         DO ji = 1, nidx
-            IF( idxice(ji) /= 0 ) THEN
-               nidx2 = nidx2 + 1
-               idxice2(nidx2) = idxice(ji)
-               zhbnew(nidx2,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
+         ipti = 0 ; iptidx(:) = 0
+         DO ji = 1, npti
+            IF( nptidx(ji) /= 0 ) THEN
+               ipti = ipti + 1
+               iptidx(ipti)   = nptidx(ji)
+               zhbnew(ipti,:) = zhbnew(ji,:)  ! adjust zhbnew to new indices
             ENDIF
          END DO
-         idxice(:) = idxice2(:)
-         nidx      = nidx2
+         nptidx(:) = iptidx(:)
+         npti      = ipti
          !
       ENDIF
@@ -181,5 +181,5 @@
       !  4) Compute g(h) 
       !-----------------------------------------------------------------------------------------------
-      IF( nidx > 0 ) THEN
+      IF( npti > 0 ) THEN
          !
          zhb0(:) = hi_max(0)   ;   zhb1(:) = hi_max(1)
@@ -189,17 +189,17 @@
          DO jl = 1, jpl
             !
-            CALL tab_2d_1d( nidx, idxice(1:nidx), h_ib_1d(1:nidx), h_i_b(:,:,jl) )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl)   )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_ib_1d(1:npti), h_i_b(:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
+            CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
             !
             IF( jl == 1 ) THEN
                !  
                ! --- g(h) for category 1 --- !
-               CALL ice_itd_glinear( zhb0(1:nidx)  , zhb1(1:nidx)  , h_ib_1d(1:nidx)  , a_i_1d(1:nidx)  ,  &   ! in
-                  &                  g0  (1:nidx,1), g1  (1:nidx,1), hL      (1:nidx,1), hR    (1:nidx,1)   )   ! out
+               CALL ice_itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
+                  &                  g0  (1:npti,1), g1  (1:npti,1), hL      (1:npti,1), hR    (1:npti,1)   )   ! out
                   !
                ! Area lost due to melting of thin ice
-               DO ji = 1, nidx
+               DO ji = 1, npti
                   !
                   IF( a_i_1d(ji) > epsi10 ) THEN
@@ -233,13 +233,13 @@
                END DO
                !
-               CALL tab_1d_2d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl)   )
-               CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
-               CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
+               CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
+               CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
+               CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
                !
             ENDIF ! jl=1
             !
             ! --- g(h) for each thickness category --- !  
-            CALL ice_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), h_i_1d(1:nidx)   , a_i_1d(1:nidx)   ,  &   ! in
-               &                  g0    (1:nidx,jl  ), g1    (1:nidx,jl), hL     (1:nidx,jl), hR    (1:nidx,jl)   )   ! out
+            CALL ice_itd_glinear( zhbnew(1:npti,jl-1), zhbnew(1:npti,jl), h_i_1d(1:npti)   , a_i_1d(1:npti)   ,  &   ! in
+               &                  g0    (1:npti,jl  ), g1    (1:npti,jl), hL     (1:npti,jl), hR    (1:npti,jl)   )   ! out
             !
          END DO
@@ -250,5 +250,5 @@
          DO jl = 1, jpl - 1
             !
-            DO ji = 1, nidx
+            DO ji = 1, npti
                !
                ! left and right integration limits in eta space
@@ -281,14 +281,14 @@
          ! 6) Shift ice between categories
          !----------------------------------------------------------------------------------------------
-         CALL ice_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )
+         CALL ice_itd_shiftice ( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )
          
          !----------------------------------------------------------------------------------------------
          ! 7) Make sure h_i >= minimum ice thickness hi_min
          !----------------------------------------------------------------------------------------------
-         CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,1)   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
          
-         DO ji = 1, nidx
+         DO ji = 1, npti
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
                a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
@@ -302,7 +302,7 @@
          END DO
          !
-         CALL tab_1d_2d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,1)   )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)   )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
          !
       ENDIF
@@ -338,5 +338,5 @@
       z2_3 = 2._wp / 3._wp
       !
-      DO ji = 1, nidx
+      DO ji = 1, npti
          !
          IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
@@ -392,13 +392,13 @@
       !!------------------------------------------------------------------
          
-      CALL tab_3d_2d( nidx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i  )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i  )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i  )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s  )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip )
-      CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d(1:nidx,1:jpl), t_su )
+      CALL tab_3d_2d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
+      CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
+      CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
+      CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
+      CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
+      CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
+      CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
+      CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
+      CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
 
       !----------------------------------------------------------------------------------------------
@@ -406,5 +406,5 @@
       !----------------------------------------------------------------------------------------------
       DO jl = 1, jpl
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl)
          END DO
@@ -415,5 +415,5 @@
       !-------------------------------------------------------------------------------
       DO jl = 1, jpl - 1
-         DO ji = 1, nidx
+         DO ji = 1, npti
             !
             jl1 = kdonor(ji,jl)
@@ -475,7 +475,7 @@
          DO jk = 1, nlay_s         !--- Snow heat content
             !
-            DO ji = 1, nidx
-               ii = MOD( idxice(ji) - 1, jpi ) + 1
-               ij = ( idxice(ji) - 1 ) / jpi + 1
+            DO ji = 1, npti
+               ii = MOD( nptidx(ji) - 1, jpi ) + 1
+               ij = ( nptidx(ji) - 1 ) / jpi + 1
                !
                jl1 = kdonor(ji,jl)
@@ -494,7 +494,7 @@
 
          DO jk = 1, nlay_i         !--- Ice heat content
-            DO ji = 1, nidx
-               ii = MOD( idxice(ji) - 1, jpi ) + 1
-               ij = ( idxice(ji) - 1 ) / jpi + 1
+            DO ji = 1, npti
+               ii = MOD( nptidx(ji) - 1, jpi ) + 1
+               ij = ( nptidx(ji) - 1 ) / jpi + 1
                !
                jl1 = kdonor(ji,jl)
@@ -517,21 +517,21 @@
       ! 3) Update ice thickness and temperature
       !-------------------------------------------------------------------------------
-      WHERE( a_i_2d(1:nidx,:) >= epsi20 )
-         h_i_2d(1:nidx,:)  =  v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:) 
-         t_su_2d(1:nidx,:)  =  zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:) 
+      WHERE( a_i_2d(1:npti,:) >= epsi20 )
+         h_i_2d(1:npti,:)  =  v_i_2d(1:npti,:) / a_i_2d(1:npti,:) 
+         t_su_2d(1:npti,:)  =  zaTsfn(1:npti,:) / a_i_2d(1:npti,:) 
       ELSEWHERE
-         h_i_2d(1:nidx,:)  = 0._wp
-         t_su_2d(1:nidx,:)  = rt0
+         h_i_2d(1:npti,:)  = 0._wp
+         t_su_2d(1:npti,:)  = rt0
       END WHERE
       !
-      CALL tab_2d_3d( nidx, idxice(1:nidx), h_i_2d (1:nidx,1:jpl), h_i  )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i  )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i  )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s  )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip )
-      CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d(1:nidx,1:jpl), t_su )
+      CALL tab_2d_3d( npti, nptidx(1:npti), h_i_2d (1:npti,1:jpl), h_i  )
+      CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i  )
+      CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i  )
+      CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s  )
+      CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i )
+      CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i )
+      CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
+      CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
+      CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       !
    END SUBROUTINE ice_itd_shiftice
@@ -564,19 +564,19 @@
       DO jl = 1, jpl-1        ! identify thicknesses that are too big
          !                    !---------------------------------------
-         nidx = 0   ;   idxice(:) = 0
+         npti = 0   ;   nptidx(:) = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN
-                  nidx = nidx + 1
-                  idxice( nidx ) = (jj - 1) * jpi + ji                  
+                  npti = npti + 1
+                  nptidx( npti ) = (jj - 1) * jpi + ji                  
                ENDIF
             END DO
          END DO
          !
-!!clem   CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d (1:nidx), h_i(:,:,jl)   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl)    )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl)    )
-         !
-         DO ji = 1, nidx
+!!clem   CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,jl)   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl)    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl)    )
+         !
+         DO ji = 1, npti
             jdonor(ji,jl)  = jl 
             ! how much of a_i you send in cat sup is somewhat arbitrary
@@ -592,10 +592,10 @@
          END DO
          !
-         IF( nidx > 0 ) THEN
-            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl=>jl+1
+         IF( npti > 0 ) THEN
+            CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl=>jl+1
             ! Reset shift parameters
-            jdonor(1:nidx,jl) = 0
-            zdaice(1:nidx,jl) = 0._wp
-            zdvice(1:nidx,jl) = 0._wp
+            jdonor(1:npti,jl) = 0
+            zdaice(1:npti,jl) = 0._wp
+            zdvice(1:npti,jl) = 0._wp
          ENDIF
          !
@@ -605,17 +605,17 @@
       DO jl = jpl-1, 1, -1    ! Identify thicknesses that are too small
          !                    !-----------------------------------------
-         nidx = 0 ; idxice(:) = 0
+         npti = 0 ; nptidx(:) = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN
-                  nidx = nidx + 1
-                  idxice( nidx ) = (jj - 1) * jpi + ji                  
+                  npti = npti + 1
+                  nptidx( npti ) = (jj - 1) * jpi + ji                  
                ENDIF
             END DO
          END DO
          !
-         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,jl+1)    ) ! jl+1 is ok
-         CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d  (1:nidx), v_i(:,:,jl+1)    ) ! jl+1 is ok
-         DO ji = 1, nidx
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,jl+1)    ) ! jl+1 is ok
+         CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d  (1:npti), v_i(:,:,jl+1)    ) ! jl+1 is ok
+         DO ji = 1, npti
             jdonor(ji,jl) = jl + 1
             zdaice(ji,jl) = a_i_1d(ji) 
@@ -623,10 +623,10 @@
          END DO
          !
-         IF( nidx > 0 ) THEN
-            CALL ice_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) )  ! Shift jl+1=>jl
+         IF( npti > 0 ) THEN
+            CALL ice_itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) )  ! Shift jl+1=>jl
             ! Reset shift parameters
-            jdonor(1:nidx,jl) = 0
-            zdaice(1:nidx,jl) = 0._wp
-            zdvice(1:nidx,jl) = 0._wp
+            jdonor(1:npti,jl) = 0
+            zdaice(1:npti,jl) = 0._wp
+            zdvice(1:npti,jl) = 0._wp
          ENDIF
          !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8565)
@@ -211,29 +211,29 @@
 
          ! select ice covered grid points
-         nidx = 0 ; idxice(:) = 0
+         npti = 0 ; nptidx(:) = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF ( a_i(ji,jj,jl) > epsi10 ) THEN     
-                  nidx         = nidx  + 1
-                  idxice(nidx) = (jj - 1) * jpi + ji
+                  npti         = npti  + 1
+                  nptidx(npti) = (jj - 1) * jpi + ji
                ENDIF
             END DO
          END DO
 
-         IF( lk_mpp )         CALL mpp_ini_ice( nidx , numout )
-
-         IF( nidx > 0 ) THEN  ! If there is no ice, do nothing.
+         IF( lk_mpp )         CALL mpp_ini_ice( npti , numout )
+
+         IF( npti > 0 ) THEN  ! If there is no ice, do nothing.
             !                                                                
                               CALL ice_thd_1d2d( jl, 1 )            ! --- Move to 1D arrays --- !
             !                                                       ! --- & Change units of e_i, e_s from J/m2 to J/m3 --- !
             !
-            s_i_new   (1:nidx) = 0._wp ; dh_s_tot (1:nidx) = 0._wp  ! --- some init --- !  (important to have them here) 
-            dh_i_surf (1:nidx) = 0._wp ; dh_i_bott(1:nidx) = 0._wp
-            dh_snowice(1:nidx) = 0._wp ; dh_i_sub (1:nidx) = 0._wp
+            s_i_new   (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp  ! --- some init --- !  (important to have them here) 
+            dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp
+            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp
             !
             IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
                               CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
                               CALL ice_thd_dh                              ! Ice/Snow thickness   
-                              CALL ice_thd_ent( e_i_1d(1:nidx,:) )         ! Ice enthalpy remapping
+                              CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping
             ENDIF
             !
@@ -294,5 +294,5 @@
       ! Recover ice temperature
       DO jk = 1, nlay_i
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ztmelts       = -tmut * sz_i_1d(ji,jk)
             ! Conversion q(S,T) -> T (second order equation)
@@ -323,5 +323,5 @@
       !!-----------------------------------------------------------------------
       !
-      DO ji = 1, nidx
+      DO ji = 1, npti
          zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) + dh_i_sub(ji) )
          IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
@@ -360,84 +360,84 @@
       CASE( 1 )            !==  from 2D to 1D  ==!
          !                 !---------------------!
-         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i             )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,kl)     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), h_i_1d(1:nidx), h_i(:,:,kl)     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), h_s_1d(1:nidx), h_s(:,:,kl)     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_su_1d(1:nidx), t_su(:,:,kl)     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     )
          DO jk = 1, nlay_s
-            CALL tab_2d_1d( nidx, idxice(1:nidx), t_s_1d(1:nidx,jk), t_s(:,:,jk,kl)   )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,kl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)   )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_2d_1d( nidx, idxice(1:nidx), t_i_1d(1:nidx,jk), t_i(:,:,jk,kl)   )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,kl)   )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), sz_i_1d(1:nidx,jk), sz_i(:,:,jk,kl)   )
-         END DO
-         !
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qprec_ice_1d(1:nidx), qprec_ice        )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qsr_ice_1d  (1:nidx), qsr_ice (:,:,kl) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), fr1_i0_1d   (1:nidx), fr1_i0           )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), fr2_i0_1d   (1:nidx), fr2_i0           )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qns_ice_1d  (1:nidx), qns_ice (:,:,kl) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), ftr_ice_1d  (1:nidx), ftr_ice (:,:,kl) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), evap_ice_1d (1:nidx), evap_ice(:,:,kl) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), dqns_ice_1d (1:nidx), dqns_ice(:,:,kl) )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d     (1:nidx), t_bo             )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sprecip_1d  (1:nidx), sprecip          ) 
-         CALL tab_2d_1d( nidx, idxice(1:nidx), fhtur_1d    (1:nidx), fhtur            )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), fhld_1d     (1:nidx), fhld             )
-         !
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_sni_1d(1:nidx), wfx_snw_sni   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_sum_1d(1:nidx), wfx_snw_sum   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_sub_1d    (1:nidx), wfx_sub       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_sub_1d(1:nidx), wfx_snw_sub   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_ice_sub_1d(1:nidx), wfx_ice_sub   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_err_sub_1d(1:nidx), wfx_err_sub   )
-         !
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_bog_1d (1:nidx), wfx_bog          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_bom_1d (1:nidx), wfx_bom          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_sum_1d (1:nidx), wfx_sum          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_sni_1d (1:nidx), wfx_sni          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam          )
-         !
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bom_1d (1:nidx), sfx_bom          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sum_1d (1:nidx), sfx_sum          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sni_1d (1:nidx), sfx_sni          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam          )
-         !
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sum_1d (1:nidx), hfx_sum          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_bom_1d (1:nidx), hfx_bom          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_bog_1d (1:nidx), hfx_bog          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dif_1d (1:nidx), hfx_dif          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d (1:nidx), hfx_opw          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out          )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)   )
+         END DO
+         !
+         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d  (1:npti), qsr_ice (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d   (1:npti), fr1_i0           )
+         CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d   (1:npti), fr2_i0           )
+         CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d  (1:npti), qns_ice (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d  (1:npti), ftr_ice (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d (1:npti), evap_ice(:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), dqns_ice_1d (1:npti), dqns_ice(:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d     (1:npti), t_bo             )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sprecip_1d  (1:npti), sprecip          ) 
+         CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d    (1:npti), fhtur            )
+         CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d     (1:npti), fhld             )
+         !
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sub_1d    (1:npti), wfx_sub       )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub   )
+         !
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          )
+         !
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam          )
+         !
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out          )
          !
          ! SIMIP diagnostics
-         CALL tab_2d_1d( nidx, idxice(1:nidx), diag_fc_bo_1d(1:nidx), diag_fc_bo   )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), diag_fc_su_1d(1:nidx), diag_fc_su   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), diag_fc_bo_1d(1:npti), diag_fc_bo   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), diag_fc_su_1d(1:npti), diag_fc_su   )
          ! ocean surface fields
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m )
 
          ! --- Change units of e_i, e_s from J/m2 to J/m3 --- !
          DO jk = 1, nlay_i
-            WHERE( h_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (h_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i
+            WHERE( h_i_1d(1:npti)>0._wp ) e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) / (h_i_1d(1:npti) * a_i_1d(1:npti)) * nlay_i
          END DO
          DO jk = 1, nlay_s
-            WHERE( h_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (h_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s
+            WHERE( h_s_1d(1:npti)>0._wp ) e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) / (h_s_1d(1:npti) * a_i_1d(1:npti)) * nlay_s
          END DO
          !
@@ -447,80 +447,80 @@
          ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
          DO jk = 1, nlay_i
-            e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * h_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i
+            e_i_1d(1:npti,jk) = e_i_1d(1:npti,jk) * h_i_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_i
          END DO
          DO jk = 1, nlay_s
-            e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * h_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s
+            e_s_1d(1:npti,jk) = e_s_1d(1:npti,jk) * h_s_1d(1:npti) * a_i_1d(1:npti) * r1_nlay_s
          END DO
          !
          ! Change thickness to volume (replaces routine ice_var_eqv2glo)
-         v_i_1d(1:nidx)  = h_i_1d(1:nidx) * a_i_1d(1:nidx)
-         v_s_1d(1:nidx)  = h_s_1d(1:nidx) * a_i_1d(1:nidx)
-         sv_i_1d(1:nidx) = s_i_1d(1:nidx) * v_i_1d(1:nidx)
+         v_i_1d(1:npti)  = h_i_1d(1:npti) * a_i_1d(1:npti)
+         v_s_1d(1:npti)  = h_s_1d(1:npti) * a_i_1d(1:npti)
+         sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti)
          
-         CALL tab_1d_2d( nidx, idxice(1:nidx), at_i_1d(1:nidx), at_i             )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,kl)     )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), h_i_1d(1:nidx), h_i(:,:,kl)     )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), h_s_1d(1:nidx), h_s(:,:,kl)     )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), t_su_1d(1:nidx), t_su(:,:,kl)     )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     )
          DO jk = 1, nlay_s
-            CALL tab_1d_2d( nidx, idxice(1:nidx), t_s_1d(1:nidx,jk), t_s(:,:,jk,kl) )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), e_s_1d(1:nidx,jk), e_s(:,:,jk,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_1d_2d( nidx, idxice(1:nidx), t_i_1d(1:nidx,jk), t_i(:,:,jk,kl) )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), e_i_1d(1:nidx,jk), e_i(:,:,jk,kl) )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), sz_i_1d(1:nidx,jk), sz_i(:,:,jk,kl) )
-         END DO
-         !
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_sni_1d(1:nidx), wfx_snw_sni )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_sum_1d(1:nidx), wfx_snw_sum )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_sub_1d    (1:nidx), wfx_sub     )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_sub_1d(1:nidx), wfx_snw_sub )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_ice_sub_1d(1:nidx), wfx_ice_sub )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_err_sub_1d(1:nidx), wfx_err_sub )
-         !
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_bog_1d (1:nidx), wfx_bog        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_bom_1d (1:nidx), wfx_bom        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_sum_1d (1:nidx), wfx_sum        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_sni_1d (1:nidx), wfx_sni        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_res_1d (1:nidx), wfx_res        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_spr_1d (1:nidx), wfx_spr        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_lam_1d (1:nidx), wfx_lam        )
-         !
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bog_1d (1:nidx), sfx_bog        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bom_1d (1:nidx), sfx_bom        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sum_1d (1:nidx), sfx_sum        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sni_1d (1:nidx), sfx_sni        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_res_1d (1:nidx), sfx_res        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_sub_1d (1:nidx), sfx_sub        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_lam_1d (1:nidx), sfx_lam        )
-         !
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d (1:nidx), hfx_thd        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_spr_1d (1:nidx), hfx_spr        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sum_1d (1:nidx), hfx_sum        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_bom_1d (1:nidx), hfx_bom        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_bog_1d (1:nidx), hfx_bog        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dif_1d (1:nidx), hfx_dif        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d (1:nidx), hfx_opw        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_rem_1d(1:nidx), hfx_err_rem )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_out_1d (1:nidx), hfx_out        )
-         !
-         CALL tab_1d_2d( nidx, idxice(1:nidx), qns_ice_1d(1:nidx), qns_ice(:,:,kl) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), ftr_ice_1d(1:nidx), ftr_ice(:,:,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
+         END DO
+         !
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sub_1d    (1:npti), wfx_sub     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sub_1d(1:npti), wfx_snw_sub )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_ice_sub_1d(1:npti), wfx_ice_sub )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_err_sub_1d(1:npti), wfx_err_sub )
+         !
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bog_1d (1:npti), wfx_bog        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_bom_1d (1:npti), wfx_bom        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sni_1d (1:npti), wfx_sni        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_res_1d (1:npti), wfx_res        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        )
+         !
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bom_1d (1:npti), sfx_bom        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sum_1d (1:npti), sfx_sum        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sni_1d (1:npti), sfx_sni        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_res_1d (1:npti), sfx_res        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        )
+         !
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d (1:npti), hfx_thd        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_spr_1d (1:npti), hfx_spr        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sum_1d (1:npti), hfx_sum        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bom_1d (1:npti), hfx_bom        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_bog_1d (1:npti), hfx_bog        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dif_1d (1:npti), hfx_dif        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d (1:npti), hfx_opw        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_snw_1d (1:npti), hfx_snw        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_sub_1d (1:npti), hfx_sub        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_res_1d (1:npti), hfx_res        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_dif_1d(1:npti), hfx_err_dif )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_err_rem_1d(1:npti), hfx_err_rem )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_out_1d (1:npti), hfx_out        )
+         !
+         CALL tab_1d_2d( npti, nptidx(1:npti), qns_ice_1d(1:npti), qns_ice(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), ftr_ice_1d(1:npti), ftr_ice(:,:,kl) )
          !
          ! SIMIP diagnostics         
-         CALL tab_1d_2d( nidx, idxice(1:nidx), t_si_1d      (1:nidx), t_si(:,:,kl) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), diag_fc_bo_1d(1:nidx), diag_fc_bo   )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), diag_fc_su_1d(1:nidx), diag_fc_su   )
+         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d      (1:npti), t_si(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), diag_fc_bo_1d(1:npti), diag_fc_bo   )
+         CALL tab_1d_2d( npti, nptidx(1:npti), diag_fc_su_1d(1:npti), diag_fc_su   )
          ! extensive variables
-         CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i (:,:,kl) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), v_s_1d (1:nidx), v_s (:,:,kl) )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sv_i_1d(1:nidx), sv_i(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
          !
       END SELECT
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90	(revision 8565)
@@ -118,5 +118,5 @@
       !
       zastar = 1._wp / ( 1._wp - (rn_dmin / zdmax)**(1._wp/rn_beta) )
-      DO ji = 1, nidx   
+      DO ji = 1, npti   
          ! --- Calculate reduction of total sea ice concentration --- !
          zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta         ! Mean floe caliper diameter [m]
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8565)
@@ -117,8 +117,8 @@
 
       ! initialize layer thicknesses and enthalpies
-      h_i_old (1:nidx,0:nlay_i+1) = 0._wp
-      eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
+      h_i_old (1:npti,0:nlay_i+1) = 0._wp
+      eh_i_old(1:npti,0:nlay_i+1) = 0._wp
       DO jk = 1, nlay_i
-         DO ji = 1, nidx
+         DO ji = 1, npti
             h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
             eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
@@ -130,5 +130,5 @@
       !------------------------------------------------------------------------------!
       !
-      DO ji = 1, nidx
+      DO ji = 1, npti
          zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
@@ -142,5 +142,5 @@
       ! (should not happen but sometimes it does)
       !------------------------------------------------------------------------------!
-      DO ji = 1, nidx
+      DO ji = 1, npti
          IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting
             ! Contribution to heat flux to the ocean [W.m-2], < 0  
@@ -159,5 +159,5 @@
       !------------------------------------------------------------!
       DO jk = 1, nlay_i
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
             zeh_i(ji)   = zeh_i(ji) + e_i_1d(ji,jk) * zh_i(ji,jk)
@@ -183,8 +183,8 @@
       ! Martin Vancoppenolle, December 2006
 
-      CALL ice_thd_snwblow( 1. - at_i_1d(1:nidx), zsnw(1:nidx) ) ! snow distribution over ice after wind blowing
-
-      zdeltah(1:nidx,:) = 0._wp
-      DO ji = 1, nidx
+      CALL ice_thd_snwblow( 1. - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
+
+      zdeltah(1:npti,:) = 0._wp
+      DO ji = 1, npti
          !-----------
          ! Snow fall
@@ -220,8 +220,8 @@
 
       ! If heat still available (zq_su > 0), then melt more snow
-      zdeltah(1:nidx,:) = 0._wp
-      zdh_s_mel(1:nidx) = 0._wp
+      zdeltah(1:npti,:) = 0._wp
+      zdh_s_mel(1:npti) = 0._wp
       DO jk = 1, nlay_s
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ! thickness change
             rswitch          = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) 
@@ -245,6 +245,6 @@
       ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
       ! clem comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
-      zdeltah(1:nidx,:) = 0._wp
-      DO ji = 1, nidx
+      zdeltah(1:npti,:) = 0._wp
+      DO ji = 1, npti
          zdh_s_sub(ji)  = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
          ! remaining evap in kg.m-2 (used for ice melting later on)
@@ -265,5 +265,5 @@
       
       ! --- Update snow diags --- !
-      DO ji = 1, nidx
+      DO ji = 1, npti
          dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
       END DO
@@ -274,5 +274,5 @@
       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
       DO jk = 1, nlay_s
-         DO ji = 1,nidx
+         DO ji = 1,npti
             rswitch       = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
             e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) *           &
@@ -285,7 +285,7 @@
       ! 3.4 Surface ice ablation 
       !--------------------------
-      zdeltah(1:nidx,:) = 0._wp ! important
+      zdeltah(1:npti,:) = 0._wp ! important
       DO jk = 1, nlay_i
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ztmelts           = - tmut * sz_i_1d(ji,jk)          ! Melting point of layer k [C]
             
@@ -375,10 +375,10 @@
       END DO
       ! update ice thickness
-      DO ji = 1, nidx
+      DO ji = 1, npti
          h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
       END DO
 
       ! remaining "potential" evap is sent to ocean
-      DO ji = 1, nidx
+      DO ji = 1, npti
          wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
       END DO
@@ -407,5 +407,5 @@
       
       ! Iterative procedure
-      DO ji = 1, nidx
+      DO ji = 1, npti
          IF(  zf_tt(ji) < 0._wp  ) THEN
             DO iter = 1, num_iter_max
@@ -478,7 +478,7 @@
       ! 4.2 Basal melt
       !----------------
-      zdeltah(1:nidx,:) = 0._wp ! important
+      zdeltah(1:npti,:) = 0._wp ! important
       DO jk = nlay_i, 1, -1
-         DO ji = 1, nidx
+         DO ji = 1, npti
             IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting 
 
@@ -554,5 +554,5 @@
       ! Update temperature, energy
       !-------------------------------------------
-      DO ji = 1, nidx
+      DO ji = 1, npti
          h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_bott(ji) )
       END DO  
@@ -563,6 +563,6 @@
       ! If heat still available for melting and snow remains, then melt more snow
       !-------------------------------------------
-      zdeltah(1:nidx,:) = 0._wp ! important
-      DO ji = 1, nidx
+      zdeltah(1:npti,:) = 0._wp ! important
+      DO ji = 1, npti
          zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
          rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) )   ! =1 if snow
@@ -592,5 +592,5 @@
       ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
       z1_rho = 1._wp / ( rhosn+rau0-rhoic )
-      DO ji = 1, nidx
+      DO ji = 1, npti
          !
          dh_snowice(ji) = MAX(  0._wp , ( rhosn * h_s_1d(ji) + (rhoic-rau0) * h_i_1d(ji) ) * z1_rho )
@@ -631,5 +631,5 @@
       ! Update temperature, energy
       !-------------------------------------------
-      DO ji = 1, nidx
+      DO ji = 1, npti
          rswitch     =  1.0 - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 
          t_su_1d(ji) =  rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0
@@ -637,5 +637,5 @@
 
       DO jk = 1, nlay_s
-         DO ji = 1,nidx
+         DO ji = 1,npti
             ! mask enthalpy
             rswitch       = 1._wp - MAX(  0._wp , SIGN( 1._wp, - h_s_1d(ji) )  )
@@ -647,5 +647,5 @@
 
       ! --- ensure that a_i = 0 where h_i = 0 ---
-      WHERE( h_i_1d(1:nidx) == 0._wp )   a_i_1d(1:nidx) = 0._wp
+      WHERE( h_i_1d(1:npti) == 0._wp )   a_i_1d(1:npti) = 0._wp
       !
    END SUBROUTINE ice_thd_dh
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8565)
@@ -222,10 +222,10 @@
       !-------------------------------------
       ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
-      nidx = 0 ; idxice(:) = 0
+      npti = 0 ; nptidx(:) = 0
       DO jj = 1, jpj
          DO ji = 1, jpi
             IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
-               nidx = nidx + 1
-               idxice( nidx ) = (jj - 1) * jpi + ji
+               npti = npti + 1
+               nptidx( npti ) = (jj - 1) * jpi + ji
             ENDIF
          END DO
@@ -237,26 +237,26 @@
       ! If ocean gains heat do nothing. Otherwise compute new ice formation
 
-      IF ( nidx > 0 ) THEN
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d(1:nidx)      , at_i        )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )
-         CALL tab_3d_2d( nidx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )
+      IF ( npti > 0 ) THEN
+
+         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti)      , at_i        )
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
-            END DO
-         END DO
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx) , qlead       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx) , t_bo        )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , ht_i_new    )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx) , zvrel       )
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d  )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx) , sss_m       )
+               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead       )
+         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo        )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel       )
+
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d  )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m       )
 
          !------------------------------------------------------------------------------|
@@ -265,8 +265,8 @@
          DO jl = 1, jpl
             DO jk = 1, nlay_i               
-               WHERE( v_i_2d(1:nidx,jl) > 0._wp )
-                  ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) / v_i_2d(1:nidx,jl) * REAL( nlay_i )
+               WHERE( v_i_2d(1:npti,jl) > 0._wp )
+                  ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) / v_i_2d(1:npti,jl) * REAL( nlay_i )
                ELSEWHERE
-                  ze_i_2d(1:nidx,jk,jl) = 0._wp
+                  ze_i_2d(1:npti,jk,jl) = 0._wp
                END WHERE
             END DO
@@ -279,6 +279,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_b(1:nidx,:) = v_i_2d(1:nidx,:) 
-         za_b(1:nidx,:) = a_i_2d(1:nidx,:)
+         zv_b(1:npti,:) = v_i_2d(1:npti,:) 
+         za_b(1:npti,:) = a_i_2d(1:npti,:)
 
          !----------------------
@@ -287,11 +287,11 @@
          SELECT CASE ( nn_icesal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(1:nidx) = rn_icesal
+            zs_newice(1:npti) = rn_icesal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(1:nidx) =   2.3
+            zs_newice(1:npti) =   2.3
          END SELECT
 
@@ -300,5 +300,5 @@
          !-------------------------
          ! We assume that new ice is formed at the seawater freezing point
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C)
             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
@@ -310,10 +310,10 @@
          ! Age of new ice
          !----------------
-         zo_newice(1:nidx) = 0._wp
+         zo_newice(1:npti) = 0._wp
 
          !-------------------
          ! Volume of new ice
          !-------------------
-         DO ji = 1, nidx
+         DO ji = 1, npti
 
             zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
@@ -340,8 +340,8 @@
          END DO
          
-         zv_frazb(1:nidx) = 0._wp
+         zv_frazb(1:npti) = 0._wp
          IF( ln_frazil ) THEN
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            DO ji = 1, nidx
+            DO ji = 1, npti
                rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
                zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
@@ -354,5 +354,5 @@
          ! Area of new ice
          !-----------------
-         DO ji = 1, nidx
+         DO ji = 1, npti
             za_newice(ji) = zv_newice(ji) / zh_newice(ji)
          END DO
@@ -367,5 +367,5 @@
          ! If lateral ice growth gives an ice concentration gt 1, then
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
-         DO ji = 1, nidx
+         DO ji = 1, npti
             IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN
                zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) )
@@ -381,5 +381,5 @@
          ! find which category to fill
          DO jl = 1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
                   a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
@@ -389,8 +389,8 @@
             END DO
          END DO
-         at_i_1d(1:nidx) = SUM( a_i_2d(1:nidx,:), dim=2 )
+         at_i_1d(1:npti) = SUM( a_i_2d(1:npti,:), dim=2 )
 
          ! Heat content
-         DO ji = 1, nidx
+         DO ji = 1, npti
             jl = jcat(ji)                                                    ! categroy in which new ice is put
             zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
@@ -398,5 +398,5 @@
 
          DO jk = 1, nlay_i
-            DO ji = 1, nidx
+            DO ji = 1, npti
                jl = jcat(ji)
                rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
@@ -413,8 +413,8 @@
 
             ! for remapping
-            h_i_old (1:nidx,0:nlay_i+1) = 0._wp
-            eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
+            h_i_old (1:npti,0:nlay_i+1) = 0._wp
+            eh_i_old(1:npti,0:nlay_i+1) = 0._wp
             DO jk = 1, nlay_i
-               DO ji = 1, nidx
+               DO ji = 1, npti
                   h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
                   eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
@@ -423,5 +423,5 @@
 
             ! new volumes including lateral/bottom accretion + residual
-            DO ji = 1, nidx
+            DO ji = 1, npti
                rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
                zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
@@ -433,5 +433,5 @@
             ENDDO
             ! --- Ice enthalpy remapping --- !
-            CALL ice_thd_ent( ze_i_2d(1:nidx,:,jl) ) 
+            CALL ice_thd_ent( ze_i_2d(1:npti,:,jl) ) 
          ENDDO
 
@@ -440,5 +440,5 @@
          !-----------------
          DO jl = 1, jpl
-            DO ji = 1, nidx
+            DO ji = 1, npti
                sv_i_2d(ji,jl) = sv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
             END DO
@@ -450,5 +450,5 @@
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) * v_i_2d(1:nidx,jl) * r1_nlay_i 
+               ze_i_2d(1:npti,jk,jl) = ze_i_2d(1:npti,jk,jl) * v_i_2d(1:npti,jl) * r1_nlay_i 
             END DO
          END DO
@@ -456,18 +456,18 @@
          ! 7) Change 2D vectors to 1D vectors 
          !------------------------------------------------------------------------------!
-         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )
-         CALL tab_2d_3d( nidx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
           DO jl = 1, jpl
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
-            END DO
-         END DO
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw )
+               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_opw_1d(1:npti), sfx_opw )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_opw_1d(1:npti), wfx_opw )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d(1:npti), hfx_thd )
+         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_opw_1d(1:npti), hfx_opw )
          !
-      ENDIF ! nidx > 0
+      ENDIF ! npti > 0
       !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8565)
@@ -78,8 +78,8 @@
       !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
       !--------------------------------------------------------------------------
-      zeh_cum0(1:nidx,0) = 0._wp 
-      zh_cum0 (1:nidx,0) = 0._wp
+      zeh_cum0(1:npti,0) = 0._wp 
+      zh_cum0 (1:npti,0) = 0._wp
       DO jk0 = 1, nlay_i+2
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + eh_i_old(ji,jk0-1)
             zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1)
@@ -91,21 +91,21 @@
       !------------------------------------
       ! new layer thickesses
-      DO ji = 1, nidx
+      DO ji = 1, npti
          zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i  
       END DO
 
       ! new layers interfaces
-      zh_cum1(1:nidx,0) = 0._wp
+      zh_cum1(1:npti,0) = 0._wp
       DO jk1 = 1, nlay_i
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji)
          END DO
       END DO
 
-      zeh_cum1(1:nidx,0) = 0._wp 
+      zeh_cum1(1:npti,0) = 0._wp 
       ! new cumulative q*h => linear interpolation
       DO jk0 = 1, nlay_i+1
          DO jk1 = 1, nlay_i-1
-            DO ji = 1, nidx
+            DO ji = 1, npti
                IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN
                   zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1  ) ) +  &
@@ -117,9 +117,9 @@
       END DO
       ! to ensure that total heat content is strictly conserved, set:
-      zeh_cum1(1:nidx,nlay_i) = zeh_cum0(1:nidx,nlay_i+2) 
+      zeh_cum1(1:npti,nlay_i) = zeh_cum0(1:npti,nlay_i+2) 
 
       ! new enthalpies
       DO jk1 = 1, nlay_i
-         DO ji = 1, nidx
+         DO ji = 1, npti
             rswitch      = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 
             qnew(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
@@ -130,5 +130,5 @@
       ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do), 
       ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
-      DO ji = 1, nidx
+      DO ji = 1, npti
          hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice *  &
             &               ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( eh_i_old(ji,0:nlay_i+1) ) ) 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8565)
@@ -70,5 +70,5 @@
          z1_time_fl = 1._wp / rn_time_fl * rdt_ice
          !
-         DO ji = 1, nidx
+         DO ji = 1, npti
 
             !---------------------------------------------------------
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8565)
@@ -144,5 +144,5 @@
 
       ! --- diag error on heat diffusion - PART 1 --- !
-      DO ji = 1, nidx
+      DO ji = 1, npti
          zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
             &           SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
@@ -152,5 +152,5 @@
       ! 1) Initialization
       !------------------
-      DO ji = 1, nidx
+      DO ji = 1, npti
          isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - h_s_1d(ji) ) )  ! is there snow or not
          ! layer thickness
@@ -159,17 +159,17 @@
       END DO
       !
-      WHERE( zh_i(1:nidx) >= epsi10 )   ;   z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)
-      ELSEWHERE                         ;   z1_h_i(1:nidx) = 0._wp
+      WHERE( zh_i(1:npti) >= epsi10 )   ;   z1_h_i(1:npti) = 1._wp / zh_i(1:npti)
+      ELSEWHERE                         ;   z1_h_i(1:npti) = 0._wp
       END WHERE
 
-      WHERE( zh_s(1:nidx) >= epsi10 )   ;   z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)
-      ELSEWHERE                         ;   z1_h_s(1:nidx) = 0._wp
+      WHERE( zh_s(1:npti) >= epsi10 )   ;   z1_h_s(1:npti) = 1._wp / zh_s(1:npti)
+      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp
       END WHERE
       !
       ! temperatures
-      ztsub  (1:nidx)   = t_su_1d(1:nidx)  ! temperature at the previous iteration
-      ztsold (1:nidx,:) = t_s_1d(1:nidx,:)   ! Old snow temperature
-      ztiold (1:nidx,:) = t_i_1d(1:nidx,:)   ! Old ice temperature
-      t_su_1d(1:nidx)   = MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
+      ztsub  (1:npti)   = t_su_1d(1:npti)  ! temperature at the previous iteration
+      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature
+      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature
+      t_su_1d(1:npti)   = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! necessary
       !
       !-------------
@@ -177,5 +177,5 @@
       !-------------
       z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
-      DO ji = 1, nidx
+      DO ji = 1, npti
          ! --- Computation of i0 --- !
          ! i0 describes the fraction of solar radiation which does not contribute
@@ -199,7 +199,7 @@
 
       ! --- Transmission/absorption of solar radiation in the ice --- !
-      zradtr_s(1:nidx,0) = zftrice(1:nidx)
+      zradtr_s(1:npti,0) = zftrice(1:npti)
       DO jk = 1, nlay_s
-         DO ji = 1, nidx
+         DO ji = 1, npti
             !                             ! radiation transmitted below the layer-th snow layer
             zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
@@ -209,7 +209,7 @@
       END DO
          
-      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
+      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) )
       DO jk = 1, nlay_i 
-         DO ji = 1, nidx
+         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) )
@@ -219,5 +219,5 @@
       END DO
 
-      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
+      ftr_ice_1d(1:npti) = zradtr_i(1:npti,nlay_i)   ! record radiation transmitted below the ice
       !
       iconv    =  0          ! number of iterations
@@ -228,6 +228,6 @@
          iconv = iconv + 1
          !
-         ztib(1:nidx,:) = t_i_1d(1:nidx,:)
-         ztsb(1:nidx,:) = t_s_1d(1:nidx,:)
+         ztib(1:npti,:) = t_i_1d(1:npti,:)
+         ztsb(1:npti,:) = t_s_1d(1:npti,:)
          !
          !--------------------------------
@@ -236,10 +236,10 @@
          IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
             !
-            DO ji = 1, nidx
+            DO ji = 1, npti
                ztcond_i(ji,0)      = rcdic + zbeta * sz_i_1d(ji,1)      / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
                ztcond_i(ji,nlay_i) = rcdic + 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, nidx
+               DO ji = 1, npti
                   ztcond_i(ji,jk) = rcdic + 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 )
@@ -249,5 +249,5 @@
          ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
             !
-            DO ji = 1, nidx
+            DO ji = 1, npti
                ztcond_i(ji,0)      = rcdic + 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 )
@@ -256,5 +256,5 @@
             END DO
             DO jk = 1, nlay_i-1
-               DO ji = 1, nidx
+               DO ji = 1, npti
                   ztcond_i(ji,jk) = rcdic + 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 )   &
@@ -264,5 +264,5 @@
             !
          ENDIF
-         ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )        
+         ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) )        
          !
          !--- G(he) : enhancement of thermal conductivity in mono-category case
@@ -270,5 +270,5 @@
          ! Used in mono-category case only to simulate an ITD implicitly
          ! Fichefet and Morales Maqueda, JGR 1997
-         zghe(1:nidx) = 1._wp
+         zghe(1:npti) = 1._wp
          !
          SELECT CASE ( nn_monocat )
@@ -277,5 +277,5 @@
 
             zepsilon = 0.1_wp
-            DO ji = 1, nidx
+            DO ji = 1, npti
                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)
@@ -292,9 +292,9 @@
          !--- Snow
          DO jk = 0, nlay_s-1
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
             END DO
          END DO
-         DO ji = 1, nidx   ! Snow-ice interface
+         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
@@ -307,9 +307,9 @@
          !--- Ice
          DO jk = 0, nlay_i
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
             END DO
          END DO
-         DO ji = 1, nidx   ! Snow-ice interface
+         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
@@ -319,5 +319,5 @@
          !--------------------------------------
          DO jk = 1, nlay_i
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zcpi = cpic + 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_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
@@ -326,5 +326,5 @@
 
          DO jk = 1, nlay_s
-            DO ji = 1, nidx
+            DO ji = 1, npti
                zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
             END DO
@@ -335,5 +335,5 @@
          !----------------------------
          IF ( ln_dqns_i ) THEN 
-            DO ji = 1, nidx
+            DO ji = 1, npti
                ! 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) - ztsub(ji) )
@@ -341,5 +341,5 @@
          ENDIF
 
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
          END DO
@@ -354,11 +354,11 @@
 
          !!ice interior terms (top equation has the same form as the others)
-         ztrid   (1:nidx,:,:) = 0._wp
-         zindterm(1:nidx,:)   = 0._wp
-         zindtbis(1:nidx,:)   = 0._wp
-         zdiagbis(1:nidx,:)   = 0._wp
+         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, nidx
+            DO ji = 1, npti
                jk = jm - nlay_s - 1
                ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
@@ -370,5 +370,5 @@
 
          jm =  nlay_s + nlay_i + 1
-         DO ji = 1, nidx
+         DO ji = 1, npti
             !!ice bottom term
             ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
@@ -380,5 +380,5 @@
 
 
-         DO ji = 1, nidx
+         DO ji = 1, npti
             !                               !---------------------!
             IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
@@ -497,5 +497,5 @@
          jm_maxt = 0
          jm_mint = nlay_i+5
-         DO ji = 1, nidx
+         DO ji = 1, npti
             jm_mint = MIN(jm_min(ji),jm_mint)
             jm_maxt = MAX(jm_max(ji),jm_maxt)
@@ -503,5 +503,5 @@
 
          DO jk = jm_mint+1, jm_maxt
-            DO ji = 1, nidx
+            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)
@@ -510,5 +510,5 @@
          END DO
 
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ! ice temperatures
             t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
@@ -516,5 +516,5 @@
 
          DO jm = nlay_i + nlay_s, nlay_s + 2, -1
-            DO ji = 1, nidx
+            DO ji = 1, npti
                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)
@@ -522,5 +522,5 @@
          END DO
 
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ! snow temperatures      
             IF( h_s_1d(ji) > 0._wp ) THEN
@@ -542,5 +542,5 @@
          ! zdti_max is a measure of error, it has to be under zdti_bnd
          zdti_max = 0._wp
-         DO ji = 1, nidx
+         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) ) )     
@@ -548,5 +548,5 @@
 
          DO jk  =  1, nlay_s
-            DO ji = 1, nidx
+            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) ) )
@@ -555,5 +555,5 @@
 
          DO jk  =  1, nlay_i
-            DO ji = 1, nidx
+            DO ji = 1, npti
                ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0 
                t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
@@ -576,5 +576,5 @@
       ! 10) Fluxes at the interfaces
       !-----------------------------
-      DO ji = 1, nidx
+      DO ji = 1, npti
          !                                ! surface ice conduction flux
          fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
@@ -589,5 +589,5 @@
       ! --- diagnose the change in non-solar flux due to surface temperature change --- !
       IF ( ln_dqns_i ) THEN
-         DO ji = 1, nidx
+         DO ji = 1, npti
             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
          END DO
@@ -597,5 +597,5 @@
       !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
       !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-      DO ji = 1, nidx
+      DO ji = 1, npti
          zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  &
             &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s )
@@ -614,5 +614,5 @@
 
       ! --- Diagnostics SIMIP --- !
-      DO ji = 1, nidx
+      DO ji = 1, npti
          !--- Conduction fluxes (positive downwards)
          diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
@@ -645,5 +645,5 @@
       !
       DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = 1, nidx
+         DO ji = 1, npti
             ztmelts      = - tmut  * sz_i_1d(ji,jk)
             t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
@@ -655,5 +655,5 @@
       END DO
       DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = 1, nidx
+         DO ji = 1, npti
             e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
          END DO
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8564)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8565)
@@ -378,5 +378,5 @@
       CASE( 1 )       !  constant salinity in time and space  !
          !            !---------------------------------------!
-         sz_i_1d(1:nidx,:) = rn_icesal
+         sz_i_1d(1:npti,:) = rn_icesal
          !
          !            !---------------------------------------------!
@@ -387,10 +387,10 @@
          !
          !                                      ! Slope of the linear profile 
-         WHERE( h_i_1d(1:nidx) > epsi20 )   ;   z_slope_s(1:nidx) = 2._wp * s_i_1d(1:nidx) / h_i_1d(1:nidx)
-         ELSEWHERE                           ;   z_slope_s(1:nidx) = 0._wp
+         WHERE( h_i_1d(1:npti) > epsi20 )   ;   z_slope_s(1:npti) = 2._wp * s_i_1d(1:npti) / h_i_1d(1:npti)
+         ELSEWHERE                           ;   z_slope_s(1:npti) = 0._wp
          END WHERE
          
          z1_dS = 1._wp / ( zsi1 - zsi0 )
-         DO ji = 1, nidx
+         DO ji = 1, npti
             zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - s_i_1d(ji) ) * z1_dS , 1._wp  )  )
             !                             ! force a constant profile when SSS too low (Baltic Sea)
@@ -400,5 +400,5 @@
          ! Computation of the profile
          DO jk = 1, nlay_i
-            DO ji = 1, nidx
+            DO ji = 1, npti
                !                          ! linear profile with 0 surface value
                zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * h_i_1d(ji) * r1_nlay_i
@@ -414,5 +414,5 @@
          !            !-------------------------------------------!                                   (mean = 2.30)
          !
-         s_i_1d(1:nidx) = 2.30_wp
+         s_i_1d(1:npti) = 2.30_wp
          !
 !!gm cf remark in ice_var_salprof routine, CASE( 3 )
@@ -420,5 +420,5 @@
             zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
             zsal =  1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) )
-            DO ji = 1, nidx
+            DO ji = 1, npti
                sz_i_1d(ji,jk) = zsal
             END DO
