Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4902)
@@ -105,5 +105,5 @@
    !! ** Global variables                                                 |
    !!-------------|-------------|---------------------------------|-------|
-   !! a_i         | a_i_b       |    Ice concentration            |       |
+   !! a_i         | a_i_1d      |    Ice concentration            |       |
    !! v_i         |      -      |    Ice volume per unit area     | m     |
    !! v_s         |      -      |    Snow volume per unit area    | m     |
@@ -111,7 +111,7 @@
    !! oa_i        !      -      !    Sea ice areal age content    | day   |
    !! e_i         !      -      !    Ice enthalpy                 | 10^9 J| 
-   !!      -      ! q_i_b       !    Ice enthalpy per unit vol.   | J/m3  | 
+   !!      -      ! q_i_1d      !    Ice enthalpy per unit vol.   | J/m3  | 
    !! e_s         !      -      !    Snow enthalpy                | 10^9 J| 
-   !!      -      ! q_s_b       !    Snow enthalpy per unit vol.  | J/m3  | 
+   !!      -      ! q_s_1d      !    Snow enthalpy per unit vol.  | J/m3  | 
    !!                                                                     |
    !!-------------|-------------|---------------------------------|-------|
@@ -120,12 +120,12 @@
    !!-------------|-------------|---------------------------------|-------|
    !!                                                                     |
-   !! ht_i        | ht_i_b      |    Ice thickness                | m     |
-   !! ht_s        ! ht_s_b      |    Snow depth                   | m     |
-   !! sm_i        ! sm_i_b      |    Sea ice bulk salinity        ! ppt   |
-   !! s_i         ! s_i_b       |    Sea ice salinity profile     ! ppt   |
+   !! ht_i        | ht_i_1d     |    Ice thickness                | m     |
+   !! ht_s        ! ht_s_1d     |    Snow depth                   | m     |
+   !! sm_i        ! sm_i_1d     |    Sea ice bulk salinity        ! ppt   |
+   !! s_i         ! s_i_1d      |    Sea ice salinity profile     ! ppt   |
    !! o_i         !      -      |    Sea ice Age                  ! days  |
-   !! t_i         ! t_i_b       |    Sea ice temperature          ! K     |
-   !! t_s         ! t_s_b       |    Snow temperature             ! K     |
-   !! t_su        ! t_su_b      |    Sea ice surface temperature  ! K     |
+   !! t_i         ! t_i_1d      |    Sea ice temperature          ! K     |
+   !! t_s         ! t_s_1d      |    Snow temperature             ! K     |
+   !! t_su        ! t_su_1d     |    Sea ice surface temperature  ! K     |
    !!                                                                     |
    !! notes: the ice model only sees a bulk (i.e., vertically averaged)   |
@@ -142,5 +142,5 @@
    !! ***         Category-summed state variables (diagnostic)        *** |
    !! ******************************************************************* |
-   !! at_i        | at_i_b      |    Total ice concentration      |       |
+   !! at_i        | at_i_1d     |    Total ice concentration      |       |
    !! vt_i        |      -      |    Total ice vol. per unit area | m     |
    !! vt_s        |      -      |    Total snow vol. per unit ar. | m     |
@@ -176,6 +176,7 @@
    REAL(wp), PUBLIC ::   ecc              !: eccentricity of the elliptical yield curve
    REAL(wp), PUBLIC ::   ahi0             !: sea-ice hor. eddy diffusivity coeff. (m2/s)
-   REAL(wp), PUBLIC ::   telast           !: timescale for elastic waves (s) !SB
-   REAL(wp), PUBLIC ::   alphaevp         !: coeficient of the internal stresses !SB
+   REAL(wp), PUBLIC ::   telast           !: timescale for elastic waves (s)
+   REAL(wp), PUBLIC ::   relast           !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
+   REAL(wp), PUBLIC ::   alphaevp         !: coeficient of the internal stresses 
    REAL(wp), PUBLIC ::   unit_fac = 1.e+09_wp  !: conversion factor for ice / snow enthalpy
    REAL(wp), PUBLIC ::   hminrhg          !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity
@@ -247,16 +248,16 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: Variation of snow mass over 1 time step     [Kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice     !: Variation of ice mass over 1 time step      [Kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: Variation of snow mass over 1 time step due to sublimation [Kg/m2]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni     !: snow ice growth 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw     !: lateral ice growth 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog     !: bottom ice growth 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn     !: dynamical ice growth 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom     !: vertical bottom melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum     !: vertical surface melt
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res     !: production (growth+melt) due to limupdate
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange over 1 time step [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice over 1 time step [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow sublimation over 1 time step [kg/m2]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange over 1 time step [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -323,11 +324,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i          !: mean sea ice salinity averaged over all categories [PSU]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   at_i_typ     !: total area   contained in each ice type [m^2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   vt_i_typ     !: total volume contained in each ice type [m^3]
-
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s        !: Snow temperatures [K]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s        !: Snow ...      
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   e_i_cat    !: ! go to trash
       
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i        !: ice temperatures          [K]
@@ -350,9 +346,9 @@
    !! * Old values of global variables
    !!--------------------------------------------------------------------------
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   old_v_s, old_v_i               !: snow and ice volumes
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   old_a_i, old_smv_i, old_oa_i   !: ???
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   old_e_s                        !: snow heat content
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   old_e_i                        !: ice temperatures
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   old_u_ice, old_v_ice           !: ice velocity (gv6 and gv7)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b               !: snow and ice volumes
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, smv_i_b, oa_i_b     !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                      !: snow heat content
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                      !: ice temperatures
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b           !: ice velocity
       
 
@@ -377,14 +373,6 @@
    !! * Ice thickness distribution variables
    !!--------------------------------------------------------------------------
-   ! REMOVE
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   ice_types      !: Vector connecting types and categories
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ice_cat_bounds !: Matrix containing the integer upper and 
-   !                                                                       !  lower boundaries of ice thickness categories
-   ! REMOVE
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   ice_ncat_types !: nb of thickness categories in each ice type
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_max         !: Boundary of ice thickness categories in thickness space
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_mean        !: Mean ice thickness in catgories 
-   ! REMOVE
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hi_max_typ     !: Boundary of ice thickness categories in thickness space
 
    !!--------------------------------------------------------------------------
@@ -475,8 +463,8 @@
          &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) )
       ii = ii + 1
-      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , at_i_typ(jpi,jpj,jpm) ,                            &
-         &      e_s(jpi,jpj,nlay_s,jpl) , vt_i_typ(jpi,jpj,jpm) , e_i_cat(jpi,jpj,jpl) , STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( t_i(jpi,jpj,jkmax,jpl) , e_i(jpi,jpj,jkmax,jpl) , s_i(jpi,jpj,jkmax,jpl) , STAT=ierr(ii) )
+      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) ,                            &
+         &      e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
+      ii = ii + 1
+      ALLOCATE( t_i(jpi,jpj,nlay_i+1,jpl) , e_i(jpi,jpj,nlay_i+1,jpl) , s_i(jpi,jpj,nlay_i+1,jpl) , STAT=ierr(ii) )
 
       ! * Moments for advection
@@ -494,13 +482,13 @@
          &      STAT=ierr(ii) )
       ii = ii + 1
-      ALLOCATE( sxe (jpi,jpj,jkmax,jpl) , sye (jpi,jpj,jkmax,jpl) , sxxe(jpi,jpj,jkmax,jpl) ,     &
-         &      syye(jpi,jpj,jkmax,jpl) , sxye(jpi,jpj,jkmax,jpl)                           , STAT=ierr(ii) )
+      ALLOCATE( sxe (jpi,jpj,nlay_i+1,jpl) , sye (jpi,jpj,nlay_i+1,jpl) , sxxe(jpi,jpj,nlay_i+1,jpl) ,     &
+         &      syye(jpi,jpj,nlay_i+1,jpl) , sxye(jpi,jpj,nlay_i+1,jpl)                           , STAT=ierr(ii) )
 
       ! * Old values of global variables
       ii = ii + 1
-      ALLOCATE( old_v_s  (jpi,jpj,jpl) , old_v_i  (jpi,jpj,jpl) , old_e_s(jpi,jpj,nlay_s,jpl) ,     &
-         &      old_a_i  (jpi,jpj,jpl) , old_smv_i(jpi,jpj,jpl) , old_e_i(jpi,jpj,jkmax ,jpl) ,     &
-         &      old_oa_i (jpi,jpj,jpl)                                                        ,     &
-         &      old_u_ice(jpi,jpj)     , old_v_ice(jpi,jpj)                                   , STAT=ierr(ii) )
+      ALLOCATE( v_s_b  (jpi,jpj,jpl) , v_i_b  (jpi,jpj,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,     &
+         &      a_i_b  (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i+1 ,jpl) ,     &
+         &      oa_i_b (jpi,jpj,jpl)                                                        ,     &
+         &      u_ice_b(jpi,jpj)     , v_ice_b(jpi,jpj)                                   , STAT=ierr(ii) )
 
       ! * Increment of global variables
@@ -512,17 +500,16 @@
          &     STAT=ierr(ii) )
       ii = ii + 1
-      ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,jkmax,jpl) , d_u_ice_dyn(jpi,jpj) ,     &
-         &      d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,jkmax,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) )
+      ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i+1,jpl) , d_u_ice_dyn(jpi,jpj) ,     &
+         &      d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i+1,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) )
       
       ! * Ice thickness distribution variables
       ii = ii + 1
-      ALLOCATE( ice_types(jpl) , ice_cat_bounds(jpm,2) , ice_ncat_types  (jpm) ,     &
-         &      hi_max (0:jpl) , hi_mean(jpl)          , hi_max_typ(0:jpl,jpm) , STAT=ierr(ii) )
+      ALLOCATE( hi_max(0:jpl), hi_mean(jpl),  STAT=ierr(ii) )
 
       ! * Ice diagnostics
       ii = ii + 1
-      ALLOCATE( dv_dt_thd(jpi,jpj,jpl) ,     &
-         &      izero    (jpi,jpj,jpl)  , diag_trp_vi(jpi,jpj) , diag_trp_vs(jpi,jpj), diag_trp_ei(jpi,jpj), diag_trp_es(jpi,jpj),     & 
-         &      diag_heat_dhc(jpi,jpj) ,  STAT=ierr(ii) )
+      ALLOCATE( dv_dt_thd(jpi,jpj,jpl), izero (jpi,jpj,jpl),    &
+         &      diag_trp_vi(jpi,jpj), diag_trp_vs  (jpi,jpj), diag_trp_ei(jpi,jpj),   & 
+         &      diag_trp_es(jpi,jpj), diag_heat_dhc(jpi,jpj),  STAT=ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4902)
@@ -66,6 +66,6 @@
       !
       !                                ! adequation jpk versus ice/snow layers/categories
-      IF( jpl   > jpk  .OR.  jpm    > jpk .OR.                                    &
-          jkmax > jpk  .OR.  nlay_s > jpk      )   CALL ctl_stop( 'STOP',         &
+      IF( jpl > jpk .OR. (nlay_i+1) > jpk .OR. nlay_s > jpk )   &
+         &      CALL ctl_stop( 'STOP',                     &
          &     'ice_init: the 3rd dimension of workspace arrays is too small.',   &
          &     'use more ocean levels or less ice/snow layers/categories.' )
@@ -174,5 +174,5 @@
       !!              limistate (only) and is changed to 99 m in ice_init
       !!------------------------------------------------------------------
-      INTEGER  ::   jl, jm               ! dummy loop index
+      INTEGER  ::   jl                   ! dummy loop index
       REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
       !!------------------------------------------------------------------
@@ -185,33 +185,6 @@
       ! 1) Ice thickness distribution parameters initialization    
       !------------------------------------------------------------------------------!
