Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8514)
@@ -40,8 +40,8 @@
 !------------------------------------------------------------------------------
           ! -- ice_rdgrft_strength -- !
-   ln_str_Hib       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)                      
+   ln_str_H79       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)                      
       rn_pstar      =   2.0e+04       !     ice strength thickness parameter (N/m2) 
       rn_crhg       =   20.0          !     ice strength conc. parameter (-)
-   ln_str_Rot       = .false.         !  ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2)    
+   ln_str_R75       = .false.         !  ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2)    
       rn_perdg      =   17.0          !     ridging work divided by pot. energy change in ridging
                    ! -- ice_rdgrft -- !
@@ -94,92 +94,95 @@
    ln_adv_Pra       = .false.         !  Advection scheme (Prather)
    ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho)
-      nn_UMx        =   5             !     order of the scheme for UMx (1-5)
+      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order)
 /
 !------------------------------------------------------------------------------
 &namice_thd     !   Ice thermodynamics
 !------------------------------------------------------------------------------
-   ln_limthd      =  .true.           !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-                   ! -- limthd_dif -- !
-   rn_kappa_i     = 1.0               !  radiation attenuation coefficient in sea ice (m-1)
-   nn_ice_thcon   = 1                 !  sea ice thermal conductivity
-                                      !     0: k = k0 + beta.S/T            (Untersteiner, 1964)
-                                      !     1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007)
-   ln_dqnsice     = .true.            !  change the surface non-solar flux with surface temperature (T) or not (F)
-   rn_cdsn        = 0.31              !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
-                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
-                    ! -- limthd_dh -- !
-   ln_limdH       = .true.            !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   rn_betas       = 0.66              !  exponent in lead-ice repratition of snow precipitation
-                                      !     betas = 1 -> equipartition, betas < 1 -> more on leads
-                    ! -- limthd_da -- !
-   ln_limdA       = .true.            !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   rn_beta        = 1.0               !     (ln_latmelt=T) coef. beta for lateral melting param. Recommended range=[0.8-1.2]
-                                      !      => decrease = more melt and melt peaks toward higher concentration (A~0.5 for beta=1 ; A~0.8 for beta=0.2)
-                                      !         0.3 = best fit for western Fram Strait and Antarctica
-                                      !         1.4 = best fit for eastern Fram Strait      
-   rn_dmin        = 8.                !     (ln_latmelt=T) minimum floe diameter for lateral melting param. Recommended range=[6-10]
-                                      !      => 6  vs 8m = +40% melting at the peak (A~0.5)
-                                      !         10 vs 8m = -20% melting
-                   ! -- limthd_lac -- !
-   ln_limdO       = .true.            !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   rn_hnewice     = 0.1               !  thickness for new ice formation in open water (m)
-   ln_frazil      = .false.           !  Frazil ice parameterization (ice collection as a function of wind)
-   rn_maxfrazb    = 1.0               !     (ln_frazil=T) maximum fraction of frazil ice collecting at the ice base
-   rn_vfrazb      = 0.417             !     (ln_frazil=T) thresold drift speed for frazil ice collecting at the ice bottom (m/s)
-   rn_Cfrazb      = 5.0               !     (ln_frazil=T) squeezing coefficient for frazil ice collecting at the ice bottom
+   ln_icethd        = .true.            !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+                     ! -- limthd_dif -- !
+   rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice (m-1)
+   ln_cndi_U64      = .false.           !  sea ice thermal conductivity: k = k0 + beta.S/T            (Untersteiner, 1964)
+   ln_cndi_P07      = .true.            !  sea ice thermal conductivity: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007)
+   ln_dqns_i        = .true.            !  change the surface non-solar flux with surface temperature (T) or not (F)
+   rn_cnd_s         =   0.31            !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
+                                        !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
+                      ! -- limthd_dh -- !
+   ln_icedH         = .true.            !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+   rn_blow_s        =   0.66            !     mesure of snow blowing into the leads
+                                        !     = 1 => no snow blowing, < 1 => some snow blowing
+                      ! -- limthd_da -- !
+   ln_icedA         = .true.            !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+      rn_beta       =   1.0             !     coef. beta for lateral melting param. Recommended range=[0.8-1.2]
+                                        !      => decrease = more melt and melt peaks toward higher concentration (A~0.5 for beta=1 ; A~0.8 for beta=0.2)
+                                        !         0.3 = best fit for western Fram Strait and Antarctica
+                                        !         1.4 = best fit for eastern Fram Strait      
+      rn_dmin       =   8.              !     minimum floe diameter for lateral melting param. Recommended range=[6-10]
+                                        !      => 6  vs 8m = +40% melting at the peak (A~0.5)
+                                        !         10 vs 8m = -20% melting
+                     ! -- limthd_lac -- !
+   ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+      rn_hinew      =   0.1             !     thickness for new ice formation in open water (m)
+      ln_frazil     = .false.           !     Frazil ice parameterization (ice collection as a function of wind)
+         rn_maxfraz =   1.0             !        maximum fraction of frazil ice collecting at the ice base
+         rn_vfraz   =   0.417           !        thresold drift speed for frazil ice collecting at the ice bottom (m/s)
+         rn_Cfraz   =   5.0             !        squeezing coefficient for frazil ice collecting at the ice bottom
                     ! -- limitd_th -- !
-   rn_himin       = 0.1               !  minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice
-                       ! -- icestp -- !
-   nn_limflx      = -1                !  LIM3 Multi-category heat flux formulation (use -1 if LIM3 is not used)
-                                      !  =-1  Use per-category fluxes, bypass redistributor, forced mode only, not yet implemented coupled
-                                      !  = 0  Average per-category fluxes (forced and coupled mode)
-                                      !  = 1  Average and redistribute per-category fluxes, forced mode only, not yet implemented coupled
-                                      !  = 2  Redistribute a single flux over categories (coupled mode only)
+   rn_himin         =   0.1             !  minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice
+                         ! -- icestp -- !
+   nn_iceflx        =  -1               !  Redistribute heat flux over ice categories
+                                        !   =-1  Do nothing (needs N(cat) fluxes)
+                                        !        ==> forced mode only
+                                        !   = 0  Average N(cat) fluxes then apply the average over the N(cat) ice
+                                        !        ==> forced and coupled modes
+                                        !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity
+                                        !        ==> forced mode only
+                                        !   = 2  Redistribute a single flux over categories
+                                        !        ==> coupled mode only
 /
 !------------------------------------------------------------------------------
 &namice_sal     !   Ice salinity
 !------------------------------------------------------------------------------
-                   ! -- limthd_sal -- !
-   ln_limdS       = .true.            !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-   nn_icesal      =  2                !  ice salinity option
-                                      !     1: constant ice salinity (S=rn_icesal)
-                                      !     2: varying salinity parameterization S(z,t)
-                                      !     3: prescribed salinity profile S(z), Schwarzacher, 1959
-   rn_icesal      =  4.               !    (nn_icesal=1) ice salinity (g/kg)
-   rn_sal_gd      =  5.               !  restoring ice salinity, gravity drainage (g/kg)
-   rn_time_gd     =  1.73e+6          !  restoring time scale, gravity drainage  (s)
-   rn_sal_fl      =  2.               !  restoring ice salinity, flushing (g/kg)
-   rn_time_fl     =  8.64e+5          !  restoring time scale, flushing (s)
-   rn_simax       = 20.               !  maximum tolerated ice salinity (g/kg)
-   rn_simin       =  0.1              !  minimum tolerated ice salinity (g/kg)
+                    ! -- limthd_sal -- !
+   ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+   nn_icesal      =   2                !  ice salinity option
+                                       !     1: constant ice salinity (S=rn_icesal)
+                                       !     2: varying salinity parameterization S(z,t)
+                                       !     3: prescribed salinity profile S(z), Schwarzacher, 1959
+   rn_icesal      =   4.               !    (nn_icesal=1) ice salinity (g/kg)
+   rn_sal_gd      =   5.               !  restoring ice salinity, gravity drainage (g/kg)
+   rn_time_gd     =   1.73e+6          !  restoring time scale, gravity drainage  (s)
+   rn_sal_fl      =   2.               !  restoring ice salinity, flushing (g/kg)
+   rn_time_fl     =   8.64e+5          !  restoring time scale, flushing (s)
+   rn_simax       =  20.               !  maximum tolerated ice salinity (g/kg)
+   rn_simin       =   0.1              !  minimum tolerated ice salinity (g/kg)
 /
 !------------------------------------------------------------------------------
 &namicemp      !   Melt ponds
 !------------------------------------------------------------------------------
-                 ! -- limmp      -- !
-   ln_pnd         = .false.         !  active melt ponds
-   ln_pnd_rad     = .false.         !  active melt ponds radiative coupling
-   ln_pnd_fw      = .false.         !  active melt ponds freshwater coupling
-   nn_pnd_scheme  = 0               !  type of melt pond scheme  : =0 prescribed ( Tsu=0 ), =1 empirical, =2 topographic
-   rn_apnd        = 0.2             !  prescribed pond fraction, at Tsu=0  : (0<rn_apnd<1, nn_pnd_scheme = 0)
-   rn_hpnd        = 0.05            !  prescribed pond depth, at Tsu=0     : (0<rn_apnd<1, nn_pnd_scheme = 0)
+                    ! -- limmp      -- !
+   ln_pnd         = .false.            !  active melt ponds
+   ln_pnd_rad     = .false.            !  active melt ponds radiative coupling
+   ln_pnd_fw      = .false.            !  active melt ponds freshwater coupling
+   nn_pnd_scheme  =   0                !  type of melt pond scheme  : =0 prescribed ( Tsu=0 ), =1 empirical, =2 topographic
+   rn_apnd        =   0.2              !  prescribed pond fraction, at Tsu=0  : (0<rn_apnd<1, nn_pnd_scheme = 0)
+   rn_hpnd        =   0.05             !  prescribed pond depth, at Tsu=0     : (0<rn_apnd<1, nn_pnd_scheme = 0)
 /
 !------------------------------------------------------------------------------
 &namice_ini     !   Ice initialization
 !------------------------------------------------------------------------------
-                  ! -- limistate -- !
-   ln_iceini      = .true.          !  activate ice initialization (T) or not (F)
-   ln_iceini_file = .false.         !  netcdf file provided for initialization (T) or not (F)
-   rn_thres_sst   =  2.0            !  maximum water temperature with initial ice (degC)
-   rn_hts_ini_n   =  0.3            !  initial real snow thickness (m), North
-   rn_hts_ini_s   =  0.3            !        "            "             South
-   rn_hti_ini_n   =  3.0            !  initial real ice thickness  (m), North
-   rn_hti_ini_s   =  1.0            !        "            "             South
-   rn_ati_ini_n   =  0.9            !  initial ice concentration   (-), North
-   rn_ati_ini_s   =  0.9            !        "            "             South
-   rn_smi_ini_n   =  6.3            !  initial ice salinity     (g/kg), North
-   rn_smi_ini_s   =  6.3            !        "            "             South
-   rn_tmi_ini_n   =  270.           !  initial ice/snw temperature (K), North
-   rn_tmi_ini_s   =  270.           !        "            "             South
+                     ! -- limistate -- !
+   ln_iceini      = .true.             !  activate ice initialization (T) or not (F)
+   ln_iceini_file = .false.            !  netcdf file provided for initialization (T) or not (F)
+   rn_thres_sst   =   2.0              !  maximum water temperature with initial ice (degC)
+   rn_hts_ini_n   =   0.3              !  initial real snow thickness (m), North
+   rn_hts_ini_s   =   0.3              !        "            "             South
+   rn_hti_ini_n   =   3.0              !  initial real ice thickness  (m), North
+   rn_hti_ini_s   =   1.0              !        "            "             South
+   rn_ati_ini_n   =   0.9              !  initial ice concentration   (-), North
+   rn_ati_ini_s   =   0.9              !        "            "             South
+   rn_smi_ini_n   =   6.3              !  initial ice salinity     (g/kg), North
+   rn_smi_ini_s   =   6.3              !        "            "             South
+   rn_tmi_ini_n   = 270.               !  initial ice/snw temperature (K), North
+   rn_tmi_ini_s   = 270.               !        "            "             South
 
    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', ''
