Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5462)
@@ -11,5 +11,4 @@
    !!   'key_lim3'                                      LIM-3 sea-ice model
    !!----------------------------------------------------------------------
-   USE par_ice        ! LIM sea-ice parameters
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
@@ -18,5 +17,5 @@
    PRIVATE
 
-   PUBLIC    ice_alloc  !  Called in iceini.F90
+   PUBLIC    ice_alloc  !  Called in sbc_lim_init
 
    !!======================================================================
@@ -110,7 +109,7 @@
    !! smv_i       |      -      |    Sea ice salt content         | ppt.m |
    !! oa_i        !      -      !    Sea ice areal age content    | day   |
-   !! e_i         !      -      !    Ice enthalpy                 | 10^9 J| 
+   !! e_i         !      -      !    Ice enthalpy                 | J/m2  | 
    !!      -      ! q_i_1d      !    Ice enthalpy per unit vol.   | J/m3  | 
-   !! e_s         !      -      !    Snow enthalpy                | 10^9 J| 
+   !! e_s         !      -      !    Snow enthalpy                | J/m2  | 
    !!      -      ! q_s_1d      !    Snow enthalpy per unit vol.  | J/m3  | 
    !!                                                                     |
@@ -148,6 +147,6 @@
    !! tm_i        |      -      |    Mean sea ice temperature     | K     |
    !! ot_i        !      -      !    Sea ice areal age content    | day   |
-   !! et_i        !      -      !    Total ice enthalpy           | 10^9 J| 
-   !! et_s        !      -      !    Total snow enthalpy          | 10^9 J| 
+   !! et_i        !      -      !    Total ice enthalpy           | J/m2  | 
+   !! et_s        !      -      !    Total snow enthalpy          | J/m2  | 
    !! bv_i        !      -      !    Mean relative brine volume   | ???   | 
    !!=====================================================================
@@ -165,67 +164,68 @@
    REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice
 
-   !                                     !!** ice-dynamic namelist (namicedyn) **
-   INTEGER , PUBLIC ::   nevp             !: number of iterations for subcycling
-   REAL(wp), PUBLIC ::   epsd             !: tolerance parameter for dynamic
-   REAL(wp), PUBLIC ::   om               !: relaxation constant
-   REAL(wp), PUBLIC ::   cw               !: drag coefficient for oceanic stress
-   REAL(wp), PUBLIC ::   pstar            !: determines ice strength (N/M), Hibler JPO79
-   REAL(wp), PUBLIC ::   c_rhg            !: determines changes in ice strength
-   REAL(wp), PUBLIC ::   creepl           !: creep limit : has to be under 1.0e-9
-   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)
-   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 ::   hminrhg          !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity
+   !                                     !!** ice-thickness distribution namelist (namiceitd) **
+   INTEGER , PUBLIC ::   nn_catbnd        !: categories distribution following: tanh function (1), or h^(-alpha) function (2)
+   REAL(wp), PUBLIC ::   rn_himean        !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only)
+
+   !                                     !!** ice-dynamics namelist (namicedyn) **
+   LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength
+   INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75)
+   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
+   INTEGER , PUBLIC ::   nn_ahi0          !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation)
+   REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_icestr = 1
+   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
+   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength (N/M), Hibler JPO79
+   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength
+   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
+   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
+   REAL(wp), PUBLIC ::   rn_ahi0_ref      !: sea-ice hor. eddy diffusivity coeff. (m2/s)
+   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
 
    !                                     !!** ice-salinity namelist (namicesal) **
-   REAL(wp), PUBLIC ::   s_i_max          !: maximum ice salinity [PSU]
-   REAL(wp), PUBLIC ::   s_i_min          !: minimum ice salinity [PSU]
-   REAL(wp), PUBLIC ::   s_i_0            !: 1st sal. value for the computation of sal .prof. [PSU]
-   REAL(wp), PUBLIC ::   s_i_1            !: 2nd sal. value for the computation of sal .prof. [PSU]
-   REAL(wp), PUBLIC ::   sal_G            !: restoring salinity for gravity drainage [PSU]
-   REAL(wp), PUBLIC ::   sal_F            !: restoring salinity for flushing [PSU]
-   REAL(wp), PUBLIC ::   time_G           !: restoring time constant for gravity drainage (= 20 days) [s]
-   REAL(wp), PUBLIC ::   time_F           !: restoring time constant for gravity drainage (= 10 days) [s]
-   REAL(wp), PUBLIC ::   bulk_sal         !: bulk salinity (ppt) in case of constant salinity
+   REAL(wp), PUBLIC ::   rn_simax         !: maximum ice salinity [PSU]
+   REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU]
+   REAL(wp), PUBLIC ::   rn_sal_gd        !: restoring salinity for gravity drainage [PSU]
+   REAL(wp), PUBLIC ::   rn_sal_fl        !: restoring salinity for flushing [PSU]
+   REAL(wp), PUBLIC ::   rn_time_gd       !: restoring time constant for gravity drainage (= 20 days) [s]
+   REAL(wp), PUBLIC ::   rn_time_fl       !: restoring time constant for gravity drainage (= 10 days) [s]
+   REAL(wp), PUBLIC ::   rn_icesal        !: bulk salinity (ppt) in case of constant salinity
 
    !                                     !!** ice-salinity namelist (namicesal) **