-
-      !- Types boundaries (integer)
-      !----------------------------
-      ice_cat_bounds(1,1) = 1
-      ice_cat_bounds(1,2) = jpl
-
-      !- Number of ice thickness categories in each ice type
-      DO jm = 1, jpm
-         ice_ncat_types(jm) = ice_cat_bounds(jm,2) - ice_cat_bounds(jm,1) + 1 
-      END DO
-
-      !- Make the correspondence between thickness categories and ice types
-      !---------------------------------------------------------------------
-      DO jm = 1, jpm       !over types
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) !over thickness categories
-            ice_types(jl) = jm
-         END DO
-      END DO
-
       IF(lwp) THEN  
-         WRITE(numout,*) ' Number of ice types jpm =      ', jpm
          WRITE(numout,*) ' Number of ice categories jpl = ', jpl
-         DO jm = 1, jpm
-            WRITE(numout,*) ' Ice type ', jm
-            WRITE(numout,*) ' Number of thickness categories ', ice_ncat_types(jm)
-            WRITE(numout,*) ' Thickness category boundaries  ', ice_cat_bounds(jm,1:2)
-         END DO
-         WRITE(numout,*) 'Ice type vector', ice_types(1:jpl)
-         WRITE(numout,*)
       ENDIF
 
@@ -219,33 +192,17 @@
       !----------------------------------
       hi_max(:) = 0._wp
-      hi_max_typ(:,:) = 0._wp
-
-      !- Type 1 - undeformed ice
-      zc1 =  3._wp / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp )
+
+      zc1 =  3._wp / REAL( jpl, wp )
       zc2 = 10._wp * zc1
       zc3 =  3._wp
 
-      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-         zx1 = REAL( jl-1 , wp ) / REAL( ice_cat_bounds(1,2) - ice_cat_bounds(1,1) + 1 , wp )
+      DO jl = 1, jpl
+         zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
          hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
       END DO
 
-      !- Fill in the hi_max_typ vector, useful in other circumstances
-      ! Tricky trick: hi_max_typ is actually not used in the code and will be removed in a
-      ! next flyspray at this time, the tricky trick will also be removed (Martin, march 08)
-      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-         hi_max_typ(jl,1) = hi_max(jl)
-      END DO
-
-      IF(lwp) WRITE(numout,*) ' Thickness category boundaries independently of ice type '
+      IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
       IF(lwp) WRITE(numout,*) ' hi_max ', hi_max(0:jpl)
 
-      IF(lwp) WRITE(numout,*) ' Thickness category boundaries inside ice types '
-      IF(lwp) THEN 
-         DO jm = 1, jpm
-            WRITE(numout,*) ' Type number ', jm
-            WRITE(numout,*) ' hi_max_typ : ', hi_max_typ(0:ice_ncat_types(jm),jm)
-         END DO
-      ENDIF
       !
       DO jl = 1, jpl
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 4902)
@@ -73,8 +73,8 @@
       !! ** Method  : Arithmetics
       !!---------------------------------------------------------------------
-      INTEGER                               , INTENT(in   ) ::   ksum   !: number of categories
-      INTEGER                               , INTENT(in   ) ::   klay   !: number of vertical layers
-      REAL(wp), DIMENSION(jpi,jpj,jkmax,jpl), INTENT(in   ) ::   pin   !: input field
-      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) ::   pout   !: output field
+      INTEGER                                  , INTENT(in   ) ::   ksum   !: number of categories
+      INTEGER                                  , INTENT(in   ) ::   klay   !: number of vertical layers
+      REAL(wp), DIMENSION(jpi,jpj,nlay_i+1,jpl), INTENT(in   ) ::   pin   !: input field
+      REAL(wp), DIMENSION(jpi,jpj)             , INTENT(  out) ::   pout   !: output field
       !
       INTEGER ::   jk, jl   ! dummy loop indices
@@ -175,18 +175,24 @@
          zei_b  = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
          zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
-            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) )
+            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
+            &             ) * area(:,:) * tms(:,:) )
          zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) )
+            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
+            &                 ) * area(:,:) * tms(:,:) )
          zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) )
+            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &                  ) * area(:,:) / unit_fac * tms(:,:) )
 
       ELSEIF( icount == 1 ) THEN
 
          zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zfs_b
+            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
+            &                ) * area(:,:) * tms(:,:) ) - zfs_b
          zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
-            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) - zfw_b
+            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
+            &                ) * area(:,:) * tms(:,:) ) - zfw_b
          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(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b
+            &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &                ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b
  
          zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 4902)
@@ -59,5 +59,6 @@
       !!
       REAL(dp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
-      REAL(dp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni, zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn 
+      REAL(dp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   &
+      &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn 
       REAL(dp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn
       REAL(dp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub  
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 4902)
@@ -83,6 +83,6 @@
          IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:)
-         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:)
+         u_ice_b(:,:) = u_ice(:,:) * tmu(:,:)
+         v_ice_b(:,:) = v_ice(:,:) * tmv(:,:)
 
          ! Rheology (ice dynamics)
@@ -243,5 +243,5 @@
       NAMELIST/namicedyn/ epsd, om, cw, angvg, pstar,   &
          &                c_rhg, creepl, ecc, ahi0,     &
-         &                nevp, telast, alphaevp, hminrhg
+         &                nevp, relast, alphaevp, hminrhg
       !!-------------------------------------------------------------------
 
@@ -269,5 +269,5 @@
          WRITE(numout,*) '   horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
          WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp
-         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
+         WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast
          WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
          WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg
@@ -287,4 +287,7 @@
       pstarh = pstar * 0.5_wp
 
+      ! elastic damping
+      telast = relast * rdt_ice
+
       !  Diffusion coefficients.
       ahiu(:,:) = ahi0 * umask(:,:,1)
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4902)
@@ -128,6 +128,8 @@
       DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
          DO ji = 1, jpi
-            IF( ( tsn(ji,jj,1,jp_tem)  - ( t_bo(ji,jj) - rt0 ) ) * tms(ji,jj) >= thres_sst ) THEN  ; zswitch(ji,jj) = 0._wp * tms(ji,jj)    ! no ice
-            ELSE                                                                                   ; zswitch(ji,jj) = 1._wp * tms(ji,jj)    !    ice
+            IF( ( tsn(ji,jj,1,jp_tem)  - ( t_bo(ji,jj) - rt0 ) ) * tms(ji,jj) >= thres_sst ) THEN 
+               zswitch(ji,jj) = 0._wp * tms(ji,jj)    ! no ice
+            ELSE                                                                                   
+               zswitch(ji,jj) = 1._wp * tms(ji,jj)    !    ice
             ENDIF
          END DO
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 4902)
@@ -692,5 +692,5 @@
 
       IF( partfun_swi == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
-         DO jl = 0, ice_cat_bounds(1,2)       ! only undeformed ice participates
+         DO jl = 0, jpl    
             DO jj = 1, jpj 
                DO ji = 1, jpi
@@ -715,5 +715,5 @@
             Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
          END DO !jl
-         DO jl = 0, ice_cat_bounds(1,2)
+         DO jl = 0, jpl
              athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
          END DO
@@ -897,6 +897,6 @@
       CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
       CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
-      CALL wrk_alloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
-      CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )
+      CALL wrk_alloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw )
+      CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
 
       ! Conservation check
@@ -1193,5 +1193,5 @@
          !-------------------------------------------------------------------------------
          !        jl1 looping 1-jpl
-         DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
+         DO jl2  = 1, jpl 
             ! over categories to which ridged ice is transferred
 !CDIR NODEP
@@ -1238,5 +1238,5 @@
          END DO                 ! jl2 (new ridges)            
 
-         DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
+         DO jl2 = 1, jpl 
 
 !CDIR NODEP
@@ -1302,6 +1302,6 @@
       CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
       CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
-      CALL wrk_dealloc( jpi, jpj, jkmax,      eirft, erdg1, erdg2, ersw )
-      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, eicen_init )
+      CALL wrk_dealloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw )
+      CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init )
       !
    END SUBROUTINE lim_itd_me_ridgeshift
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 4902)
@@ -6,5 +6,5 @@
    !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
    !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3
-   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age and types
+   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age
    !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked
    !!----------------------------------------------------------------------
@@ -66,5 +66,5 @@
       INTEGER, INTENT(in) ::   kt   ! time step index
       !
-      INTEGER ::   ji,jj, jk, jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index         
+      INTEGER ::   ji, jj, jk, jl   ! dummy loop index         
       !
       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
@@ -86,9 +86,5 @@
       ! Given thermodynamic growth rates, transport ice between
       ! thickness categories.
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF( ice_ncat_types(jm) > 1 )   CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt )
-      END DO
+      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt )
       !
       CALL lim_var_glo2eqv    ! only for info
@@ -123,10 +119,10 @@
             CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ')
             CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ')
-            DO ja = 1, nlay_i
+            DO jk = 1, nlay_i
                CALL prt_ctl_info(' ')
-               CALL prt_ctl_info(' - Layer : ', ivar1=ja)
+               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
                CALL prt_ctl_info('   ~~~~~~~')
-               CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
-               CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
+               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ')
+               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ')
             END DO
          END DO
@@ -140,5 +136,5 @@
    !
 
-   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt )
+   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_th_rem ***
@@ -153,5 +149,4 @@
       INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
       INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
-      INTEGER , INTENT (in) ::   ntyp    ! Number of the type used
       INTEGER , INTENT (in) ::   kt      ! Ocean time step 
       !
@@ -171,5 +166,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hL          ! left boundary for the ITD for each thickness
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_o     ! old ice thickness
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume
@@ -189,5 +184,5 @@
       CALL wrk_alloc( jpi,jpj, zremap_flag )    ! integer
       CALL wrk_alloc( jpi,jpj,jpl-1, zdonor )   ! integer
-      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_o, dummy_es )
+      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
       CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )   
       CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
@@ -218,5 +213,4 @@
          WRITE(numout,*) ' klbnd :       ', klbnd
          WRITE(numout,*) ' kubnd :       ', kubnd
-         WRITE(numout,*) ' ntyp  :       ', ntyp 
       ENDIF
 
@@ -227,7 +221,7 @@
                zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes
                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb
-               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes
-               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb
-               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 
+               zindb             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes
+               zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * zindb
+               IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 
             END DO
          END DO
@@ -274,11 +268,11 @@
             !
             zhbnew(ii,ij,jl) = hi_max(jl)
-            IF ( old_a_i(ii,ij,jl) > epsi10 .AND. old_a_i(ii,ij,jl+1) > epsi10 ) THEN
+            IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN
                !interpolate between adjacent category growth rates
-               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) )
-               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) )
-            ELSEIF ( old_a_i(ii,ij,jl) > epsi10) THEN
+               zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) )
+               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) )
+            ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN
                zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)
-            ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN
+            ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN
                zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)
             ENDIF
@@ -321,6 +315,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme
-            zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er
+            zhb0(ji,jj) = hi_max(0) ! 0eme
+            zhb1(ji,jj) = hi_max(1) ! 1er
 
             zhbnew(ji,jj,klbnd-1) = 0._wp