@@ -194,21 +197,22 @@
 &namice_alb     !   albedo parameters
 !------------------------------------------------------------------------------
-   nn_ice_alb   =    1   !  parameterization of ice/snow albedo
-                         !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
-                         !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004),
-                         !        giving cloud-sky albedo
-   rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
-   rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
-   rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
-   rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
-   rn_alb_dpnd  =  0.27  !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' )
+   nn_ice_alb     =    1               !  parameterization of ice/snow albedo
+                                       !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
+                                       !     1: "home made" based on Brandt et al. (JClim 2005)
+                                       !                         and Grenfell & Perovich (JGR 2004), giving cloud-sky albedo
+                                       !     2: as 1 with melt ponds
+   rn_alb_sdry    =   0.85             !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
+   rn_alb_smlt    =   0.75             !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
+   rn_alb_idry    =   0.60             !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
+   rn_alb_imlt    =   0.50             !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
+   rn_alb_dpnd    =   0.27             !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' )
 /
 !------------------------------------------------------------------------------
 &namice_dia     !   Diagnostics
 !------------------------------------------------------------------------------
-   ln_limdiachk   =  .false.        !  check online the heat, mass & salt budgets (T) or not (F)
-   ln_limdiahsb   =  .false.        !  output the heat, mass & salt budgets (T) or not (F)
-   ln_limctl      =  .false.        !  ice points output for debug (T or F)
-   iiceprt        =    10           !  i-index for debug
-   jiceprt        =    10           !  j-index for debug
-/
+   ln_icediachk   = .false.            !  check online the heat, mass & salt budgets (T) or not (F)
+   ln_icediahsb   = .false.            !  output the heat, mass & salt budgets (T) or not (F)
+   ln_icectl      = .false.            !  ice points output for debug (T or F)
+   iiceprt        =  10                !  i-index for debug
+   jiceprt        =  10                !  j-index for debug
+/
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8514)
@@ -165,7 +165,7 @@
    
    !                                     !!** ice-diagnostics namelist (namice_dia) **
-   LOGICAL , PUBLIC ::   ln_limdiachk     !: flag for ice diag (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_limdiahsb     !: flag for ice diag (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_limctl        !: flag for sea-ice points output (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_icectl        !: flag for sea-ice points output (T) or not (F)
    INTEGER , PUBLIC ::   iiceprt          !: debug i-point
    INTEGER , PUBLIC ::   jiceprt          !: debug j-point
@@ -181,8 +181,8 @@
    INTEGER , PUBLIC ::   nn_UMx           !: order of the UMx advection scheme   
                                           ! -- icerdgrft -- !
-   LOGICAL , PUBLIC ::   ln_str_Hib       !: ice strength parameterization (Hibler79)
+   LOGICAL , PUBLIC ::   ln_str_H79       !: ice strength parameterization (Hibler79)
    REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79
    REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength
-   LOGICAL , PUBLIC ::   ln_str_Rot       !: ice strength parameterization (Rothrock75)
+   LOGICAL , PUBLIC ::   ln_str_R75       !: ice strength parameterization (Rothrock75)
    REAL(wp), PUBLIC ::   rn_perdg         !: ridging work divided by pot. energy change in ridging
                                           ! -- icerhg -- !
@@ -199,36 +199,38 @@
 
    !                                     !!** ice-thermodynamics namelist (namice_thd) **
-   LOGICAL , PUBLIC ::   ln_limthd        !: flag for ice thermo (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_icethd        !: flag for ice thermo (T) or not (F)
                                           ! -- icethd_dif -- !
    REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   INTEGER , PUBLIC ::   nn_ice_thcon     !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)
-   LOGICAL , PUBLIC ::   ln_dqnsice       !: change non-solar surface flux with changing surface temperature (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_cndi_U64      !: thermal conductivity: Untersteiner (1964)
+   LOGICAL , PUBLIC ::   ln_cndi_P07      !: thermal conductivity: Pringle et al (2007)
+   LOGICAL , PUBLIC ::   ln_dqns_i        !: change non-solar surface flux with changing surface temperature (T) or not (F)
    INTEGER , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1-4) or not (0)
-   REAL(wp), PUBLIC ::   rn_cdsn          !: thermal conductivity of the snow [W/m/K]
+   REAL(wp), PUBLIC ::   rn_cnd_s         !: thermal conductivity of the snow [W/m/K]
                                           ! -- icethd_dh -- !
-   LOGICAL , PUBLIC ::   ln_limdH         !: activate ice thickness change from growing/melting (T) or not (F)
-   REAL(wp), PUBLIC ::   rn_betas         !: coef. for partitioning of snowfall between leads and sea ice
+   LOGICAL , PUBLIC ::   ln_icedH         !: activate ice thickness change from growing/melting (T) or not (F)
+   REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice
                                           ! -- icethd_da -- !
-   LOGICAL , PUBLIC ::   ln_limdA         !: activate lateral melting param. (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_icedA         !: activate lateral melting param. (T) or not (F)
    REAL(wp), PUBLIC ::   rn_beta          !: coef. beta for lateral melting param.
    REAL(wp), PUBLIC ::   rn_dmin          !: minimum floe diameter for lateral melting param.
                                           ! -- icethd_lac -- !
-   LOGICAL , PUBLIC ::   ln_limdO         !: activate ice growth in open-water (T) or not (F)
-   REAL(wp), PUBLIC ::   rn_hnewice       !: thickness for new ice formation (m)
+   LOGICAL , PUBLIC ::   ln_icedO         !: activate ice growth in open-water (T) or not (F)
+   REAL(wp), PUBLIC ::   rn_hinew         !: thickness for new ice formation (m)
    LOGICAL , PUBLIC ::   ln_frazil        !: use of frazil ice collection as function of wind (T) or not (F)
-   REAL(wp), PUBLIC ::   rn_maxfrazb      !: maximum portion of frazil ice collecting at the ice bottom
-   REAL(wp), PUBLIC ::   rn_vfrazb        !: threshold drift speed for collection of bottom frazil ice
-   REAL(wp), PUBLIC ::   rn_Cfrazb        !: squeezing coefficient for collection of bottom frazil ice
+   REAL(wp), PUBLIC ::   rn_maxfraz       !: maximum portion of frazil ice collecting at the ice bottom
+   REAL(wp), PUBLIC ::   rn_vfraz         !: threshold drift speed for collection of bottom frazil ice
+   REAL(wp), PUBLIC ::   rn_Cfraz         !: squeezing coefficient for collection of bottom frazil ice
                                           ! -- iceitd_th -- !
    REAL(wp), PUBLIC ::   rn_himin         !: minimum ice thickness
                                           ! --  -- !
-   INTEGER , PUBLIC ::   nn_limflx        !: LIM3 Multi-category heat flux formulation
-   !                                      !   =-1  Use of per-category fluxes
-   !                                      !   = 0  Average per-category fluxes
-   !                                      !   = 1  Average then redistribute per-category fluxes
+   INTEGER , PUBLIC ::   nn_iceflx        !: Redistribute heat flux over ice categories
+   !                                      !   =-1  Do nothing (needs N(cat) fluxes)
+   !                                      !   = 0  Average N(cat) fluxes then apply the average over the N(cat) ice 
+   !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice
+   !                                      !                                   using T-ice and albedo sensitivity
    !                                      !   = 2  Redistribute a single flux over categories
 
    !                                     !!** ice-salinity namelist (namice_sal) **
-   LOGICAL , PUBLIC ::   ln_limdS         !: activate gravity drainage and flushing (T) or not (F)
+   LOGICAL , PUBLIC ::   ln_icedS         !: activate gravity drainage and flushing (T) or not (F)
    INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model
    !                                      ! 1 - constant salinity in both space and time
@@ -354,5 +356,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sm_i      !: Sea-Ice Bulk salinity (ppt)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   smv_i     !: Sea-Ice Bulk salinity times volume per area (ppt.m)
-   !                                                                    !  this is an extensive variable that has to be transported
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age (s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area (s)
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8514)
@@ -37,10 +37,10 @@
    PUBLIC   ice_adv_init   ! called by icestp
 
-   INTEGER, PUBLIC ::              nice_dyn   ! choice of the type of advection scheme
+   INTEGER, PUBLIC ::              nice_dyn   !: choice of the type of advection scheme
    !                                  ! associated indices:
-   INTEGER, PUBLIC, PARAMETER ::   np_dynNO   = 0   ! no ice dynamics and ice advection
-   INTEGER, PUBLIC, PARAMETER ::   np_dynFULL = 1   ! full ice dynamics  (rheology + advection + ridging/rafting + correction)
-   INTEGER, PUBLIC, PARAMETER ::   np_dyn     = 2   ! no ridging/rafting (rheology + advection                   + correction)
-   INTEGER, PUBLIC, PARAMETER ::   np_dynPURE = 3   ! pure dynamics      (rheology + advection) 
+   INTEGER, PUBLIC, PARAMETER ::   np_dynNO   = 0   !: no ice dynamics and ice advection
+   INTEGER, PUBLIC, PARAMETER ::   np_dynFULL = 1   !: full ice dynamics  (rheology + advection + ridging/rafting + correction)
+   INTEGER, PUBLIC, PARAMETER ::   np_dyn     = 2   !: no ridging/rafting (rheology + advection                   + correction)
+   INTEGER, PUBLIC, PARAMETER ::   np_dynPURE = 3   !: pure dynamics      (rheology + advection) 
 
    !! * Substitution
@@ -90,5 +90,5 @@
 
       ! conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       
       ! store old values for diag
@@ -224,10 +224,10 @@
       
       ! conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
         
       ! --------------
       ! control prints
       ! --------------
-      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
+      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
       !
       IF( nn_timing == 1 )  CALL timing_stop('iceadv')
@@ -264,5 +264,5 @@
       IF(lwp) THEN                     ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_adv_init : ice parameters for ice dynamics '
+         WRITE(numout,*) 'ice_adv_init: ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_adv'
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8514)
@@ -328,6 +328,6 @@
       IF(lwp) THEN                      ! Control print
          WRITE(numout,*)
-         WRITE(numout,*) 'albedo : set albedo parameters'
-         WRITE(numout,*) '~~~~~~~'
+         WRITE(numout,*) 'ice_alb_init: set albedo parameters'
+         WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_alb : albedo '
          WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8514)
@@ -64,5 +64,5 @@
       ENDIF
       !                             !--- conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
       !
@@ -86,5 +86,5 @@
     
       !                             !-----------------------------------------------------
-      IF (  nn_icesal == 2  ) THEN  !  salinity stays in bounds [Simin,Simax]            !
+      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        !
       !                             !-----------------------------------------------------
          zzc = rhoic * r1_rdtice
@@ -176,9 +176,9 @@
       !
       !                                !--- conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
       !                                !--- control prints
       IF( ln_ctl )                    CALL ice_prt3D( 'icecor' )
-      IF( ln_limctl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
+      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
       !
       IF( nn_timing == 1 )   CALL timing_stop('icecor')
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8514)
@@ -54,5 +54,5 @@
       !!                     + test if ice concentration and volume are > 0
       !!
-      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiachk=true
+      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
       !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to
@@ -167,5 +167,5 @@
       !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
       !!
-      !! ** Method  : This is an online diagnostics which can be activated with ln_limdiachk=true
+      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
       !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8514)
@@ -164,5 +164,5 @@
       INTEGER            ::   ios, ierror   ! local integer
       !!
-      NAMELIST/namice_dia/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt  
+      NAMELIST/namice_dia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt  
       !!----------------------------------------------------------------------
       !
@@ -181,11 +181,11 @@
          WRITE(numout,*) ' ~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_dia : '
-         WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_limdiachk = ', ln_limdiachk
-         WRITE(numout,*) '      Output          heat/mass/salt budget      ln_limdiahsb = ', ln_limdiahsb
-         WRITE(numout,*) '      control prints for a given grid point         ln_limctl = ', ln_limctl
+         WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_icediachk = ', ln_icediachk
+         WRITE(numout,*) '      Output          heat/mass/salt budget      ln_icediahsb = ', ln_icediahsb
+         WRITE(numout,*) '      control prints for a given grid point         ln_icectl = ', ln_icectl
          WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')'
       ENDIF
       !      
-      IF( ln_limdiahsb ) THEN
+      IF( ln_icediahsb ) THEN
          ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror )
          IF( ierror > 0 )  THEN