-   INTEGER , PUBLIC ::   num_sal             !: salinity configuration used in the model
+   INTEGER , PUBLIC ::   nn_icesal           !: salinity configuration used in the model
    !                                         ! 1 - constant salinity in both space and time
    !                                         ! 2 - prognostic salinity (s(z,t))
    !                                         ! 3 - salinity profile, constant in time
-   INTEGER , PUBLIC ::   thcon_i_swi         !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)
+   INTEGER , PUBLIC ::   nn_ice_thcon        !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)
+   INTEGER , PUBLIC ::   nn_monocat          !: virtual ITD mono-category parameterizations (1) or not (0)
+   LOGICAL , PUBLIC ::   ln_it_qnsice        !: iterate surface flux with changing surface temperature or not (F)
 
    !                                     !!** ice-mechanical redistribution namelist (namiceitdme)
-   REAL(wp), PUBLIC ::   Cs               !: fraction of shearing energy contributing to ridging            
-   REAL(wp), PUBLIC ::   Cf               !: ratio of ridging work to PE loss
-   REAL(wp), PUBLIC ::   fsnowrdg         !: fractional snow loss to the ocean during ridging
-   REAL(wp), PUBLIC ::   fsnowrft         !: fractional snow loss to the ocean during ridging
-   REAL(wp), PUBLIC ::   Gstar            !: fractional area of young ice contributing to ridging
-   REAL(wp), PUBLIC ::   astar            !: equivalent of G* for an exponential participation function
-   REAL(wp), PUBLIC ::   Hstar            !: thickness that determines the maximal thickness of ridged ice
-   REAL(wp), PUBLIC ::   hparmeter        !: threshold thickness (m) for rafting / ridging 
-   REAL(wp), PUBLIC ::   Craft            !: coefficient for smoothness of the hyperbolic tangent in rafting
-   REAL(wp), PUBLIC ::   ridge_por        !: initial porosity of ridges (0.3 regular value)
-   REAL(wp), PUBLIC ::   betas            !: coef. for partitioning of snowfall between leads and sea ice
-   REAL(wp), PUBLIC ::   kappa_i          !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   REAL(wp), PUBLIC ::   nconv_i_thd      !: maximal number of iterations for heat diffusion
-   REAL(wp), PUBLIC ::   maxer_i_thd      !: maximal tolerated error (C) for heat diffusion
+   REAL(wp), PUBLIC ::   rn_cs            !: fraction of shearing energy contributing to ridging            
+   REAL(wp), PUBLIC ::   rn_fsnowrdg      !: fractional snow loss to the ocean during ridging
+   REAL(wp), PUBLIC ::   rn_fsnowrft      !: fractional snow loss to the ocean during ridging
+   REAL(wp), PUBLIC ::   rn_gstar         !: fractional area of young ice contributing to ridging
+   REAL(wp), PUBLIC ::   rn_astar         !: equivalent of G* for an exponential participation function
+   REAL(wp), PUBLIC ::   rn_hstar         !: thickness that determines the maximal thickness of ridged ice
+   REAL(wp), PUBLIC ::   rn_hraft         !: threshold thickness (m) for rafting / ridging 
+   REAL(wp), PUBLIC ::   rn_craft         !: coefficient for smoothness of the hyperbolic tangent in rafting
+   REAL(wp), PUBLIC ::   rn_por_rdg       !: initial porosity of ridges (0.3 regular value)
+   REAL(wp), PUBLIC ::   rn_betas         !: coef. for partitioning of snowfall between leads and sea ice
+   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
+   REAL(wp), PUBLIC ::   nn_conv_dif      !: maximal number of iterations for heat diffusion
+   REAL(wp), PUBLIC ::   rn_terr_dif      !: maximal tolerated error (C) for heat diffusion
 
    !                                     !!** ice-mechanical redistribution namelist (namiceitdme)