@@ -343,5 +337,5 @@
       !-----------------------------------------------------------------------------------------------
       !- 7.1 g(h) for category 1 at start of time step
-      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd),         &
+      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd),         &
          &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   &
          &                  hR(:,:,klbnd), zremap_flag )
@@ -368,5 +362,5 @@
                   ! Constrain new thickness <= ht_i
                   zdamax = a_i(ii,ij,klbnd) * & 
-                     (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0
+                     (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0
                   !ice area lost due to melting of thin ice
                   zda0   = MIN(zda0, zdamax)
@@ -382,10 +376,7 @@
             ELSE ! if ice accretion
                ! ji, a_i > epsi10; zdh0 > 0
-               IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 
+               zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 
                ! zhbnew was 0, and is shifted to the right to account for thin ice
                ! growth in openwater (F0 = f1)
-               IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0 
-               ! in other types there is
-               ! no open water growth (F0 = 0)
             ENDIF ! zdh0 
 
@@ -493,5 +484,5 @@
       CALL wrk_dealloc( jpi,jpj, zremap_flag )    ! integer
       CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )   ! integer
-      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_o, dummy_es )
+      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es )
       CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )   
       CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )   
@@ -839,5 +830,5 @@
    
 
-   SUBROUTINE lim_itd_th_reb( klbnd, kubnd, ntyp )
+   SUBROUTINE lim_itd_th_reb( klbnd, kubnd )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_th_reb ***
@@ -849,5 +840,4 @@
       INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point
       INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied
-      INTEGER , INTENT (in) ::   ntyp    ! number of the ice type involved in the rebinning process
       !
       INTEGER ::   ji,jj, jl   ! dummy loop indices
@@ -889,12 +879,12 @@
 
       !------------------------------------------------------------------------------
-      ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd)
+      ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd)
       !------------------------------------------------------------------------------
       DO jj = 1, jpj 
          DO ji = 1, jpi 
             IF( a_i(ji,jj,klbnd) > epsi10 ) THEN
-               IF( ht_i(ji,jj,klbnd) <= hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN
-                  a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp) 
-                  ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp)
+               IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN
+                  a_i(ji,jj,klbnd)  = v_i(ji,jj,klbnd) / hi_max(0) 
+                  ht_i(ji,jj,klbnd) = hi_max(0)
                ENDIF
             ENDIF
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4902)
@@ -118,9 +118,9 @@
       CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface
       CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! non-solar flux at ocean surface
-      CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux at ice surface
-      CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  ! non-solar flux at ice surface
-      CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux transmitted thru ice
+      CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )  !     solar flux at ice surface
+      CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )  ! non-solar flux at ice surface
+      CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )  !     solar flux transmitted thru ice
       CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )  
-      CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) )
+      CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * a_i_b(:,:,:), dim=3 ) )
 
       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
@@ -139,5 +139,5 @@
                zfcm1 = qsr_tot(ji,jj)
                DO jl = 1, jpl
-                  zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * old_a_i(ji,jj,jl)
+                  zfcm1 = zfcm1 + ( ftr_ice(ji,jj,jl) - qsr_ice(ji,jj,jl) ) * a_i_b(ji,jj,jl)
                END DO
             ELSE
@@ -145,5 +145,5 @@
                zfcm1   = pfrld(ji,jj) * qsr(ji,jj)
                DO jl = 1, jpl
-                  zfcm1   = zfcm1 + old_a_i(ji,jj,jl) * ftr_ice(ji,jj,jl)
+                  zfcm1   = zfcm1 + a_i_b(ji,jj,jl) * ftr_ice(ji,jj,jl)
                END DO
             ENDIF
@@ -181,9 +181,10 @@
 
             ! mass flux from ice/ocean
-            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj)
+            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
+                           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj)
 
             ! mass flux at the ocean/ice interface
             fmmflx(ji,jj) = - wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model
-            emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
+            emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj)       ! mass flux + F/M mass flux (always ice/ocean mass exchange)
             
          END DO
@@ -193,5 +194,6 @@
       !      salt flux at the ocean surface      !
       !------------------------------------------!
-      sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:)
+      sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:)   &
+         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:)
 
       !-------------------------------------------------------------!
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4902)
@@ -22,5 +22,5 @@
    USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain variables
-   USE oce     , ONLY : fraqsr_1lev
+   USE oce     , ONLY : fraqsr_1lev 
    USE ice            ! LIM: sea-ice variables
    USE par_ice        ! LIM: sea-ice parameters
@@ -153,6 +153,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * old_a_i(ji,jj,jl)
-                  zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * old_a_i(ji,jj,jl)
+                  zqsr(ji,jj) = zqsr(ji,jj) - qsr_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
+                  zqns(ji,jj) = zqns(ji,jj) - qns_ice(ji,jj,jl) * a_i_b(ji,jj,jl)
                END DO
             END DO
@@ -282,19 +282,19 @@
             !-------------------------
 
-            CALL tab_2d_1d( nbpb, at_i_b     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, a_i_b      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, ht_i_b     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, ht_s_b     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-
-            CALL tab_2d_1d( nbpb, t_su_b     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, sm_i_b     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, at_i_1d     (1:nbpb), at_i            , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, a_i_1d      (1:nbpb), a_i(:,:,jl)     , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, ht_i_1d     (1:nbpb), ht_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, ht_s_1d     (1:nbpb), ht_s(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
+
+            CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
             DO jk = 1, nlay_s
-               CALL tab_2d_1d( nbpb, t_s_b(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-               CALL tab_2d_1d( nbpb, q_s_b(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+               CALL tab_2d_1d( nbpb, t_s_1d(1:nbpb,jk), t_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+               CALL tab_2d_1d( nbpb, q_s_1d(1:nbpb,jk), e_s(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
             END DO
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nbpb, t_i_b(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-               CALL tab_2d_1d( nbpb, q_i_b(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
-               CALL tab_2d_1d( nbpb, s_i_b(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+               CALL tab_2d_1d( nbpb, t_i_1d(1:nbpb,jk), t_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+               CALL tab_2d_1d( nbpb, q_i_1d(1:nbpb,jk), e_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
+               CALL tab_2d_1d( nbpb, s_i_1d(1:nbpb,jk), s_i(:,:,jk,jl)  , jpi, jpj, npb(1:nbpb) )
             END DO
 
@@ -310,5 +310,5 @@
             ENDIF
             CALL tab_2d_1d( nbpb, dqns_ice_1d(1:nbpb), dqns_ice(:,:,jl), jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, t_bo_b     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, t_bo_1d     (1:nbpb), t_bo            , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, sprecip_1d (1:nbpb), sprecip         , jpi, jpj, npb(1:nbpb) ) 
             CALL tab_2d_1d( nbpb, fhtur_1d   (1:nbpb), fhtur           , jpi, jpj, npb(1:nbpb) )
@@ -361,5 +361,5 @@
 
             ! --- Ice enthalpy remapping --- !
-            CALL lim_thd_ent( 1, nbpb, q_i_b(1:nbpb,:) ) 
+            CALL lim_thd_ent( 1, nbpb, q_i_1d(1:nbpb,:) ) 
                                             
             !---------------------------------!
@@ -377,18 +377,18 @@
             !--------------------------------
 
-               CALL tab_1d_2d( nbpb, at_i          , npb, at_i_b    (1:nbpb)   , jpi, jpj )
-               CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_b    (1:nbpb)   , jpi, jpj )
-               CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_b    (1:nbpb)   , jpi, jpj )
-               CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_b     (1:nbpb)   , jpi, jpj )
-               CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_b    (1:nbpb)   , jpi, jpj )
-               CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_b    (1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, at_i          , npb, at_i_1d    (1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, ht_i(:,:,jl)  , npb, ht_i_1d    (1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, ht_s(:,:,jl)  , npb, ht_s_1d    (1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj )
+               CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj )
             DO jk = 1, nlay_s
-               CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_b     (1:nbpb,jk), jpi, jpj)
-               CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_b     (1:nbpb,jk), jpi, jpj)
+               CALL tab_1d_2d( nbpb, t_s(:,:,jk,jl), npb, t_s_1d     (1:nbpb,jk), jpi, jpj)
+               CALL tab_1d_2d( nbpb, e_s(:,:,jk,jl), npb, q_s_1d     (1:nbpb,jk), jpi, jpj)
             END DO
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_b     (1:nbpb,jk), jpi, jpj)
-               CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_b     (1:nbpb,jk), jpi, jpj)
-               CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_b     (1:nbpb,jk), jpi, jpj)
+               CALL tab_1d_2d( nbpb, t_i(:,:,jk,jl), npb, t_i_1d     (1:nbpb,jk), jpi, jpj)
+               CALL tab_1d_2d( nbpb, e_i(:,:,jk,jl), npb, q_i_1d     (1:nbpb,jk), jpi, jpj)
+               CALL tab_1d_2d( nbpb, s_i(:,:,jk,jl), npb, s_i_1d     (1:nbpb,jk), jpi, jpj)
             END DO
                CALL tab_1d_2d( nbpb, qlead         , npb, qlead_1d  (1:nbpb)   , jpi, jpj )
@@ -532,15 +532,15 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
-            ztmelts       =  -tmut * s_i_b(ji,jk) + rtt
+            ztmelts       =  -tmut * s_i_1d(ji,jk) + rtt
             ! Conversion q(S,T) -> T (second order equation)
             zaaa          =  cpic
-            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus
+            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus
             zccc          =  lfus * ( ztmelts - rtt )
             zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
-            t_i_b(ji,jk)  =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
+            t_i_1d(ji,jk)  =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa )
             
             ! mask temperature
-            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 
-            t_i_b(ji,jk) =  zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt
+            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
+            t_i_1d(ji,jk) =  zswitch * t_i_1d(ji,jk) + ( 1._wp - zswitch ) * rtt
          END DO 
       END DO 
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4902)
@@ -129,5 +129,5 @@
       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
       CALL wrk_alloc( jpij, zintermelt )
-      CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )
+      CALL wrk_alloc( jpij, nlay_i+1, zdeltah, zh_i )
       CALL wrk_alloc( jpij, icount )
       
@@ -155,6 +155,6 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
-            h_i_old (ji,jk) = ht_i_b(ji) / REAL( nlay_i )
-            qh_i_old(ji,jk) = q_i_b(ji,jk) * h_i_old(ji,jk)
+            h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i )
+            qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
          ENDDO
       ENDDO
@@ -165,11 +165,11 @@
       !
       DO ji = kideb, kiut
-         zinda      = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
-         ztmelts    = zinda * rtt + ( 1._wp - zinda ) * rtt
+         zinda         = 1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )
+         ztmelts       = zinda * rtt + ( 1._wp - zinda ) * rtt
 
          zfdum      = 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) 
 
-         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
+         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO
@@ -181,13 +181,13 @@
       !------------------------------------------------------------------------------!
       DO ji = kideb, kiut
-         IF( t_s_b(ji,1) > rtt ) THEN !!! Internal melting
+         IF( t_s_1d(ji,1) > rtt ) THEN !!! Internal melting
             ! Contribution to heat flux to the ocean [W.m-2], < 0  
-            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
+            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! Contribution to mass flux
-            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
+            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! updates
-            ht_s_b(ji)   = 0._wp
-            q_s_b (ji,1) = 0._wp
-            t_s_b (ji,1) = rtt
+            ht_s_1d(ji)   = 0._wp
+            q_s_1d (ji,1) = 0._wp
+            t_s_1d (ji,1) = rtt
          END IF
       END DO
@@ -198,10 +198,10 @@
       !
       DO ji = kideb, kiut     
-         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
+         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
       END DO
       !
       DO jk = 1, nlay_s
          DO ji = kideb, kiut
-            zqh_s(ji) =  zqh_s(ji) + q_s_b(ji,jk) * zh_s(ji)
+            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji)
          END DO
       END DO
@@ -209,6 +209,6 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
-            zh_i(ji,jk) = ht_i_b(ji) / REAL( nlay_i )
-            zqh_i(ji)   = zqh_i(ji) + q_i_b(ji,jk) * zh_i(ji,jk)
+            zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i )
+            zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
          END DO
       END DO
@@ -237,5 +237,5 @@
          !-----------
          ! thickness change
-         zcoeff = ( 1._wp - ( 1._wp - at_i_b(ji) )**betas ) / at_i_b(ji) 
+         zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) 
          zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
          ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
@@ -243,9 +243,9 @@
          IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
          ! heat flux from snow precip (>0, W.m-2)
-         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
+         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
          ! mass flux, <0
-         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_b(ji) * zdh_s_pre(ji) * r1_rdtice
+         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice
          ! update thickness
-         ht_s_b    (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )
+         ht_s_1d    (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )
 
          !---------------------
@@ -258,12 +258,12 @@
          zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 
          ! heat used to melt snow (W.m-2, >0)
-         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
+         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
          ! snow melting only = water into the ocean (then without snow precip), >0
-         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice
+         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice
          
          ! updates available heat + thickness
          zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )      