@@ -202,5 +202,5 @@
    SUBROUTINE ice_dia_rst( cdrw, kt )
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE limdia_rst  ***
+      !!                   ***  ROUTINE icedia_rst  ***
       !!                     
       !! ** Purpose :   Read or write DIA file in restart file
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8514)
@@ -143,9 +143,9 @@
                                 CALL blk_ice_flx( t_su, alb_ice )    ! 
          IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
-         IF( nn_limflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+         IF( nn_iceflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
       CASE ( jp_purecpl )           !--- coupled formulation
                                 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su )
-         IF( nn_limflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
+         IF( nn_iceflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx )
          !
       END SELECT
@@ -168,5 +168,5 @@
 
 
-   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_limflx )
+   SUBROUTINE ice_flx_dist( ptn_ice, palb_ice, pqns_ice, pqsr_ice, pdqn_ice, pevap_ice, pdevap_ice, k_iceflx )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE ice_flx_dist  ***
@@ -177,8 +177,12 @@
       !! ** Method  :   average then redistribute
       !!
-      !! ** Action  :
-      !!---------------------------------------------------------------------
-      INTEGER                   , INTENT(in   ) ::   k_limflx   ! =-1 do nothing; =0 average ;
-      !                                                         ! = 1 average and redistribute ; =2 redistribute
+      !! ** Action  :   depends on k_iceflx
+      !!                = -1  Do nothing (needs N(cat) fluxes)
+      !!                =  0  Average N(cat) fluxes then apply the average over the N(cat) ice 
+      !!                =  1  Average N(cat) fluxes then redistribute over the N(cat) ice
+      !!                                                 using T-ice and albedo sensitivity
+      !!                =  2  Redistribute a single flux over categories
+      !!---------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   k_iceflx   ! redistributor
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptn_ice    ! ice surface temperature
       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   palb_ice   ! ice albedo
@@ -208,5 +212,5 @@
       END WHERE
       
-      SELECT CASE( k_limflx )       !==  averaged on all ice categories  ==!
+      SELECT CASE( k_iceflx )       !==  averaged on all ice categories  ==!
       !
       CASE( 0 , 1 )
@@ -231,5 +235,5 @@
       END SELECT
       !
-      SELECT CASE( k_limflx )       !==  redistribution on all ice categories  ==!
+      SELECT CASE( k_iceflx )       !==  redistribution on all ice categories  ==!
       !
       CASE( 1 , 2 )
Index: anches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90	(revision 8513)
+++ 	(revision )
@@ -1,610 +1,0 @@
-MODULE iceist
-   !!======================================================================
-   !!                     ***  MODULE  iceist  ***
-   !!              Initialisation of diagnostics ice variables
-   !!======================================================================
-   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
-   !!            3.0  ! 2007    (M. Vancoppenolle)  Rewrite for ice cats
-   !!            3.0  ! 2009-11 (M. Vancoppenolle)  Enhanced version for ice cats
-   !!            3.0  ! 2011-02 (G. Madec)  dynamical allocation
-   !!             -   ! 2014    (C. Rousset)  add N/S initializations
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                       LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   ice_ist       :  initialization of diagnostics ice variables
-   !!   ice_ist_init  :  initialization of ice state and namelist read
-   !!----------------------------------------------------------------------
-   USE phycst         ! physical constant
-   USE oce            ! dynamics and tracers variables
-   USE dom_oce        ! ocean domain
-   USE sbc_oce , ONLY : sst_m, sss_m 
-   USE sbc_ice , ONLY : tn_ice
-   USE eosbn2         ! equation of state
-   USE ice            ! sea-ice variables
-   USE par_oce        ! ocean parameters
-   USE icevar         ! ice_var_salprof
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-   USE fldread        ! read input fields
-   USE iom
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_ist      ! called by icestp.F90
-
-   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
-   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
-   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
-   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
-   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
-   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
-   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
-   !
-   ! ** namelist (namice_ini) **
-   LOGICAL  ::   ln_iceini        ! initialization or not
-   LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file
-   REAL(wp) ::   rn_thres_sst     ! threshold water temperature for initial sea ice
-   REAL(wp) ::   rn_hts_ini_n     ! initial snow thickness in the north
-   REAL(wp) ::   rn_hts_ini_s     ! initial snow thickness in the south
-   REAL(wp) ::   rn_hti_ini_n     ! initial ice thickness in the north
-   REAL(wp) ::   rn_hti_ini_s     ! initial ice thickness in the south
-   REAL(wp) ::   rn_ati_ini_n     ! initial leads area in the north
-   REAL(wp) ::   rn_ati_ini_s     ! initial leads area in the south
-   REAL(wp) ::   rn_smi_ini_n     ! initial salinity 
-   REAL(wp) ::   rn_smi_ini_s     ! initial salinity
-   REAL(wp) ::   rn_tmi_ini_n     ! initial temperature
-   REAL(wp) ::   rn_tmi_ini_s     ! initial temperature
-   
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
-   !! $Id: iceist.F90 8378 2017-07-26 13:55:59Z clem $
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-!!gm  better name:  ice_istate
-   SUBROUTINE ice_ist
-      !!-------------------------------------------------------------------
-      !!                    ***  ROUTINE ice_ist  ***
-      !!
-      !! ** Purpose :   defined the sea-ice initial state
-      !!
-      !! ** Method  :   This routine will put some ice where ocean
-      !!                is at the freezing point, then fill in ice 
-      !!                state variables using prescribed initial 
-      !!                values in the namelist            
-      !!
-      !! ** Steps   :   1) Read namelist
-      !!                2) Basal temperature; ice and hemisphere masks
-      !!                3) Fill in the ice thickness distribution using gaussian
-      !!                4) Fill in space-dependent arrays for state variables
-      !!                5) Diagnostic arrays
-      !!                6) Lateral boundary conditions
-      !!
-      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
-      !!              where there is no ice (clem: I do not know why, is it mandatory?) 
-      !!--------------------------------------------------------------------
-      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
-      REAL(wp)   :: ztmelts, zdh
-      INTEGER    :: i_hemis, i_fill, jl0
-      REAL(wp)   :: zarg, zV, zconv, zdv
-      REAL(wp), DIMENSION(jpi,jpj)     :: zswitch    ! ice indicator
-      REAL(wp), DIMENSION(jpi,jpj)     :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj)     :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
-      INTEGER , DIMENSION(4)           :: itest
-      !--------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'ice_ist : sea-ice initialization '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~ '
-
-      !--------------------------------------------------------------------
-      ! 1) Read namelist
-      !--------------------------------------------------------------------
-      !
-      CALL ice_ist_init
-
-      ! init surface temperature
-      DO jl = 1, jpl
-         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
-         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
-      END DO
-
-      ! init basal temperature (considered at freezing point)
-      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
-      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
-
-
-      !--------------------------------------------------------------------
-      ! 2) Initialization of sea ice state variables
-      !--------------------------------------------------------------------
-      IF( ln_iceini ) THEN
-         !
-         IF( ln_iceini_file )THEN
-         !
-            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
-            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
-            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
-            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
-            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
-            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
-            !
-            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 
-            ELSEWHERE                       ; zswitch(:,:) = 0._wp
-            END WHERE
-            !
-         ELSE ! ln_iceini_file = F
-
-            !--------------------------------------------------------------------
-            ! 3) Basal temperature, ice mask
-            !--------------------------------------------------------------------
-            ! no ice if sst <= t-freez + ttest
-            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 
-            ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1)
-            END WHERE
-
-            !-----------------------------
-            ! 3.1) Hemisphere-dependent arrays
-            !-----------------------------
-            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( ff_t(ji,jj) >= 0._wp ) THEN
-                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
-                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
-                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
-                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
-                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
-                  ELSE
-                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
-                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
-                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
-                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
-                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
-                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
-                  ENDIF
-               END DO
-            END DO
-            !
-         ENDIF ! ln_iceini_file
-         
-         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
-         !---------------------------------------------------------------------
-         ! 3.2) Distribute ice concentration and thickness into the categories
-         !---------------------------------------------------------------------
-         ! a gaussian distribution for ice concentration is used
-         ! then we check whether the distribution fullfills
-         ! volume and area conservation, positivity and ice categories bounds
-         zh_i_ini(:,:,:) = 0._wp 
-         za_i_ini(:,:,:) = 0._wp
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !
-               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
-
-                  !--- jl0: most likely index where cc will be maximum
-                  jl0 = jpl
-                  DO jl = 1, jpl
-                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
-                        jl0 = jl
-                        CYCLE
-                     ENDIF
-                  END DO
-                  !
-                  ! initialisation of tests
-                  itest(:)  = 0
-                  
-                  i_fill = jpl + 1                                             !====================================
-                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
-                     ! iteration                                               !====================================
-                     i_fill = i_fill - 1
-
-                     ! initialisation of ice variables for each try
-                     zh_i_ini(ji,jj,:) = 0._wp 
-                     za_i_ini(ji,jj,:) = 0._wp
-                     itest(:) = 0
-                     !
-                     ! *** case very thin ice: fill only category 1
-                     IF ( i_fill == 1 ) THEN
-                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
-                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
-
-                     ! *** case ice is thicker: fill categories >1
-                     ELSE
-
-                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
-                        DO jl = 1, i_fill-1
-                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
-                        END DO
-                        !
-                        !--- Concentrations
-                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
-                        DO jl = 1, i_fill - 1
-                           IF( jl /= jl0 )THEN
-                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
-                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
-                           ENDIF
-                        END DO
-                        !
-                        ! Concentration in the last (i_fill) category
-                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
-
-                        ! Ice thickness in the last (i_fill) category
-                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
-                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
-
-                        ! clem: correction if concentration of upper cat is greater than lower cat
-                        !       (it should be a gaussian around jl0 but sometimes it is not)
-                        IF ( jl0 /= jpl ) THEN
-                           DO jl = jpl, jl0+1, -1
-                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
-                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
-                                 zh_i_ini(ji,jj,jl    ) = 0._wp
-                                 za_i_ini(ji,jj,jl    ) = 0._wp
-                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
-                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
-                              END IF
-                           ENDDO
-                        ENDIF
-                        !
-                     ENDIF ! case ice is thick or thin
-
-                     !---------------------
-                     ! Compatibility tests
-                     !---------------------
-                     ! Test 1: area conservation
-                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
-                     IF ( zconv < epsi06 ) itest(1) = 1
-                     
-                     ! Test 2: volume conservation
-                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &
-                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
-                     IF ( zconv < epsi06 ) itest(2) = 1
-                     
-                     ! Test 3: thickness of the last category is in-bounds ?
-                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
-                     
-                     ! Test 4: positivity of ice concentrations
-                     itest(4) = 1
-                     DO jl = 1, i_fill
-                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
-                     END DO
-                     !                                      !============================
-                  END DO                                    ! end iteration on categories
-                  !                                         !============================
-                  !
-                  IF( lwp .AND. SUM(itest) /= 4 ) THEN 
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
-                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
-                     WRITE(numout,*)
-                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
-                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
-                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
-                  ENDIF
-               
-               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
-               !
-            END DO   
-         END DO   
-
-         !---------------------------------------------------------------------
-         ! 3.3) Space-dependent arrays for ice state variables
-         !---------------------------------------------------------------------
-
-         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
-         DO jl = 1, jpl ! loop over categories
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
-                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
-                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
-                  o_i(ji,jj,jl)   = 0._wp                                                     ! age (0 day)
-                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
-
-                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
-                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
-                  ELSE
-                    ht_s(ji,jj,jl)= 0._wp
-                  ENDIF
-
-                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
-                  ! In case snow load is in excess that would lead to transformation from snow to ice
-                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
-                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
-                  ! recompute ht_i, ht_s avoiding out of bounds values
-                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
-                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
-
-                  ! ice volume, salt content, age content
-                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
-                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
-                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
-                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
-               END DO
-            END DO
-         END DO
-
-         ! for constant salinity in time
-         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
-            CALL ice_var_salprof
-            smv_i = sm_i * v_i
-         ENDIF
-            
-         ! Snow temperature and heat content
-         DO jk = 1, nlay_s
-            DO jl = 1, jpl ! loop over categories
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
-                     ! Snow energy of melting
-                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
-
-                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
-                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
-                  END DO
-               END DO
-            END DO
-         END DO
-
-         ! Ice salinity, temperature and heat content
-         DO jk = 1, nlay_i
-            DO jl = 1, jpl ! loop over categories
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
-                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
-                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
-
-                     ! heat content per unit volume
-                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
-                        -   rcp     * ( ztmelts - rt0 ) )
-
-                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
-                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
-                  END DO
-               END DO
-            END DO
-         END DO
-
-         tn_ice (:,:,:) = t_su (:,:,:)
-
-         ! MV MP 2016
-         ! Melt pond volume and fraction
-         IF ( ln_pnd ) THEN
-
-            DO jl = 1, jpl
-
-               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
-               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
-               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
-               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
-
-            END DO
-
-         ELSE
-
-            a_ip(:,:,:)      = 0._wp
-            v_ip(:,:,:)      = 0._wp
-            a_ip_frac(:,:,:) = 0._wp
-            h_ip     (:,:,:) = 0._wp
-
-         ENDIF
-         ! END MV MP 2016
-
-      ELSE ! if ln_iceini=false
-         a_i  (:,:,:) = 0._wp
-         v_i  (:,:,:) = 0._wp
-         v_s  (:,:,:) = 0._wp
-         smv_i(:,:,:) = 0._wp
-         oa_i (:,:,:) = 0._wp
-         ht_i (:,:,:) = 0._wp
-         ht_s (:,:,:) = 0._wp
-         sm_i (:,:,:) = 0._wp
-         o_i  (:,:,:) = 0._wp
-
-         e_i(:,:,:,:) = 0._wp
-         e_s(:,:,:,:) = 0._wp
-
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
-            END DO
-            DO jk = 1, nlay_s
-               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
-            END DO
-         END DO
-
-         a_ip(:,:,:)      = 0._wp
-         v_ip(:,:,:)      = 0._wp
-         a_ip_frac(:,:,:) = 0._wp
-         h_ip     (:,:,:) = 0._wp
-
-      ENDIF ! ln_iceini
-      
-      at_i (:,:) = 0.0_wp
-      DO jl = 1, jpl
-         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
-      END DO
-
-      !--------------------------------------------------------------------
-      ! 4) Global ice variables for output diagnostics                    | 
-      !--------------------------------------------------------------------
-      u_ice (:,:)     = 0._wp
-      v_ice (:,:)     = 0._wp
-      stress1_i(:,:)  = 0._wp
-      stress2_i(:,:)  = 0._wp
-      stress12_i(:,:) = 0._wp
-
-      !--------------------------------------------------------------------
-      ! 5) Moments for advection
-      !--------------------------------------------------------------------
-
-      sxopw (:,:) = 0._wp 
-      syopw (:,:) = 0._wp 
-      sxxopw(:,:) = 0._wp 
-      syyopw(:,:) = 0._wp 
-      sxyopw(:,:) = 0._wp
-
-      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
-      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
-      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
-      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
-      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
-
-      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
-      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
-      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
-      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
-      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
-
-      sxsal  (:,:,:)  = 0._wp
-      sysal  (:,:,:)  = 0._wp
-      sxxsal (:,:,:)  = 0._wp
-      syysal (:,:,:)  = 0._wp
-      sxysal (:,:,:)  = 0._wp
-
-      sxage  (:,:,:)  = 0._wp
-      syage  (:,:,:)  = 0._wp
-      sxxage (:,:,:)  = 0._wp
-      syyage (:,:,:)  = 0._wp
-      sxyage (:,:,:)  = 0._wp
-
-      !------------------------------------
-      ! 6) store fields at before time-step
-      !------------------------------------
-      ! it is only necessary for the 1st interpolation by Agrif
-      a_i_b  (:,:,:)   = a_i  (:,:,:)
-      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
-      v_i_b  (:,:,:)   = v_i  (:,:,:)
-      v_s_b  (:,:,:)   = v_s  (:,:,:)
-      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
-      smv_i_b(:,:,:)   = smv_i(:,:,:)
-      oa_i_b (:,:,:)   = oa_i (:,:,:)
-      u_ice_b(:,:)     = u_ice(:,:)
-      v_ice_b(:,:)     = v_ice(:,:)
-
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
-         sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
-         syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
-         sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
-         syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
-         sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
-      ENDIF
-      ! END MV MP 2016
-
-!!!clem
-!!      ! Output the initial state and forcings
-!!      CALL dia_wri_state( 'output.init', nit000 )
-!!!      
-
-   END SUBROUTINE ice_ist
-
-   SUBROUTINE ice_ist_init
-      !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_ist_init  ***
-      !!        
-      !! ** Purpose : Definition of initial state of the ice 
-      !!
-      !! ** Method : Read the namice_ini namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
-      !!
-      !! ** input : 
-      !!        Namelist namice_ini
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code 
-      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
-      !!-----------------------------------------------------------------------------
-      !
-      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
-      INTEGER :: ji,jj
-      INTEGER :: ifpr, ierror
-      !
-      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
-      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
-      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
-      !
-      NAMELIST/namice_ini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
-         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
-         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
-         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
-      !!-----------------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )              ! Namelist namice_ini in reference namelist : Ice initial state
-      READ  ( numnam_ice_ref, namice_ini, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_ini in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namice_ini in configuration namelist : Ice initial state
-      READ  ( numnam_ice_cfg, namice_ini, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_ini in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_ini )
-
-      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
-      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
-      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
-
-      ! Define the initial parameters
-      ! -------------------------
-
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_ist_init : ice parameters inititialisation '
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_ini'
-         WRITE(numout,*) '      initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
-         WRITE(numout,*) '      ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file
-         WRITE(numout,*) '      threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
-         WRITE(numout,*) '      initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
-         WRITE(numout,*) '      initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
-         WRITE(numout,*) '      initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
-         WRITE(numout,*) '      initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
-         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
-         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
-         WRITE(numout,*) '      initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
-         WRITE(numout,*) '      initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
-         WRITE(numout,*) '      initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
-         WRITE(numout,*) '      initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
-      ENDIF
-
-      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
-         !
-         ! set si structure
-         ALLOCATE( si(jpfldi), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'Ice_ini in iceist: unable to allocate si structure' )   ;   RETURN
-         ENDIF
-
-         DO ifpr = 1, jpfldi
-            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
-            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
-         END DO
-
-         ! fill si with slf_i and control print
-         CALL fld_fill( si, slf_i, cn_dir, 'ice_ist', 'lim istate ini', 'numnam_ice' )
-
-         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
-
-      ENDIF
-
-   END SUBROUTINE ice_ist_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :         Empty module          NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE iceist
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8514)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8514)
@@ -0,0 +1,609 @@
+MODULE iceistate
+   !!======================================================================
+   !!                     ***  MODULE  iceistate  ***
+   !!              Initialisation of diagnostics ice variables
+   !!======================================================================
+   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
+   !!            3.0  ! 2007    (M. Vancoppenolle)  Rewrite for ice cats
+   !!            3.0  ! 2009-11 (M. Vancoppenolle)  Enhanced version for ice cats
+   !!            3.0  ! 2011-02 (G. Madec)  dynamical allocation
+   !!             -   ! 2014    (C. Rousset)  add N/S initializations
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                       LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_istate       :  initialization of diagnostics ice variables
+   !!   ice_istate_init  :  initialization of ice state and namelist read
+   !!----------------------------------------------------------------------
+   USE phycst         ! physical constant
+   USE oce            ! dynamics and tracers variables
+   USE dom_oce        ! ocean domain
+   USE sbc_oce , ONLY : sst_m, sss_m 
+   USE sbc_ice , ONLY : tn_ice
+   USE eosbn2         ! equation of state
+   USE ice            ! sea-ice variables
+   USE par_oce        ! ocean parameters
+   USE icevar         ! ice_var_salprof
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE fldread        ! read input fields
+   USE iom
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_istate      ! called by icestp.F90
+
+   INTEGER , PARAMETER ::   jpfldi = 6           ! maximum number of files to read
+   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness (m)    at T-point
+   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snow thicknes (m)    at T-point
+   INTEGER , PARAMETER ::   jp_ati = 3           ! index of ice fraction (%) at T-point
+   INTEGER , PARAMETER ::   jp_tsu = 4           ! index of ice surface temp (K)    at T-point
+   INTEGER , PARAMETER ::   jp_tmi = 5           ! index of ice temp at T-point
+   INTEGER , PARAMETER ::   jp_smi = 6           ! index of ice sali at T-point
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
+   !
+   ! ** namelist (namice_ini) **
+   LOGICAL  ::   ln_iceini        ! initialization or not
+   LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file
+   REAL(wp) ::   rn_thres_sst     ! threshold water temperature for initial sea ice
+   REAL(wp) ::   rn_hts_ini_n     ! initial snow thickness in the north
+   REAL(wp) ::   rn_hts_ini_s     ! initial snow thickness in the south
+   REAL(wp) ::   rn_hti_ini_n     ! initial ice thickness in the north
+   REAL(wp) ::   rn_hti_ini_s     ! initial ice thickness in the south
+   REAL(wp) ::   rn_ati_ini_n     ! initial leads area in the north
+   REAL(wp) ::   rn_ati_ini_s     ! initial leads area in the south
+   REAL(wp) ::   rn_smi_ini_n     ! initial salinity 
+   REAL(wp) ::   rn_smi_ini_s     ! initial salinity
+   REAL(wp) ::   rn_tmi_ini_n     ! initial temperature
+   REAL(wp) ::   rn_tmi_ini_s     ! initial temperature
+   
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
+   !! $Id: iceistate.F90 8378 2017-07-26 13:55:59Z clem $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_istate
+      !!-------------------------------------------------------------------
+      !!                    ***  ROUTINE ice_istate  ***
+      !!
+      !! ** Purpose :   defined the sea-ice initial state
+      !!
+      !! ** Method  :   This routine will put some ice where ocean
+      !!                is at the freezing point, then fill in ice 
+      !!                state variables using prescribed initial 
+      !!                values in the namelist            
+      !!
+      !! ** Steps   :   1) Read namelist
+      !!                2) Basal temperature; ice and hemisphere masks
+      !!                3) Fill in the ice thickness distribution using gaussian
+      !!                4) Fill in space-dependent arrays for state variables
+      !!                5) Diagnostic arrays
+      !!                6) Lateral boundary conditions
+      !!
+      !! ** Notes   : o_i, t_su, t_s, t_i, s_i must be filled everywhere, even
+      !!              where there is no ice (clem: I do not know why, is it mandatory?) 
+      !!--------------------------------------------------------------------
+      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
+      REAL(wp)   :: ztmelts, zdh
+      INTEGER    :: i_hemis, i_fill, jl0
+      REAL(wp)   :: zarg, zV, zconv, zdv
+      REAL(wp), DIMENSION(jpi,jpj)     :: zswitch    ! ice indicator
+      REAL(wp), DIMENSION(jpi,jpj)     :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj)     :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
+      INTEGER , DIMENSION(4)           :: itest
+      !--------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'ice_istate : sea-ice initialization '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~ '
+
+      !--------------------------------------------------------------------
+      ! 1) Read namelist
+      !--------------------------------------------------------------------
+      !
+      CALL ice_istate_init
+
+      ! init surface temperature
+      DO jl = 1, jpl
+         t_su  (:,:,jl) = rt0 * tmask(:,:,1)
+         tn_ice(:,:,jl) = rt0 * tmask(:,:,1)
+      END DO
+
+      ! init basal temperature (considered at freezing point)
+      CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
+      t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
+
+
+      !--------------------------------------------------------------------
+      ! 2) Initialization of sea ice state variables
+      !--------------------------------------------------------------------
+      IF( ln_iceini ) THEN
+         !
+         IF( ln_iceini_file )THEN
+         !
+            zht_i_ini(:,:)  = si(jp_hti)%fnow(:,:,1)
+            zht_s_ini(:,:)  = si(jp_hts)%fnow(:,:,1)
+            zat_i_ini(:,:)  = si(jp_ati)%fnow(:,:,1)
+            zts_u_ini(:,:)  = si(jp_tsu)%fnow(:,:,1)
+            ztm_i_ini(:,:)  = si(jp_tmi)%fnow(:,:,1)
+            zsm_i_ini(:,:)  = si(jp_smi)%fnow(:,:,1)
+            !
+            WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 
+            ELSEWHERE                       ; zswitch(:,:) = 0._wp
+            END WHERE
+            !
+         ELSE ! ln_iceini_file = F
+
+            !--------------------------------------------------------------------
+            ! 3) Basal temperature, ice mask
+            !--------------------------------------------------------------------
+            ! no ice if sst <= t-freez + ttest
+            WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 
+            ELSEWHERE                                                                  ; zswitch(:,:) = tmask(:,:,1)
+            END WHERE
+
+            !-----------------------------
+            ! 3.1) Hemisphere-dependent arrays
+            !-----------------------------
+            ! assign initial thickness, concentration, snow depth and salinity to an hemisphere-dependent array
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( ff_t(ji,jj) >= 0._wp ) THEN
+                     zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj)
+                     zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj)
+                     zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj)
+                     zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
+                     zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj)
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj)
+                  ELSE
+                     zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj)
+                     zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj)
+                     zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj)
+                     zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
+                     zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj)
+                     ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            !
+         ENDIF ! ln_iceini_file
+         
+         zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:)   ! ice volume
+         !---------------------------------------------------------------------
+         ! 3.2) Distribute ice concentration and thickness into the categories
+         !---------------------------------------------------------------------
+         ! a gaussian distribution for ice concentration is used
+         ! then we check whether the distribution fullfills
+         ! volume and area conservation, positivity and ice categories bounds
+         zh_i_ini(:,:,:) = 0._wp 
+         za_i_ini(:,:,:) = 0._wp
+         !
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN
+
+                  !--- jl0: most likely index where cc will be maximum
+                  jl0 = jpl
+                  DO jl = 1, jpl
+                     IF ( ( zht_i_ini(ji,jj) >  hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN
+                        jl0 = jl
+                        CYCLE
+                     ENDIF
+                  END DO
+                  !
+                  ! initialisation of tests
+                  itest(:)  = 0
+                  
+                  i_fill = jpl + 1                                             !====================================
+                  DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
+                     ! iteration                                               !====================================
+                     i_fill = i_fill - 1
+
+                     ! initialisation of ice variables for each try
+                     zh_i_ini(ji,jj,:) = 0._wp 
+                     za_i_ini(ji,jj,:) = 0._wp
+                     itest(:) = 0
+                     !
+                     ! *** case very thin ice: fill only category 1
+                     IF ( i_fill == 1 ) THEN
+                        zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj)
+                        za_i_ini(ji,jj,1) = zat_i_ini(ji,jj)
+
+                     ! *** case ice is thicker: fill categories >1
+                     ELSE
+
+                        ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
+                        DO jl = 1, i_fill-1
+                           zh_i_ini(ji,jj,jl) = hi_mean(jl)
+                        END DO
+                        !
+                        !--- Concentrations
+                        za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl))
+                        DO jl = 1, i_fill - 1
+                           IF( jl /= jl0 )THEN
+                              zarg               = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) )
+                              za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2)
+                           ENDIF
+                        END DO
+                        !
+                        ! Concentration in the last (i_fill) category
+                        za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) )
+
+                        ! Ice thickness in the last (i_fill) category
+                        zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) )
+                        zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 
+
+                        ! clem: correction if concentration of upper cat is greater than lower cat
+                        !       (it should be a gaussian around jl0 but sometimes it is not)
+                        IF ( jl0 /= jpl ) THEN
+                           DO jl = jpl, jl0+1, -1
+                              IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN
+                                 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl)
+                                 zh_i_ini(ji,jj,jl    ) = 0._wp
+                                 za_i_ini(ji,jj,jl    ) = 0._wp
+                                 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1)  &
+                                    &                     + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 )
+                              END IF
+                           ENDDO
+                        ENDIF
+                        !
+                     ENDIF ! case ice is thick or thin
+
+                     !---------------------
+                     ! Compatibility tests
+                     !---------------------
+                     ! Test 1: area conservation
+                     zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) )
+                     IF ( zconv < epsi06 ) itest(1) = 1
+                     
+                     ! Test 2: volume conservation
+                     zconv = ABS(       zat_i_ini(ji,jj)       * zht_i_ini(ji,jj)   &
+                        &        - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) )
+                     IF ( zconv < epsi06 ) itest(2) = 1
+                     
+                     ! Test 3: thickness of the last category is in-bounds ?
+                     IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
+                     
+                     ! Test 4: positivity of ice concentrations
+                     itest(4) = 1
+                     DO jl = 1, i_fill
+                        IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0
+                     END DO
+                     !                                      !============================
+                  END DO                                    ! end iteration on categories
+                  !                                         !============================
+                  !
+                  IF( lwp .AND. SUM(itest) /= 4 ) THEN 
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' !!!! ALERT itest is not equal to 4      !!! '
+                     WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure '
+                     WRITE(numout,*)
+                     WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:)
+                     WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj)
+                     WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj)
+                  ENDIF
+               
+               ENDIF !  zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp
+               !
+            END DO   
+         END DO   
+
+         !---------------------------------------------------------------------
+         ! 3.3) Space-dependent arrays for ice state variables
+         !---------------------------------------------------------------------
+
+         ! Ice concentration, thickness and volume, ice salinity, ice age, surface temperature
+         DO jl = 1, jpl ! loop over categories
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  a_i(ji,jj,jl)   = zswitch(ji,jj) * za_i_ini(ji,jj,jl)                       ! concentration
+                  ht_i(ji,jj,jl)  = zswitch(ji,jj) * zh_i_ini(ji,jj,jl)                       ! ice thickness
+                  sm_i(ji,jj,jl)  = zswitch(ji,jj) * zsm_i_ini(ji,jj)                         ! salinity
+                  o_i(ji,jj,jl)   = 0._wp                                                     ! age (0 day)
+                  t_su(ji,jj,jl)  = zswitch(ji,jj) * zts_u_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 ! surf temp
+
+                  IF( zht_i_ini(ji,jj) > 0._wp )THEN
+                    ht_s(ji,jj,jl)= ht_i(ji,jj,jl) * ( zht_s_ini(ji,jj) / zht_i_ini(ji,jj) )  ! snow depth
+                  ELSE
+                    ht_s(ji,jj,jl)= 0._wp
+                  ENDIF
+
+                  ! This case below should not be used if (ht_s/ht_i) is ok in namelist
+                  ! In case snow load is in excess that would lead to transformation from snow to ice
+                  ! Then, transfer the snow excess into the ice (different from icethd_dh)
+                  zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 ) 
+                  ! recompute ht_i, ht_s avoiding out of bounds values
+                  ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
+                  ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
+
+                  ! ice volume, salt content, age content
+                  v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
+                  v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
+                  smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
+                  oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
+               END DO
+            END DO
+         END DO
+
+         ! for constant salinity in time
+         IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
+            CALL ice_var_salprof
+            smv_i = sm_i * v_i
+         ENDIF
+            
+         ! Snow temperature and heat content
+         DO jk = 1, nlay_s
+            DO jl = 1, jpl ! loop over categories
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     t_s(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0
+                     ! Snow energy of melting
+                     e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+
+                     ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
+                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
+                  END DO
+               END DO
+            END DO
+         END DO
+
+         ! Ice salinity, temperature and heat content
+         DO jk = 1, nlay_i
+            DO jl = 1, jpl ! loop over categories
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     t_i(ji,jj,jk,jl) = zswitch(ji,jj) * ztm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rt0 
+                     s_i(ji,jj,jk,jl) = zswitch(ji,jj) * zsm_i_ini(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * rn_simin
+                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K
+
+                     ! heat content per unit volume
+                     e_i(ji,jj,jk,jl) = zswitch(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                        +   lfus    * ( 1._wp - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
+                        -   rcp     * ( ztmelts - rt0 ) )
+
+                     ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
+                  END DO
+               END DO
+            END DO
+         END DO
+
+         tn_ice (:,:,:) = t_su (:,:,:)
+
+         ! MV MP 2016
+         ! Melt pond volume and fraction
+         IF ( ln_pnd ) THEN
+
+            DO jl = 1, jpl
+
+               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
+               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
+               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
+               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
+
+            END DO
+
+         ELSE
+
+            a_ip(:,:,:)      = 0._wp
+            v_ip(:,:,:)      = 0._wp
+            a_ip_frac(:,:,:) = 0._wp
+            h_ip     (:,:,:) = 0._wp
+
+         ENDIF
+         ! END MV MP 2016
+
+      ELSE ! if ln_iceini=false
+         a_i  (:,:,:) = 0._wp
+         v_i  (:,:,:) = 0._wp
+         v_s  (:,:,:) = 0._wp
+         smv_i(:,:,:) = 0._wp
+         oa_i (:,:,:) = 0._wp
+         ht_i (:,:,:) = 0._wp
+         ht_s (:,:,:) = 0._wp
+         sm_i (:,:,:) = 0._wp
+         o_i  (:,:,:) = 0._wp
+
+         e_i(:,:,:,:) = 0._wp
+         e_s(:,:,:,:) = 0._wp
+
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               t_i(:,:,jk,jl) = rt0 * tmask(:,:,1)
+            END DO
+            DO jk = 1, nlay_s
+               t_s(:,:,jk,jl) = rt0 * tmask(:,:,1)
+            END DO
+         END DO
+
+         a_ip(:,:,:)      = 0._wp
+         v_ip(:,:,:)      = 0._wp
+         a_ip_frac(:,:,:) = 0._wp
+         h_ip     (:,:,:) = 0._wp
+
+      ENDIF ! ln_iceini
+      
+      at_i (:,:) = 0.0_wp
+      DO jl = 1, jpl
+         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
+      END DO
+
+      !--------------------------------------------------------------------
+      ! 4) Global ice variables for output diagnostics                    | 
+      !--------------------------------------------------------------------
+      u_ice (:,:)     = 0._wp
+      v_ice (:,:)     = 0._wp
+      stress1_i(:,:)  = 0._wp
+      stress2_i(:,:)  = 0._wp
+      stress12_i(:,:) = 0._wp
+
+      !--------------------------------------------------------------------
+      ! 5) Moments for advection
+      !--------------------------------------------------------------------
+
+      sxopw (:,:) = 0._wp 
+      syopw (:,:) = 0._wp 
+      sxxopw(:,:) = 0._wp 
+      syyopw(:,:) = 0._wp 
+      sxyopw(:,:) = 0._wp
+
+      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
+      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
+      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
+      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
+      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
+
+      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
+      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
+      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
+      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
+      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
+
+      sxsal  (:,:,:)  = 0._wp
+      sysal  (:,:,:)  = 0._wp
+      sxxsal (:,:,:)  = 0._wp
+      syysal (:,:,:)  = 0._wp
+      sxysal (:,:,:)  = 0._wp
+
+      sxage  (:,:,:)  = 0._wp
+      syage  (:,:,:)  = 0._wp
+      sxxage (:,:,:)  = 0._wp
+      syyage (:,:,:)  = 0._wp
+      sxyage (:,:,:)  = 0._wp
+
+      !------------------------------------
+      ! 6) store fields at before time-step
+      !------------------------------------
+      ! it is only necessary for the 1st interpolation by Agrif
+      a_i_b  (:,:,:)   = a_i  (:,:,:)
+      e_i_b  (:,:,:,:) = e_i  (:,:,:,:)
+      v_i_b  (:,:,:)   = v_i  (:,:,:)
+      v_s_b  (:,:,:)   = v_s  (:,:,:)
+      e_s_b  (:,:,:,:) = e_s  (:,:,:,:)
+      smv_i_b(:,:,:)   = smv_i(:,:,:)
+      oa_i_b (:,:,:)   = oa_i (:,:,:)
+      u_ice_b(:,:)     = u_ice(:,:)
+      v_ice_b(:,:)     = v_ice(:,:)
+
+      ! MV MP 2016
+      IF ( nn_pnd_scheme > 0 ) THEN
+         sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
+         syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
+         sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
+         syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
+         sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
+      ENDIF
+      ! END MV MP 2016
+
+!!!clem
+!!      ! Output the initial state and forcings
+!!      CALL dia_wri_state( 'output.init', nit000 )
+!!!      
+
+   END SUBROUTINE ice_istate
+
+   SUBROUTINE ice_istate_init
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_istate_init  ***
+      !!        
+      !! ** Purpose : Definition of initial state of the ice 
+      !!
+      !! ** Method : Read the namice_ini namelist and check the parameter 
+      !!       values called at the first timestep (nit000)
+      !!
+      !! ** input : 
+      !!        Namelist namice_ini
+      !!
+      !! history :
+      !!  8.5  ! 03-08 (C. Ethe) original code 
+      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
+      !!-----------------------------------------------------------------------------
+      !
+      INTEGER :: ios,ierr,inum_ice                 ! Local integer output status for namelist read
+      INTEGER :: ji,jj
+      INTEGER :: ifpr, ierror
+      !
+      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files
+      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi
+      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read
+      !
+      NAMELIST/namice_ini/ ln_iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s,  &
+         &                rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, &
+         &                rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s,                             &
+         &                sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, cn_dir
+      !!-----------------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namice_ini in reference namelist : Ice initial state
+      READ  ( numnam_ice_ref, namice_ini, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_ini in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namice_ini in configuration namelist : Ice initial state
+      READ  ( numnam_ice_cfg, namice_ini, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_ini in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namice_ini )
+
+      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
+      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
+      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
+
+      ! Define the initial parameters
+      ! -------------------------
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namice_ini'
+         WRITE(numout,*) '      initialization with ice (T) or not (F)       ln_iceini     = ', ln_iceini
+         WRITE(numout,*) '      ice initialization from a netcdf file      ln_iceini_file  = ', ln_iceini_file
+         WRITE(numout,*) '      threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst
+         WRITE(numout,*) '      initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n
+         WRITE(numout,*) '      initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s 
+         WRITE(numout,*) '      initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n
+         WRITE(numout,*) '      initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s
+         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n
+         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s
+         WRITE(numout,*) '      initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n
+         WRITE(numout,*) '      initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s
+         WRITE(numout,*) '      initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n
+         WRITE(numout,*) '      initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s
+      ENDIF
+
+      IF( ln_iceini_file ) THEN                      ! Ice initialization using input file
+         !
+         ! set si structure
+         ALLOCATE( si(jpfldi), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'Ice_ini in iceistate: unable to allocate si structure' )   ;   RETURN
+         ENDIF
+
+         DO ifpr = 1, jpfldi
+            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
+            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
+         END DO
+
+         ! fill si with slf_i and control print
+         CALL fld_fill( si, slf_i, cn_dir, 'ice_istate', 'ice istate ini', 'numnam_ice' )
+
+         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
+
+      ENDIF
+
+   END SUBROUTINE ice_istate_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :         Empty module          NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE iceistate
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8514)
@@ -81,5 +81,5 @@
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
 