-   INTEGER , PUBLIC ::   ridge_scheme_swi !: scheme used for ice ridging
-   INTEGER , PUBLIC ::   raft_swi         !: rafting of ice or not                        
-   INTEGER , PUBLIC ::   partfun_swi      !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007)
-   INTEGER , PUBLIC ::   brinstren_swi    !: use brine volume to diminish ice strength
-
-   REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( ecc * ecc )
-   REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cw
-
+   LOGICAL , PUBLIC ::   ln_rafting      !: rafting of ice or not                        
+   INTEGER , PUBLIC ::   nn_partfun      !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007)
+
+   REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( rn_ecc * rn_ecc )
+   REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cio
+   REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i
+   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s 
+   !
    !                                     !!** switch for presence of ice or not 
    REAL(wp), PUBLIC ::   rswitch
-
+   !
    !                                     !!** define some parameters 
-   REAL(wp), PUBLIC, PARAMETER ::   unit_fac = 1.e+09_wp  !: conversion factor for ice / snow enthalpy
    REAL(wp), PUBLIC, PARAMETER ::   epsi06   = 1.e-06_wp  !: small number 
    REAL(wp), PUBLIC, PARAMETER ::   epsi10   = 1.e-10_wp  !: small number 
@@ -266,4 +266,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2]
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics) [s-1]
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bom     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -282,4 +286,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations 
@@ -296,7 +301,4 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
-
-   ! temporary arrays for dummy version of the code
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dh_i_surf2D, dh_i_bott2D, q_s
 
    !!--------------------------------------------------------------------------
@@ -333,5 +335,5 @@
       
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i        !: ice temperatures          [K]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i        !: ice thermal contents [Giga J]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i        !: ice thermal contents    [J/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   s_i        !: ice salinities          [PSU]
 
@@ -356,57 +358,47 @@
    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
-      
-
-   !!--------------------------------------------------------------------------
-   !! * Increment of global variables
-   !!--------------------------------------------------------------------------
+            
+   !!--------------------------------------------------------------------------
+   !! * Ice thickness distribution variables
+   !!--------------------------------------------------------------------------
+   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 
+
+   !!--------------------------------------------------------------------------
+   !! * Ice Run
+   !!--------------------------------------------------------------------------
+   !                                                  !!: ** Namelist namicerun read in sbc_lim_init **
+   INTEGER          , PUBLIC ::   jpl             !: number of ice  categories 
+   INTEGER          , PUBLIC ::   nlay_i          !: number of ice  layers 
+   INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers 
+   CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory
+   CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory
+   LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
+   LOGICAL          , PUBLIC ::   ln_icectl       !: flag for sea-ice points output (T) or not (F)
+   REAL(wp)         , PUBLIC ::   rn_amax         !: maximum ice concentration
+   INTEGER          , PUBLIC ::   iiceprt         !: debug i-point
+   INTEGER          , PUBLIC ::   jiceprt         !: debug j-point
+   !
+   !!--------------------------------------------------------------------------
+   !! * Ice diagnostics
+   !!--------------------------------------------------------------------------
+   ! Increment of global variables
    ! thd refers to changes induced by thermodynamics
    ! trp   ''         ''     ''       advection (transport of ice)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_a_i_thd  , d_a_i_trp                 !: icefractions                  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_v_s_thd  , d_v_s_trp                 !: snow volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_v_i_thd  , d_v_i_trp                 !: ice  volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_smv_i_thd, d_smv_i_trp               !:     
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_sm_i_fl  , d_sm_i_gd                 !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_sm_i_se  , d_sm_i_si  , d_sm_i_la    !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   d_oa_i_thd , d_oa_i_trp                !:
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_s_thd  , d_e_s_trp     !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   d_e_i_thd  , d_e_i_trp     !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   d_u_ice_dyn, d_v_ice_dyn   !: ice velocity 
-      
-   !!--------------------------------------------------------------------------
-   !! * Ice thickness distribution variables
-   !!--------------------------------------------------------------------------
-   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 
-
-   !!--------------------------------------------------------------------------
-   !! * Ice Run
-   !!--------------------------------------------------------------------------
-   !                                                  !!: ** Namelist namicerun read in iceini **
-   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
-   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
-   LOGICAL               , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
-   LOGICAL               , PUBLIC ::   ln_nicep        !: flag for sea-ice points output (T) or not (F)
-   REAL(wp)              , PUBLIC ::   cai             !: atmospheric drag over sea ice
-   REAL(wp)              , PUBLIC ::   cao             !: atmospheric drag over ocean
-   REAL(wp)              , PUBLIC ::   amax            !: maximum ice concentration
+   LOGICAL , PUBLIC                                        ::   ln_limdiahsb  !: flag for ice diag (T) or not (F)
+   LOGICAL , PUBLIC                                        ::   ln_limdiaout  !: flag for ice diag (T) or not (F)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vi   !: transport of ice volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_vs   !: transport of snw volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_ei   !: transport of ice enthalpy (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_es   !: transport of snw enthalpy (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_trp_smv  !: transport of salt content
    !
-   !!--------------------------------------------------------------------------
-   !! * Ice diagnostics
-   !!--------------------------------------------------------------------------
-   !! Check if everything down here is necessary
-   LOGICAL , PUBLIC                                      ::   ln_limdiahsb  !: flag for ice diag (T) or not (F)
-   LOGICAL , PUBLIC                                      ::   ln_limdiaout  !: flag for ice diag (T) or not (F)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dv_dt_thd     !: thermodynamic growth rates 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_vi   !: transport of ice volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_vs   !: transport of snw volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_ei   !: transport of ice enthalpy (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_es   !: transport of snw enthalpy (W/m2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_heat     !: snw/ice heat content variation   [W/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_smvi     !: ice salt content variation   [] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vice     !: ice volume variation   [m/s] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vsnw     !: snw volume variation   [m/s] 
    !
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_heat_dhc !: snw/ice heat content variation   [W/m2] 
-   !
-   INTEGER , PUBLIC ::   jiindx, jjindx        !: indexes of the debugging point
-
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
@@ -422,5 +414,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(19), ii
+      INTEGER :: ierr(17), ii
       !!-----------------------------------------------------------------
 
@@ -439,15 +431,17 @@
 
       ii = ii + 1
-      ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,      &
-         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,      &
-         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,    &
+      ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,                        &
+         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        &
+         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                        &
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
-         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     &
-         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      &
-         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      &
-         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,   &
-         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , hfx_err_rem(jpi,jpj), &
-         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,  &
-         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) , &
+         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
+         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                        &
+         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     &
+         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                        &
+         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    &
+         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
+         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) ,                                   &
+         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,                           &
+         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) ,    &
          &      hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,  STAT=ierr(ii) )
 
@@ -464,8 +458,7 @@
          &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) )
       ii = ii + 1