-         ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )
-         zh_s  (ji) = ht_s_b(ji) / REAL( nlay_s )
+         ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) )
+         zh_s  (ji) = ht_s_1d(ji) / REAL( nlay_s )
 
          ENDIF
@@ -275,17 +275,17 @@
          DO ji = kideb, kiut
             ! thickness change
-            zindh            = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) ) 
-            zindq            = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_b(ji,jk) + epsi20 ) ) 
-            zdeltah  (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_b(ji,jk), epsi20 )
+            zindh            = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 
+            zindq            = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) + epsi20 ) ) 
+            zdeltah  (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 )
             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
             zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)    
             ! heat used to melt snow(W.m-2, >0)
-            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice 
+            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
 
             ! updates available heat + thickness
-            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_b(ji,jk) )
-            ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )
+            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
+            ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
 
          END DO
@@ -304,15 +304,15 @@
          ! forced  mode: snow thickness change due to sublimation
          DO ji = kideb, kiut
-            zdh_s_sub(ji)  =  MAX( - ht_s_b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
+            zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
             ! Heat flux by sublimation [W.m-2], < 0
             !      sublimate first snow that had fallen, then pre-existing snow
             zcoeff         =      ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) )   * zqprec(ji) +   &
-               &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) )  &
-               &  * a_i_b(ji) * r1_rdtice
+               &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) )  &
+               &  * a_i_1d(ji) * r1_rdtice
             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff
             ! Mass flux by sublimation
-            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
+            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice
             ! new snow thickness
-            ht_s_b(ji)     =  MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
+            ht_s_1d(ji)     =  MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) )
          END DO
       ENDIF
@@ -321,5 +321,5 @@
       DO ji = kideb, kiut
          dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
-         zh_s(ji)       = ht_s_b(ji) / REAL( nlay_s )
+         zh_s(ji)       = ht_s_1d(ji) / REAL( nlay_s )
       END DO ! ji
 
@@ -331,9 +331,9 @@
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
-            zindh  =  MAX(  0._wp , SIGN( 1._wp, - ht_s_b(ji) + epsi20 )  )
-            q_s_b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) *             &
+            zindh  =  MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) + epsi20 )  )
+            q_s_1d(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_1d(ji), epsi20 ) *             &
               &            ( (   MAX( 0._wp, dh_s_tot(ji) )              ) * zqprec(ji) +  &
-              &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )
-            zq_s(ji)     =  zq_s(ji) + q_s_b(ji,jk)
+              &              ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) )
+            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
          END DO
       END DO
@@ -345,7 +345,7 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut 
-            zEi            = - q_i_b(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]
-
-            ztmelts        = - tmut * s_i_b(ji,jk) + rtt           ! Melting point of layer k [K]
+            zEi            = - q_i_1d(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]
+
+            ztmelts        = - tmut * s_i_1d(ji,jk) + rtt           ! Melting point of layer k [K]
 
             zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
@@ -367,15 +367,15 @@
             zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
 
-            ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
-            sfx_sum_1d(ji)   = sfx_sum_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
+            ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
+            sfx_sum_1d(ji)   = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
 
             ! Contribution to heat flux [W.m-2], < 0
-            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
 
             ! Total heat flux used in this process [W.m-2], > 0  
-            hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
+            hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
 
             ! Contribution to mass flux
-            wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+            wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
            
             ! record which layers have disappeared (for bottom melting) 
@@ -387,5 +387,5 @@
 
             ! update heat content (J.m-2) and layer thickness
-            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
+            qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
             h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
          END DO
@@ -393,5 +393,5 @@
       ! update ice thickness
       DO ji = kideb, kiut
-         ht_i_b(ji) =  MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )
+         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
       END DO
 
@@ -423,5 +423,5 @@
       !clem debug. Just to be sure that enthalpy at nlay_i+1 is null
       DO ji = kideb, kiut
-         q_i_b(ji,nlay_i+1) = 0._wp
+         q_i_1d(ji,nlay_i+1) = 0._wp
       END DO
 
@@ -445,9 +445,9 @@
 
                s_i_new(ji)        = zswitch_sal * zfracs * sss_m(ii,ij)  &  ! New ice salinity
-                                  + ( 1. - zswitch_sal ) * sm_i_b(ji) 
+                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
                ! New ice growth
                ztmelts            = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
 
-               zt_i_new           = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
+               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
                
                zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
@@ -455,5 +455,5 @@
                   &               + rcp  * ( ztmelts-rtt )          
 
-               zEw                = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
+               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
 
                zdE                = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
@@ -461,5 +461,5 @@
                dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
 
-               q_i_b(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
+               q_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
                
             ENDIF ! fc_bo_i
@@ -476,5 +476,5 @@
             ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K)
             
-            zt_i_new       = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
+            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
             
             zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
@@ -482,22 +482,22 @@
                &               + rcp  * ( ztmelts-rtt )          
             
-            zEw            = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
+            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
             
             zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0)
             
             ! Contribution to heat flux to the ocean [W.m-2], >0  
-            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
 
             ! Total heat flux used in this process [W.m-2], <0  
-            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
+            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
             
             ! Contribution to salt flux, <0
-            sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
+            sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt * r1_rdtice
 
             ! Contribution to mass flux, <0
-            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
+            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * r1_rdtice
 
             ! update heat content (J.m-2) and layer thickness
-            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
+            qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_1d(ji,nlay_i+1)
             h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
          ENDIF
@@ -512,12 +512,12 @@
             IF(  zf_tt(ji)  >=  0._wp  .AND. jk > icount(ji) ) THEN   ! do not calculate where layer has already disappeared from surface melting 
 
-               ztmelts = - tmut * s_i_b(ji,jk) + rtt  ! Melting point of layer jk (K)
-
-               IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
+               ztmelts = - tmut * s_i_1d(ji,jk) + rtt  ! Melting point of layer jk (K)
+
+               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
                   zintermelt(ji)    = 1._wp
 
-                  zEi               = - q_i_b(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0)
-
-                  !!zEw               = rcp * ( t_i_b(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_b (J/kg, <0)
+                  zEi               = - q_i_1d(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0)
+
+                  !!zEw               = rcp * ( t_i_1d(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0)
 
                   zdE               = 0._wp                         ! Specific enthalpy difference   (J/kg, <0)
@@ -532,19 +532,19 @@
 
                   ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 
-                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice
-
-                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
-                  sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
+                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice
+
+                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
+                  sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
                                     
                   ! Contribution to mass flux
-                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
 
                   ! update heat content (J.m-2) and layer thickness
-                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
+                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
 
                ELSE                               !!! Basal melting
 
-                  zEi               = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEi               = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
 
                   zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
@@ -567,17 +567,17 @@
 
                   ! Contribution to heat flux to the ocean [W.m-2], <0  
-                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
-
-                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
-                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
+                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice
+
+                  ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
+                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
                   
                   ! Total heat flux used in this process [W.m-2], >0  
-                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
+                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice
                   
                   ! Contribution to mass flux
-                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
+                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
 
                   ! update heat content (J.m-2) and layer thickness
-                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
+                  qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
                   h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
                ENDIF
@@ -602,10 +602,10 @@
 !
 !               ! excessive energy is sent to lateral ablation
-!               zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
-!               zq_1cat(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
+!               zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) )
+!               zq_1cat(ji) =  zinda * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0
 !
 !               ! correct salt and mass fluxes
-!               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
-!               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdvres * r1_rdtice
+!               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
+!               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice
 !            ENDIF
 !         END DO
@@ -616,5 +616,5 @@
       !-------------------------------------------
       DO ji = kideb, kiut
-         ht_i_b(ji) =  MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )
+         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_bott(ji) )
       END DO  
 
@@ -627,21 +627,21 @@
       DO ji = kideb, kiut
          zq_rema(ji)     = zq_su(ji) + zq_bo(ji) 
-!         zindh           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) )   ! =1 if snow
+!         zindh           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) )   ! =1 if snow
 !         zindq           = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) )
 !         zdeltah  (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
-!         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_b(ji) ) ) ! bound melting
+!         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting
 !         zdh_s_mel(ji)   = zdh_s_mel(ji) + zdeltah(ji,1)    
 !         dh_s_tot (ji)   = dh_s_tot(ji) + zdeltah(ji,1)
-!         ht_s_b   (ji)   = ht_s_b(ji)   + zdeltah(ji,1)
+!         ht_s_1d   (ji)   = ht_s_1d(ji)   + zdeltah(ji,1)
 !        
 !         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2)
 !         ! heat used to melt snow
-!         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
+!         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
 !         ! Contribution to mass flux
-!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
+!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
 !    
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
-         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_b(ji) ) * r1_rdtice
+         hfx_out(ii,ij)  = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
 
          IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
@@ -656,22 +656,22 @@
       DO ji = kideb, kiut
          !
-         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic )  )
-
-         ht_i_b(ji)     = ht_i_b(ji) + dh_snowice(ji)
-         ht_s_b(ji)     = ht_s_b(ji) - dh_snowice(ji)
+         dh_snowice(ji) = MAX(  0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic )  )
+
+         ht_i_1d(ji)     = ht_i_1d(ji) + dh_snowice(ji)
+         ht_s_1d(ji)     = ht_s_1d(ji) - dh_snowice(ji)
 
          ! Salinity of snow ice
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
-         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_b(ji)
+         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
 
          ! entrapment during snow ice formation
          ! new salinity difference stored (to be used in limthd_ent.F90)
          IF (  num_sal == 2  ) THEN
-            zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
+            zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi10 ) )
             ! salinity dif due to snow-ice formation
-            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch     
+            dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch     
             ! salinity dif due to bottom growth 
             IF (  zf_tt(ji)  < 0._wp ) THEN
-               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch
+               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi10 ) * zswitch
             ENDIF
          ENDIF
@@ -685,20 +685,20 @@
          
          ! Contribution to heat flux
-         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice 
+         hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 
 
          ! Contribution to salt flux