-      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !-----------------------------------------------------------------------------------------------
@@ -291,5 +291,5 @@
          CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
          CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)   )
          
          DO ji = 1, nidx
@@ -307,9 +307,9 @@
          CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1)   )
          CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d  (1:nidx), a_i(:,:,1)    )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)    )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1)   )
          !
       ENDIF
       !
-      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
    END SUBROUTINE ice_itd_rem
@@ -664,5 +664,5 @@
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_itd_init : Initialization of ice cat distribution '
+         WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_itd : '
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8514)
@@ -141,5 +141,5 @@
       ENDIF      
       !                    ! conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !-----------------------------------------------------------------------------!
@@ -288,5 +288,5 @@
       !-----------------------------------------------------------------------------!
       !                    ! conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !                    ! control prints
@@ -774,5 +774,5 @@
 
       !                              !--------------------------------------------------!
-      IF( ln_str_Rot ) THEN          ! Ice strength => Rothrock (1975) method           !
+      IF( ln_str_R75 ) THEN          ! Ice strength => Rothrock (1975) method           !
       !                              !--------------------------------------------------!
          z1_3 = 1._wp / 3._wp
@@ -793,5 +793,5 @@
          ismooth = 1
       !                              !--------------------------------------------------!
-      ELSEIF( ln_str_Hib ) THEN      ! Ice strength => Hibler (1979) method             !
+      ELSEIF( ln_str_H79 ) THEN      ! Ice strength => Hibler (1979) method             !
       !                              !--------------------------------------------------!
          strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