-      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) )
+      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,jpl) , e_i(jpi,jpj,nlay_i,jpl) , s_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
 
       ! * Moments for advection
@@ -483,24 +476,12 @@
          &      STAT=ierr(ii) )
       ii = ii + 1
-      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) )
+      ALLOCATE( sxe (jpi,jpj,nlay_i,jpl) , sye (jpi,jpj,nlay_i,jpl) , sxxe(jpi,jpj,nlay_i,jpl) ,     &
+         &      syye(jpi,jpj,nlay_i,jpl) , sxye(jpi,jpj,nlay_i,jpl)                            , STAT=ierr(ii) )
 
       ! * Old values of global variables
       ii = ii + 1
       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
-      ii = ii + 1
-      ALLOCATE( d_a_i_thd(jpi,jpj,jpl) , d_a_i_trp (jpi,jpj,jpl) , d_v_s_thd  (jpi,jpj,jpl) , d_v_s_trp  (jpi,jpj,jpl) ,   &
-         &      d_v_i_thd(jpi,jpj,jpl) , d_v_i_trp (jpi,jpj,jpl) , d_smv_i_thd(jpi,jpj,jpl) , d_smv_i_trp(jpi,jpj,jpl) ,   &     
-         &      d_sm_i_fl(jpi,jpj,jpl) , d_sm_i_gd (jpi,jpj,jpl) , d_sm_i_se  (jpi,jpj,jpl) , d_sm_i_si  (jpi,jpj,jpl) ,   &
-         &      d_sm_i_la(jpi,jpj,jpl) , d_oa_i_thd(jpi,jpj,jpl) , d_oa_i_trp (jpi,jpj,jpl) ,   &
-         &     STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( d_e_s_thd(jpi,jpj,nlay_s,jpl) , d_e_i_thd(jpi,jpj,nlay_i,jpl) , d_u_ice_dyn(jpi,jpj) ,     &
-         &      d_e_s_trp(jpi,jpj,nlay_s,jpl) , d_e_i_trp(jpi,jpj,nlay_i,jpl) , d_v_ice_dyn(jpi,jpj) , STAT=ierr(ii) )
+         &      a_i_b  (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) ,     &
+         &      oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj)     , v_ice_b(jpi,jpj)          , STAT=ierr(ii) )
       
       ! * Ice thickness distribution variables
@@ -510,7 +491,7 @@
       ! * Ice diagnostics
       ii = ii + 1
-      ALLOCATE( dv_dt_thd(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( diag_trp_vi(jpi,jpj), diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj),   & 
+         &      diag_trp_es(jpi,jpj), diag_trp_smv(jpi,jpj), diag_heat  (jpi,jpj),   &
+         &      diag_smvi  (jpi,jpj), diag_vice   (jpi,jpj), diag_vsnw  (jpi,jpj), STAT=ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