-         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice 
+         sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_1d(ji) * zfmdt * r1_rdtice 
           
          ! Contribution to mass flux
          ! All snow is thrown in the ocean, and seawater is taken to replace the volume
-         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
-         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
+         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice
+         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhosn * r1_rdtice
 
          ! update heat content (J.m-2) and layer thickness
-         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
+         qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_1d(ji,1) + zfmdt * zEw
          h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
          
          ! Total ablation (to debug)
-         IF( ht_i_b(ji) <= 0._wp )   a_i_b(ji) = 0._wp
+         IF( ht_i_1d(ji) <= 0._wp )   a_i_1d(ji) = 0._wp
 
       END DO !ji
@@ -710,6 +710,6 @@
       !clem bug: we should take snow into account here
       DO ji = kideb, kiut
-         zindh    =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) ) 
-         t_su_b(ji) =  zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt
+         zindh    =  1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 
+         t_su_1d(ji) =  zindh * t_su_1d(ji) + ( 1.0 - zindh ) * rtt
       END DO  ! ji
 
@@ -717,8 +717,8 @@
          DO ji = kideb,kiut
             ! mask enthalpy
-            zinda        =  MAX(  0._wp , SIGN( 1._wp, - ht_s_b(ji) )  )
-            q_s_b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)
-            ! recalculate t_s_b from q_s_b
-            t_s_b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
+            zinda        =  MAX(  0._wp , SIGN( 1._wp, - ht_s_1d(ji) )  )
+            q_s_1d(ji,jk) = ( 1.0 - zinda ) * q_s_1d(ji,jk)
+            ! recalculate t_s_1d from q_s_1d
+            t_s_1d(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
          END DO
       END DO
@@ -727,5 +727,5 @@
       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
       CALL wrk_dealloc( jpij, zintermelt )
-      CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )
+      CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i )
       CALL wrk_dealloc( jpij, icount )
       !
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4902)
@@ -75,10 +75,10 @@
       !!
       !! ** Inputs / Ouputs : (global commons)
-      !!           surface temperature : t_su_b
-      !!           ice/snow temperatures   : t_i_b, t_s_b
-      !!           ice salinities          : s_i_b
+      !!           surface temperature : t_su_1d
+      !!           ice/snow temperatures   : t_i_1d, t_s_1d
+      !!           ice salinities          : s_i_1d
       !!           number of layers in the ice/snow: nlay_i, nlay_s
       !!           profile of the ice/snow layers : z_i, z_s
-      !!           total ice/snow thickness : ht_i_b, ht_s_b
+      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
       !!
       !! ** External : 
@@ -98,5 +98,5 @@
       INTEGER ::   ii, ij      ! temporary dummy loop index
       INTEGER ::   numeq       ! current reference number of equation
-      INTEGER ::   layer       ! vertical dummy loop index 
+      INTEGER ::   jk       ! vertical dummy loop index 
       INTEGER ::   nconv       ! number of iterations in iterative procedure
       INTEGER ::   minnumeqmin, maxnumeqmax
@@ -114,6 +114,6 @@
       REAL(wp) ::   zerritmax   ! current maximal error on temperature 
       REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point
-      REAL(wp), POINTER, DIMENSION(:) ::   ztsuold     ! old surface temperature (before the iterative procedure )
-      REAL(wp), POINTER, DIMENSION(:) ::   ztsuoldit   ! surface temperature at previous iteration
+      REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure )
+      REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration
       REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
@@ -129,5 +129,5 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice
       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztiold      ! Old temperature in the ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice
       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice
       REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
@@ -137,12 +137,12 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow
       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s       ! Eta factor in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp      ! Temporary temperature in the snow to check the convergence
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsold       ! Temporary temperature in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s          ! Vertical cotes of the layers in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zindterm   ! Independent term
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zindtbis   ! temporary independent term
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zindterm    ! Independent term
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zindtbis    ! temporary independent term
       REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid     ! tridiagonal system terms
       ! diag errors on heat
       REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err
@@ -150,10 +150,10 @@
       ! 
       CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow )
-      CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
+      CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
       CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
-      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
-      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0)
-      CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  )
-      CALL wrk_alloc( jpij, jkmax+2, 3, ztrid )
+      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
+      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)
+      CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis  )
+      CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid )
 
       CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err )
@@ -162,6 +162,6 @@
       zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
       DO ji = kideb, kiut
-         zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  &
-            &           SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 
+         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
+            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
       END DO
 
@@ -172,10 +172,10 @@
       DO ji = kideb , kiut
          ! is there snow or not
-         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  )
+         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  )
          ! surface temperature of fusion
          ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt
          ! layer thickness
-         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i )
-         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
+         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i )
+         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
       END DO
 
@@ -187,13 +187,13 @@
       z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
 
-      DO layer = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
-         DO ji = kideb , kiut
-            z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )
-         END DO
-      END DO
-
-      DO layer = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
-         DO ji = kideb , kiut
-            z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )
+      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
+         DO ji = kideb , kiut
+            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s )
+         END DO
+      END DO
+
+      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
+         DO ji = kideb , kiut
+            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i )
          END DO
       END DO
@@ -216,8 +216,8 @@
       DO ji = kideb , kiut
          ! switches
-         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  ) 
+         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )  ) 
          ! hs > 0, isnow = 1
          zhsu (ji) = hnzst  ! threshold for the computation of i0
-         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )     
+         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) )     
 
          i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
@@ -226,5 +226,5 @@
          !            a function of the cloud cover
          !
-         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0)
+         !i0(ji)     =  (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0)
          !formula used in Cice
       END DO
@@ -248,10 +248,10 @@
       END DO
 
-      DO layer = 1, nlay_s          ! Radiation through snow
+      DO jk = 1, nlay_s          ! Radiation through snow
          DO ji = kideb, kiut
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
             !                             ! radiation absorbed by the layer-th snow layer
-            zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)
+            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
          END DO
       END DO
@@ -261,10 +261,10 @@
       END DO
 
-      DO layer = 1, nlay_i          ! Radiation through ice
+      DO jk = 1, nlay_i          ! Radiation through ice
          DO ji = kideb, kiut
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
             !                             ! radiation absorbed by the layer-th ice layer
-            zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)
+            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
          END DO
       END DO
@@ -280,19 +280,19 @@
       !
       DO ji = kideb, kiut        ! Old surface temperature
-         ztsuold  (ji) =  t_su_b(ji)                              ! temperature at the beg of iter pr.
-         ztsuoldit(ji) =  t_su_b(ji)                              ! temperature at the previous iter
-         t_su_b   (ji) =  MIN( t_su_b(ji), ztfs(ji) - ztsu_err )  ! necessary
+         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
+         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
+         t_su_1d   (ji) =  MIN( t_su_1d(ji), ztfs(ji) - ztsu_err )  ! necessary
          zerrit   (ji) =  1000._wp                                ! initial value of error
       END DO
 
-      DO layer = 1, nlay_s       ! Old snow temperature
-         DO ji = kideb , kiut
-            ztsold(ji,layer) =  t_s_b(ji,layer)
-         END DO
-      END DO
-
-      DO layer = 1, nlay_i       ! Old ice temperature
-         DO ji = kideb , kiut
-            ztiold(ji,layer) =  t_i_b(ji,layer)
+      DO jk = 1, nlay_s       ! Old snow temperature
+         DO ji = kideb , kiut
+            ztsb(ji,jk) =  t_s_1d(ji,jk)
+         END DO
+      END DO
+
+      DO jk = 1, nlay_i       ! Old ice temperature
+         DO ji = kideb , kiut
+            ztib(ji,jk) =  t_i_1d(ji,jk)
          END DO
       END DO
@@ -311,12 +311,12 @@
          IF( thcon_i_swi == 0 ) THEN      ! Untersteiner (1964) formula
             DO ji = kideb , kiut
-               ztcond_i(ji,0)        = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)
+               ztcond_i(ji,0)        = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt)
                ztcond_i(ji,0)        = MAX(ztcond_i(ji,0),zkimin)
             END DO
-            DO layer = 1, nlay_i-1
+            DO jk = 1, nlay_i-1
                DO ji = kideb , kiut
-                  ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) /  &
-                     MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)
-                  ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin)
+                  ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
+                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)
+                  ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin)
                END DO
             END DO
@@ -325,19 +325,20 @@
          IF( thcon_i_swi == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
             DO ji = kideb , kiut
-               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt )   &
-                  &                   - 0.011_wp * ( t_i_b(ji,1) - rtt )  
+               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt )   &
+                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rtt )  
                ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
             END DO
-            DO layer = 1, nlay_i-1
+            DO jk = 1, nlay_i-1
                DO ji = kideb , kiut
-                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   &
-                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   &
-                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )  
-                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin )
+                  ztcond_i(ji,jk) = rcdic +                                                                     & 
+                     &                 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                          &
+                     &                 / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt)   &
+                     &               - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt )  
+                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
                END DO
             END DO
             DO ji = kideb , kiut
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   &
-                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )  
+               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt)   &
+                  &                        - 0.011_wp * ( t_bo_1d(ji) - rtt )  
                ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
             END DO
@@ -355,15 +356,15 @@
          END DO
 
-         DO layer = 1, nlay_s-1
-            DO ji = kideb , kiut
-               zkappa_s(ji,layer)  = 2.0 * rcdsn / &
+         DO jk = 1, nlay_s-1
+            DO ji = kideb , kiut
+               zkappa_s(ji,jk)  = 2.0 * rcdsn / &
                   MAX(epsi10,2.0*zh_s(ji))
             END DO
          END DO
 
-         DO layer = 1, nlay_i-1
+         DO jk = 1, nlay_i-1
             DO ji = kideb , kiut
                !-- Ice kappa factors
-               zkappa_i(ji,layer)  = 2.0*ztcond_i(ji,layer)/ &
+               zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ &
                   MAX(epsi10,2.0*zh_i(ji)) 
             END DO
@@ -384,18 +385,18 @@
          !------------------------------------------------------------------------------|
          !
-         DO layer = 1, nlay_i
-            DO ji = kideb , kiut
-               ztitemp(ji,layer)   = t_i_b(ji,layer)
-               zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ &
-                  MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)
-               zeta_i(ji,layer)    = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &
+         DO jk = 1, nlay_i
+            DO ji = kideb , kiut
+               ztitemp(ji,jk)   = t_i_1d(ji,jk)
+               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ &
+                  MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10)
+               zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), &
                   epsi10)
             END DO
          END DO
 
-         DO layer = 1, nlay_s
-            DO ji = kideb , kiut
-               ztstemp(ji,layer) = t_s_b(ji,layer)
-               zeta_s(ji,layer)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
+         DO jk = 1, nlay_s
+            DO ji = kideb , kiut
+               ztstemp(ji,jk) = t_s_1d(ji,jk)
+               zeta_s(ji,jk)  = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)
             END DO
          END DO
@@ -408,5 +409,5 @@
             DO ji = kideb , kiut
                ! 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_b(ji) - ztsuoldit(ji) )
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
             END DO
          ENDIF
@@ -432,5 +433,5 @@
          !!ice interior terms (top equation has the same form as the others)
 
-         DO numeq=1,jkmax+2
+         DO numeq=1,nlay_i+3
             DO ji = kideb , kiut
                ztrid(ji,numeq,1) = 0.
@@ -445,11 +446,11 @@
          DO numeq = nlay_s + 2, nlay_s + nlay_i 
             DO ji = kideb , kiut