@@ -865,6 +865,6 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_rdgrft/ ln_str_Hib, rn_pstar, rn_crhg, &
-         &                    ln_str_Rot, rn_perdg,          &
+      NAMELIST/namice_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
+         &                    ln_str_R75, rn_perdg,          &
          &                    rn_cs     ,                    &
          &                    ln_partf_lin, rn_gstar,        &
@@ -885,11 +885,11 @@
       IF (lwp) THEN                          ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for ridging/rafting '
+         WRITE(numout,*) 'ice_rdgrft_init: ice parameters for ridging/rafting '
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_rdgrft'
-         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_Hib   = ', ln_str_Hib 
+         WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
          WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
          WRITE(numout,*) '            2nd bulk-rhelogy parameter                         rn_crhg      = ', rn_crhg
-         WRITE(numout,*) '      ice strength parameterization Rothrock (1975)            ln_str_Rot   = ', ln_str_Rot 
+         WRITE(numout,*) '      ice strength parameterization Rothrock (1975)            ln_str_R75   = ', ln_str_R75 
          WRITE(numout,*) '            Ratio of ridging work to PotEner change in ridging rn_perdg     = ', rn_perdg 
          WRITE(numout,*) '      Fraction of shear energy contributing to ridging         rn_cs        = ', rn_cs 
@@ -910,6 +910,6 @@
       ENDIF
       !
-      IF ( ( ln_str_Hib .AND. ln_str_Rot ) .OR. ( .NOT.ln_str_Hib .AND. .NOT.ln_str_Rot ) ) THEN
-         CALL ctl_stop( 'ice_rdgrft_init: choose one and only one formulation for ice strength (ln_str_Hib or ln_str_Rot)' )
+      IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN
+         CALL ctl_stop( 'ice_rdgrft_init: choose one and only one formulation for ice strength (ln_str_H79 or ln_str_R75)' )
       ENDIF
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8514)
@@ -70,5 +70,5 @@
       !
       !                             ! -- conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !                      
       IF( ln_landfast ) THEN        ! -- Landfast ice parameterization: define max bottom friction
@@ -96,5 +96,5 @@
       !
       !                                                   !- conservation test
-      IF( ln_limdiachk   )   CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       IF( ln_ctl         )   CALL ice_prt3D  ('icerhg')   !- Control prints
       IF( nn_timing == 1 )   CALL timing_stop('icerhg')   !- timing
@@ -132,5 +132,5 @@
       IF(lwp) THEN                     ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_rhg_init : ice parameters for ice dynamics '
+         WRITE(numout,*) 'ice_rhg_init: ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_rhg'
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8514)
@@ -90,5 +90,5 @@
       ENDIF
       !
-      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - Beginning the time step - ' )   ! control print
+      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - Beginning the time step - ' )   ! control print
    END SUBROUTINE ice_rst_opn
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8514)
@@ -49,5 +49,5 @@
    USE limmp          ! sea-ice: melt ponds
    ! END MV MP 2016
-   USE iceist         ! sea-ice: initial state
+   USE iceistate      ! sea-ice: initial state
    USE icethd_sal     ! sea-ice: thermodynamics and salinity
    USE iceitd         ! sea-ice: remapping thickness distribution
@@ -163,5 +163,5 @@
          IF( .NOT. Agrif_Root()     )   CALL agrif_interp_lim3('T')   ! -- AGRIF 
 #endif