-               layer              = numeq - nlay_s - 1
-               ztrid(ji,numeq,1)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer-1)
-               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + &
-                  zkappa_i(ji,layer))
-               ztrid(ji,numeq,3)  =  - zeta_i(ji,layer)*zkappa_i(ji,layer)
-               zindterm(ji,numeq) =  ztiold(ji,layer) + zeta_i(ji,layer)* &
-                  zradab_i(ji,layer)
+               jk              = numeq - nlay_s - 1
+               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk-1)
+               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + &
+                  zkappa_i(ji,jk))
+               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk)
+               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* &
+                  zradab_i(ji,jk)
             END DO
          ENDDO
@@ -462,12 +463,12 @@
                +  zkappa_i(ji,nlay_i-1) )
             ztrid(ji,numeq,3)  =  0.0
-            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &
+            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &
                ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
-               *  t_bo_b(ji) ) 
+               *  t_bo_1d(ji) ) 
          ENDDO
 
 
          DO ji = kideb , kiut
-            IF ( ht_s_b(ji).gt.0.0 ) THEN
+            IF ( ht_s_1d(ji).gt.0.0 ) THEN
                !
                !------------------------------------------------------------------------------|
@@ -477,11 +478,11 @@
                !!snow interior terms (bottom equation has the same form as the others)
                DO numeq = 3, nlay_s + 1
-                  layer =  numeq - 1
-                  ztrid(ji,numeq,1)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer-1)
-                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + &
-                     zkappa_s(ji,layer) )
-                  ztrid(ji,numeq,3)   =  - zeta_s(ji,layer)*zkappa_s(ji,layer)
-                  zindterm(ji,numeq)  =  ztsold(ji,layer) + zeta_s(ji,layer)* &
-                     zradab_s(ji,layer)
+                  jk =  numeq - 1
+                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk-1)
+                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + &
+                     zkappa_s(ji,jk) )
+                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
+                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* &
+                     zradab_s(ji,jk)
                END DO
 
@@ -490,8 +491,8 @@
                   ztrid(ji,nlay_s+2,3)    =  0.0
                   zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
-                     t_bo_b(ji) 
+                     t_bo_1d(ji) 
                ENDIF
 
-               IF ( t_su_b(ji) .LT. rtt ) THEN
+               IF ( t_su_1d(ji) .LT. rtt ) THEN
 
                   !------------------------------------------------------------------------------|
@@ -506,5 +507,5 @@
                   ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
                   ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
-                  zindterm(ji,1) = dzf(ji)*t_su_b(ji)   - zf(ji)
+                  zindterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji)
 
                   !!first layer of snow equation
@@ -512,5 +513,5 @@
                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)
                   ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
+                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
 
                ELSE 
@@ -529,7 +530,7 @@
                      zkappa_s(ji,0) * zg1s )
                   ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *            &
+                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            &
                      ( zradab_s(ji,1) +                         &
-                     zkappa_s(ji,0) * zg1s * t_su_b(ji) ) 
+                     zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
                ENDIF
             ELSE
@@ -539,5 +540,5 @@
                !------------------------------------------------------------------------------|
                !
-               IF (t_su_b(ji) .LT. rtt) THEN
+               IF (t_su_1d(ji) .LT. rtt) THEN
                   !
                   !------------------------------------------------------------------------------|
@@ -553,5 +554,5 @@
                   ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1    
                   ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_b(ji) - zf(ji)
+                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
 
                   !!first layer of ice equation
@@ -560,5 +561,5 @@
                      + zkappa_i(ji,0) * zg1 )
                   ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)  
-                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)  
+                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)  
 
                   !!case of only one layer in the ice (surface & ice equations are altered)
@@ -573,6 +574,6 @@
                      ztrid(ji,numeqmin(ji)+1,3)  =  0.0
 
-                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1)* &
-                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )
+                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* &
+                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) )
                   ENDIF
 
@@ -593,6 +594,6 @@
                      zg1)  
                   ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
-                     zkappa_i(ji,0) * zg1 * t_su_b(ji) ) 
+                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
+                     zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
 
                   !!case of only one layer in the ice (surface & ice equations are altered)
@@ -602,7 +603,7 @@
                         zkappa_i(ji,1))
                      ztrid(ji,numeqmin(ji),3)  =  0.0
-                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1)* &
-                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) &
-                        + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
+                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* &
+                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) &
+                        + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
                   ENDIF
 
@@ -623,5 +624,5 @@
 
          maxnumeqmax = 0
-         minnumeqmin = jkmax+4
+         minnumeqmin = nlay_i+5
 
          DO ji = kideb , kiut
@@ -632,7 +633,7 @@
          END DO
 
-         DO layer = minnumeqmin+1, maxnumeqmax
-            DO ji = kideb , kiut
-               numeq               =  min(max(numeqmin(ji)+1,layer),numeqmax(ji))
+         DO jk = minnumeqmin+1, maxnumeqmax
+            DO ji = kideb , kiut
+               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
                   ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
@@ -644,12 +645,12 @@
          DO ji = kideb , kiut
             ! ice temperatures
-            t_i_b(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
+            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
          END DO
 
          DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1
             DO ji = kideb , kiut
-               layer    =  numeq - nlay_s - 1
-               t_i_b(ji,layer)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
-                  t_i_b(ji,layer+1))/zdiagbis(ji,numeq)
+               jk    =  numeq - nlay_s - 1
+               t_i_1d(ji,jk)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
+                  t_i_1d(ji,jk+1))/zdiagbis(ji,numeq)
             END DO
          END DO
@@ -657,15 +658,15 @@
          DO ji = kideb , kiut
             ! snow temperatures      
-            IF (ht_s_b(ji).GT.0._wp) &
-               t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
-               *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
-               *        MAX(0.0,SIGN(1.0,ht_s_b(ji))) 
+            IF (ht_s_1d(ji).GT.0._wp) &
+               t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
+               *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) &
+               *        MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 
 
             ! surface temperature
-            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  )
-            ztsuoldit(ji) = t_su_b(ji)
-            IF( t_su_b(ji) < ztfs(ji) ) &
-               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1)   &
-               &          + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
+            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) )  )  )
+            ztsubit(ji) = t_su_1d(ji)
+            IF( t_su_1d(ji) < ztfs(ji) ) &
+               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   &
+               &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
          END DO
          !
@@ -677,20 +678,20 @@
          ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd
          DO ji = kideb , kiut
-            t_su_b(ji) =  MAX(  MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp  )
-            zerrit(ji) =  ABS( t_su_b(ji) - ztsuoldit(ji) )     
-         END DO
-
-         DO layer  =  1, nlay_s
-            DO ji = kideb , kiut
-               t_s_b(ji,layer) = MAX(  MIN( t_s_b(ji,layer), rtt ), 190._wp  )
-               zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))
-            END DO
-         END DO
-
-         DO layer  =  1, nlay_i
-            DO ji = kideb , kiut
-               ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt 
-               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
-               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
+            t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp  )
+            zerrit(ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
+         END DO
+
+         DO jk  =  1, nlay_s
+            DO ji = kideb , kiut
+               t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rtt ), 190._wp  )
+               zerrit(ji)      = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk)))
+            END DO
+         END DO
+
+         DO jk  =  1, nlay_i
+            DO ji = kideb , kiut
+               ztmelt_i        = -tmut * s_i_1d(ji,jk) + rtt 
+               t_i_1d(ji,jk) =  MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp)
+               zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk)))
             END DO
          END DO
@@ -717,11 +718,11 @@
       DO ji = kideb, kiut
          ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 
-         IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )
+         IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) )
          !                                ! surface ice conduction flux
-         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  )
-         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   &
-            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji))
+         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) )  )
+         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
+            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
          !                                ! bottom ice conduction flux
-         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) )
+         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
       END DO
 
@@ -730,8 +731,10 @@
       !-----------------------------------------
       DO ji = kideb, kiut
-         IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC
-            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
+         IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC
+            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
+               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
          ELSE                         ! case T_su = 0degC
-            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji)
+            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
+               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
          ENDIF
       END DO
@@ -742,8 +745,8 @@
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
       DO ji = kideb, kiut
-         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) +  &
-            &                              SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) )
+         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
+            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) )
          zhfx_err(ji)   = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 
-         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_b(ji)
+         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
       END DO 
 
@@ -767,15 +770,16 @@
       DO ji = kideb, kiut
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
-         hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
+         hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) )
       END DO
    
       !
       CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
-      CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )
+      CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
       CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu )
-      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
-      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 )
-      CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )
-      CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid )
+      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   &
+         &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
+      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
+      CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis )
+      CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid )
       CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
 
@@ -798,8 +802,8 @@
       DO jk = 1, nlay_i             ! Sea ice energy of melting
          DO ji = kideb, kiut
-            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt 
-            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) )
-            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             &
-               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   &
+            ztmelts      = - tmut  * s_i_1d(ji,jk) + rtt 
+            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) )
+            q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) )                                             &
+               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) )   &
                &                   - rcp  *                 ( ztmelts-rtt )  ) 
          END DO
@@ -807,5 +811,5 @@
       DO jk = 1, nlay_s             ! Snow energy of melting
          DO ji = kideb, kiut
-            q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )
+            q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus )
          END DO
       END DO
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 4902)
@@ -146,5 +146,5 @@
       ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
       DO ji = kideb, kiut
-         hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_b(ji) * r1_rdtice *  &
+         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( qh_i_old(ji,0:nlay_i+1) ) ) 
       END DO
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4902)
@@ -72,9 +72,9 @@
       !!             - Computation of variation of ice volume and mass
       !!             - Computation of frldb after lateral accretion and 
-      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)      
+      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
       !!------------------------------------------------------------------------
-      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
-      INTEGER ::   layer, nbpac     ! local integers 
-      INTEGER ::   ii, ij, iter   !   -       -
+      INTEGER ::   ji,jj,jk,jl      ! dummy loop indices
+      INTEGER ::   nbpac            ! local integers 
+      INTEGER ::   ii, ij, iter     !   -       -
       REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde  ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
@@ -90,5 +90,5 @@
       REAL(wp) ::   zv_newfra
   
-      INTEGER , POINTER, DIMENSION(:) ::   jcat      ! indexes of categories where new ice grows
+      INTEGER , POINTER, DIMENSION(:) ::   jcat        ! indexes of categories where new ice grows
       REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
 
@@ -102,16 +102,15 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
       REAL(wp), POINTER, DIMENSION(:) ::   zat_i_1d    ! total ice fraction    
-      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
-      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb   ! accretion of frazil ice at the ice bottom
+      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
       REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
 
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d     ! 1-D version of a_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d     ! 1-D version of v_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d    ! 1-D version of oa_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d   ! 1-D version of smv_i
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_b      ! old volume of ice in category jl
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_b      ! old area of ice in category jl
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d  ! 1-D version of oa_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i
+
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d !: 1-D version of e_i
 
       REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
@@ -120,7 +119,7 @@
       CALL wrk_alloc( jpij, jcat )   ! integer
       CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d )
-      CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
-      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d )
+      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
+      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
+      CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d )
       CALL wrk_alloc( jpi,jpj, zvrel )
 
@@ -303,9 +302,8 @@
 
          CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw, jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
          CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
 