-         IF( ln_limthd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
+         IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
          !
          !
@@ -191,5 +191,5 @@
          ! --- ice thermodynamics --- !
          !----------------------------!
-         IF( ln_limthd )            CALL ice_thd( kt )          ! -- Ice thermodynamics      
+         IF( ln_icethd )            CALL ice_thd( kt )          ! -- Ice thermodynamics      
 
          ! MV MP 2016
@@ -197,5 +197,5 @@
          ! END MV MP 2016
 
-         IF( ln_limthd )            CALL ice_cor( kt , 2 )      ! -- Corrections
+         IF( ln_icethd )            CALL ice_cor( kt , 2 )      ! -- Corrections
          ! ---
 # if defined key_agrif
@@ -214,5 +214,5 @@
 !!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
 !!# endif
-         IF( ln_limdiahsb )         CALL ice_dia( kt )        ! -- Diagnostics and outputs 
+         IF( ln_icediahsb )         CALL ice_dia( kt )        ! -- Diagnostics and outputs 
          !
                                     CALL ice_wri( 1 )         ! -- Ice outputs 
@@ -223,5 +223,5 @@
          IF( lrst_ice )             CALL ice_rst_write( kt )  ! -- Ice restart file 
          !
-         IF( ln_limctl )            CALL ice_ctl( kt )        ! alerts in case of model crash
+         IF( ln_icectl )            CALL ice_ctl( kt )        ! alerts in case of model crash
          !
       ENDIF   ! End sea-ice time step only
@@ -249,6 +249,6 @@
       !!----------------------------------------------------------------------
       IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'ice_init : update ocean surface boudary condition' 
-      IF(lwp) WRITE(numout,*) '~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping'
+      IF(lwp) WRITE(numout,*) 'ice_init: Arrays allocation & Initialization off all routines & init state' 
+      IF(lwp) WRITE(numout,*) '~~~~~~~~'
       !
       !                                ! Open the reference and configuration namelist files and namelist output file
@@ -285,5 +285,5 @@
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
-         CALL ice_ist
+         CALL ice_istate
       ELSE                                    ! start from a restart file
          CALL ice_rst_read
@@ -360,5 +360,5 @@
       ENDIF
       !
-      IF( ln_bdy .AND. ln_limdiachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY')
+      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY')
       !
       rdt_ice   = REAL(nn_fsbc) * rdt           !--- sea-ice timestep and inverse
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8514)
@@ -97,5 +97,5 @@
       
       ! conservation test
-      IF( ln_limdiachk )   CALL ice_cons_hsm( 0, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
+      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
 
       CALL ice_var_glo2eqv
@@ -179,7 +179,7 @@
       
       ! In case we bypass open-water ice formation
-      IF( .NOT. ln_limdO )  qlead(:,:) = 0._wp
+      IF( .NOT. ln_icedO )  qlead(:,:) = 0._wp
       ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
-      IF( .NOT. ln_limdH ) THEN
+      IF( .NOT. ln_icedH ) THEN
          hfx_in(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
          fhtur (:,:) = 0._wp
@@ -225,5 +225,5 @@
             dh_snowice(1:nidx) = 0._wp ; dh_i_sub (1:nidx) = 0._wp
             !
-            IF( ln_limdH ) THEN                                     ! --- growing/melting --- !
+            IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
                               CALL ice_thd_dif                             ! Ice/Snow Temperature profile
                               CALL ice_thd_dh                              ! Ice/Snow thickness   
@@ -236,8 +236,8 @@
             !
 !!gm please create a new logical (l_thd_lam or a better explicit name) set one for all in icestp.F90 module
-!!gm        l_thd_lam = ln_limdH .AND. ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 )
+!!gm        l_thd_lam = ln_icedH .AND. ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 )
 !!gm        by the way, the different options associated with nn_monocat =1 to 4  are quite impossible to identify
 !!gm        more comment to add when ready the namelist, with an explicit print in the ocean.output
-            IF( ln_limdH ) THEN
+            IF( ln_icedH ) THEN
                IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN
                               CALL ice_thd_lam                      ! --- extra lateral melting if monocat --- !
@@ -245,5 +245,5 @@
             END IF
             !
-            IF( ln_limdA )    CALL ice_thd_da                       ! --- lateral melting --- !
+            IF( ln_icedA )    CALL ice_thd_da                       ! --- lateral melting --- !
             !
                               CALL ice_thd_1d2d( jl, 2 )            ! --- Change units of e_i, e_s from J/m3 to J/m2 --- !
@@ -257,5 +257,5 @@
       oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
 
-      IF( ln_limdiachk )   CALL ice_cons_hsm( 1, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
+      IF( ln_icediachk )   CALL ice_cons_hsm( 1, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
       !
                            CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
@@ -264,7 +264,7 @@
       IF( jpl > 1 )        CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- !
       !
-      IF( ln_limdO )       CALL ice_thd_lac                ! --- frazil ice growing in leads --- !
-      !
-      IF( ln_limctl )      CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print
+      IF( ln_icedO )       CALL ice_thd_lac                ! --- frazil ice growing in leads --- !
+      !
+      IF( ln_icectl )      CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print
       IF( ln_ctl )         CALL ice_prt3D( 'icethd' )      ! Control print
       !
@@ -539,9 +539,9 @@
       INTEGER  ::   ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_thd/ ln_limthd, rn_kappa_i, nn_ice_thcon, ln_dqnsice, rn_cdsn,                       &
-         &                ln_limdH, rn_betas,                                                             &
-         &                ln_limdA, rn_beta, rn_dmin,                                                     &
-         &                ln_limdO, rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, rn_himin,   &
-         &                nn_limflx
+      NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,          &
+         &                ln_icedH, rn_blow_s,                                                             &
+         &                ln_icedA, rn_beta, rn_dmin,                                                     &
+         &                ln_icedO, rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz, rn_himin,   &
+         &                nn_iceflx
       !!-------------------------------------------------------------------
       !
@@ -560,47 +560,52 @@
          WRITE(numout,*) '~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_thd'
-         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                       ln_limthd    = ', ln_limthd
+         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                       ln_icethd    = ', ln_icethd
          WRITE(numout,*) '   -- icethd_dif --'
          WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-         WRITE(numout,*) '      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon
-         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqnsice   = ', ln_dqnsice
-         WRITE(numout,*) '      thermal conductivity of the snow                        rn_cdsn      = ', rn_cdsn
+         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64
+         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07
+         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
+         WRITE(numout,*) '      thermal conductivity of the snow                        rn_cnd_s     = ', rn_cnd_s
          WRITE(numout,*) '   -- icethd_dh --'
-         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_limdH     = ', ln_limdH
-         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas
+         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_icedH     = ', ln_icedH
+         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall          rn_blow_s    = ', rn_blow_s
          WRITE(numout,*) '   -- icethd_da --'
-         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_limdA     = ', ln_limdA
+         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_icedA     = ', ln_icedA
          WRITE(numout,*) '      Coef. beta for lateral melting param.                   rn_beta      = ', rn_beta
          WRITE(numout,*) '      Minimum floe diameter for lateral melting param.        rn_dmin      = ', rn_dmin
          WRITE(numout,*) '   -- icethd_lac --'
-         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_limdO     = ', ln_limdO
-         WRITE(numout,*) '      ice thickness for lateral accretion                     rn_hnewice   = ', rn_hnewice
+         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_icedO     = ', ln_icedO
+         WRITE(numout,*) '      ice thickness for lateral accretion                     rn_hinew     = ', rn_hinew
          WRITE(numout,*) '      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil
-         WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb  = ', rn_maxfrazb
-         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil rn_vfrazb    = ', rn_vfrazb
-         WRITE(numout,*) '      Squeezing coefficient for collection of frazil          rn_Cfrazb    = ', rn_Cfrazb
+         WRITE(numout,*) '      Maximum proportion of frazil ice collecting at bottom   rn_maxfraz   = ', rn_maxfraz
+         WRITE(numout,*) '      Threshold relative drift speed for collection of frazil rn_vfraz     = ', rn_vfraz
+         WRITE(numout,*) '      Squeezing coefficient for collection of frazil          rn_Cfraz     = ', rn_Cfraz
          WRITE(numout,*) '   -- iceitd --'
          WRITE(numout,*) '      minimum ice thickness                                   rn_himin     = ', rn_himin 
          WRITE(numout,*) '   -- icestp --'
-         WRITE(numout,*) '      Multicategory heat flux formulation                     nn_limflx    = ', nn_limflx
+         WRITE(numout,*) '      Multicategory heat flux formulation                     nn_iceflx    = ', nn_iceflx
       ENDIF
       !
-      IF ( rn_hnewice < rn_himin )   CALL ctl_stop( 'STOP', 'ice_thd_init : rn_hnewice should be >= rn_himin' )
+      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
+         CALL ctl_stop( 'ice_thd_init: choose one and only one formulation for thermal conductivity (ln_cndi_U64 or ln_cndi_P07)' )
+      ENDIF
+      !
+      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'STOP', 'ice_thd_init : rn_hinew should be >= rn_himin' )
       !
       IF(lwp) WRITE(numout,*)
-      SELECT CASE( nn_limflx )         ! LIM3 Multi-category heat flux formulation
+      SELECT CASE( nn_iceflx )         ! LIM3 Multi-category heat flux formulation
       CASE( -1  )
-         IF(lwp) WRITE(numout,*) '   LIM3: use per-category fluxes (nn_limflx = -1) '
-         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
+         IF(lwp) WRITE(numout,*) '   LIM3: use per-category fluxes (nn_iceflx = -1) '
+         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in coupled mode must be 0 or 2' )
       CASE(  0  )
-         IF(lwp) WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_limflx = 0) '
+         IF(lwp) WRITE(numout,*) '   LIM3: use average per-category fluxes (nn_iceflx = 0) '
       CASE(  1  )
-         IF(lwp) WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_limflx = 1) '
-         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_limflx for LIM3 in coupled mode must be 0 or 2' )
+         IF(lwp) WRITE(numout,*) '   LIM3: use average then redistribute per-category fluxes (nn_iceflx = 1) '
+         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in coupled mode must be 0 or 2' )
       CASE(  2  )
-         IF(lwp) WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_limflx = 2) '
-         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_limflx for LIM3 in forced mode cannot be 2' )
+         IF(lwp) WRITE(numout,*) '   LIM3: Redistribute a single flux over categories (nn_iceflx = 2) '
+         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for LIM3 in forced mode cannot be 2' )
       CASE DEFAULT
-         CALL ctl_stop( 'ice_thd_init: LIM3 option, nn_limflx, should be between -1 and 2' )
+         CALL ctl_stop( 'ice_thd_init: LIM3 option, nn_iceflx, should be between -1 and 2' )
       END SELECT
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8514)
@@ -20,4 +20,5 @@
    USE ice            ! LIM variables
    USE ice1D          ! LIM thermodynamics
+   USE icethd_sal     ! sea-ice: salinity profiles
    !
    USE in_out_manager ! I/O manager
@@ -107,8 +108,8 @@
       !!------------------------------------------------------------------
 
-      ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values)
+      ! Discriminate between time varying salinity and constant
       SELECT CASE( nn_icesal )                  ! varying salinity or not
-         CASE( 1, 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
-         CASE( 2 )      ;   zswitch_sal = 1._wp   ! varying salinity profile
+         CASE( 1 , 3 )   ;   zswitch_sal = 0._wp   ! prescribed salinity profile
+         CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
       END SELECT
 
@@ -406,6 +407,6 @@
 
       num_iter_max = 1
-      IF( nn_icesal == 2 ) num_iter_max = 5
-
+      IF( nn_icesal == 2 )   num_iter_max = 5  ! salinity varying in time
+      
       ! Iterative procedure
       DO ji = 1, nidx
@@ -586,5 +587,5 @@
          hfx_out_1d(ji)  = hfx_out_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice
 
-         IF( ln_limctl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
+         IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
       END DO
 
@@ -613,5 +614,5 @@
          sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_rdtice 
 
-         ! virtual salt flux to keep salinity constant
+         ! Case constant salinity in time: virtual salt flux to keep salinity constant
          IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
             sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean
@@ -665,5 +666,5 @@
       REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
       REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
-      pout = ( 1._wp - ( pin )**rn_betas )
+      pout = ( 1._wp - ( pin )**rn_blow_s )
    END SUBROUTINE ice_thd_snwblow_2d
 
@@ -671,5 +672,5 @@
       REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
       REAL(wp), DIMENSION(:), INTENT(inout) :: pout
-      pout = ( 1._wp - ( pin )**rn_betas )
+      pout = ( 1._wp - ( pin )**rn_blow_s )
    END SUBROUTINE ice_thd_snwblow_1d
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8514)
@@ -287,5 +287,5 @@
          !------------------------------------------------------------------------------|
          !
-         IF( nn_ice_thcon == 0 ) THEN      ! Untersteiner (1964) formula
+         IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
             DO ji = 1 , nidx
                ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
@@ -299,7 +299,6 @@
                END DO
             END DO
-         ENDIF
-
-         IF( nn_ice_thcon == 1 ) THEN      ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T
+
+         ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
             DO ji = 1 , nidx
                ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
@@ -347,6 +346,6 @@
 
                ! Effective thickness he (zhe)
-               zfac     = 1._wp / ( rn_cdsn + zkimean )
-               zratio_s = rn_cdsn   * zfac
+               zfac     = 1._wp / ( rn_cnd_s + zkimean )
+               zratio_s = rn_cnd_s   * zfac
                zratio_i = zkimean * zfac
                zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
@@ -371,11 +370,11 @@
          DO ji = 1, nidx
             zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
-            zkappa_s(ji,0)        = zghe(ji) * rn_cdsn * zfac
-            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cdsn * zfac
+            zkappa_s(ji,0)        = zghe(ji) * rn_cnd_s * zfac
+            zkappa_s(ji,nlay_s)   = zghe(ji) * rn_cnd_s * zfac
          END DO
 
          DO jk = 1, nlay_s-1
             DO ji = 1 , nidx
-               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cdsn / MAX( epsi10, 2.0 * zh_s(ji) )
+               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rn_cnd_s / MAX( epsi10, 2.0 * zh_s(ji) )
             END DO
          END DO
@@ -393,6 +392,6 @@
             zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
             zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
-            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cdsn * ztcond_i(ji,0) / & 
-           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cdsn * zh_i(ji) ) )
+            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rn_cnd_s * ztcond_i(ji,0) / & 
+           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rn_cnd_s * zh_i(ji) ) )
             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
          END DO