@@ -320,6 +318,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_old(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
-         za_old(1:nbpac,:) = za_i_1d(1:nbpac,:)
+         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
+         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
          !----------------------
          ! Thickness of new ice
@@ -328,5 +326,5 @@
             zh_newice(ji) = hiccrit
          END DO
-         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_b(1:nbpac)
+         IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
@@ -352,6 +350,6 @@
          DO ji = 1, nbpac
             ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
-               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) )   &
+            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                             &
+               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) )   &
                &                       - rcp  *         ( ztmelts - rtt )  )
          END DO ! ji
@@ -371,5 +369,5 @@
             zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
 
-            zEw           = rcp * ( t_bo_b(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_b [J/kg]
+            zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg]
                                                                    ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 
                                                                    
@@ -442,5 +440,5 @@
          DO ji = 1, nbpac
             jl = jcat(ji)                                                    ! categroy in which new ice is put
-            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice
+            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
          END DO
 
@@ -450,5 +448,5 @@
                zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
                ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      &
-                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) )  &
+                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
                   &        * zinda / MAX( zv_i_1d(ji,jl), epsi20 )
             END DO
@@ -492,5 +490,5 @@
             DO ji = 1, nbpac
                zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes
-               zoa_i_1d(ji,jl)  = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb   
+               zoa_i_1d(ji,jl)  = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb   
             END DO 
          END DO   
@@ -501,5 +499,5 @@
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl)
+               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
                zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
             END DO
@@ -520,5 +518,4 @@
          CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj )
          CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
-         CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
 
          CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
@@ -544,7 +541,7 @@
       CALL wrk_dealloc( jpij, jcat )   ! integer
       CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d )
-      CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
-      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d )
+      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
+      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )
+      CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d )
       CALL wrk_dealloc( jpi,jpj, zvrel )
       !
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 4902)
@@ -60,5 +60,5 @@
       !---------------------------------------------------------
       DO ji = kideb, kiut
-         sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)
+         sm_i_1d(ji) = sm_i_1d(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji)
       END DO
  
@@ -66,9 +66,9 @@
       ! 1) Constant salinity, constant in time                                       |
       !------------------------------------------------------------------------------|
-!!gm comment: if num_sal = 1 s_i_new, s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
+!!gm comment: if num_sal = 1 s_i_new, s_i_1d and sm_i_1d can be set to bulk_sal one for all in the initialisation phase !!
 !!gm           ===>>>   simplification of almost all test on num_sal value
       IF(  num_sal == 1  ) THEN
-            s_i_b  (kideb:kiut,1:nlay_i) =  bulk_sal
-            sm_i_b (kideb:kiut)          =  bulk_sal 
+            s_i_1d (kideb:kiut,1:nlay_i) =  bulk_sal
+            sm_i_1d(kideb:kiut)          =  bulk_sal 
             s_i_new(kideb:kiut)          =  bulk_sal
       ENDIF
@@ -83,6 +83,6 @@
             ! Switches 
             !----------
-            iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt )        )    ! =1 if summer 
-            igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo
+            iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rtt )        )     ! =1 if summer 
+            igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) )    ! =1 if t_su < t_bo
 
             !---------------------
@@ -90,7 +90,7 @@
             !---------------------
             ! drainage by gravity drainage
-            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
+            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_1d(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
             ! drainage by flushing  
-            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
+            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_1d(ji) - sal_F , 0._wp ) / time_F * rdt_ice
 
             !-----------------
@@ -99,10 +99,10 @@
             ! only drainage terms ( gravity drainage and flushing )
             ! snow ice / bottom sources are added in lim_thd_ent to conserve energy
-            sm_i_b(ji) = sm_i_b(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
+            sm_i_1d(ji) = sm_i_1d(ji) + dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji)
 
             !----------------------------
             ! Salt flux - brine drainage
             !----------------------------
-            sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_b(ji) * ht_i_b(ji) * ( dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) ) * r1_rdtice
+            sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoic * a_i_1d(ji) * ht_i_1d(ji) * ( dsm_i_fl_1d(ji) + dsm_i_gd_1d(ji) ) * r1_rdtice
 
          END DO
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4902)
@@ -63,5 +63,5 @@
       INTEGER, INTENT(in) ::   kt   ! number of iteration
       !
-      INTEGER  ::   ji, jj, jk, jl, layer   ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jl, jn      ! dummy loop indices
       INTEGER  ::   initad                  ! number of sub-timestep for the advection
       INTEGER  ::   ierr                    ! error status
@@ -85,5 +85,5 @@
       CALL wrk_alloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
       CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
-      CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e )
+      CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, zs0e )
 
       CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem
@@ -167,5 +167,5 @@
 
          IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
-            DO jk = 1,initad
+            DO jn = 1,initad
                CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
@@ -197,16 +197,16 @@
                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
-                     CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
-                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
-                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
-                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
-                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
-                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
+                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
+                     CALL lim_adv_x( zusnit, u_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                     CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                   END DO
                END DO
             END DO
          ELSE
-            DO jk = 1, initad
+            DO jn = 1, initad
                CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0ow (:,:), sxopw(:,:),   &             !--- ice open water area
                   &                                       sxxopw(:,:), syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
@@ -239,11 +239,11 @@
                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0c0 (:,:,jl), sxc0 (:,:,jl),   &
                      &                                       sxxc0 (:,:,jl), syc0 (:,:,jl), syyc0 (:,:,jl), sxyc0 (:,:,jl)  )
-                  DO layer = 1, nlay_i                                                           !--- ice heat contents ---
-                     CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
-                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
-                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
-                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,layer,jl), sxe (:,:,layer,jl),   & 
-                        &                                       sxxe(:,:,layer,jl), sye (:,:,layer,jl),   &
-                        &                                       syye(:,:,layer,jl), sxye(:,:,layer,jl) )
+                  DO jk = 1, nlay_i                                                           !--- ice heat contents ---
+                     CALL lim_adv_y( zusnit, v_ice, 1._wp , zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
+                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, zs0e(:,:,jk,jl), sxe (:,:,jk,jl),   & 
+                        &                                       sxxe(:,:,jk,jl), sye (:,:,jk,jl),   &
+                        &                                       syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                   END DO
                END DO
@@ -506,5 +506,5 @@
       CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold )
       CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi )
-      CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e )
+      CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, zs0e )
 
       CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zaiold, zhimax )   ! clem
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4902)
@@ -69,6 +69,5 @@
       !!                
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl, jm    ! dummy loop indices
-      INTEGER  ::   jbnd1, jbnd2
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
       INTEGER  ::   i_ice_switch
       REAL(wp) ::   zsal
@@ -93,9 +92,5 @@
       ! Rebin categories with thickness out of bounds
       !----------------------------------------------------
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
-      END DO
+      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
 
       at_i(:,:) = 0._wp
@@ -126,11 +121,5 @@
       ! Final thickness distribution rebinning
       ! --------------------------------------
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
-         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
-         ENDIF
-      END DO
+      IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl)
 
       !-----------------
@@ -161,14 +150,14 @@
       ! Diagnostics
       ! -------------------------------------------------
-      d_u_ice_dyn(:,:)     = u_ice(:,:)     - old_u_ice(:,:)
-      d_v_ice_dyn(:,:)     = v_ice(:,:)     - old_v_ice(:,:)
-      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - old_a_i  (:,:,:)
-      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - old_v_s  (:,:,:)  
-      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - old_v_i  (:,:,:)   
-      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - old_e_s  (:,:,:,:)  
-      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)
-      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
+      d_u_ice_dyn(:,:)     = u_ice(:,:)     - u_ice_b(:,:)
+      d_v_ice_dyn(:,:)     = v_ice(:,:)     - v_ice_b(:,:)
+      d_a_i_trp  (:,:,:)   = a_i  (:,:,:)   - a_i_b  (:,:,:)
+      d_v_s_trp  (:,:,:)   = v_s  (:,:,:)   - v_s_b  (:,:,:)  
+      d_v_i_trp  (:,:,:)   = v_i  (:,:,:)   - v_i_b  (:,:,:)   
+      d_e_s_trp  (:,:,:,:) = e_s  (:,:,:,:) - e_s_b  (:,:,:,:)  
+      d_e_i_trp  (:,:,1:nlay_i,:) = e_i  (:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:)
+      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - oa_i_b (:,:,:)
       d_smv_i_trp(:,:,:)   = 0._wp
-      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
+      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:)
 
       ! conservation test
@@ -186,5 +175,5 @@
          CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update1  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
          CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update1  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
-         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update1  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
+         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update1  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
 
          DO jl = 1, jpl
@@ -199,26 +188,26 @@
             CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update1  : o_i         : ')
             CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update1  : a_i         : ')
-            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_a_i     : ')
+            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : a_i_b       : ')
             CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_a_i_trp   : ')
             CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update1  : v_i         : ')
-            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_i     : ')
+            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_i_b       : ')
             CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_i_trp   : ')
             CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
-            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update1  : old_v_s     : ')
+            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ')
             CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update1  : d_v_s_trp   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
-            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i1    : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1_b      : ')
             CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : de_i1_trp   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
-            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : old_e_i2    : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2_b      : ')
             CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : de_i2_trp   : ')
             CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
-            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update1  : old_e_snow  : ')
+            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ')
             CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : d_e_s_trp   : ')
             CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update1  : smv_i       : ')
-            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update1  : old_smv_i   : ')
+            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update1  : smv_i_b     : ')
             CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update1  : d_smv_i_trp : ')
             CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update1  : oa_i        : ')
-            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update1  : old_oa_i    : ')
+            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update1  : oa_i_b      : ')
             CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update1  : d_oa_i_trp  : ')
 
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4902)
@@ -67,6 +67,5 @@
       !!
       !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl, jm    ! dummy loop indices
-      INTEGER  ::   jbnd1, jbnd2
+      INTEGER  ::   ji, jj, jk, jl    ! dummy loop indices
       INTEGER  ::   i_ice_switch
       REAL(wp) ::   zh, zsal
@@ -89,26 +88,19 @@
       ! Rebin categories with thickness out of bounds
       !----------------------------------------------------
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
-      END DO
+      IF ( jpl > 1 )   CALL lim_itd_th_reb(1, jpl)
 
       !----------------------------------------------------------------------
       ! Constrain the thickness of the smallest category above hiclim
       !----------------------------------------------------------------------
-      DO jm = 1, jpm
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               jl = ice_cat_bounds(jm,1)
-               IF( v_i(ji,jj,jl) > 0._wp .AND. ht_i(ji,jj,jl) < hiclim ) THEN
-                  zh             = hiclim / ht_i(ji,jj,jl)
-                  ht_s(ji,jj,jl) = ht_s(ji,jj,jl) * zh
-                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) * zh
-                  a_i (ji,jj,jl) = a_i(ji,jj,jl)  / zh
-               ENDIF
-            END DO !ji
-         END DO !jj
-      END DO !jm
+      DO jj = 1, jpj 
+         DO ji = 1, jpi
+            IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < hiclim ) THEN
+               zh             = hiclim / ht_i(ji,jj,1)
+               ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh
+               ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh
+               a_i (ji,jj,1) = a_i(ji,jj,1)  / zh
+            ENDIF
+         END DO
+      END DO
       
       !-----------------------------------------------------
@@ -139,11 +131,5 @@
       ! Final thickness distribution rebinning
       ! --------------------------------------
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
-         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
-         ENDIF
-      END DO
+      IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl )
 
       !-----------------
@@ -196,12 +182,12 @@
       ! Diagnostics
       ! -------------------------------------------------