@@ -423,5 +422,5 @@
          !------------------------------------------------------------------------------|
          !
-         IF ( ln_dqnsice ) THEN 
+         IF ( ln_dqns_i ) THEN 
             DO ji = 1 , nidx
                ! update of the non solar flux according to the update in T_su
@@ -708,7 +707,7 @@
       !--- Snow-ice interfacial temperature (diagnostic SIMIP)
       DO ji = 1, nidx
-         zfac        = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
+         zfac        = 1. / MAX( epsi10 , rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) )
          IF( zh_s(ji) >= 1.e-3 ) THEN
-            t_si_1d(ji) = ( rn_cdsn        * zh_i(ji) * t_s_1d(ji,1) + &
+            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + &
                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac
          ELSE
@@ -718,5 +717,5 @@
       ! END MV SIMIP 2016
 
-      IF( ln_limctl .AND. lwp ) THEN
+      IF( ln_icectl .AND. lwp ) THEN
          WRITE(numout,*) ' zdti_max : ', zdti_max
          WRITE(numout,*) ' iconv    : ', iconv
@@ -748,5 +747,5 @@
 
       ! --- diagnose the change in non-solar flux due to surface temperature change --- !
-      IF ( ln_dqnsice ) THEN
+      IF ( ln_dqns_i ) THEN
          DO ji = 1, nidx
             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8514)
@@ -26,4 +26,5 @@
    USE icethd_ent     ! sea-ice: thermodynamics, enthalpy
    USE icevar         ! sea-ice: operations
+   USE icethd_sal     ! sea-ice: salinity profiles
    !
    USE in_out_manager ! I/O manager
@@ -112,5 +113,5 @@
       !!-----------------------------------------------------------------------!
 
-      IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       CALL ice_var_agg(1)
@@ -133,5 +134,5 @@
 
       ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp )   ;   hicol(:,:) = rn_hnewice
+      WHERE( qlead(:,:) < 0._wp )   ;   hicol(:,:) = rn_hinew
       ELSEWHERE                     ;   hicol(:,:) = 0._wp
       END WHERE
@@ -348,5 +349,5 @@
             DO ji = 1, nidx
                rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
-               zfrazb        = rswitch * ( TANH( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb
+               zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
                zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
                zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
@@ -478,5 +479,5 @@
       ENDIF ! nidx > 0
       !
-      IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
    END SUBROUTINE ice_thd_lac
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8514)
@@ -55,7 +55,7 @@
       SELECT CASE ( nn_icesal )
       !
-      !              !---------------------------------------------!
-      CASE( 2 )      !  time varying salinity with linear profile  !
-      !              !---------------------------------------------!
+      !               !---------------------------------------------!
+      CASE( 2 )       !  time varying salinity with linear profile  !
+      !               !---------------------------------------------!
          DO ji = 1, nidx
 
@@ -71,5 +71,5 @@
             sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si
 
-            IF( ln_limdS ) THEN
+            IF( ln_icedS ) THEN
                !---------------------------------------------------------
                !  Update ice salinity from brine drainage and flushing
@@ -91,7 +91,7 @@
          CALL ice_var_salprof1d
          !
-      !                 !---------------------------------------------!
-      CASE( 3 )         ! constant salinity with a fix profile        ! (Schwarzacher (1959) multiyear salinity profile(mean = 2.30)
-      !                 !---------------------------------------------!
+      !         !---------------------------------------------!
+      CASE( 3 ) ! constant salinity with a fixed profile      ! (Schwarzacher (1959) multiyear salinity profile(mean = 2.30)
+      !         !---------------------------------------------!
          CALL ice_var_salprof1d
       !
@@ -114,5 +114,5 @@
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_sal/ ln_limdS , nn_icesal , rn_icesal, rn_sal_gd, rn_time_gd,   &
+      NAMELIST/namice_sal/ ln_icedS , nn_icesal , rn_icesal, rn_sal_gd, rn_time_gd,   &
          &                rn_sal_fl, rn_time_fl, rn_simax , rn_simin 
       !!-------------------------------------------------------------------
@@ -132,5 +132,5 @@
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_sal'
-         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_limdS   = ', ln_limdS
+         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_icedS   = ', ln_icedS
          WRITE(numout,*) '      switch for salinity                                     nn_icesal  = ', nn_icesal
          WRITE(numout,*) '      bulk salinity value if nn_icesal = 1                    rn_icesal  = ', rn_icesal
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8514)
@@ -121,5 +121,5 @@
 
       ! --- case we bypass ice thermodynamics --- !
-      IF( .NOT. ln_limthd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
+      IF( .NOT. ln_icethd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
          hfx_in   (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
          hfx_out  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
@@ -208,7 +208,7 @@
 
       !                    ! conservation test
-      IF( ln_limdiachk .AND. .NOT. ln_bdy)   CALL ice_cons_final( 'iceupdate' )
+      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final( 'iceupdate' )
       !                    ! control prints
-      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' )
+      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' )
       IF( ln_ctl    )   CALL ice_prt3D( 'iceupdate' )
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8514)
@@ -191,5 +191,5 @@
       o_i(:,:,:)  = oa_i(:,:,:) * z1_a_i(:,:,:)    !--- ice age
 
-      IF( nn_icesal == 2 )THEN                     !--- salinity (with a minimum value imposed everywhere)
+      IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere)
          WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )
          ELSEWHERE                      ;   sm_i(:,:,:) = rn_simin
@@ -289,20 +289,17 @@
       !!------------------------------------------------------------------
 
-!!gm  much much more secure to defined when reading nn_icesal in the namelist integers =1, 2, 3  with explicit names
-!!       for example np_Scst_noProfile = 1  ;  np_Svar_linProfile = 2   ;   np_Scst_fixProfile
-
 !!gm Question: Remove the option 3 ?  How many years since it last use ? 
 
       SELECT CASE ( nn_icesal )
       !
-      !              !---------------------------------------!
-      CASE( 1 )      !  constant salinity in time and space  !
-         !           !---------------------------------------!
+      !               !---------------------------------------!
+      CASE( 1 )       !  constant salinity in time and space  !
+         !            !---------------------------------------!
          s_i (:,:,:,:) = rn_icesal
          sm_i(:,:,:)   = rn_icesal
          !
-         !           !---------------------------------------------!
-      CASE( 2 )      !  time varying salinity with linear profile  !
-         !           !---------------------------------------------!
+         !            !---------------------------------------------!
+      CASE( 2 )       !  time varying salinity with linear profile  !
+         !            !---------------------------------------------!
          !
          ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
@@ -343,7 +340,7 @@
          DEALLOCATE( z_slope_s , zalpha )
          !
-         !           !-------------------------------------------!
-      CASE( 3 )      ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
-         !           !-------------------------------------------!                                   (mean = 2.30)
+         !            !-------------------------------------------!
+      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
+         !            !-------------------------------------------!                                   (mean = 2.30)
          !
          sm_i(:,:,:) = 2.30_wp
@@ -415,12 +412,12 @@
       SELECT CASE ( nn_icesal )
       !
-      !              !---------------------------------------!
-      CASE( 1 )      !  constant salinity in time and space  !
-         !           !---------------------------------------!
+      !               !---------------------------------------!
+      CASE( 1 )       !  constant salinity in time and space  !
+         !            !---------------------------------------!
          s_i_1d(:,:) = rn_icesal
          !
-         !           !---------------------------------------------!
-      CASE( 2 )      !  time varying salinity with linear profile  !
-         !           !---------------------------------------------!
+         !            !---------------------------------------------!
+      CASE( 2 )       !  time varying salinity with linear profile  !
+         !            !---------------------------------------------!
          !
          ALLOCATE( z_slope_s(jpij) )
@@ -445,7 +442,7 @@
          DEALLOCATE( z_slope_s )
 
-         !           !-------------------------------------------!
-      CASE( 3 )      ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
-         !           !-------------------------------------------!                                   (mean = 2.30)
+         !            !-------------------------------------------!
+      CASE( 3 )       ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
+         !            !-------------------------------------------!                                   (mean = 2.30)
          !
          sm_i_1d(:) = 2.30_wp
@@ -572,5 +569,5 @@
       !!                  ice thickness distribution follows a gaussian law
       !!               around the concentration of the most likely ice thickness
-      !!                           (similar as iceist.F90)
+      !!                           (similar as iceistate.F90)
       !!
       !! ** Method:   Iterative procedure
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8514)
@@ -73,5 +73,5 @@
                         CALL ice_var_zapsmall
                         CALL ice_var_agg(1)
-      IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
+      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
       !
    END SUBROUTINE bdy_ice
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/SETTE/sette.sh
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/SETTE/sette.sh	(revision 8513)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/SETTE/sette.sh	(revision 8514)
@@ -143,5 +143,5 @@
 #                     15 & 16 
 
-for config in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
+for config in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 
 
 do
@@ -312,5 +312,6 @@
     set_namelist namelist_cfg jpnj 8
     set_namelist namelist_cfg jpnij 32
-    set_namelist namelist_ice_cfg ln_limdiachk .true.
+    set_namelist namelist_ice_cfg ln_icediachk .true.
+##    set_namelist namelist_ice_cfg ln_icediahsb .true.
     set_namelist namelist_top_cfg ln_trcdta .false.
     # put ln_ironsed, ln_river, ln_ndepo, ln_dust to false
@@ -351,5 +352,6 @@
     set_namelist namelist_cfg jpnj 8
     set_namelist namelist_cfg jpnij 32
-    set_namelist namelist_ice_cfg ln_limdiachk .true.
+    set_namelist namelist_ice_cfg ln_icediachk .true.
+##    set_namelist namelist_ice_cfg ln_icediahsb .true.
     set_namelist namelist_top_cfg ln_rsttr .true.
     set_namelist namelist_top_cfg nn_rsttr 2
@@ -413,5 +415,6 @@
     set_namelist namelist_cfg jpnj 8
     set_namelist namelist_cfg jpnij 32
-    set_namelist namelist_ice_cfg ln_limdiachk .true.
+    set_namelist namelist_ice_cfg ln_icediachk .true.
+    set_namelist namelist_ice_cfg ln_icediahsb .true.
     set_namelist namelist_top_cfg ln_trcdta .false.
     # put ln_ironsed, ln_river, ln_ndepo, ln_dust to false
@@ -455,5 +458,6 @@
     set_namelist namelist_cfg jpnj 4
     set_namelist namelist_cfg jpnij 32
-    set_namelist namelist_ice_cfg ln_limdiachk .true.
+    set_namelist namelist_ice_cfg ln_icediachk .true.
+    set_namelist namelist_ice_cfg ln_icediahsb .true.
     set_namelist namelist_top_cfg ln_trcdta .false.
     # put ln_ironsed, ln_river, ln_ndepo, ln_dust to false
@@ -814,5 +818,5 @@
     set_namelist namelist_cfg jpnj 8
     set_namelist namelist_cfg jpnij 32
-    set_namelist namelist_ice_cfg ln_limdiachk .true.
+    set_namelist namelist_ice_cfg ln_icediachk .true.
     if [ ${USING_MPMD} == "yes" ] ; then
        set_xio_using_server iodef.xml true
@@ -840,5 +844,5 @@
     set_namelist namelist_cfg ln_rstart .true.
     set_namelist namelist_cfg nn_rstctl 2
-    set_namelist namelist_ice_cfg ln_limdiachk .true.
+    set_namelist namelist_ice_cfg ln_icediachk .true.
     set_namelist namelist_ice_cfg cn_icerst_in \"SAS_00000050_restart_ice\"
     if [ ${USING_MPMD} == "yes" ] ; then