-      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:) 
-      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:)
-      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)  
-      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 
-      d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:)
-      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:)
+      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - a_i_b(:,:,:) 
+      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - v_s_b(:,:,:)
+      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - v_i_b(:,:,:)  
+      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - e_s_b(:,:,:,:) 
+      d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:)
+      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - oa_i_b (:,:,:)
       d_smv_i_thd(:,:,:) = 0._wp
-      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
+      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:)
       ! diag only (clem)
       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
@@ -211,5 +197,6 @@
          DO ji = 1, jpi            
             diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  & 
-               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj)   
+               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) )    &
+               &                   ) * unit_fac * r1_rdtice / area(ji,jj)   
          END DO
       END DO
@@ -228,5 +215,5 @@
          CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update2  : strength    :')
          CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update2  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
-         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update2  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
+         CALL prt_ctl(tab2d_1=u_ice_b    , clinfo1=' lim_update2  : u_ice_b     :', tab2d_2=v_ice_b    , clinfo2=' v_ice_b     :')
 
          DO jl = 1, jpl
@@ -241,26 +228,26 @@
             CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update2  : o_i         : ')
             CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update2  : a_i         : ')
-            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_a_i     : ')
+            CALL prt_ctl(tab2d_1=a_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : a_i_b       : ')
             CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_a_i_thd   : ')
             CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update2  : v_i         : ')
-            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_i     : ')
+            CALL prt_ctl(tab2d_1=v_i_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_i_b       : ')
             CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_i_thd   : ')
             CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
-            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update2  : old_v_s     : ')
+            CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
             CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update2  : d_v_s_thd   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
-            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i1    : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1_b      : ')
             CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : de_i1_thd   : ')
             CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
-            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : old_e_i2    : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2_b      : ')
             CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : de_i2_thd   : ')
             CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
-            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update2  : old_e_snow  : ')
+            CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
             CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : d_e_s_thd   : ')
             CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update2  : smv_i       : ')
-            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update2  : old_smv_i   : ')
+            CALL prt_ctl(tab2d_1=smv_i_b    (:,:,jl)        , clinfo1= ' lim_update2  : smv_i_b     : ')
             CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update2  : d_smv_i_thd : ')
             CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update2  : oa_i        : ')
-            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update2  : old_oa_i    : ')
+            CALL prt_ctl(tab2d_1=oa_i_b     (:,:,jl)        , clinfo1= ' lim_update2  : oa_i_b      : ')
             CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update2  : d_oa_i_thd  : ')
 
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4902)
@@ -464,5 +464,5 @@
       ! Vertically constant, constant in time
       !---------------------------------------
-      IF( num_sal == 1 )   s_i_b(:,:) = bulk_sal
+      IF( num_sal == 1 )   s_i_1d(:,:) = bulk_sal
 
       !------------------------------------------------------
@@ -473,5 +473,5 @@
          !
          DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
-            z_slope_s(ji) = 2._wp * sm_i_b(ji) / MAX( epsi10 , ht_i_b(ji) )
+            z_slope_s(ji) = 2._wp * sm_i_1d(ji) / MAX( epsi10 , ht_i_1d(ji) )
          END DO
 
@@ -489,16 +489,16 @@
                ij =     ( npb(ji) - 1 ) / jpi + 1
                ! zind0 = 1 if sm_i le s_i_0 and 0 otherwise
-               zind0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_b(ji) ) ) 
+               zind0  = MAX( 0._wp , SIGN( 1._wp  , s_i_0 - sm_i_1d(ji) ) ) 
                ! zind01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws 
-               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_b(ji) ) ) 
+               zind01 = ( 1._wp - zind0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) ) 
                ! if 2.sm_i GE sss_m then zindbal = 1
                ! this is to force a constant salinity profile in the Baltic Sea
-               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_b(ji) - sss_m(ii,ij) ) )
+               zindbal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
                !
-               zalpha = (  zind0 + zind01 * ( sm_i_b(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zindbal )
+               zalpha = (  zind0 + zind01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zindbal )
                !
-               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_b(ji) * dummy_fac2
+               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2
                ! weighting the profile
-               s_i_b(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_b(ji)
+               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
             END DO ! ji
          END DO ! jk
@@ -512,5 +512,5 @@
       IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
          !
-         sm_i_b(:) = 2.30_wp
+         sm_i_1d(:) = 2.30_wp
          !
 !CDIR NOVERRCHK
@@ -519,5 +519,5 @@
             zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
             DO ji = kideb, kiut
-               s_i_b(ji,jk) = zsal
+               s_i_1d(ji,jk) = zsal
             END DO
          END DO
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 4902)
@@ -298,26 +298,48 @@
       !!----------------------------------------------------------------------
 
-      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2"   , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
+      CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   &
+      &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
 
       CALL histend( kid, snc4set )   ! end of the file definition
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 4902)
@@ -12,5 +12,4 @@
 
    !                                             !!! ice thermodynamics
-   INTEGER, PUBLIC, PARAMETER ::   jkmax    = 6   !: maximum number of ice layers
    INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 5   !: number of ice layers
    INTEGER, PUBLIC, PARAMETER ::   nlay_s   = 1   !: number of snow layers
@@ -18,5 +17,4 @@
    !                                             !!! ice mechanical redistribution
    INTEGER, PUBLIC, PARAMETER ::   jpl      = 5   !: number of ice categories
-   INTEGER, PUBLIC, PARAMETER ::   jpm      = 1   !: number of ice types
 
    !!----------------------------------------------------------------------
Index: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4901)
+++ branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4902)
@@ -34,5 +34,5 @@
    !!-----------------------------
    !: In ice thermodynamics, to spare memory, the vectors are folded
-   !: from 1D to 2D vectors. The following variables, with ending _1d (or _b)
+   !: from 1D to 2D vectors. The following variables, with ending _1d
    !: are the variables corresponding to 2d vectors
 
@@ -40,11 +40,11 @@
    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: correspondance between points (lateral accretion)
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlead_1d      !: <==> the 2D  qlead
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftr_ice_1d    !: <==> the 2D  ftr_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d    !: <==> the 2D  qsr_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr1_i0_1d     !: <==> the 2D  fr1_i0
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr2_i0_1d     !: <==> the 2D  fr2_i0
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d    !: <==> the 2D  qns_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_b        !: <==> the 2D  t_bo
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlead_1d     
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ftr_ice_1d   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qsr_ice_1d  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr1_i0_1d   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fr2_i0_1d   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d     
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d
@@ -65,28 +65,27 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_res_1d
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_1d    !: <==> the 2D  wfx_snw
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d    !: <==> the 2D  wfx_sub
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bog_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bom_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sum_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sni_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_opw_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_res_1d    !: <==> the 2D  wfx_ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_spr_1d    !: <==> the 2D  wfx_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bog_1d    
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_bom_1d   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sum_1d  
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sni_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_opw_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_res_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_spr_1d
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bri_1d    !: <==> the 2D sfx_bri
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bog_1d    !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bom_1d    !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sum_1d    !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sni_1d    !: 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_opw_1d    !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d    !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bri_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bog_1d    
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bom_1d    
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sum_1d    
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sni_1d    
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_opw_1d   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d  
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  at_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d       !: <==> the 2D  fhtur
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d        !: <==> the 2D  at_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d      !: <==> the 2D  fhtur
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhld_1d       !: <==> the 2D  fhld
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dqns_ice_1d   !: <==> the 2D  dqns_ice
@@ -100,27 +99,27 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_se_1d   !: Ice salinity variations due to basal salt entrapment
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dsm_i_si_1d   !: Ice salinity variations due to lateral accretion    
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hicol_b       !: Ice collection thickness accumulated in fleads
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hicol_1d      !: Ice collection thickness accumulated in leads
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_b      !: <==> the 2D  t_su
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_b       !: <==> the 2D  a_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_b      !: <==> the 2D  ht_s
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_b      !: <==> the 2D  ht_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su       !: Surface Conduction flux 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i     !: Bottom  Conduction flux 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot    !: Snow accretion/ablation        [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_surf   !: Ice surface accretion/ablation [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott   !: Ice bottom accretion/ablation  [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice  !: Snow ice formation             [m of ice]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_b      !: Ice bulk salinity [ppt]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new     !: Salinity of new ice at the bottom
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_s
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_s_1d       !: <==> the 2D  ht_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_surf     !: Ice surface accretion/ablation [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_1d       !: Ice bulk salinity [ppt]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_b   !: corresponding to the 2D var  t_s
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_b   !: corresponding to the 2D var  t_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   s_i_b   !: profiled ice salinity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_i_b   !:    Ice  enthalpy per unit volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_s_b   !:    Snow enthalpy per unit volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d   !: corresponding to the 2D var  t_s
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_1d   !: corresponding to the 2D var  t_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   s_i_1d   !: profiled ice salinity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_i_1d   !:    Ice  enthalpy per unit volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_s_1d   !:    Snow enthalpy per unit volume
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qh_i_old  !: ice heat content (q*h, J.m-2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old   !: ice thickness layer (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qh_i_old !: ice heat content (q*h, J.m-2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old  !: ice thickness layer (m)
 
    INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point
@@ -146,27 +145,31 @@
          &      qsr_ice_1d (jpij) ,     &
          &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) ,     &
-         &      t_bo_b   (jpij) ,                                          &
-         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , &
+         &      t_bo_1d   (jpij) ,                                          &
+         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,     & 
+         &      hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , &
-         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) )
+         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , &
+         &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) )
       !
-      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     &
-         &      fhtur_1d   (jpij) , wfx_ice_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,     &
-         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  wfx_res_1d (jpij) ,  &
+      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_1d     (jpij) ,     &
+         &      fhtur_1d   (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,     &
+         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , &
+         &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  wfx_res_1d (jpij) ,  &
          &      dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,     &
          &      tatm_ice_1d(jpij) ,      &   
          &      i0         (jpij) ,     &  
-         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , &
+         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,   &
+         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , &
          &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,     &     
-         &      dsm_i_si_1d(jpij) , hicol_b    (jpij)                     , STAT=ierr(2) )
+         &      dsm_i_si_1d(jpij) , hicol_1d    (jpij)                     , STAT=ierr(2) )
       !
-      ALLOCATE( t_su_b    (jpij) , a_i_b    (jpij) , ht_i_b   (jpij) ,    &   
-         &      ht_s_b    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &    
+      ALLOCATE( t_su_1d    (jpij) , a_i_1d    (jpij) , ht_i_1d   (jpij) ,    &   
+         &      ht_s_1d    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &    
          &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &    
-         &      dh_snowice(jpij) , sm_i_b   (jpij) , s_i_new  (jpij) ,    &
-         &      t_s_b(jpij,nlay_s),                                       &
-         &      t_i_b(jpij,jkmax), s_i_b(jpij,jkmax)                ,     &            
-         &      q_i_b(jpij,jkmax), q_s_b(jpij,jkmax)                ,     &
-         &      qh_i_old(jpij,0:jkmax), h_i_old(jpij,0:jkmax) , STAT=ierr(3))
+         &      dh_snowice(jpij) , sm_i_1d   (jpij) , s_i_new  (jpij) ,    &
+         &      t_s_1d(jpij,nlay_s),                                       &
+         &      t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1)                ,     &            
+         &      q_i_1d(jpij,nlay_i+1), q_s_1d(jpij,nlay_i+1)                ,     &
+         &      qh_i_old(jpij,0:nlay_i+1), h_i_old(jpij,0:nlay_i+1) , STAT=ierr(3))
       !
       thd_ice_alloc = MAXVAL( ierr )
