Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_ice_cfg
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_ice_cfg	(revision 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_ice_cfg	(revision 8531)
@@ -1,64 +1,84 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !! ESIM configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref
-!!              1 - Generic parameters                 (namice_run)
-!!              2 - Ice thickness discretization       (namice_itd)
-!!              3 - Ice dynamics                       (namice_dyn)
-!!              4 - Ice ridging/rafting                (namice_rdgrft)
-!!              5 - Ice rheology                       (namice_rhg)
-!!              6 - Ice advection                      (namice_adv)
-!!              7 - Ice thermodynamics                 (namice_thd)
-!!              8 - Ice salinity                       (namice_sal)
-!!              9 - Ice melt ponds                     (namice_mp)
-!!             10 - Ice initialization                 (namice_ini)
-!!             11 - Ice/snow albedos                   (namice_alb)
-!!             12 - Ice diagnostics                    (namice_dia)
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface forcing                (namforcing)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (nammp)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !
 !------------------------------------------------------------------------------
-&namice_run     !   Generic parameters
+&nampar     !   Generic parameters
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_itd     !   Ice discretization
+&namitd     !   Ice discretization
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_dyn     !   Ice dynamics
+&namdyn     !   Ice dynamics
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_rdgrft  !   Ice ridging/rafting
+&namdyn_rdgrft  !   Ice ridging/rafting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_rhg     !   Ice rheology
+&namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_adv     !   Ice advection
+&namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_thd     !   Ice thermodynamics
+&namforcing     !   Ice surface forcing
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_sal     !   Ice salinity
+&namthd     !   Ice thermodynamics
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicemp      !   Melt ponds
+&namthd_zdf     !   Ice heat diffusion
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_ini     !   Ice initialization
+&namthd_da     !   Ice lateral melting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_alb     !   albedo parameters
+&namthd_do     !   Ice growth in open water
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namice_dia     !   Diagnostics
+&namthd_sal     !   Ice salinity
 !------------------------------------------------------------------------------
 /
+!------------------------------------------------------------------------------
+&nammp      !   Melt ponds
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namini     !   Ice initialization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namalb     !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namdia     !   Diagnostics
+!------------------------------------------------------------------------------
+/
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8531)
@@ -1,20 +1,24 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !! ESIM namelist:  
-!!              1 - Generic parameters                 (namice_run)
-!!              2 - Ice thickness discretization       (namice_itd)
-!!              3 - Ice dynamics                       (namice_dyn)
-!!              4 - Ice ridging/rafting                (namice_rdgrft)
-!!              5 - Ice rheology                       (namice_rhg)
-!!              6 - Ice advection                      (namice_adv)
-!!              7 - Ice thermodynamics                 (namice_thd)
-!!              8 - Ice salinity                       (namice_sal)
-!!              9 - Ice melt ponds                     (namice_mp)
-!!             10 - Ice initialization                 (namice_ini)
-!!             11 - Ice/snow albedos                   (namice_alb)
-!!             12 - Ice diagnostics                    (namice_dia)
+!!              1 - Generic parameters                 (nampar)
+!!              2 - Ice thickness discretization       (namitd)
+!!              3 - Ice dynamics                       (namdyn)
+!!              4 - Ice ridging/rafting                (namdyn_rdgrft)
+!!              5 - Ice rheology                       (namdyn_rhg)
+!!              6 - Ice advection                      (namdyn_adv)
+!!              7 - Ice surface forcing                (namforcing)
+!!              8 - Ice thermodynamics                 (namthd)
+!!              9 - Ice heat diffusion                 (namthd_zdf)
+!!             10 - Ice lateral melting                (namthd_da)
+!!             11 - Ice growth in open water           (namthd_do)
+!!             12 - Ice salinity                       (namthd_sal)
+!!             13 - Ice melt ponds                     (nammp)
+!!             14 - Ice initialization                 (namini)
+!!             15 - Ice/snow albedos                   (namalb)
+!!             16 - Ice diagnostics                    (namdia)
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 !
 !------------------------------------------------------------------------------
-&namice_run     !   Generic parameters
+&nampar     !   Generic parameters
 !------------------------------------------------------------------------------
    jpl              =   5             !  number of ice  categories
@@ -35,5 +39,5 @@
 /
 !------------------------------------------------------------------------------
-&namice_itd     !   Ice discretization
+&namitd     !   Ice discretization
 !------------------------------------------------------------------------------
    rn_himean        =   2.0           !  expected domain-average ice thickness (m)
@@ -41,5 +45,5 @@
 /
 !------------------------------------------------------------------------------
-&namice_dyn     !   Ice dynamics
+&namdyn     !   Ice dynamics
 !------------------------------------------------------------------------------
    ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
@@ -49,5 +53,4 @@
       rn_vice       =   0.            !        prescribed ice v-velocity
    rn_ishlat        =   2.            !  free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
-   rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-)
    ln_landfast      = .false.         !  landfast ice parameterization (T or F)                           
       rn_gamma      =   0.15          !     fraction of ocean depth that ice must reach to initiate landfast
@@ -59,5 +62,5 @@
 /
 !------------------------------------------------------------------------------
-&namice_rdgrft  !   Ice ridging/rafting
+&namdyn_rdgrft  !   Ice ridging/rafting
 !------------------------------------------------------------------------------
           ! -- ice_rdgrft_strength -- !
@@ -86,5 +89,5 @@
 /
 !------------------------------------------------------------------------------
-&namice_rhg     !   Ice rheology
+&namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
    ln_rhg_EVP       = .true.          !  EVP rheology
@@ -96,5 +99,5 @@
 /
 !------------------------------------------------------------------------------
-&namice_adv     !   Ice advection
+&namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
    ln_adv_Pra       = .false.         !  Advection scheme (Prather)
@@ -103,89 +106,99 @@
 /
 !------------------------------------------------------------------------------
-&namice_thd     !   Ice thermodynamics
-!------------------------------------------------------------------------------
-                     ! -- icethd_dif -- !
-   rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice [1/m]
-   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)
-                      ! -- icethd_dh -- !
-   ln_icedH         = .true.            !  activate ice thickness change from growing/melting (T) or not (F)
-   rn_blow_s        =   0.66            !     mesure of snow blowing into the leads
-                                        !     = 1 => no snow blowing, < 1 => some snow blowing
-                      ! -- icethd_da -- !
-   ln_icedA         = .true.            !  activate lateral melting param. (T) or not (F)
-      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
-                     ! -- icethd_lac -- !
-   ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F)
-      rn_hinew      =   0.1             !     thickness for new ice formation in open water (m), must be larger than rn_hnewice
-      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
-                         ! -- 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
-!------------------------------------------------------------------------------
-                    ! -- icethd_sal -- !
-   ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F)
-   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)
-/
-!------------------------------------------------------------------------------
-&namice_ini     !   Ice initialization
-!------------------------------------------------------------------------------
-                     ! -- iceistate -- !
-   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              !  max delta temp. above Tfreeze with initial ice = (sst - tfreeze)
-   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
+&namforcing     !   Ice surface forcing
+!------------------------------------------------------------------------------
+   rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-)
+   rn_blow_s        =   0.66          !  mesure of snow blowing into the leads
+                                      !     = 1 => no snow blowing, < 1 => some snow blowing
+   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
+/
+!------------------------------------------------------------------------------
+&namthd     !   Ice thermodynamics
+!------------------------------------------------------------------------------
+   ln_icedH         = .true.          !  activate ice thickness change from growing/melting (T) or not (F)
+   ln_icedA         = .true.          !  activate lateral melting param. (T) or not (F)
+   ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
+   ln_icedS         = .true.          !  activate gravity drainage and flushing (T) or not (F)
+/
+!------------------------------------------------------------------------------
+&namthd_zdf     !   Ice heat diffusion
+!------------------------------------------------------------------------------
+   ln_zdf_Beer      = .true.          !  Heat diffusion follows a Beer law
+   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)
+   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)
+   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m]
+   ln_dqns_i        = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F)
+/
+!------------------------------------------------------------------------------
+&namthd_da     !   Ice lateral melting
+!------------------------------------------------------------------------------
+   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
+/
+!------------------------------------------------------------------------------
+&namthd_do     !   Ice growth in open water
+!------------------------------------------------------------------------------
+   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice
+   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
+/
+!------------------------------------------------------------------------------
+&namthd_sal     !   Ice salinity
+!------------------------------------------------------------------------------
+   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)
+/
+!------------------------------------------------------------------------------
+&nammp      !   Melt ponds
+!------------------------------------------------------------------------------
+   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)
+/
+!------------------------------------------------------------------------------
+&namini     !   Ice initialization
+!------------------------------------------------------------------------------
+   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           !  max delta temp. above Tfreeze with initial ice = (sst - tfreeze)
+   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'  , '' , '', ''
@@ -198,24 +211,24 @@
 /
 !------------------------------------------------------------------------------
-&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
-                                       !     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_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
-/
+&namalb     !   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
+                                      !     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 ( '' )
+/
+!------------------------------------------------------------------------------
+&namdia     !   Diagnostics
+!------------------------------------------------------------------------------
+   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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8531)
@@ -153,8 +153,9 @@
    !! * Share Module variables
    !!--------------------------------------------------------------------------
-   !                                     !!** ice-generic parameters namelist (namice_run) **
+   !                                     !!** ice-generic parameters namelist (nampar) **
    INTEGER           , PUBLIC ::   jpl             !: number of ice  categories 
    INTEGER           , PUBLIC ::   nlay_i          !: number of ice  layers 
    INTEGER           , PUBLIC ::   nlay_s          !: number of snow layers 
+   INTEGER           , PUBLIC ::   nn_monocat      !: virtual ITD mono-category parameterizations (1-4) or not (0)
    LOGICAL           , PUBLIC ::   ln_icedyn       !: flag for ice dynamics (T) or not (F)
    LOGICAL           , PUBLIC ::   ln_icethd       !: flag for ice thermo   (T) or not (F)
@@ -166,8 +167,8 @@
    CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory
 
-   !                                     !!** ice-itd namelist (namice_itd) **
+   !                                     !!** ice-itd namelist (namitd) **
    REAL(wp), PUBLIC ::   rn_himin         !: minimum ice thickness
    
-   !                                     !!** ice-dynamics namelist (namice_dyn) **
+   !                                     !!** ice-dynamics namelist (namdyn) **
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
@@ -177,5 +178,5 @@
    REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice) 
    !
-   !                                     !!** ice-rheology namelist (namice_rhg) **
+   !                                     !!** ice-rheology namelist (namrhg) **
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
    REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
@@ -183,26 +184,7 @@
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
    !
-   !                                     !!** ice-thermodynamics namelist (namice_thd) **
-                                          ! -- icethd_dif -- !
-   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   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_cnd_s         !: thermal conductivity of the snow [W/m/K]
+   !                                     !!** ice-thermodynamics namelist (namthd) **
                                           ! -- icethd_dh -- !
-   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_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_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_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
                                           ! -- icethd -- !
    INTEGER , PUBLIC ::   nn_iceflx        !: Redistribute heat flux over ice categories
@@ -213,5 +195,5 @@
    !                                      !   = 2  Redistribute a single flux over categories
 
-   !                                     !!** ice-salinity namelist (namice_sal) **
+   !                                     !!** ice-salinity namelist (namthd_sal) **
    INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model
    !                                      ! 1 - constant salinity in both space and time
@@ -223,5 +205,5 @@
 
    ! MV MP 2016
-   !                                     !!** melt pond namelist (namicemp)
+   !                                     !!** melt pond namelist (nammp)
    LOGICAL , PUBLIC ::   ln_pnd           !: activate ponds or not
    LOGICAL , PUBLIC ::   ln_pnd_rad       !: ponds radiatively active or not
@@ -231,5 +213,5 @@
    REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1), only if nn_pnd_scheme = 0
    ! END MV MP 2016
-   !                                     !!** ice-diagnostics namelist (namice_dia) **
+   !                                     !!** ice-diagnostics namelist (namdia) **
    LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)
    LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F)
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8531)
@@ -42,5 +42,5 @@
    INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme
    !
-   ! ** namelist (namice_adv) **
+   ! ** namelist (namdyn_adv) **
    LOGICAL ::   ln_adv_Pra   ! Prather        advection scheme
    LOGICAL ::   ln_adv_UMx   ! Ultimate-Macho advection scheme
@@ -132,22 +132,22 @@
       !!      dynamics
       !!
-      !! ** Method  :  Read the namice_adv namelist and check the ice-dynamic
+      !! ** Method  :  Read the namdyn_adv namelist and check the ice-dynamic
       !!       parameter values called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namice_adv
+      !! ** input   :   Namelist namdyn_adv
       !!-------------------------------------------------------------------
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx
+      NAMELIST/namdyn_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx
       !!-------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )         ! Namelist namice_adv in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namice_adv, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_adv in reference namelist', lwp )
+      REWIND( numnam_ice_ref )         ! Namelist namdyn_adv in reference namelist : Ice dynamics
+      READ  ( numnam_ice_ref, namdyn_adv, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_adv in reference namelist', lwp )
       !
-      REWIND( numnam_ice_cfg )         ! Namelist namice_adv in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namice_adv, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_adv in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_adv )
+      REWIND( numnam_ice_cfg )         ! Namelist namdyn_adv in configuration namelist : Ice dynamics
+      READ  ( numnam_ice_cfg, namdyn_adv, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_adv in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namdyn_adv )
       !
       IF(lwp) THEN                     ! control print
@@ -155,8 +155,8 @@
          WRITE(numout,*) 'ice_adv_init: ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_adv'
-         WRITE(numout,*) '      type of advection scheme (Prather)                     ln_adv_Pra = ', ln_adv_Pra 
-         WRITE(numout,*) '      type of advection scheme (Ulimate-Macho)               ln_adv_UMx = ', ln_adv_UMx 
-         WRITE(numout,*) '         order of the Ultimate-Macho scheme                      nn_UMx = ', nn_UMx
+         WRITE(numout,*) '   Namelist namdyn_adv:'
+         WRITE(numout,*) '      type of advection scheme (Prather)             ln_adv_Pra = ', ln_adv_Pra 
+         WRITE(numout,*) '      type of advection scheme (Ulimate-Macho)       ln_adv_UMx = ', ln_adv_UMx 
+         WRITE(numout,*) '         order of the Ultimate-Macho scheme          nn_UMx     = ', nn_UMx
       ENDIF
       !
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8531)
@@ -34,5 +34,5 @@
    REAL(wp) , PARAMETER ::   r1_c2 = 1. / rc2
    !
-   ! ** albedo namelist (namice_alb)
+   ! ** albedo namelist (namalb)
    INTEGER  ::   nn_ice_alb       ! type of albedo scheme: 0: Shine & Henderson-Sellers (JGR 1985)
    !                                      !                         1: "home made" based on Brandt et al. (JClim 2005)
@@ -310,19 +310,19 @@
       !! ** Purpose :   initializations for the albedo parameters
       !!
-      !! ** Method  :   Read the namelist namice_alb
+      !! ** Method  :   Read the namelist namalb
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_alb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
-      !!----------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )              ! Namelist namice_alb in reference namelist : Albedo parameters
-      READ  ( numnam_ice_ref, namice_alb, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namice_alb in configuration namelist : Albedo parameters
-      READ  ( numnam_ice_cfg, namice_alb, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_alb )
+      NAMELIST/namalb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namalb in reference namelist : Albedo parameters
+      READ  ( numnam_ice_ref, namalb, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namalb in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namalb in configuration namelist : Albedo parameters
+      READ  ( numnam_ice_cfg, namalb, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namalb in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namalb )
       !
       IF(lwp) THEN                      ! Control print
@@ -330,5 +330,5 @@
          WRITE(numout,*) 'ice_alb_init: set albedo parameters'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_alb : albedo '
+         WRITE(numout,*) '   Namelist namalb:'
          WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb
          WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8531)
@@ -48,5 +48,5 @@
    INTEGER FUNCTION ice_dia_alloc()
       !!---------------------------------------------------------------------!
-      !!                ***  ROUTINE ice_rdgrft_alloc ***
+      !!                ***  ROUTINE ice_dia_alloc ***
       !!---------------------------------------------------------------------!
       ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc )
@@ -175,15 +175,15 @@
       INTEGER            ::   ios, ierror   ! local integer
       !!
-      NAMELIST/namice_dia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt  
+      NAMELIST/namdia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt  
       !!----------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )      ! Namelist namice_dia in reference namelist : Parameters for ice
-      READ  ( numnam_ice_ref, namice_dia, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dia in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )      ! Namelist namice_dia in configuration namelist : Parameters for ice
-      READ  ( numnam_ice_cfg, namice_dia, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dia in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_dia )
+      REWIND( numnam_ice_ref )      ! Namelist namdia in reference namelist : Parameters for ice
+      READ  ( numnam_ice_ref, namdia, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdia in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )      ! Namelist namdia in configuration namelist : Parameters for ice
+      READ  ( numnam_ice_cfg, namdia, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdia in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namdia )
       !
       IF(lwp) THEN                  ! control print
@@ -191,8 +191,8 @@
          WRITE(numout,*) 'ice_dia_init: ice diagnostics'
          WRITE(numout,*) ' ~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_dia : '
+         WRITE(numout,*) '   Namelist namdia:'
          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,*) '      control prints for a given grid point      ln_icectl    = ', ln_icectl
          WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')'
       ENDIF
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8531)
@@ -25,4 +25,5 @@
    USE lib_mpp        ! MPP library
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager
    USE lib_fortran    ! glob_sum
    USE timing         ! Timing
@@ -209,24 +210,24 @@
       !!      dynamics
       !!
-      !! ** Method  :  Read the namice_dyn namelist and check the ice-dynamic
+      !! ** Method  :  Read the namdyn namelist and check the ice-dynamic
       !!       parameter values called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namice_dyn
+      !! ** input   :   Namelist namdyn
       !!-------------------------------------------------------------------
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_dyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
-         &                 rn_ishlat  , rn_cio   ,                         &
-         &                 ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax
-      !!-------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )         ! Namelist namice_dyn in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namice_dyn, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dyn in reference namelist', lwp )
-      !
-      REWIND( numnam_ice_cfg )         ! Namelist namice_dyn in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namice_dyn, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_dyn in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_dyn )
+      NAMELIST/namdyn/ ln_dynFULL, ln_dynRHGADV, ln_dynADV, rn_uice, rn_vice,  &
+         &             rn_ishlat  ,                                            &
+         &             ln_landfast, rn_gamma , rn_icebfr, rn_lfrelax
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )         ! Namelist namdyn in reference namelist : Ice dynamics
+      READ  ( numnam_ice_ref, namdyn, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn in reference namelist', lwp )
+      !
+      REWIND( numnam_ice_cfg )         ! Namelist namdyn in configuration namelist : Ice dynamics
+      READ  ( numnam_ice_cfg, namdyn, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namdyn )
       !
       IF(lwp) THEN                     ! control print
@@ -234,5 +235,5 @@
          WRITE(numout,*) 'ice_dyn_init: ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_dyn'
+         WRITE(numout,*) '   Namelist namdyn:'
          WRITE(numout,*) '      Full ice dynamics      (rhg + adv + ridge/raft + corr)  ln_dynFULL   = ', ln_dynFULL
          WRITE(numout,*) '      No ridge/raft & No cor (rhg + adv)                      ln_dynRHGADV = ', ln_dynRHGADV
@@ -240,10 +241,9 @@
          WRITE(numout,*) '           with prescribed velocity given by '
          WRITE(numout,*) '               a uniform field               (u,v)_ice = (rn_uice,rn_vice) = (', rn_uice,',', rn_vice,')'
-         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat     = ', rn_ishlat
-         WRITE(numout,*) '      drag coefficient for oceanic stress                    rn_cio        = ', rn_cio
-         WRITE(numout,*) '      Landfast: param (T or F)                               ln_landfast   = ', ln_landfast
-         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_gamma      = ', rn_gamma
-         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr     = ', rn_icebfr
-         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax    = ', rn_lfrelax
+         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics         rn_ishlat    = ', rn_ishlat
+         WRITE(numout,*) '      Landfast: param (T or F)                                ln_landfast  = ', ln_landfast
+         WRITE(numout,*) '         fraction of ocean depth that ice must reach          rn_gamma     = ', rn_gamma
+         WRITE(numout,*) '         maximum bottom stress per unit area of contact       rn_icebfr    = ', rn_icebfr
+         WRITE(numout,*) '         relax time scale (s-1) to reach static friction      rn_lfrelax   = ', rn_lfrelax
       ENDIF
       !                             !== set the choice of ice dynamics ==!
@@ -267,5 +267,9 @@
       ENDIF
       !                                      !--- NO Landfast ice : set to zero once for all
-      IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp 
+      IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp
+      !
+      CALL ice_rdgrft_init          ! set ice ridging/rafting parameters
+      CALL ice_rhg_init             ! set ice rheology parameters
+      CALL ice_adv_init             ! set ice advection parameters
       !
    END SUBROUTINE ice_dyn_init
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8531)
@@ -30,6 +30,7 @@
    PRIVATE
 
-   PUBLIC ice_forcing_tau  ! routine called by icestp.F90
-   PUBLIC ice_forcing_flx  ! routine called by icestp.F90
+   PUBLIC ice_forcing_tau   ! called by icestp.F90
+   PUBLIC ice_forcing_flx   ! called by icestp.F90
+   PUBLIC ice_forcing_init  ! called by icestp.F90
 
    !! * Substitutions
@@ -257,4 +258,59 @@
    END SUBROUTINE ice_flx_dist
 
+   SUBROUTINE ice_forcing_init
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_forcing_init  ***
+      !!
+      !! ** Purpose : Physical constants and parameters linked to the ice
+      !!      dynamics
+      !!
+      !! ** Method  :  Read the namforcing namelist and check the ice-dynamic
+      !!       parameter values called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namforcing
+      !!-------------------------------------------------------------------
+      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )         ! Namelist namforcing in reference namelist : Ice dynamics
+      READ  ( numnam_ice_ref, namforcing, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namforcing in reference namelist', lwp )
+      !
+      REWIND( numnam_ice_cfg )         ! Namelist namforcing in configuration namelist : Ice dynamics
+      READ  ( numnam_ice_cfg, namforcing, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namforcing in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namforcing )
+      !
+      IF(lwp) THEN                     ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_forcing_init: ice parameters for ice dynamics '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namforcing:'
+         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio    = ', rn_cio
+         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s = ', rn_blow_s
+         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx = ', nn_iceflx
+      ENDIF
+      !
+      IF(lwp) WRITE(numout,*)
+      SELECT CASE( nn_iceflx )         ! ESIM Multi-category heat flux formulation
+      CASE( -1  )
+         IF(lwp) WRITE(numout,*) '   ESIM: use per-category fluxes (nn_iceflx = -1) '
+         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for ESIM in coupled mode must be 0 or 2' )
+      CASE(  0  )
+         IF(lwp) WRITE(numout,*) '   ESIM: use average per-category fluxes (nn_iceflx = 0) '
+      CASE(  1  )
+         IF(lwp) WRITE(numout,*) '   ESIM: use average then redistribute per-category fluxes (nn_iceflx = 1) '
+         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for ESIM in coupled mode must be 0 or 2' )
+      CASE(  2  )
+         IF(lwp) WRITE(numout,*) '   ESIM: Redistribute a single flux over categories (nn_iceflx = 2) '
+         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_iceflx for ESIM in forced mode cannot be 2' )
+      CASE DEFAULT
+         CALL ctl_stop( 'ice_thd_init: ESIM option, nn_iceflx, should be between -1 and 2' )
+      END SELECT
+      !
+   END SUBROUTINE ice_forcing_init
+
 #else
    !!----------------------------------------------------------------------
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8531)
@@ -49,5 +49,5 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read)
    !
-   ! ** namelist (namice_ini) **
+   ! ** namelist (namini) **
    LOGICAL  ::   ln_iceini        ! initialization or not
    LOGICAL  ::   ln_iceini_file   ! Ice initialization state from 2D netcdf file
@@ -512,9 +512,9 @@
       !! ** Purpose : Definition of initial state of the ice 
       !!
-      !! ** Method : Read the namice_ini namelist and check the parameter 
+      !! ** Method : Read the namini namelist and check the parameter 
       !!       values called at the first timestep (nit000)
       !!
       !! ** input : 
-      !!        Namelist namice_ini
+      !!        Namelist namini
       !!
       !! history :
@@ -531,18 +531,18 @@
       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
+      NAMELIST/namini/ 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 )
+      REWIND( numnam_ice_ref )              ! Namelist namini in reference namelist : Ice initial state
+      READ  ( numnam_ice_ref, namini, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namini in configuration namelist : Ice initial state
+      READ  ( numnam_ice_cfg, namini, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namini in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namini )
 
       slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
@@ -555,5 +555,5 @@
          WRITE(numout,*) 'ice_istate_init: ice parameters inititialisation '
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_ini'
+         WRITE(numout,*) '   Namelist namini:'
          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
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8531)
@@ -38,5 +38,5 @@
    PUBLIC   ice_itd_reb   ! called in iceerr
 
-   ! ** ice-thickness distribution namelist (namice_itd) **
+   ! ** ice-thickness distribution namelist (namitd) **
    REAL(wp) ::   rn_himean        ! mean thickness of the domain (used to compute the distribution)
 
@@ -643,5 +643,5 @@
       !! ** Purpose :   Initializes the ice thickness distribution
       !! ** Method  :   ...
-      !! ** input   :   Namelist namice_itd
+      !! ** input   :   Namelist namitd
       !!-------------------------------------------------------------------
       INTEGER  ::   jl    ! dummy loop index
@@ -649,15 +649,15 @@
       REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      -
       !!
-      NAMELIST/namice_itd/ rn_himean, rn_himin
-      !!------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )      ! Namelist namice_itd in reference namelist : Parameters for ice
-      READ  ( numnam_ice_ref, namice_itd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_itd in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )      ! Namelist namice_itd in configuration namelist : Parameters for ice
-      READ  ( numnam_ice_cfg, namice_itd, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_itd in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_itd )
+      NAMELIST/namitd/ rn_himean, rn_himin
+      !!------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )      ! Namelist namitd in reference namelist : Parameters for ice
+      READ  ( numnam_ice_ref, namitd, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )      ! Namelist namitd in configuration namelist : Parameters for ice
+      READ  ( numnam_ice_cfg, namitd, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namitd in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namitd )
       !
       IF(lwp) THEN                  ! control print
@@ -665,5 +665,5 @@
          WRITE(numout,*) 'ice_itd_init: Initialization of ice cat distribution '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_itd : '
+         WRITE(numout,*) '   Namelist namitd: '
          WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean
          WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin 
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8531)
@@ -12,5 +12,4 @@
    !!   'key_lim3'                                      LIM-3 sea-ice model
    !!----------------------------------------------------------------------
-   USE par_oce        ! ocean parameters
    USE dom_oce        ! ocean domain
    USE phycst         ! physical constants (ocean directory) 
@@ -31,24 +30,24 @@
    PRIVATE
 
-   PUBLIC   ice_rdgrft               ! called by ice_stp
-   PUBLIC   ice_rdgrft_strength      ! called by icerhg_evp
-   PUBLIC   ice_rdgrft_init          ! called by ice_stp
+   PUBLIC   ice_rdgrft        ! called by icestp
+   PUBLIC   ice_strength      ! called by icerhg_evp
+   PUBLIC   ice_rdgrft_init   ! called by icedyn
 
    ! Variables shared among ridging subroutines
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apartf   ! participation function; fraction of ridging/closing associated w/ category n
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hi_hrdg  ! thickness of ridging ice / mean ridge thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
    !
-   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
-   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
-   REAL(wp)            ::   zdrho               ! 
+   REAL(wp), PARAMETER ::   hrdg_hi_min = 1.1_wp    ! min ridge thickness multiplier: min(hrdg/hi)
+   REAL(wp), PARAMETER ::   hi_hrft     = 0.5_wp    ! rafting multipliyer: (hi/hraft)
+   REAL(wp)            ::   zdrho                   ! 
    !
-   ! ** namelist (namice_rdgrft) **
+   ! ** namelist (namdyn_rdgrft) **
    LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
    REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
@@ -83,7 +82,7 @@
       !!                ***  ROUTINE ice_rdgrft_alloc ***
       !!---------------------------------------------------------------------!
-      ALLOCATE( asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     &
-         &      hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     &
-         &      hrmax(jpi,jpj,jpl) , krdg  (jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
+      ALLOCATE( asum (jpi,jpj)     , apartf (jpi,jpj,0:jpl) , aksum (jpi,jpj)     ,     &
+         &      hrmin(jpi,jpj,jpl) , hraft  (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     &
+         &      hrmax(jpi,jpj,jpl) , hi_hrdg(jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
 
       IF( lk_mpp                )   CALL mpp_sum ( ice_rdgrft_alloc )
@@ -124,5 +123,5 @@
       INTEGER  ::   niter              ! local integer 
       INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging
-      REAL(wp) ::   za                 ! local scalar
+      REAL(wp) ::   zfac               ! local scalar
       REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s)
       !                                                 ! (ridging ice area - area of new ridges) / dt
@@ -150,5 +149,5 @@
       !-----------------------------------------------------------------------------!
       !
-      CALL ice_rdgrft_prep                             ! prepare ridging
+      CALL rdgrft_prep                             ! prepare ridging
       !
       DO jj = 1, jpj                                        ! Initialize arrays.
@@ -203,5 +202,5 @@
       DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
 
-         ! 3.2 closing_gross
+         ! 3.1 closing_gross
          !-----------------------------------------------------------------------------!
          ! Based on the ITD of ridging and ridged ice, convert the net
@@ -216,9 +215,9 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
-               IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
-                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
-               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
-                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
+               zfac   = ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
+               IF    ( zfac < 0._wp .AND. zfac > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i
+                  opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
+               ELSEIF( zfac > 0._wp .AND. zfac > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
+                  opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
                ENDIF
             END DO
@@ -232,7 +231,7 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
-                  IF( za  >  a_i(ji,jj,jl) ) THEN
-                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za
+                  zfac = apartf(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
+                  IF( zfac  >  a_i(ji,jj,jl) ) THEN
+                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / zfac
                   ENDIF
                END DO
@@ -240,14 +239,14 @@
          END DO
 
-         ! 3.3 Redistribute area, volume, and energy.
+         ! 3.2 Redistribute area, volume, and energy.
          !-----------------------------------------------------------------------------!
-         CALL ice_rdgrft_ridgeshift( opning, closing_gross )
+         CALL rdgrft_shift( opning, closing_gross )
          
-         ! 3.4 Compute total area of ice plus open water after ridging.
+         ! 3.3 Compute total area of ice plus open water after ridging.
          !-----------------------------------------------------------------------------!
          ! This is in general not equal to one because of divergence during transport
          asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
 
-         ! 3.5 Do we keep on iterating?
+         ! 3.4 Do we keep on iterating?
          !-----------------------------------------------------------------------------!
          ! Check whether asum = 1.  If not (because the closing and opening
@@ -276,5 +275,5 @@
          !
          IF( iterate_ridging == 1 ) THEN
-            CALL ice_rdgrft_prep
+            CALL rdgrft_prep
             IF( niter  >  nitermax ) THEN
                WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
@@ -295,7 +294,7 @@
 
 
-   SUBROUTINE ice_rdgrft_prep
+   SUBROUTINE rdgrft_prep
       !!---------------------------------------------------------------------!
-      !!                ***  ROUTINE ice_rdgrft_prep ***
+      !!                ***  ROUTINE rdgrft_prep ***
       !!
       !! ** Purpose :   preparation for ridging and strength calculations
@@ -319,7 +318,17 @@
       END WHERE
 
-      !------------------------------------------------------------------------------!
-      ! 1) Participation function 
-      !------------------------------------------------------------------------------!
+      !-----------------------------------------------------------------
+      ! 1) Participation function: a(h) = b(h).g(h) (apartf)
+      !-----------------------------------------------------------------
+      ! Compute the participation function apartf; this is analogous to
+      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
+      ! area lost from category n due to ridging/closing
+      ! apartf(n)   = total area lost due to ridging/closing
+      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 
+      !
+      ! The expressions for apartf are found by integrating b(h)g(h) between
+      ! the category boundaries.
+      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1
+      !-----------------------------------------------------------------
       !
       ! Compute total area of ice plus open water.
@@ -337,16 +346,4 @@
       END DO
 
-      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
-      !--------------------------------------------------------------------------------------------------
-      ! Compute the participation function athorn; this is analogous to
-      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
-      ! area lost from category n due to ridging/closing
-      ! athorn(n)   = total area lost due to ridging/closing
-      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 
-      !
-      ! The expressions for athorn are found by integrating b(h)g(h) between
-      ! the category boundaries.
-      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1
-      !-----------------------------------------------------------------
       !
       IF( ln_partf_lin ) THEN          !--- Linear formulation (Thorndike et al., 1975)
@@ -355,11 +352,11 @@
                DO ji = 1, jpi
                   IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN
-                     athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
+                     apartf(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * &
                         &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar )
                   ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN
-                     athorn(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
+                     apartf(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  &
                         &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar )
                   ELSE
-                     athorn(ji,jj,jl) = 0._wp
+                     apartf(ji,jj,jl) = 0._wp
                   ENDIF
                END DO
@@ -374,5 +371,5 @@
          END DO
          DO jl = 0, jpl
-            athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
+            apartf(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl)
          END DO
          !
@@ -384,15 +381,15 @@
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl)
-                  araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl)
+                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jj,jl)
+                  araft (ji,jj,jl) = apartf(ji,jj,jl) - aridge(ji,jj,jl)
                END DO
             END DO
          END DO
       ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone
-         aridge(:,:,:) = athorn(:,:,1:jpl)
+         aridge(:,:,:) = apartf(:,:,1:jpl)
          araft (:,:,:) = 0._wp
       ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone   
         aridge(:,:,:) = 0._wp
-         araft (:,:,:) = athorn(:,:,1:jpl)
+         araft (:,:,:) = apartf(:,:,1:jpl)
       ELSE                                               !- no ridging & no rafting
          aridge(:,:,:) = 0._wp
@@ -408,5 +405,5 @@
       ! This parameterization is a modified version of Hibler (1980).
       ! The mean ridging thickness, zhmean, is proportional to hi^(0.5)
-      !  and for very thick ridging ice must be >= krdgmin*hi
+      !  and for very thick ridging ice must be >= hrdg_hi_min*hi
       !
       ! The minimum ridging thickness, hrmin, is equal to 2*hi 
@@ -426,24 +423,25 @@
       !-----------------------------------------------------------------
 
-      aksum(:,:) = athorn(:,:,0)
+      aksum(:,:) = apartf(:,:,0)
+      zdummy = 1._wp / hi_hrft
       ! Transfer function
       DO jl = 1, jpl !all categories have a specific transfer function
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( athorn(ji,jj,jl) > 0._wp ) THEN
-                  zhmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
-                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
-                  hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
-                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
-                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
+               IF ( apartf(ji,jj,jl) > 0._wp ) THEN
+                  zhmean            = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * hrdg_hi_min )
+                  hrmin  (ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) )
+                  hrmax  (ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl)
+                  hraft  (ji,jj,jl) = ht_i(ji,jj,jl) * zdummy
+                  hi_hrdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 )
                   !
                   ! Normalization factor : aksum, ensures mass conservation
-                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) )    &
-                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          )
+                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - hi_hrdg(ji,jj,jl) )    &
+                     &                        + araft (ji,jj,jl) * ( 1._wp - hi_hrft )
                ELSE
-                  hrmin(ji,jj,jl)  = 0._wp 
-                  hrmax(ji,jj,jl)  = 0._wp 
-                  hraft(ji,jj,jl)  = 0._wp 
-                  krdg (ji,jj,jl)  = 1._wp
+                  hrmin  (ji,jj,jl) = 0._wp 
+                  hrmax  (ji,jj,jl) = 0._wp 
+                  hraft  (ji,jj,jl) = 0._wp 
+                  hi_hrdg(ji,jj,jl) = 1._wp
                ENDIF
             END DO
@@ -451,10 +449,10 @@
       END DO
       !
-   END SUBROUTINE ice_rdgrft_prep
-
-
-   SUBROUTINE ice_rdgrft_ridgeshift( opning, closing_gross )
+   END SUBROUTINE rdgrft_prep
+
+
+   SUBROUTINE rdgrft_shift( opning, closing_gross )
       !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ice_rdgrft_strength ***
+      !!                ***  ROUTINE rdgrft_shift ***
       !!
       !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
@@ -463,14 +461,14 @@
       !!              and add to thicker ice categories.
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges
-      !
-      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
-      INTEGER ::   ij                ! horizontal index, combines i and j loops
-      INTEGER ::   icells            ! number of cells with a_i > puny
-      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration and new area going to jl2
-
-      INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices
-      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to jl2
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   opning         ! rate of opening due to divergence/shear
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges
+      !
+      INTEGER  ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
+      INTEGER  ::   ij                         ! horizontal index, combines i and j loops
+      INTEGER  ::   icells                     ! number of cells with a_i > puny
+      REAL(wp) ::   hL, hR, farea              ! left and right limits of integration and new area going to jl2
+
+      INTEGER , DIMENSION(jpij) ::   indxi, indxj     ! compressed indices
+      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol    ! new ridge volume going to jl2
 
       REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged 
@@ -482,12 +480,12 @@
       REAL(wp), DIMENSION(jpij) ::   aprdg2           ! pond area of ridging ice
       ! END MV MP 2016
-      REAL(wp), DIMENSION(jpij) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
-
-      REAL(wp), DIMENSION(jpij) ::   vrdg1   ! volume of ice ridged
-      REAL(wp), DIMENSION(jpij) ::   vrdg2   ! volume of new ridges
-      REAL(wp), DIMENSION(jpij) ::   vsw     ! volume of seawater trapped into ridges
-      REAL(wp), DIMENSION(jpij) ::   srdg1   ! sal*volume of ice ridged
-      REAL(wp), DIMENSION(jpij) ::   srdg2   ! sal*volume of new ridges
-      REAL(wp), DIMENSION(jpij) ::   smsw    ! sal*volume of water trapped into ridges
+      REAL(wp), DIMENSION(jpij) ::   dhr, dhr2        ! hrmax - hrmin  &  hrmax^2 - hrmin^2
+
+      REAL(wp), DIMENSION(jpij) ::   vrdg1            ! volume of ice ridged
+      REAL(wp), DIMENSION(jpij) ::   vrdg2            ! volume of new ridges
+      REAL(wp), DIMENSION(jpij) ::   vsw              ! volume of seawater trapped into ridges
+      REAL(wp), DIMENSION(jpij) ::   srdg1            ! sal*volume of ice ridged
+      REAL(wp), DIMENSION(jpij) ::   srdg2            ! sal*volume of new ridges
+      REAL(wp), DIMENSION(jpij) ::   smsw             ! sal*volume of water trapped into ridges
       REAL(wp), DIMENSION(jpij) ::   oirdg1, oirdg2   ! ice age of ice ridged
 
@@ -503,8 +501,8 @@
       REAL(wp), DIMENSION(jpij) ::   oirft1, oirft2   ! ice age of ice rafted
 
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft      ! ice energy of rafting ice
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1      ! enth*volume of ice ridged
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2      ! enth*volume of new ridges
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw       ! enth of water trapped into ridges
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirft     ! ice energy of rafting ice
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg1     ! enth*volume of ice ridged
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   erdg2     ! enth*volume of new ridges
+      REAL(wp), DIMENSION(jpij,nlay_i) ::   ersw      ! enth of water trapped into ridges
       !!----------------------------------------------------------------------
 
@@ -515,10 +513,10 @@
          DO ji = 1, jpi
             ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
-               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
+               &                     ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
          END DO
       END DO
 
       !-----------------------------------------------------------------
-      ! 3) Pump everything from ice which is being ridged / rafted
+      ! 2) Pump everything from ice which is being ridged / rafted
       !-----------------------------------------------------------------
       ! Compute the area, volume, and energy of ice ridging in each
@@ -528,10 +526,10 @@
 
          !------------------------------------------------
-         ! 3.1) Identify grid cells with nonzero ridging
+         ! 2.1) Identify grid cells with nonzero ridging
          !------------------------------------------------
          icells = 0
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
+               IF( apartf(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
                   icells = icells + 1
                   indxi(icells) = ji
@@ -545,5 +543,5 @@
 
             !--------------------------------------------------------------------
-            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
+            ! 2.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
             !--------------------------------------------------------------------
             ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
@@ -551,13 +549,13 @@
 
             !---------------------------------------------------------------
-            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 
+            ! 2.3) Compute ridging /rafting fractions, make sure afrac <=1 
             !---------------------------------------------------------------
             afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
             afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
-            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
-            arft2(ij) = arft1(ij) * kraft
+            ardg2(ij) = ardg1(ij) * hi_hrdg(ji,jj,jl1)
+            arft2(ij) = arft1(ij) * hi_hrft
 
             !--------------------------------------------------------------------------
-            ! 3.4) Substract area, volume, and energy from ridging 
+            ! 2.4) Substract area, volume, and energy from ridging 
             !     / rafting category n1.
             !--------------------------------------------------------------------------
@@ -571,5 +569,5 @@
             IF ( nn_pnd_scheme > 0 ) THEN
                aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij)
-               aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
+               aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1)
                vprdg(ij)  = v_ip(ji,jj, jl1) * afrac(ij)
             ENDIF
@@ -577,5 +575,5 @@
             srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
             oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
-            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1) 
+            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1) 
 
             ! rafting volumes, heat contents ...
@@ -585,5 +583,5 @@
             IF ( nn_pnd_scheme > 0 ) THEN
                aprft1(ij) = a_ip (ji,jj,  jl1) * afrft(ij)
-               aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * kraft
+               aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * hi_hrft
                vprft(ij)  = v_ip(ji,jj,jl1)    * afrft(ij)
             ENDIF
@@ -593,8 +591,8 @@
             smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
             oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
-            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft 
+            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * hi_hrft 
 
             !-----------------------------------------------------------------
-            ! 3.5) Compute properties of new ridges
+            ! 2.5) Compute properties of new ridges
             !-----------------------------------------------------------------
             smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
@@ -612,5 +610,5 @@
 
             !------------------------------------------            
-            ! 3.7 Put the snow somewhere in the ocean
+            ! 2.6 Put the snow somewhere in the ocean
             !------------------------------------------            
             !  Place part of the snow lost by ridging into the ocean. 
@@ -627,5 +625,5 @@
             ! MV MP 2016
             !------------------------------------------            
-            ! 3.X Put the melt pond water in the ocean
+            ! 2.7 Put the melt pond water in the ocean
             !------------------------------------------            
             !  Place part of the melt pond volume into the ocean. 
@@ -637,5 +635,5 @@
 
             !-----------------------------------------------------------------
-            ! 3.8 Compute quantities used to apportion ice among categories
+            ! 2.8 Compute quantities used to apportion ice among categories
             ! in the n2 loop below
             !-----------------------------------------------------------------
@@ -662,5 +660,5 @@
 
          !--------------------------------------------------------------------
-         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
+         ! 2.9 Compute ridging ice enthalpy, remove it from ridging ice and
          !      compute ridged ice enthalpy 
          !--------------------------------------------------------------------
@@ -689,5 +687,5 @@
 
          !-------------------------------------------------------------------------------
-         ! 4) Add area, volume, and energy of new ridge to each category jl2
+         ! 3) Add area, volume, and energy of new ridge to each category jl2
          !-------------------------------------------------------------------------------
          DO jl2  = 1, jpl 
@@ -710,5 +708,5 @@
 !!gm see above               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
                IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ij) = 1._wp
-               ELSE                                                                                 ;   zswitch(ij) = 0._wp                  
+               ELSE                                                                                 ;   zswitch(ij) = 0._wp
                ENDIF
                !
@@ -743,10 +741,10 @@
       END DO ! jl1 (deforming categories)
       !
-   END SUBROUTINE ice_rdgrft_ridgeshift
-
-
-   SUBROUTINE ice_rdgrft_strength
+   END SUBROUTINE rdgrft_shift
+
+
+   SUBROUTINE ice_strength
       !!----------------------------------------------------------------------
-      !!                ***  ROUTINE ice_rdgrft_strength ***
+      !!                ***  ROUTINE ice_strength ***
       !!
       !! ** Purpose :   computes ice strength used in dynamics routines of ice thickness
@@ -767,5 +765,5 @@
 
       !                              !--------------------------------------------------!
-      CALL ice_rdgrft_prep           ! Thickness distribution of ridging and ridged ice !
+      CALL rdgrft_prep               ! Thickness distribution of ridging and ridged ice !
       !                              !--------------------------------------------------!
 
@@ -775,11 +773,11 @@
          z1_3 = 1._wp / 3._wp
          DO jl = 1, jpl
-            WHERE( athorn(:,:,jl) > 0._wp )
-               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
-                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
-                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice
-                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         &
-                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         &
-                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) )
+            WHERE( apartf(:,:,jl) > 0._wp )
+               strength(:,:) =  -         apartf(:,:,jl) * ht_i   (:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice
+                  &             + 2._wp * araft (:,:,jl) * ht_i   (:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice
+                  &             +         aridge(:,:,jl) * hi_hrdg(:,:,jl) * z1_3 *         &  ! PE gain from ridging ice
+                  &                      ( hrmax(:,:,jl) * hrmax  (:,:,jl) +                &
+                  &                        hrmin(:,:,jl) * hrmin  (:,:,jl) +                &
+                  &                        hrmax(:,:,jl) * hrmin  (:,:,jl) )
             ELSEWHERE
                strength(:,:) = 0._wp
@@ -844,5 +842,5 @@
       END SELECT
       !
-   END SUBROUTINE ice_rdgrft_strength
+   END SUBROUTINE ice_strength
 
 
@@ -854,13 +852,13 @@
       !!                to the mechanical ice redistribution
       !!
-      !! ** Method  :   Read the namice_rdgrft namelist 
+      !! ** Method  :   Read the namdyn_rdgrft namelist 
       !!                and check the parameters values 
       !!                called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namice_rdgrft
+      !! ** input   :   Namelist namdyn_rdgrft
       !!-------------------------------------------------------------------
       INTEGER :: ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
+      NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, &
          &                    ln_str_R75, rn_perdg,          &
          &                    rn_csrdg  ,                    &
@@ -872,11 +870,11 @@
       !
       REWIND( numnam_ice_ref )              ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution
-      READ  ( numnam_ice_ref, namice_rdgrft, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in reference namelist', lwp )
-      !
-      REWIND( numnam_ice_cfg )              ! Namelist namice_rdgrft in configuration namelist : Ice mechanical ice redistribution
-      READ  ( numnam_ice_cfg, namice_rdgrft, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rdgrft in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_rdgrft )
+      READ  ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp )
+      !
+      REWIND( numnam_ice_cfg )              ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution
+      READ  ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namdyn_rdgrft )
       !
       IF (lwp) THEN                          ! control print
@@ -884,5 +882,5 @@
          WRITE(numout,*) 'ice_rdgrft_init: ice parameters for ridging/rafting '
          WRITE(numout,*) '~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_rdgrft'
+         WRITE(numout,*) '   Namelist namdyn_rdgrft:'
          WRITE(numout,*) '      ice strength parameterization Hibler (1979)              ln_str_H79   = ', ln_str_H79 
          WRITE(numout,*) '            1st bulk-rheology parameter                        rn_pstar     = ', rn_pstar
@@ -915,5 +913,5 @@
       ENDIF
       !                              ! allocate tke arrays
-      IF( ice_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_rdgrft_init : unable to allocate arrays' )
+      IF( ice_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_rdgrft_init: unable to allocate arrays' )
       !
   END SUBROUTINE ice_rdgrft_init
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8531)
@@ -101,22 +101,22 @@
       !!      dynamics
       !!
-      !! ** Method  :  Read the namice_rhg namelist and check the ice-dynamic
+      !! ** Method  :  Read the namdyn_rhg namelist and check the ice-dynamic
       !!       parameter values called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namice_rhg
+      !! ** input   :   Namelist namdyn_rhg
       !!-------------------------------------------------------------------
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
       !!-------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )         ! Namelist namice_rhg in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namice_rhg, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rhg in reference namelist', lwp )
+      REWIND( numnam_ice_ref )         ! Namelist namdyn_rhg in reference namelist : Ice dynamics
+      READ  ( numnam_ice_ref, namdyn_rhg, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rhg in reference namelist', lwp )
       !
-      REWIND( numnam_ice_cfg )         ! Namelist namice_rhg in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namice_rhg, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_rhg in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_rhg )
+      REWIND( numnam_ice_cfg )         ! Namelist namdyn_rhg in configuration namelist : Ice dynamics
+      READ  ( numnam_ice_cfg, namdyn_rhg, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rhg in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namdyn_rhg )
       !
       IF(lwp) THEN                     ! control print
@@ -124,10 +124,10 @@
          WRITE(numout,*) 'ice_rhg_init: ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_rhg'
-         WRITE(numout,*) '      rheology EVP (icerhg_evp)                                   ln_rhg_EVP    = ', ln_rhg_EVP
-         WRITE(numout,*) '         creep limit                                              rn_creepl     = ', rn_creepl
-         WRITE(numout,*) '         eccentricity of the elliptical yield curve               rn_ecc        = ', rn_ecc
-         WRITE(numout,*) '         number of iterations for subcycling                      nn_nevp       = ', nn_nevp
-         WRITE(numout,*) '         ratio of elastic timescale over ice time step            rn_relast     = ', rn_relast
+         WRITE(numout,*) '   Namelist namdyn_rhg:'
+         WRITE(numout,*) '      rheology EVP (icerhg_evp)                            ln_rhg_EVP = ', ln_rhg_EVP
+         WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl
+         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc
+         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp    = ', nn_nevp
+         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast  = ', rn_relast
       ENDIF
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8531)
@@ -253,5 +253,5 @@
 
       ! Ice strength
-      CALL ice_rdgrft_strength
+      CALL ice_strength
 
       ! scale factors
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8531)
@@ -142,5 +142,5 @@
          ! --- ice dynamics and advection  --- !
          !-------------------------------------!
-         CALL ice_diag0             ! set diag of mass, heat and salt fluxes to 0
+         CALL diag_set0             ! set diag of mass, heat and salt fluxes to 0
          CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary) 
          !
@@ -243,5 +243,5 @@
       IF(lwm) CALL ctl_opn( numoni, 'output.namelist.ice', 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, 1 )
       !
-      CALL ice_run_init                ! set some ice run parameters
+      CALL par_init                ! set some ice run parameters
       !
       !                                ! Allocate the ice arrays (sbc_ice already allocated in sbc_init)
@@ -268,14 +268,12 @@
       CALL ice_var_glo2eqv
       !
+      CALL ice_forcing_init            ! set ice-ocean and ice-atm. coupling parameters
+      !
       IF( ln_icedyn ) THEN
          CALL ice_dyn_init             ! set ice dynamics parameters
-         CALL ice_rdgrft_init          ! set ice ridging/rafting parameters
-         CALL ice_rhg_init             ! set ice rheology parameters
-         CALL ice_adv_init             ! set ice advection parameters
       ENDIF
-
+      !
       IF( ln_icethd ) THEN
          CALL ice_thd_init             ! set ice thermodynics parameters
-         CALL ice_thd_sal_init         ! set ice salinity parameters
       ENDIF   
       !
@@ -299,40 +297,40 @@
 
 
-   SUBROUTINE ice_run_init
+   SUBROUTINE par_init
       !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_run_init ***
+      !!                  ***  ROUTINE par_init ***
       !!
       !! ** Purpose :   Definition some run parameter for ice model
       !!
-      !! ** Method  :   Read the namice_run namelist and check the parameter
+      !! ** Method  :   Read the nampar namelist and check the parameter
       !!                values called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namice_run
+      !! ** input   :   Namelist nampar
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  &
-         &                 cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
+      NAMELIST/nampar/ jpl, nlay_i, nlay_s, nn_monocat, ln_icedyn, ln_icethd, rn_amax_n, rn_amax_s,  &
+         &             cn_icerst_in, cn_icerst_indir, cn_icerst_out, cn_icerst_outdir
       !!-------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )      ! Namelist namice_run in reference namelist : Parameters for ice
-      READ  ( numnam_ice_ref, namice_run, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_run in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )      ! Namelist namice_run in configuration namelist : Parameters for ice
-      READ  ( numnam_ice_cfg, namice_run, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_run in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_run )
+      REWIND( numnam_ice_ref )      ! Namelist nampar in reference namelist : Parameters for ice
+      READ  ( numnam_ice_ref, nampar, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampar in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )      ! Namelist nampar in configuration namelist : Parameters for ice
+      READ  ( numnam_ice_cfg, nampar, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampar in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, nampar )
       !
       IF(lwp) THEN                  ! control print
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
-         WRITE(numout,*) ' ~~~~~~'
-         WRITE(numout,*) '   Namelist namice_run : '
+         WRITE(numout,*) 'par_init: ice parameters shared among all the routines'
+         WRITE(numout,*) ' ~~~~~~~'
+         WRITE(numout,*) '   Namelist nampar: '
          WRITE(numout,*) '      number of ice  categories                              jpl    = ', jpl
          WRITE(numout,*) '      number of ice  layers                                  nlay_i = ', nlay_i
          WRITE(numout,*) '      number of snow layers                                  nlay_s = ', nlay_s
          WRITE(numout,*) '      virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat
-         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                  ln_icedyn  = ', ln_icedyn
-         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                  ln_icethd  = ', ln_icethd
+         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                   ln_icedyn = ', ln_icedyn
+         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                   ln_icethd = ', ln_icethd
          WRITE(numout,*) '      maximum ice concentration for NH                              = ', rn_amax_n 
          WRITE(numout,*) '      maximum ice concentration for SH                              = ', rn_amax_s
@@ -346,8 +344,8 @@
       ENDIF
       IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN
-         CALL ctl_stop( 'STOP', 'ice_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
+         CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' )
       ENDIF
       !
-      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY')
+      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
       !
       rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and inverse
@@ -358,5 +356,5 @@
       r1_nlay_s = 1._wp / REAL( nlay_s, wp )
       !
-   END SUBROUTINE ice_run_init
+   END SUBROUTINE par_init
 
 
@@ -393,7 +391,7 @@
 
 
-   SUBROUTINE ice_diag0
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_diag0  ***
+   SUBROUTINE diag_set0
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE diag_set0  ***
       !!
       !! ** purpose :  set ice-ocean and ice-atm. fluxes to zeros at the beggining
@@ -442,5 +440,5 @@
       tau_icebfr(:,:) = 0._wp; ! landfast ice param only (clem: important to keep the init here)
       
-   END SUBROUTINE ice_diag0
+   END SUBROUTINE diag_set0
 
 #else
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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8531)
@@ -28,10 +28,10 @@
       &                 fr1_i0, fr2_i0
    USE ice1D          ! thermodynamic sea-ice variables
-   USE icethd_dif     ! vertical diffusion
+   USE icethd_zdf     ! vertical diffusion
    USE icethd_dh      ! ice-snow growth and melt
    USE icethd_da      ! lateral melting
    USE icethd_sal     ! ice salinity
    USE icethd_ent     ! ice enthalpy redistribution
-   USE icethd_lac     ! lateral accretion
+   USE icethd_do      ! lateral accretion
    USE iceitd         ! remapping thickness distribution
    USE icetab         ! 1D <==> 2D transformation
@@ -51,4 +51,10 @@
    PUBLIC   ice_thd_init    ! called by ice_init
 
+   !!** namelist (namthd) **
+   LOGICAL ::   ln_icedH         ! activate ice thickness change from growing/melting (T) or not (F)
+   LOGICAL ::   ln_icedA         ! activate lateral melting param. (T) or not (F)
+   LOGICAL ::   ln_icedO         ! activate ice growth in open-water (T) or not (F)
+   LOGICAL ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -70,9 +76,11 @@
       !!               at the ice base, snow acc.,heat budget of the leads)
       !!             - selection of the icy points and put them in an array
-      !!             - call ice_thd_dif  for vertical heat diffusion
+      !!             - call ice_thd_zdf  for vertical heat diffusion
       !!             - call ice_thd_dh   for vertical ice growth and melt
       !!             - call ice_thd_ent  for enthalpy remapping
       !!             - call ice_thd_sal  for ice desalination
       !!             - call ice_thd_temp to  retrieve temperature from ice enthalpy
+      !!             - call ice_thd_lam  for extra lateral ice melt if monocat
+      !!             - call ice_thd_da   for lateral ice melt
       !!             - back to the geographic grid
       !!---------------------------------------------------------------------
@@ -223,10 +231,10 @@
             !
             IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
-                              CALL ice_thd_dif                             ! Ice/Snow Temperature profile
+                              CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
                               CALL ice_thd_dh                              ! Ice/Snow thickness   
                               CALL ice_thd_ent( e_i_1d(1:nidx,:) )         ! Ice enthalpy remapping
             ENDIF
             !
-                              CALL ice_thd_sal                      ! --- Ice salinity --- !    
+                              CALL ice_thd_sal( ln_icedS )          ! --- Ice salinity --- !    
             !
                               CALL ice_thd_temp                     ! --- temperature update --- !
@@ -261,5 +269,5 @@
       IF( jpl > 1 )        CALL ice_itd_rem( kt )          ! --- Transport ice between thickness categories --- !
       !
-      IF( ln_icedO )       CALL ice_thd_lac                ! --- frazil ice growing in leads --- !
+      IF( ln_icedO )       CALL ice_thd_do                 ! --- frazil ice growing in leads --- !
       !
       ! controls
@@ -523,81 +531,46 @@
       !!                   ***  ROUTINE ice_thd_init *** 
       !!                 
-      !! ** Purpose :   Physical constants and parameters linked to the ice 
-      !!              thermodynamics
+      !! ** Purpose :   Physical constants and parameters associated with
+      !!                ice thermodynamics
       !!
-      !! ** Method  :   Read the namice_thd namelist and check the ice-thermo
-      !!              parameter values called at the first timestep (nit000)
+      !! ** Method  :   Read the namthd namelist and check the parameters
+      !!                called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namicether
+      !! ** input   :   Namelist namthd
       !!-------------------------------------------------------------------
       INTEGER  ::   ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_thd/ 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,          &
-         &                 nn_iceflx
+      NAMELIST/namthd/ ln_icedH, ln_icedA, ln_icedO, ln_icedS
       !!-------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )              ! Namelist namice_thd in reference namelist : Ice thermodynamics
-      READ  ( numnam_ice_ref, namice_thd, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_thd in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namice_thd in configuration namelist : Ice thermodynamics
-      READ  ( numnam_ice_cfg, namice_thd, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_thd in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_thd )
+      REWIND( numnam_ice_ref )              ! Namelist namthd in reference namelist : Ice thermodynamics
+      READ  ( numnam_ice_ref, namthd, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namthd in configuration namelist : Ice thermodynamics
+      READ  ( numnam_ice_cfg, namthd, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd )
       !
       !
       IF(lwp) THEN                          ! control print
-         WRITE(numout,*) 'ice_thd_init : Ice Thermodynamics'
-         WRITE(numout,*) '~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_thd'
-         WRITE(numout,*) '   -- icethd_dif --'
-         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-         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_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_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_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_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,*) '   -- icestp --'
-         WRITE(numout,*) '      Multicategory heat flux formulation                     nn_iceflx    = ', nn_iceflx
+         WRITE(numout,*) 'ice_thd_init: Ice Thermodynamics'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namthd:'
+         WRITE(numout,*) '      activate ice thick change from top/bot (T) or not (F)   ln_icedH  = ', ln_icedH
+         WRITE(numout,*) '      activate lateral melting (T) or not (F)                 ln_icedA  = ', ln_icedA
+         WRITE(numout,*) '      activate ice growth in open-water (T) or not (F)        ln_icedO  = ', ln_icedO
+         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_icedS  = ', ln_icedS
+     ENDIF
+      !
+                       CALL ice_thd_zdf_init   ! set ice heat diffusion parameters
+      IF( ln_icedA )   CALL ice_thd_da_init    ! set ice lateral melting parameters
+      IF( ln_icedO )   CALL ice_thd_do_init    ! set ice growth in open water parameters
+                       CALL ice_thd_sal_init   ! set ice salinity parameters
+      !
+      IF( ln_icedS .AND. nn_icesal == 1 ) THEN
+         ln_icedS = .FALSE.
+         CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)')
       ENDIF
-      !
-      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( 'ice_thd_init : rn_hinew should be >= rn_himin' )
-      !
-      IF(lwp) WRITE(numout,*)
-      SELECT CASE( nn_iceflx )         ! LIM3 Multi-category heat flux formulation
-      CASE( -1  )
-         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_iceflx = 0) '
-      CASE(  1  )
-         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_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_iceflx, should be between -1 and 2' )
-      END SELECT
       !
    END SUBROUTINE ice_thd_init
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90	(revision 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90	(revision 8531)
@@ -17,4 +17,5 @@
    USE ice1D          ! sea-ice: thermodynamic 1D variables
    !
+   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)  
@@ -23,5 +24,10 @@
    PRIVATE
 
-   PUBLIC   ice_thd_da   ! called by icethd.F90
+   PUBLIC   ice_thd_da        ! called by icethd.F90
+   PUBLIC   ice_thd_da_init   ! called by icestp.F90
+
+   ! ** namelist (namthd_da) **
+   REAL(wp) ::   rn_beta     ! coef. beta for lateral melting param.
+   REAL(wp) ::   rn_dmin     ! minimum floe diameter for lateral melting param.
 
    !!----------------------------------------------------------------------
@@ -149,5 +155,42 @@
       !
    END SUBROUTINE ice_thd_da
-   
+
+   SUBROUTINE ice_thd_da_init
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_thd_da_init *** 
+      !!                 
+      !! ** Purpose :   Physical constants and parameters associated with
+      !!                ice thermodynamics
+      !!
+      !! ** Method  :   Read the namthd_da namelist and check the parameters
+      !!                called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namthd_da
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ios   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namthd_da/ rn_beta, rn_dmin
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namthd_da in reference namelist : Ice thermodynamics
+      READ  ( numnam_ice_ref, namthd_da, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_da in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namthd_da in configuration namelist : Ice thermodynamics
+      READ  ( numnam_ice_cfg, namthd_da, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_da in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd_da )
+      !
+      !
+      IF(lwp) THEN                          ! control print
+         WRITE(numout,*) 'ice_thd_da_init: Ice lateral melting'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namthd_da:'
+         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
+      ENDIF
+      !
+   END SUBROUTINE ice_thd_da_init
+  
 #else
    !!----------------------------------------------------------------------
Index: anches/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 8530)
+++ 	(revision )
@@ -1,723 +1,0 @@
-MODULE icethd_dif
-   !!======================================================================
-   !!                       ***  MODULE icethd_dif ***
-   !!                       heat diffusion in sea ice 
-   !!                   computation of surface and inner T  
-   !!======================================================================
-   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
-   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
-   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
-   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
-   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   USE par_oce        ! ocean parameters
-   USE phycst         ! physical constants (ocean directory) 
-   USE ice            ! sea-ice: variables
-   USE ice1D          ! sea-ice: thermodynamics
-   !
-   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)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_thd_dif   ! called by ice_thd
-
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
-   !! $Id: icethd_dif.F90 8420 2017-08-08 12:18:46Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_thd_dif
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE ice_thd_dif  ***
-      !! ** Purpose :
-      !!           This routine determines the time evolution of snow and sea-ice 
-      !!           temperature profiles.
-      !! ** Method  :
-      !!           This is done by solving the heat equation diffusion with
-      !!           a Neumann boundary condition at the surface and a Dirichlet one
-      !!           at the bottom. Solar radiation is partially absorbed into the ice.
-      !!           The specific heat and thermal conductivities depend on ice salinity
-      !!           and temperature to take into account brine pocket melting. The 
-      !!           numerical
-      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid 
-      !!           in the ice and snow system.
-      !!
-      !!           The successive steps of this routine are
-      !!           1.  Thermal conductivity at the interfaces of the ice layers
-      !!           2.  Internal absorbed radiation
-      !!           3.  Scale factors due to non-uniform grid
-      !!           4.  Kappa factors
-      !!           Then iterative procedure begins
-      !!           5.  specific heat in the ice
-      !!           6.  eta factors
-      !!           7.  surface flux computation
-      !!           8.  tridiagonal system terms
-      !!           9.  solving the tridiagonal system with Gauss elimination
-      !!           Iterative procedure ends according to a criterion on evolution
-      !!           of temperature
-      !!
-      !! ** Inputs / Ouputs : (global commons)
-      !!           surface temperature : t_su_1d
-      !!           ice/snow temperatures   : t_i_1d, t_s_1d
-      !!           ice salinities          : s_i_1d
-      !!           number of layers in the ice/snow: nlay_i, nlay_s
-      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
-      !!------------------------------------------------------------------
-      INTEGER ::   ji, jk         ! spatial loop index
-      INTEGER ::   numeq          ! current reference number of equation
-      INTEGER ::   minnumeqmin, maxnumeqmax
-      INTEGER ::   iconv          ! number of iterations in iterative procedure
-      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
-      
-      INTEGER, DIMENSION(jpij) ::   numeqmin   ! reference number of top equation
-      INTEGER, DIMENSION(jpij) ::   numeqmax   ! reference number of bottom equation
-      
-      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
-      REAL(wp) ::   zg1       =  2._wp        !
-      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
-      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
-      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
-      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
-      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
-      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
-      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
-      REAL(wp) ::   z1_hsu
-      REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
-      REAL(wp) ::   zcpi                      ! Ice specific heat
-      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
-      REAL(wp) ::   zfac                      ! dummy factor
-      
-      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
-      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
-      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
-      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
-      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
-      
-      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
-      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
-      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
-      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
-      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
-      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
-      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
-      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
-      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
-      
-      ! Mono-category
-      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
-      REAL(wp) :: zhe           ! dummy factor
-      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity
-      !!------------------------------------------------------------------     
-
-      ! --- diag error on heat diffusion - PART 1 --- !
-      DO ji = 1, nidx
-         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
-            &           SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
-      END DO
-
-      !------------------------------------------------------------------------------!
-      ! 1) Initialization                                                            !
-      !------------------------------------------------------------------------------!
-      DO ji = 1, nidx
-         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
-         ! layer thickness
-         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
-         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
-      END DO
-      !
-      WHERE( zh_i(1:nidx) >= epsi10 )   ;   z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)
-      ELSEWHERE                         ;   z1_h_i(1:nidx) = 0._wp
-      END WHERE
-
-      WHERE( zh_s(1:nidx) >= epsi10 )   ;   z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)
-      ELSEWHERE                         ;   z1_h_s(1:nidx) = 0._wp
-      END WHERE
-      !
-      ! temperatures
-      ztsub  (1:nidx)   = t_su_1d(1:nidx)  ! temperature at the previous iteration
-      ztsold (1:nidx,:) = t_s_1d(1:nidx,:)   ! Old snow temperature
-      ztiold (1:nidx,:) = t_i_1d(1:nidx,:)   ! Old ice temperature
-      t_su_1d(1:nidx)   = MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
-      !
-      !------------------------------------------------------------------------------|
-      ! 2) Radiation	                                                             |
-      !------------------------------------------------------------------------------|
-      !
-      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
-      DO ji = 1, nidx
-         !-------------------
-         ! Computation of i0
-         !-------------------
-         ! i0 describes the fraction of solar radiation which does not contribute
-         ! to the surface energy budget but rather penetrates inside the ice.
-         ! We assume that no radiation is transmitted through the snow
-         ! If there is no no snow
-         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
-         ! zftrice = io.qsr_ice       is below the surface 
-         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
-         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
-         zfac = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )     
-         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )
-
-         !-------------------------------------------------------
-         ! Solar radiation absorbed / transmitted at the surface
-         ! Derivative of the non solar flux
-         !-------------------------------------------------------
-         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
-         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
-         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
-         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
-      END DO
-
-      !---------------------------------------------------------
-      ! Transmission - absorption of solar radiation in the ice
-      !---------------------------------------------------------
-      zradtr_s(1:nidx,0) = zftrice(1:nidx)
-      DO jk = 1, nlay_s
-         DO ji = 1, nidx
-            !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
-            !                             ! radiation absorbed by the layer-th snow layer
-            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
-         END DO
-      END DO
-         
-      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
-      DO jk = 1, nlay_i 
-         DO ji = 1, nidx
-            !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
-            !                             ! radiation absorbed by the layer-th ice layer
-            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
-         END DO
-      END DO
-
-      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
-
-      !------------------------------------------------------------------------------|
-      !  3) Iterative procedure begins                                               |
-      !------------------------------------------------------------------------------|
-      !
-      iconv    =  0          ! number of iterations
-      zdti_max =  1000._wp   ! maximal value of error on all points
-      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )
-         !
-         iconv = iconv + 1
-         !
-         ztib(1:nidx,:) = t_i_1d(1:nidx,:)
-         ztsb(1:nidx,:) = t_s_1d(1:nidx,:)
-         !
-         !------------------------------------------------------------------------------|
-         ! 4) Sea ice thermal conductivity                                              |
-         !------------------------------------------------------------------------------|
-         !
-         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 )
-               ztcond_i(ji,nlay_i) = rcdic + zbeta * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
-            END DO
-            DO jk = 1, nlay_i-1
-               DO ji = 1, nidx
-                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
-                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
-               END DO
-            END DO
-            !
-         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.09_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
-                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  &
-                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
-            END DO
-            DO jk = 1, nlay_i-1
-               DO ji = 1, nidx
-                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /          &
-                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   &
-                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
-               END DO
-            END DO
-            !
-         ENDIF
-         ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )        
-         !
-         !------------------------------------------------------------------------------|
-         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
-         !------------------------------------------------------------------------------|
-         !
-         ! Computation of effective thermal conductivity G(h)
-         ! Used in mono-category case only to simulate an ITD implicitly
-         ! Fichefet and Morales Maqueda, JGR 1997
-
-         zghe(1:nidx) = 1._wp
-         
-         SELECT CASE ( nn_monocat )
-
-         CASE ( 1 , 3 )
-
-            zepsilon = 0.1_wp
-            DO ji = 1, nidx
-   
-               ! Mean sea ice thermal conductivity
-               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )
-
-               ! Effective thickness he (zhe)
-               zhe  = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )
-
-               ! G(he)
-               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
-                  zghe(ji) =  MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )
-               ENDIF
-   
-            END DO
-
-         END SELECT
-         !
-         !------------------------------------------------------------------------------|
-         !  6) kappa factors                                                            |
-         !------------------------------------------------------------------------------|
-         !
-         !--- Snow
-         DO jk = 0, nlay_s-1
-            DO ji = 1, nidx
-               zkappa_s(ji,jk)  = zghe(ji) * rn_cnd_s * z1_h_s(ji)
-            END DO
-         END DO
-         DO ji = 1, nidx   ! Snow-ice interface
-            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
-            IF( zfac > epsi10 ) THEN
-               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
-            ELSE
-               zkappa_s(ji,nlay_s) = 0._wp
-            ENDIF
-         END DO
-
-         !--- Ice
-         DO jk = 0, nlay_i
-            DO ji = 1, nidx
-               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
-            END DO
-         END DO
-         DO ji = 1, nidx   ! Snow-ice interface
-            zkappa_i(ji,0)     = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
-         END DO
-         !
-         !------------------------------------------------------------------------------|
-         ! 7) Sea ice specific heat, eta factors                                        |
-         !------------------------------------------------------------------------------|
-         !
-         DO jk = 1, nlay_i
-            DO ji = 1, nidx
-               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
-            END DO
-         END DO
-
-         DO jk = 1, nlay_s
-            DO ji = 1, nidx
-               zeta_s(ji,jk)  = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
-            END DO
-         END DO
-         !
-         !------------------------------------------------------------------------------|
-         ! 8) surface flux computation                                                  |
-         !------------------------------------------------------------------------------|
-         !
-         IF ( ln_dqns_i ) THEN 
-            DO ji = 1, nidx
-               ! update of the non solar flux according to the update in T_su
-               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
-            END DO
-         ENDIF
-
-         DO ji = 1, nidx
-            zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
-         END DO
-         !
-         !------------------------------------------------------------------------------|
-         ! 9) tridiagonal system terms                                                  |
-         !------------------------------------------------------------------------------|
-         !
-         !!layer denotes the number of the layer in the snow or in the ice
-         !!numeq denotes the reference number of the equation in the tridiagonal
-         !!system, terms of tridiagonal system are indexed as following :
-         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
-
-         !!ice interior terms (top equation has the same form as the others)
-
-         DO numeq=1,nlay_i+3
-            DO ji = 1, nidx
-               ztrid(ji,numeq,1) = 0.
-               ztrid(ji,numeq,2) = 0.
-               ztrid(ji,numeq,3) = 0.
-               zindterm(ji,numeq)= 0.
-               zindtbis(ji,numeq)= 0.
-               zdiagbis(ji,numeq)= 0.
-            ENDDO
-         ENDDO
-
-         DO numeq = nlay_s + 2, nlay_s + nlay_i 
-            DO ji = 1, nidx
-               jk                 = numeq - nlay_s - 1
-               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
-               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
-               zindterm(ji,numeq) =  ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-            END DO
-         ENDDO
-
-         numeq =  nlay_s + nlay_i + 1
-         DO ji = 1, nidx
-            !!ice bottom term
-            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
-            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
-            ztrid(ji,numeq,3)  =  0.0
-            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
-         ENDDO
-
-
-         DO ji = 1, nidx
-            IF ( ht_s_1d(ji) > 0.0 ) THEN
-               !
-               !------------------------------------------------------------------------------|
-               !  snow-covered cells                                                          |
-               !------------------------------------------------------------------------------|
-               !
-               !!snow interior terms (bottom equation has the same form as the others)
-               DO numeq = 3, nlay_s + 1
-                  jk                  =  numeq - 1
-                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
-                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
-                  zindterm(ji,numeq)  =  ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
-               END DO
-
-               !!case of only one layer in the ice (ice equation is altered)
-               IF ( nlay_i == 1 ) THEN
-                  ztrid(ji,nlay_s+2,3)    =  0.0
-                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
-               ENDIF
-
-               IF ( t_su_1d(ji) < rt0 ) THEN
-
-                  !------------------------------------------------------------------------------|
-                  !  case 1 : no surface melting - snow present                                  |
-                  !------------------------------------------------------------------------------|
-                  numeqmin(ji)    =  1
-                  numeqmax(ji)    =  nlay_i + nlay_s + 1
-
-                  !!surface equation
-                  ztrid(ji,1,1)  = 0.0
-                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
-                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
-                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji)
-
-                  !!first layer of snow equation
-                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
-                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
-
-               ELSE 
-                  !
-                  !------------------------------------------------------------------------------|
-                  !  case 2 : surface is melting - snow present                                  |
-                  !------------------------------------------------------------------------------|
-                  !
-                  numeqmin(ji)    =  2
-                  numeqmax(ji)    =  nlay_i + nlay_s + 1
-
-                  !!first layer of snow equation
-                  ztrid(ji,2,1)  =  0.0
-                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
-                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
-               ENDIF
-            ELSE
-               !
-               !------------------------------------------------------------------------------|
-               !  cells without snow                                                          |
-               !------------------------------------------------------------------------------|
-               !
-               IF ( t_su_1d(ji) < rt0 ) THEN
-                  !
-                  !------------------------------------------------------------------------------|
-                  !  case 3 : no surface melting - no snow                                       |
-                  !------------------------------------------------------------------------------|
-                  !
-                  numeqmin(ji)      =  nlay_s + 1
-                  numeqmax(ji)      =  nlay_i + nlay_s + 1
-
-                  !!surface equation	
-                  ztrid(ji,numeqmin(ji),1)   =  0.0
-                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
-                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji)
-
-                  !!first layer of ice equation
-                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
-                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
-                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1)  
-                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
-
-                  !!case of only one layer in the ice (surface & ice equations are altered)
-
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,numeqmin(ji),1)    =  0.0
-                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
-                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
-                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
-                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
-
-                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1) *  &
-                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
-                  ENDIF
-
-               ELSE
-
-                  !
-                  !------------------------------------------------------------------------------|
-                  ! case 4 : surface is melting - no snow                                        |
-                  !------------------------------------------------------------------------------|
-                  !
-                  numeqmin(ji)    =  nlay_s + 2
-                  numeqmax(ji)    =  nlay_i + nlay_s + 1
-
-                  !!first layer of ice equation
-                  ztrid(ji,numeqmin(ji),1)      =  0.0
-                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
-                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1) *  &
-                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
-
-                  !!case of only one layer in the ice (surface & ice equations are altered)
-                  IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,numeqmin(ji),1)  =  0.0
-                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,numeqmin(ji),3)  =  0.0
-                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
-                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
-                  ENDIF
-
-               ENDIF
-            ENDIF
-
-         END DO
-         !
-         !------------------------------------------------------------------------------|
-         ! 10) tridiagonal system solving                                               |
-         !------------------------------------------------------------------------------|
-         !
-         ! Solve the tridiagonal system with Gauss elimination method.
-         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, 
-         ! McGraw-Hill 1984.	
-
-         maxnumeqmax = 0
-         minnumeqmin = nlay_i+5
-
-         DO ji = 1, nidx
-            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
-            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
-            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
-            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
-         END DO
-
-         DO jk = minnumeqmin+1, maxnumeqmax
-            DO ji = 1, nidx
-               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
-               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
-               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
-            END DO
-         END DO
-
-         DO ji = 1, nidx
-            ! ice temperatures
-            t_i_1d(ji,nlay_i) =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
-         END DO
-
-         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
-            DO ji = 1, nidx
-               jk    =  numeq - nlay_s - 1
-               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
-            END DO
-         END DO
-
-         DO ji = 1, nidx
-            ! snow temperatures      
-            IF( ht_s_1d(ji) > 0._wp ) THEN
-               t_s_1d(ji,nlay_s) =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
-                  &                 / zdiagbis(ji,nlay_s+1)
-            ENDIF
-            ! surface temperature
-            ztsub(ji) = t_su_1d(ji)
-            IF( t_su_1d(ji) < rt0 ) THEN
-               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
-                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
-            ENDIF
-         END DO
-         !
-         !--------------------------------------------------------------------------
-         !  11) Has the scheme converged ?, end of the iterative procedure         |
-         !--------------------------------------------------------------------------
-         !
-         ! check that nowhere it has started to melt
-         ! zdti_max is a measure of error, it has to be under zdti_bnd
-         zdti_max = 0._wp
-         DO ji = 1, nidx
-            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
-            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
-         END DO
-
-         DO jk  =  1, nlay_s
-            DO ji = 1, nidx
-               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
-               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
-            END DO
-         END DO
-
-         DO jk  =  1, nlay_i
-            DO ji = 1, nidx
-               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
-               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
-               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
-            END DO
-         END DO
-
-         ! Compute spatial maximum over all errors
-         ! note that this could be optimized substantially by iterating only the non-converging points
-         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
-
-      END DO  ! End of the do while iterative procedure
-
-      IF( ln_icectl .AND. lwp ) THEN
-         WRITE(numout,*) ' zdti_max : ', zdti_max
-         WRITE(numout,*) ' iconv    : ', iconv
-      ENDIF
-
-      !
-      !-------------------------------------------------------------------------!
-      !   12) Fluxes at the interfaces                                          !
-      !-------------------------------------------------------------------------!
-      DO ji = 1, nidx
-         !                                ! surface ice conduction flux
-         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
-            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
-         !                                ! bottom ice conduction flux
-         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
-      END DO
-
-      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
-      CALL ice_thd_enmelt
-
-      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
-      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) 
-         END DO
-      END IF
-
-      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
-      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
-      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-      DO ji = 1, nidx
-         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
-            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
-
-         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
-         ELSE                          ! case T_su = 0degC
-            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
-         ENDIF
-         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
-
-         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
-         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
-
-      END DO   
-
-      ! --- Diagnostics SIMIP --- !
-      DO ji = 1, nidx
-         !--- Conduction fluxes (positive downwards)
-         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
-         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
-
-         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
-         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
-         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
-            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
-            t_si_1d(ji) = t_su_1d(ji)
-         ENDIF
-      END DO
-      !
-   END SUBROUTINE ice_thd_dif
-
-
-   SUBROUTINE ice_thd_enmelt
-      !!-----------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_thd_enmelt *** 
-      !!                 
-      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
-      !!
-      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts  ! local scalar 
-      !!-------------------------------------------------------------------
-      !
-      DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = 1, nidx
-            ztmelts      = - tmut  * s_i_1d(ji,jk)
-            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
-                                                                !   (sometimes dif scheme produces abnormally high temperatures)   
-            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
-               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
-               &                    - rcp  *   ztmelts )
-         END DO
-      END DO
-      DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = 1, nidx
-            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
-         END DO
-      END DO
-      !
-   END SUBROUTINE ice_thd_enmelt
-
-#else
-   !!----------------------------------------------------------------------
-   !!                   Dummy Module                 No LIM-3 sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE icethd_dif
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8531)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8531)
@@ -0,0 +1,541 @@
+MODULE icethd_do
+   !!======================================================================
+   !!                       ***  MODULE icethd_do   ***
+   !!                lateral thermodynamic growth of the ice 
+   !!======================================================================
+   !! History :  LIM  ! 2005-12 (M. Vancoppenolle)  Original code
+   !!             -   ! 2006-01 (M. Vancoppenolle)  add ITD
+   !!            3.0  ! 2007-07 (M. Vancoppenolle)  Mass and energy conservation tested
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                       LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_thd_do   : ice growth in open water (=lateral accretion of ice)
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
+   USE dom_oce        ! domain variables
+   USE phycst         ! physical constants
+   USE sbc_oce , ONLY : sss_m
+   USE sbc_ice , ONLY : utau_ice, vtau_ice
+   USE ice1D          ! sea-ice: thermodynamics
+   USE ice            ! sea-ice: variables
+   USE icetab         ! sea-ice: 2D <==> 1D
+   USE icectl         ! sea-ice: conservation
+   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
+   USE lib_mpp        ! MPP library
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_thd_do        ! called by ice_thd
+   PUBLIC   ice_thd_do_init   ! called by ice_stp
+
+   ! ** namelist (namthd_do) **
+   REAL(wp) ::   rn_hinew         ! thickness for new ice formation (m)
+   LOGICAL  ::   ln_frazil        ! use of frazil ice collection as function of wind (T) or not (F)
+   REAL(wp) ::   rn_maxfraz       ! maximum portion of frazil ice collecting at the ice bottom
+   REAL(wp) ::   rn_vfraz         ! threshold drift speed for collection of bottom frazil ice
+   REAL(wp) ::   rn_Cfraz         ! squeezing coefficient for collection of bottom frazil ice
+
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
+   !! $Id: icethd_do.F90 8420 2017-08-08 12:18:46Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_thd_do
+      !!-------------------------------------------------------------------
+      !!               ***   ROUTINE ice_thd_do  ***
+      !!  
+      !! ** Purpose : Computation of the evolution of the ice thickness and 
+      !!      concentration as a function of the heat balance in the leads.
+      !!      It is only used for lateral accretion
+      !!       
+      !! ** Method  : Ice is formed in the open water when ocean lose heat
+      !!      (heat budget of open water Bl is negative) .
+      !!      Computation of the increase of 1-A (ice concentration) fol-
+      !!      lowing the law :
+      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
+      !!       where - h0 is the thickness of ice created in the lead
+      !!             - a is a minimum fraction for leads
+      !!             - F is a monotonic non-increasing function defined as:
+      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
+      !!             - exld is the exponent closure rate (=2 default val.)
+      !! 
+      !! ** Action : - Adjustment of snow and ice thicknesses and heat
+      !!                content in brine pockets
+      !!             - Updating ice internal temperature
+      !!             - Computation of variation of ice volume and mass
+      !!             - Computation of a_i after lateral accretion and 
+      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
+      !!------------------------------------------------------------------------
+      INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices
+      INTEGER  ::   iter     !   -       -
+      REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars
+      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
+      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
+
+      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
+      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
+      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
+      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
+     
+      REAL(wp) ::   zv_newfra
+  
+      INTEGER , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows
+      REAL(wp), DIMENSION(jpij) ::   zswinew     ! switch for new ice or not
+
+      REAL(wp), DIMENSION(jpij) ::   zv_newice   ! volume of accreted ice
+      REAL(wp), DIMENSION(jpij) ::   za_newice   ! fractional area of accreted ice
+      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! thickness of accreted ice
+      REAL(wp), DIMENSION(jpij) ::   ze_newice   ! heat content of accreted ice
+      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of accreted ice
+      REAL(wp), DIMENSION(jpij) ::   zo_newice   ! age of accreted ice
+      REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
+      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
+      REAL(wp), DIMENSION(jpij) ::   zat_i_1d    ! total ice fraction    
+      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
+      REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
+
+      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b      ! old volume of ice in category jl
+      REAL(wp), DIMENSION(jpij,jpl) ::   za_b      ! old area of ice in category jl
+      REAL(wp), DIMENSION(jpij,jpl) ::   za_i_1d   ! 1-D version of a_i
+      REAL(wp), DIMENSION(jpij,jpl) ::   zv_i_1d   ! 1-D version of v_i
+      REAL(wp), DIMENSION(jpij,jpl) ::   zsmv_i_1d ! 1-D version of smv_i
+
+      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_1d !: 1-D version of e_i
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   zvrel     ! relative ice / frazil velocity
+
+      REAL(wp) :: zcai = 1.4e-3_wp                     ! ice-air drag (clem: should be dependent on coupling/forcing used)
+      !!-----------------------------------------------------------------------!
+
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+
+      CALL ice_var_agg(1)
+      CALL ice_var_glo2eqv
+
+      !------------------------------------------------------------------------------!
+      ! 3) Collection thickness of ice formed in leads and polynyas
+      !------------------------------------------------------------------------------!    
+      ! hicol is the thickness of new ice formed in open water
+      ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1)
+      ! Frazil ice forms in open water, is transported by wind
+      ! accumulates at the edge of the consolidated ice edge
+      ! where it forms aggregates of a specific thickness called
+      ! collection thickness.
+
+      ! Note : the following algorithm currently breaks vectorization
+      ! 
+
+      zvrel(:,:) = 0._wp
+
+      ! Default new ice thickness
+      WHERE( qlead(:,:) < 0._wp )   ;   hicol(:,:) = rn_hinew
+      ELSEWHERE                     ;   hicol(:,:) = 0._wp
+      END WHERE
+
+      IF( ln_frazil ) THEN
+
+         !--------------------
+         ! Physical constants
+         !--------------------
+         hicol(:,:) = 0._wp
+
+         zhicrit = 0.04 ! frazil ice thickness
+         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
+         zsqcd   = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)
+         zgamafr = 0.03
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
+                  !-------------
+                  ! Wind stress
+                  !-------------
+                  ! C-grid wind stress components
+                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
+                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
+                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
+                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
+                  ! Square root of wind stress
+                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
+
+                  !---------------------
+                  ! Frazil ice velocity
+                  !---------------------
+                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
+                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
+                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
+
+                  !-------------------
+                  ! Pack ice velocity
+                  !-------------------
+                  ! C-grid ice velocity
+                  zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
+                  zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
+
+                  !-----------------------------------
+                  ! Relative frazil/pack ice velocity
+                  !-----------------------------------
+                  ! absolute relative velocity
+                  rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
+                  zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
+                     &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
+                  zvrel(ji,jj) = SQRT( zvrel2 )
+
+                  !---------------------
+                  ! Iterative procedure
+                  !---------------------
+                  hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
+                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
+
+                  iter = 1
+                  DO WHILE ( iter < 20 ) 
+                     zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   &
+                        &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2
+                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
+
+                     hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 )
+                     iter = iter + 1
+                  END DO
+
+               ENDIF ! end of selection of pixels where ice forms
+
+            END DO 
+         END DO 
+         ! 
+         CALL lbc_lnk_multi( zvrel, 'T', 1., hicol, 'T', 1.  )
+
+      ENDIF ! End of computation of frazil ice collection thickness
+
+      !------------------------------------------------------------------------------!
+      ! 4) Identify grid points where new ice forms
+      !------------------------------------------------------------------------------!
+
+      !-------------------------------------
+      ! Select points for new ice formation
+      !-------------------------------------
+      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
+      nidx = 0 ; idxice(:) = 0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
+               nidx = nidx + 1
+               idxice( nidx ) = (jj - 1) * jpi + ji
+            ENDIF
+         END DO
+      END DO
+
+      !------------------------------
+      ! Move from 2-D to 1-D vectors
+      !------------------------------
+      ! If ocean gains heat do nothing. Otherwise compute new ice formation
+
+      IF ( nidx > 0 ) THEN
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zat_i_1d  (1:nidx)     , at_i           )
+         DO jl = 1, jpl
+            CALL tab_2d_1d( nidx, idxice(1:nidx), za_i_1d  (1:nidx,jl), a_i  (:,:,jl)  )
+            CALL tab_2d_1d( nidx, idxice(1:nidx), zv_i_1d  (1:nidx,jl), v_i  (:,:,jl)  )
+            CALL tab_2d_1d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl)  )
+            DO jk = 1, nlay_i
+               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)   )
+            END DO
+         END DO
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx)     , qlead       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx)     , t_bo        )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx)     , sfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx)     , wfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hicol_1d  (1:nidx)     , hicol       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx)     , zvrel       )
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx)     , hfx_thd     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx)     , hfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx)     , rn_amax_2d  )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx)     , sss_m       )
+
+         !------------------------------------------------------------------------------|
+         ! 2) Convert units for ice internal energy
+         !------------------------------------------------------------------------------|
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               DO ji = 1, nidx
+                  IF( zv_i_1d(ji,jl) > 0._wp )   ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i )
+               END DO
+            END DO
+         END DO
+         !------------------------------------------------------------------------------!
+         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
+         !------------------------------------------------------------------------------!
+
+         !-----------------------------------------
+         ! Keep old ice areas and volume in memory
+         !-----------------------------------------
+         zv_b(1:nidx,:) = zv_i_1d(1:nidx,:) 
+         za_b(1:nidx,:) = za_i_1d(1:nidx,:)
+
+         !----------------------
+         ! Thickness of new ice
+         !----------------------
+         zh_newice(1:nidx) = hicol_1d(1:nidx)
+
+         !----------------------
+         ! Salinity of new ice 
+         !----------------------
+         SELECT CASE ( nn_icesal )
+         CASE ( 1 )                    ! Sice = constant 
+            zs_newice(1:nidx) = rn_icesal
+         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
+            DO ji = 1, nidx
+               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
+            END DO
+         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
+            zs_newice(1:nidx) =   2.3
+         END SELECT
+
+         !-------------------------
+         ! Heat content of new ice
+         !-------------------------
+         ! We assume that new ice is formed at the seawater freezing point
+         DO ji = 1, nidx
+            ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K)
+            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                                         &
+               &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+               &                       - rcp  *         ( ztmelts - rt0 )  )
+         END DO
+
+         !----------------
+         ! Age of new ice
+         !----------------
+         DO ji = 1, nidx
+            zo_newice(ji) = 0._wp
+         END DO
+
+         !-------------------
+         ! Volume of new ice
+         !-------------------
+         DO ji = 1, nidx
+
+            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
+
+            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
+                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 
+                                                                   
+            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
+                                              
+            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0) 
+                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
+            zv_newice(ji) = - zfmdt * r1_rhoic
+
+            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
+
+            ! Contribution to heat flux to the ocean [W.m-2], >0  
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
+            ! Total heat flux used in this process [W.m-2]  
+            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
+            ! mass flux
+            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice
+            ! salt flux
+            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
+         END DO
+         
+         zv_frazb(:) = 0._wp
+         IF( ln_frazil ) THEN
+            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
+            DO ji = 1, nidx
+               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+               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)
+            END DO
+         END IF
+         
+         !-----------------
+         ! Area of new ice
+         !-----------------
+         DO ji = 1, nidx
+            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
+         END DO
+
+         !------------------------------------------------------------------------------!
+         ! 6) Redistribute new ice area and volume into ice categories                  !
+         !------------------------------------------------------------------------------!
+
+         !------------------------
+         ! 6.1) lateral ice growth
+         !------------------------
+         ! If lateral ice growth gives an ice concentration gt 1, then
+         ! we keep the excessive volume in memory and attribute it later to bottom accretion
+         DO ji = 1, nidx
+            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
+               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
+               za_newice(ji) = za_newice(ji) - zda_res  (ji)
+               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
+            ELSE
+               zda_res(ji) = 0._wp
+               zdv_res(ji) = 0._wp
+            ENDIF
+         END DO
+
+         ! find which category to fill
+         zat_i_1d(:) = 0._wp
+         DO jl = 1, jpl
+            DO ji = 1, nidx
+               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
+                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
+                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
+                  jcat    (ji)    = jl
+               ENDIF
+               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
+            END DO
+         END DO
+
+         ! Heat content
+         DO ji = 1, nidx
+            jl = jcat(ji)                                                    ! categroy in which new ice is put
+            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
+         END DO
+
+         DO jk = 1, nlay_i
+            DO ji = 1, nidx
+               jl = jcat(ji)
+               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
+               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
+                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
+                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
+            END DO
+         END DO
+
+         !------------------------------------------------
+         ! 6.2) bottom ice growth + ice enthalpy remapping
+         !------------------------------------------------
+         DO jl = 1, jpl
+
+            ! for remapping
+            h_i_old (1:nidx,0:nlay_i+1) = 0._wp
+            eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
+            DO jk = 1, nlay_i
+               DO ji = 1, nidx
+                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
+                  eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
+               END DO
+            END DO
+
+            ! new volumes including lateral/bottom accretion + residual
+            DO ji = 1, nidx
+               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
+               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
+               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)               
+               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
+               ! for remapping
+               h_i_old (ji,nlay_i+1) = zv_newfra
+               eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
+            ENDDO
+            ! --- Ice enthalpy remapping --- !
+            CALL ice_thd_ent( ze_i_1d(1:nidx,:,jl) ) 
+         ENDDO
+
+         !-----------------
+         ! Update salinity
+         !-----------------
+         DO jl = 1, jpl
+            DO ji = 1, nidx
+               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
+               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
+            END DO
+         END DO
+
+         !------------------------------------------------------------------------------!
+         ! 8) Change units for e_i
+         !------------------------------------------------------------------------------!    
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               DO ji = 1, nidx
+                  ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i 
+               END DO
+            END DO
+         END DO
+         !------------------------------------------------------------------------------!
+         ! 7) Change 2D vectors to 1D vectors 
+         !------------------------------------------------------------------------------!
+         DO jl = 1, jpl
+            CALL tab_1d_2d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl)  )
+            CALL tab_1d_2d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl)  )
+            CALL tab_1d_2d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i (:,:,jl)   )
+            DO jk = 1, nlay_i
+               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)  )
+            END DO
+         END DO
+         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw  )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw  )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd  )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw  )
+         !
+      ENDIF ! nidx > 0
+      !
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      !
+   END SUBROUTINE ice_thd_do
+
+   SUBROUTINE ice_thd_do_init
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_thd_do_init *** 
+      !!                 
+      !! ** Purpose :   Physical constants and parameters associated with
+      !!                ice thermodynamics
+      !!
+      !! ** Method  :   Read the namthd_do namelist and check the parameters
+      !!                called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namthd_do
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ios   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namthd_do/ rn_hinew, ln_frazil, rn_maxfraz, rn_vfraz, rn_Cfraz
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namthd_do in reference namelist : Ice thermodynamics
+      READ  ( numnam_ice_ref, namthd_do, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_do in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namthd_do in configuration namelist : Ice thermodynamics
+      READ  ( numnam_ice_cfg, namthd_do, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_do in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd_do )
+      !
+      !
+      IF(lwp) THEN                          ! control print
+         WRITE(numout,*) 'ice_thd_do_init: Ice growth in open water'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namthd_do:'
+         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_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
+      ENDIF
+      !
+      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_do_init : rn_hinew should be >= rn_himin' )
+      !
+   END SUBROUTINE ice_thd_do_init
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                               NO  LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE icethd_do
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8531)
@@ -33,5 +33,5 @@
    PRIVATE
 
-   PUBLIC   ice_thd_ent         ! called by icethd and icethd_lac
+   PUBLIC   ice_thd_ent         ! called by icethd and icethd_do
 
    !!----------------------------------------------------------------------
@@ -130,5 +130,5 @@
 
       ! --- diag error on heat remapping --- !
-      ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_lac), 
+      ! comment: if input h_i_old and eh_i_old are already multiplied by a_i (as in icethd_do), 
       ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0
       DO ji = 1, nidx
Index: anches/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 8530)
+++ 	(revision )
@@ -1,491 +1,0 @@
-MODULE icethd_lac
-   !!======================================================================
-   !!                       ***  MODULE icethd_lac   ***
-   !!                lateral thermodynamic growth of the ice 
-   !!======================================================================
-   !! History :  LIM  ! 2005-12 (M. Vancoppenolle)  Original code
-   !!             -   ! 2006-01 (M. Vancoppenolle)  add ITD
-   !!            3.0  ! 2007-07 (M. Vancoppenolle)  Mass and energy conservation tested
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                       LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_lat_acr   : lateral accretion of ice
-   !!----------------------------------------------------------------------
-   USE par_oce        ! ocean parameters
-   USE dom_oce        ! domain variables
-   USE phycst         ! physical constants
-   USE sbc_oce , ONLY : sss_m
-   USE sbc_ice , ONLY : utau_ice, vtau_ice
-   USE ice1D          ! sea-ice: thermodynamics
-   USE ice            ! sea-ice: variables
-   USE icetab         ! sea-ice: 2D <==> 1D
-   USE icectl         ! sea-ice: conservation
-   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
-   USE lib_mpp        ! MPP library
-   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_thd_lac   ! called by ice_thd
-
-   !!----------------------------------------------------------------------
-   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
-   !! $Id: icethd_lac.F90 8420 2017-08-08 12:18:46Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE ice_thd_lac
-      !!-------------------------------------------------------------------
-      !!               ***   ROUTINE ice_thd_lac  ***
-      !!  
-      !! ** Purpose : Computation of the evolution of the ice thickness and 
-      !!      concentration as a function of the heat balance in the leads.
-      !!      It is only used for lateral accretion
-      !!       
-      !! ** Method  : Ice is formed in the open water when ocean lose heat
-      !!      (heat budget of open water Bl is negative) .
-      !!      Computation of the increase of 1-A (ice concentration) fol-
-      !!      lowing the law :
-      !!      (dA/dt)acc = F[ (1-A)/(1-a) ] * [ Bl / (Li*h0) ]
-      !!       where - h0 is the thickness of ice created in the lead
-      !!             - a is a minimum fraction for leads
-      !!             - F is a monotonic non-increasing function defined as:
-      !!                  F(X)=( 1 - X**exld )**(1.0/exld)
-      !!             - exld is the exponent closure rate (=2 default val.)
-      !! 
-      !! ** Action : - Adjustment of snow and ice thicknesses and heat
-      !!                content in brine pockets
-      !!             - Updating ice internal temperature
-      !!             - Computation of variation of ice volume and mass
-      !!             - Computation of a_i after lateral accretion and 
-      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
-      !!------------------------------------------------------------------------
-      INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices
-      INTEGER  ::   iter     !   -       -
-      REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars
-      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
-      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
-
-      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
-      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
-      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
-      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
-     
-      REAL(wp) ::   zv_newfra
-  
-      INTEGER , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows
-      REAL(wp), DIMENSION(jpij) ::   zswinew     ! switch for new ice or not
-
-      REAL(wp), DIMENSION(jpij) ::   zv_newice   ! volume of accreted ice
-      REAL(wp), DIMENSION(jpij) ::   za_newice   ! fractional area of accreted ice
-      REAL(wp), DIMENSION(jpij) ::   zh_newice   ! thickness of accreted ice
-      REAL(wp), DIMENSION(jpij) ::   ze_newice   ! heat content of accreted ice
-      REAL(wp), DIMENSION(jpij) ::   zs_newice   ! salinity of accreted ice
-      REAL(wp), DIMENSION(jpij) ::   zo_newice   ! age of accreted ice
-      REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
-      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
-      REAL(wp), DIMENSION(jpij) ::   zat_i_1d    ! total ice fraction    
-      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
-      REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
-
-      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b      ! old volume of ice in category jl
-      REAL(wp), DIMENSION(jpij,jpl) ::   za_b      ! old area of ice in category jl
-      REAL(wp), DIMENSION(jpij,jpl) ::   za_i_1d   ! 1-D version of a_i
-      REAL(wp), DIMENSION(jpij,jpl) ::   zv_i_1d   ! 1-D version of v_i
-      REAL(wp), DIMENSION(jpij,jpl) ::   zsmv_i_1d ! 1-D version of smv_i
-
-      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_1d !: 1-D version of e_i
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   zvrel     ! relative ice / frazil velocity
-
-      REAL(wp) :: zcai = 1.4e-3_wp                     ! ice-air drag (clem: should be dependent on coupling/forcing used)
-      !!-----------------------------------------------------------------------!
-
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-
-      CALL ice_var_agg(1)
-      CALL ice_var_glo2eqv
-
-      !------------------------------------------------------------------------------!
-      ! 3) Collection thickness of ice formed in leads and polynyas
-      !------------------------------------------------------------------------------!    
-      ! hicol is the thickness of new ice formed in open water
-      ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1)
-      ! Frazil ice forms in open water, is transported by wind
-      ! accumulates at the edge of the consolidated ice edge
-      ! where it forms aggregates of a specific thickness called
-      ! collection thickness.
-
-      ! Note : the following algorithm currently breaks vectorization
-      ! 
-
-      zvrel(:,:) = 0._wp
-
-      ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp )   ;   hicol(:,:) = rn_hinew
-      ELSEWHERE                     ;   hicol(:,:) = 0._wp
-      END WHERE
-
-      IF( ln_frazil ) THEN
-
-         !--------------------
-         ! Physical constants
-         !--------------------
-         hicol(:,:) = 0._wp
-
-         zhicrit = 0.04 ! frazil ice thickness
-         ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
-         zsqcd   = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)
-         zgamafr = 0.03
-
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-               IF ( qlead(ji,jj) < 0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN ! activated if cooling and no landfast
-                  !-------------
-                  ! Wind stress
-                  !-------------
-                  ! C-grid wind stress components
-                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
-                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
-                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
-                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
-                  ! Square root of wind stress
-                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
-
-                  !---------------------
-                  ! Frazil ice velocity
-                  !---------------------
-                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
-                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
-                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
-
-                  !-------------------
-                  ! Pack ice velocity
-                  !-------------------
-                  ! C-grid ice velocity
-                  zvgx    = ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
-                  zvgy    = ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
-
-                  !-----------------------------------
-                  ! Relative frazil/pack ice velocity
-                  !-----------------------------------
-                  ! absolute relative velocity
-                  rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
-                  zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
-                     &               + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) * rswitch
-                  zvrel(ji,jj) = SQRT( zvrel2 )
-
-                  !---------------------
-                  ! Iterative procedure
-                  !---------------------
-                  hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
-                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
-
-                  iter = 1
-                  DO WHILE ( iter < 20 ) 
-                     zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   &
-                        &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2
-                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
-
-                     hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 )
-                     iter = iter + 1
-                  END DO
-
-               ENDIF ! end of selection of pixels where ice forms
-
-            END DO 
-         END DO 
-         ! 
-         CALL lbc_lnk_multi( zvrel, 'T', 1., hicol, 'T', 1.  )
-
-      ENDIF ! End of computation of frazil ice collection thickness
-
-      !------------------------------------------------------------------------------!
-      ! 4) Identify grid points where new ice forms
-      !------------------------------------------------------------------------------!
-
-      !-------------------------------------
-      ! Select points for new ice formation
-      !-------------------------------------
-      ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
-      nidx = 0 ; idxice(:) = 0
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF ( qlead(ji,jj)  <  0._wp .AND. tau_icebfr(ji,jj) == 0._wp ) THEN
-               nidx = nidx + 1
-               idxice( nidx ) = (jj - 1) * jpi + ji
-            ENDIF
-         END DO
-      END DO
-
-      !------------------------------
-      ! Move from 2-D to 1-D vectors
-      !------------------------------
-      ! If ocean gains heat do nothing. Otherwise compute new ice formation
-
-      IF ( nidx > 0 ) THEN
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zat_i_1d  (1:nidx)     , at_i           )
-         DO jl = 1, jpl
-            CALL tab_2d_1d( nidx, idxice(1:nidx), za_i_1d  (1:nidx,jl), a_i  (:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), zv_i_1d  (1:nidx,jl), v_i  (:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl)  )
-            DO jk = 1, nlay_i
-               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)   )
-            END DO
-         END DO
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx)     , qlead       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx)     , t_bo        )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx)     , sfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx)     , wfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hicol_1d  (1:nidx)     , hicol       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx)     , zvrel       )
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx)     , hfx_thd     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx)     , hfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx)     , rn_amax_2d  )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx)     , sss_m       )
-
-         !------------------------------------------------------------------------------|
-         ! 2) Convert units for ice internal energy
-         !------------------------------------------------------------------------------|
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  IF( zv_i_1d(ji,jl) > 0._wp )   ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i )
-               END DO
-            END DO
-         END DO
-         !------------------------------------------------------------------------------!
-         ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice
-         !------------------------------------------------------------------------------!
-
-         !-----------------------------------------
-         ! Keep old ice areas and volume in memory
-         !-----------------------------------------
-         zv_b(1:nidx,:) = zv_i_1d(1:nidx,:) 
-         za_b(1:nidx,:) = za_i_1d(1:nidx,:)
-
-         !----------------------
-         ! Thickness of new ice
-         !----------------------
-         zh_newice(1:nidx) = hicol_1d(1:nidx)
-
-         !----------------------
-         ! Salinity of new ice 
-         !----------------------
-         SELECT CASE ( nn_icesal )
-         CASE ( 1 )                    ! Sice = constant 
-            zs_newice(1:nidx) = rn_icesal
-         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
-            DO ji = 1, nidx
-               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_1d(ji) )
-            END DO
-         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(1:nidx) =   2.3
-         END SELECT
-
-         !-------------------------
-         ! Heat content of new ice
-         !-------------------------
-         ! We assume that new ice is formed at the seawater freezing point
-         DO ji = 1, nidx
-            ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                                         &
-               &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
-               &                       - rcp  *         ( ztmelts - rt0 )  )
-         END DO
-
-         !----------------
-         ! Age of new ice
-         !----------------
-         DO ji = 1, nidx
-            zo_newice(ji) = 0._wp
-         END DO
-
-         !-------------------
-         ! Volume of new ice
-         !-------------------
-         DO ji = 1, nidx
-
-            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
-
-            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
-                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 
-                                                                   
-            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
-                                              
-            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0) 
-                                                                   ! clem: we use qlead instead of zqld (icethd) because we suppose we are at the freezing point   
-            zv_newice(ji) = - zfmdt * r1_rhoic
-
-            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
-
-            ! Contribution to heat flux to the ocean [W.m-2], >0  
-            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
-            ! Total heat flux used in this process [W.m-2]  
-            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
-            ! mass flux
-            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice
-            ! salt flux
-            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
-         END DO
-         
-         zv_frazb(:) = 0._wp
-         IF( ln_frazil ) THEN
-            ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            DO ji = 1, nidx
-               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
-               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)
-            END DO
-         END IF
-         
-         !-----------------
-         ! Area of new ice
-         !-----------------
-         DO ji = 1, nidx
-            za_newice(ji) = zv_newice(ji) / zh_newice(ji)
-         END DO
-
-         !------------------------------------------------------------------------------!
-         ! 6) Redistribute new ice area and volume into ice categories                  !
-         !------------------------------------------------------------------------------!
-
-         !------------------------
-         ! 6.1) lateral ice growth
-         !------------------------
-         ! If lateral ice growth gives an ice concentration gt 1, then
-         ! we keep the excessive volume in memory and attribute it later to bottom accretion
-         DO ji = 1, nidx
-            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
-               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
-               za_newice(ji) = za_newice(ji) - zda_res  (ji)
-               zv_newice(ji) = zv_newice(ji) - zdv_res  (ji)
-            ELSE
-               zda_res(ji) = 0._wp
-               zdv_res(ji) = 0._wp
-            ENDIF
-         END DO
-
-         ! find which category to fill
-         zat_i_1d(:) = 0._wp
-         DO jl = 1, jpl
-            DO ji = 1, nidx
-               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
-                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
-                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
-                  jcat    (ji)    = jl
-               ENDIF
-               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
-            END DO
-         END DO
-
-         ! Heat content
-         DO ji = 1, nidx
-            jl = jcat(ji)                                                    ! categroy in which new ice is put
-            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
-         END DO
-
-         DO jk = 1, nlay_i
-            DO ji = 1, nidx
-               jl = jcat(ji)
-               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
-               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
-                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
-                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
-            END DO
-         END DO
-
-         !------------------------------------------------
-         ! 6.2) bottom ice growth + ice enthalpy remapping
-         !------------------------------------------------
-         DO jl = 1, jpl
-
-            ! for remapping
-            h_i_old (1:nidx,0:nlay_i+1) = 0._wp
-            eh_i_old(1:nidx,0:nlay_i+1) = 0._wp
-            DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
-                  eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
-               END DO
-            END DO
-
-            ! new volumes including lateral/bottom accretion + residual
-            DO ji = 1, nidx
-               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
-               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
-               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)               
-               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
-               ! for remapping
-               h_i_old (ji,nlay_i+1) = zv_newfra
-               eh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
-            ENDDO
-            ! --- Ice enthalpy remapping --- !
-            CALL ice_thd_ent( ze_i_1d(1:nidx,:,jl) ) 
-         ENDDO
-
-         !-----------------
-         ! Update salinity
-         !-----------------
-         DO jl = 1, jpl
-            DO ji = 1, nidx
-               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
-               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
-            END DO
-         END DO
-
-         !------------------------------------------------------------------------------!
-         ! 8) Change units for e_i
-         !------------------------------------------------------------------------------!    
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i 
-               END DO
-            END DO
-         END DO
-         !------------------------------------------------------------------------------!
-         ! 7) Change 2D vectors to 1D vectors 
-         !------------------------------------------------------------------------------!
-         DO jl = 1, jpl
-            CALL tab_1d_2d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl)  )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl)  )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i (:,:,jl)   )
-            DO jk = 1, nlay_i
-               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)  )
-            END DO
-         END DO
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw  )
-         !
-      ENDIF ! nidx > 0
-      !
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
-      !
-   END SUBROUTINE ice_thd_lac
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option                               NO  LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-#endif
-
-   !!======================================================================
-END MODULE icethd_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 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8531)
@@ -27,13 +27,12 @@
    PRIVATE
 
-   PUBLIC   ice_thd_sal        ! called by icethd module
+   PUBLIC   ice_thd_sal        ! called by icethd
    PUBLIC   ice_thd_sal_init   ! called by ice_init
    
-   ! ** namelist (namsal) **
-   LOGICAL  ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
-   REAL(wp) ::   rn_sal_gd        !    restoring salinity for gravity drainage [PSU]
-   REAL(wp) ::   rn_time_gd       !    restoring time constant for gravity drainage (= 20 days) [s]
-   REAL(wp) ::   rn_sal_fl        !    restoring salinity for flushing [PSU]
-   REAL(wp) ::   rn_time_fl       !    restoring time constant for gravity drainage (= 10 days) [s]
+   ! ** namelist (namthd_sal) **
+   REAL(wp) ::   rn_sal_gd        ! restoring salinity for gravity drainage [PSU]
+   REAL(wp) ::   rn_time_gd       ! restoring time constant for gravity drainage (= 20 days) [s]
+   REAL(wp) ::   rn_sal_fl        ! restoring salinity for flushing [PSU]
+   REAL(wp) ::   rn_time_fl       ! restoring time constant for gravity drainage (= 10 days) [s]
 
    !!----------------------------------------------------------------------
@@ -44,5 +43,5 @@
 CONTAINS
 
-   SUBROUTINE ice_thd_sal
+   SUBROUTINE ice_thd_sal( ld_sal )
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE ice_thd_sal  ***    
@@ -55,4 +54,5 @@
       !!               -> nn_icesal = 3 -> Sice = S(z)   [multiyear ice]
       !!---------------------------------------------------------------------
+      LOGICAL, INTENT(in) ::   ld_sal            ! gravity drainage and flushing or not 
       INTEGER  ::   ji, jk                       ! dummy loop indices 
       REAL(wp) ::   iflush, igravdr              ! local scalars
@@ -78,5 +78,5 @@
             sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si
 
-            IF( ln_icedS ) THEN
+            IF( ld_sal ) THEN
                !---------------------------------------------------------
                !  Update ice salinity from brine drainage and flushing
@@ -114,23 +114,23 @@
       !! ** Purpose :   initialization of ice salinity parameters
       !!
-      !! ** Method  :   Read the namice_sal namelist and check the parameter
-      !!              values called at the first timestep (nit000)
+      !! ** Method  :   Read the namthd_sal namelist and check the parameter
+      !!                values called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namice_sal
+      !! ** input   :   Namelist namthd_sal
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_sal/ ln_icedS , nn_icesal , rn_icesal, rn_sal_gd, rn_time_gd,   &
+      NAMELIST/namthd_sal/ nn_icesal, rn_icesal, rn_sal_gd, rn_time_gd,   &
          &                 rn_sal_fl, rn_time_fl, rn_simax , rn_simin 
       !!-------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )              ! Namelist namice_sal in reference namelist : Ice salinity
-      READ  ( numnam_ice_ref, namice_sal, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_sal in reference namelist', lwp )
+      REWIND( numnam_ice_ref )              ! Namelist namthd_sal in reference namelist : Ice salinity
+      READ  ( numnam_ice_ref, namthd_sal, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_sal in reference namelist', lwp )
       !
-      REWIND( numnam_ice_cfg )              ! Namelist namice_sal in configuration namelist : Ice salinity
-      READ  ( numnam_ice_cfg, namice_sal, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_sal in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namice_sal )
+      REWIND( numnam_ice_cfg )              ! Namelist namthd_sal in configuration namelist : Ice salinity
+      READ  ( numnam_ice_cfg, namthd_sal, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_sal in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd_sal )
       !
       IF(lwp) THEN                           ! control print
@@ -138,6 +138,5 @@
          WRITE(numout,*) 'ice_thd_sal_init : Ice parameters for salinity '
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namice_sal'
-         WRITE(numout,*) '      activate gravity drainage and flushing (T) or not (F)   ln_icedS   = ', ln_icedS
+         WRITE(numout,*) '   Namelist namthd_sal:'
          WRITE(numout,*) '      switch for salinity                                     nn_icesal  = ', nn_icesal
          WRITE(numout,*) '      bulk salinity value if nn_icesal = 1                    rn_icesal  = ', rn_icesal
@@ -148,9 +147,4 @@
          WRITE(numout,*) '      Maximum tolerated ice salinity                          rn_simax   = ', rn_simax
          WRITE(numout,*) '      Minimum tolerated ice salinity                          rn_simin   = ', rn_simin
-      ENDIF
-      !
-      IF( ln_icedS .AND. nn_icesal == 1 ) THEN
-         ln_icedS = .FALSE.
-         CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)')
       ENDIF
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8531)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8531)
@@ -0,0 +1,778 @@
+MODULE icethd_zdf
+   !!======================================================================
+   !!                       ***  MODULE icethd_zdf ***
+   !!                       heat diffusion in sea ice 
+   !!                   computation of surface and inner T  
+   !!======================================================================
+   !! History :  LIM  ! 02-2003 (M. Vancoppenolle) original 1D code
+   !!                 ! 06-2005 (M. Vancoppenolle) 3d version
+   !!                 ! 11-2006 (X Fettweis) Vectorization by Xavier
+   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation
+   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
+   USE phycst         ! physical constants (ocean directory) 
+   USE ice            ! sea-ice: variables
+   USE ice1D          ! sea-ice: thermodynamics
+   !
+   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)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_thd_zdf        ! called by ice_thd
+   PUBLIC   ice_thd_zdf_init   ! called by ice_stp
+
+   !!** namelist (namthd_zdf) **
+   LOGICAL  ::   ln_zdf_Beer      ! Heat diffusion follows a Beer Law
+   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964)
+   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007)
+   REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K]
+   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
+   LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F)
+
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
+   !! $Id: icethd_zdf.F90 8420 2017-08-08 12:18:46Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_thd_zdf
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE ice_thd_zdf  ***
+      !! ** Purpose :
+      !!           This routine determines the time evolution of snow and sea-ice 
+      !!           temperature profiles.
+      !! ** Method  :
+      !!           This is done by solving the heat equation diffusion with
+      !!           a Neumann boundary condition at the surface and a Dirichlet one
+      !!           at the bottom. Solar radiation is partially absorbed into the ice.
+      !!           The specific heat and thermal conductivities depend on ice salinity
+      !!           and temperature to take into account brine pocket melting. The 
+      !!           numerical
+      !!           scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid 
+      !!           in the ice and snow system.
+      !!
+      !!           The successive steps of this routine are
+      !!           1.  Thermal conductivity at the interfaces of the ice layers
+      !!           2.  Internal absorbed radiation
+      !!           3.  Scale factors due to non-uniform grid
+      !!           4.  Kappa factors
+      !!           Then iterative procedure begins
+      !!           5.  specific heat in the ice
+      !!           6.  eta factors
+      !!           7.  surface flux computation
+      !!           8.  tridiagonal system terms
+      !!           9.  solving the tridiagonal system with Gauss elimination
+      !!           Iterative procedure ends according to a criterion on evolution
+      !!           of temperature
+      !!
+      !! ** Inputs / Ouputs : (global commons)
+      !!           surface temperature : t_su_1d
+      !!           ice/snow temperatures   : t_i_1d, t_s_1d
+      !!           ice salinities          : s_i_1d
+      !!           number of layers in the ice/snow: nlay_i, nlay_s
+      !!           total ice/snow thickness : ht_i_1d, ht_s_1d
+      !!------------------------------------------------------------------
+      INTEGER ::   ji, jk         ! spatial loop index
+      INTEGER ::   numeq          ! current reference number of equation
+      INTEGER ::   minnumeqmin, maxnumeqmax
+      INTEGER ::   iconv          ! number of iterations in iterative procedure
+      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
+      
+      INTEGER, DIMENSION(jpij) ::   numeqmin   ! reference number of top equation
+      INTEGER, DIMENSION(jpij) ::   numeqmax   ! reference number of bottom equation
+      
+      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
+      REAL(wp) ::   zg1       =  2._wp        !
+      REAL(wp) ::   zgamma    =  18009._wp    ! for specific heat
+      REAL(wp) ::   zbeta     =  0.117_wp     ! for thermal conductivity (could be 0.13)
+      REAL(wp) ::   zraext_s  =  10._wp       ! extinction coefficient of radiation in the snow
+      REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
+      REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
+      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature 
+      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
+      REAL(wp) ::   z1_hsu
+      REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
+      REAL(wp) ::   zcpi                      ! Ice specific heat
+      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
+      REAL(wp) ::   zfac                      ! dummy factor
+      
+      REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
+      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
+      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness
+      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness
+      REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
+      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
+      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
+      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
+      REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
+      
+      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice
+      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztib        ! Temporary temperature in the ice to check the convergence
+      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsb        ! Temporary temperature in the snow to check the convergence
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
+      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
+      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
+      REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
+      REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
+      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
+      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
+      
+      ! Mono-category
+      REAL(wp) :: zepsilon      ! determines thres. above which computation of G(h) is done
+      REAL(wp) :: zhe           ! dummy factor
+      REAL(wp) :: zcnd_i        ! mean sea ice thermal conductivity
+      !!------------------------------------------------------------------     
+
+      ! --- diag error on heat diffusion - PART 1 --- !
+      DO ji = 1, nidx
+         zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
+            &           SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
+      END DO
+
+      !------------------------------------------------------------------------------!
+      ! 1) Initialization                                                            !
+      !------------------------------------------------------------------------------!
+      DO ji = 1, nidx
+         isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
+         ! layer thickness
+         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
+         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
+      END DO
+      !
+      WHERE( zh_i(1:nidx) >= epsi10 )   ;   z1_h_i(1:nidx) = 1._wp / zh_i(1:nidx)
+      ELSEWHERE                         ;   z1_h_i(1:nidx) = 0._wp
+      END WHERE
+
+      WHERE( zh_s(1:nidx) >= epsi10 )   ;   z1_h_s(1:nidx) = 1._wp / zh_s(1:nidx)
+      ELSEWHERE                         ;   z1_h_s(1:nidx) = 0._wp
+      END WHERE
+      !
+      ! temperatures
+      ztsub  (1:nidx)   = t_su_1d(1:nidx)  ! temperature at the previous iteration
+      ztsold (1:nidx,:) = t_s_1d(1:nidx,:)   ! Old snow temperature
+      ztiold (1:nidx,:) = t_i_1d(1:nidx,:)   ! Old ice temperature
+      t_su_1d(1:nidx)   = MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
+      !
+      !------------------------------------------------------------------------------|
+      ! 2) Radiation	                                                             |
+      !------------------------------------------------------------------------------|
+      !
+      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
+      DO ji = 1, nidx
+         !-------------------
+         ! Computation of i0
+         !-------------------
+         ! i0 describes the fraction of solar radiation which does not contribute
+         ! to the surface energy budget but rather penetrates inside the ice.
+         ! We assume that no radiation is transmitted through the snow
+         ! If there is no no snow
+         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
+         ! zftrice = io.qsr_ice       is below the surface 
+         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
+         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
+         zfac = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )     
+         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) )
+
+         !-------------------------------------------------------
+         ! Solar radiation absorbed / transmitted at the surface
+         ! Derivative of the non solar flux
+         !-------------------------------------------------------
+         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
+         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
+         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
+         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
+      END DO
+
+      !---------------------------------------------------------
+      ! Transmission - absorption of solar radiation in the ice
+      !---------------------------------------------------------
+      zradtr_s(1:nidx,0) = zftrice(1:nidx)
+      DO jk = 1, nlay_s
+         DO ji = 1, nidx
+            !                             ! radiation transmitted below the layer-th snow layer
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * zh_s(ji) * REAL(jk) )
+            !                             ! radiation absorbed by the layer-th snow layer
+            zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
+         END DO
+      END DO
+         
+      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
+      DO jk = 1, nlay_i 
+         DO ji = 1, nidx
+            !                             ! radiation transmitted below the layer-th ice layer
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * zh_i(ji) * REAL(jk) )
+            !                             ! radiation absorbed by the layer-th ice layer
+            zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
+         END DO
+      END DO
+
+      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
+
+      !------------------------------------------------------------------------------|
+      !  3) Iterative procedure begins                                               |
+      !------------------------------------------------------------------------------|
+      !
+      iconv    =  0          ! number of iterations
+      zdti_max =  1000._wp   ! maximal value of error on all points
+      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )
+         !
+         iconv = iconv + 1
+         !
+         ztib(1:nidx,:) = t_i_1d(1:nidx,:)
+         ztsb(1:nidx,:) = t_s_1d(1:nidx,:)
+         !
+         !------------------------------------------------------------------------------|
+         ! 4) Sea ice thermal conductivity                                              |
+         !------------------------------------------------------------------------------|
+         !
+         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 )
+               ztcond_i(ji,nlay_i) = rcdic + zbeta * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji)  - rt0 )
+            END DO
+            DO jk = 1, nlay_i-1
+               DO ji = 1, nidx
+                  ztcond_i(ji,jk) = rcdic + zbeta * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
+                     &                      MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+               END DO
+            END DO
+            !
+         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.09_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )  &
+                  &                        - 0.011_wp * ( t_i_1d(ji,1) - rt0 )
+               ztcond_i(ji,nlay_i) = rcdic + 0.09_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )  &
+                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )
+            END DO
+            DO jk = 1, nlay_i-1
+               DO ji = 1, nidx
+                  ztcond_i(ji,jk) = rcdic + 0.09_wp * 0.5_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /          &
+                     &                     MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )   &
+                     &                    - 0.011_wp * ( 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )
+               END DO
+            END DO
+            !
+         ENDIF
+         ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) )        
+         !
+         !------------------------------------------------------------------------------|
+         !  5) G(he) - enhancement of thermal conductivity in mono-category case        |
+         !------------------------------------------------------------------------------|
+         !
+         ! Computation of effective thermal conductivity G(h)
+         ! Used in mono-category case only to simulate an ITD implicitly
+         ! Fichefet and Morales Maqueda, JGR 1997
+
+         zghe(1:nidx) = 1._wp
+         
+         SELECT CASE ( nn_monocat )
+
+         CASE ( 1 , 3 )
+
+            zepsilon = 0.1_wp
+            DO ji = 1, nidx
+   
+               ! Mean sea ice thermal conductivity
+               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )
+
+               ! Effective thickness he (zhe)
+               zhe  = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )
+
+               ! G(he)
+               IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
+                  zghe(ji) =  MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )
+               ENDIF
+   
+            END DO
+
+         END SELECT
+         !
+         !------------------------------------------------------------------------------|
+         !  6) kappa factors                                                            |
+         !------------------------------------------------------------------------------|
+         !
+         !--- Snow
+         DO jk = 0, nlay_s-1
+            DO ji = 1, nidx
+               zkappa_s(ji,jk)  = zghe(ji) * rn_cnd_s * z1_h_s(ji)
+            END DO
+         END DO
+         DO ji = 1, nidx   ! Snow-ice interface
+            zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) )
+            IF( zfac > epsi10 ) THEN
+               zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac
+            ELSE
+               zkappa_s(ji,nlay_s) = 0._wp
+            ENDIF
+         END DO
+
+         !--- Ice
+         DO jk = 0, nlay_i
+            DO ji = 1, nidx
+               zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji)
+            END DO
+         END DO
+         DO ji = 1, nidx   ! Snow-ice interface
+            zkappa_i(ji,0)     = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
+         END DO
+         !
+         !------------------------------------------------------------------------------|
+         ! 7) Sea ice specific heat, eta factors                                        |
+         !------------------------------------------------------------------------------|
+         !
+         DO jk = 1, nlay_i
+            DO ji = 1, nidx
+               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztiold(ji,jk) - rt0 ), epsi10 )
+               zeta_i(ji,jk) = rdt_ice * r1_rhoic * z1_h_i(ji) / MAX( epsi10, zcpi ) 
+            END DO
+         END DO
+
+         DO jk = 1, nlay_s
+            DO ji = 1, nidx
+               zeta_s(ji,jk)  = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
+            END DO
+         END DO
+         !
+         !------------------------------------------------------------------------------|
+         ! 8) surface flux computation                                                  |
+         !------------------------------------------------------------------------------|
+         !
+         IF ( ln_dqns_i ) THEN 
+            DO ji = 1, nidx
+               ! update of the non solar flux according to the update in T_su
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
+            END DO
+         ENDIF
+
+         DO ji = 1, nidx
+            zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
+         END DO
+         !
+         !------------------------------------------------------------------------------|
+         ! 9) tridiagonal system terms                                                  |
+         !------------------------------------------------------------------------------|
+         !
+         !!layer denotes the number of the layer in the snow or in the ice
+         !!numeq denotes the reference number of the equation in the tridiagonal
+         !!system, terms of tridiagonal system are indexed as following :
+         !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
+
+         !!ice interior terms (top equation has the same form as the others)
+
+         DO numeq=1,nlay_i+3
+            DO ji = 1, nidx
+               ztrid(ji,numeq,1) = 0.
+               ztrid(ji,numeq,2) = 0.
+               ztrid(ji,numeq,3) = 0.
+               zindterm(ji,numeq)= 0.
+               zindtbis(ji,numeq)= 0.
+               zdiagbis(ji,numeq)= 0.
+            ENDDO
+         ENDDO
+
+         DO numeq = nlay_s + 2, nlay_s + nlay_i 
+            DO ji = 1, nidx
+               jk                 = numeq - nlay_s - 1
+               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
+               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
+               zindterm(ji,numeq) =  ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+            END DO
+         ENDDO
+
+         numeq =  nlay_s + nlay_i + 1
+         DO ji = 1, nidx
+            !!ice bottom term
+            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
+            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
+            ztrid(ji,numeq,3)  =  0.0
+            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
+         ENDDO
+
+
+         DO ji = 1, nidx
+            IF ( ht_s_1d(ji) > 0.0 ) THEN
+               !
+               !------------------------------------------------------------------------------|
+               !  snow-covered cells                                                          |
+               !------------------------------------------------------------------------------|
+               !
+               !!snow interior terms (bottom equation has the same form as the others)
+               DO numeq = 3, nlay_s + 1
+                  jk                  =  numeq - 1
+                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
+                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
+                  zindterm(ji,numeq)  =  ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+               END DO
+
+               !!case of only one layer in the ice (ice equation is altered)
+               IF ( nlay_i == 1 ) THEN
+                  ztrid(ji,nlay_s+2,3)    =  0.0
+                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+               ENDIF
+
+               IF ( t_su_1d(ji) < rt0 ) THEN
+
+                  !------------------------------------------------------------------------------|
+                  !  case 1 : no surface melting - snow present                                  |
+                  !------------------------------------------------------------------------------|
+                  numeqmin(ji)    =  1
+                  numeqmax(ji)    =  nlay_i + nlay_s + 1
+
+                  !!surface equation
+                  ztrid(ji,1,1)  = 0.0
+                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
+                  ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
+                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji)
+
+                  !!first layer of snow equation
+                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
+                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
+                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
+
+               ELSE 
+                  !
+                  !------------------------------------------------------------------------------|
+                  !  case 2 : surface is melting - snow present                                  |
+                  !------------------------------------------------------------------------------|
+                  !
+                  numeqmin(ji)    =  2
+                  numeqmax(ji)    =  nlay_i + nlay_s + 1
+
+                  !!first layer of snow equation
+                  ztrid(ji,2,1)  =  0.0
+                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
+                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
+                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
+               ENDIF
+            ELSE
+               !
+               !------------------------------------------------------------------------------|
+               !  cells without snow                                                          |
+               !------------------------------------------------------------------------------|
+               !
+               IF ( t_su_1d(ji) < rt0 ) THEN
+                  !
+                  !------------------------------------------------------------------------------|
+                  !  case 3 : no surface melting - no snow                                       |
+                  !------------------------------------------------------------------------------|
+                  !
+                  numeqmin(ji)      =  nlay_s + 1
+                  numeqmax(ji)      =  nlay_i + nlay_s + 1
+
+                  !!surface equation	
+                  ztrid(ji,numeqmin(ji),1)   =  0.0
+                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
+                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
+                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji)
+
+                  !!first layer of ice equation
+                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
+                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
+                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1)  
+                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
+
+                  !!case of only one layer in the ice (surface & ice equations are altered)
+
+                  IF ( nlay_i == 1 ) THEN
+                     ztrid(ji,numeqmin(ji),1)    =  0.0
+                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
+                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
+                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
+                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
+                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
+
+                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1) *  &
+                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
+                  ENDIF
+
+               ELSE
+
+                  !
+                  !------------------------------------------------------------------------------|
+                  ! case 4 : surface is melting - no snow                                        |
+                  !------------------------------------------------------------------------------|
+                  !
+                  numeqmin(ji)    =  nlay_s + 2
+                  numeqmax(ji)    =  nlay_i + nlay_s + 1
+
+                  !!first layer of ice equation
+                  ztrid(ji,numeqmin(ji),1)      =  0.0
+                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
+                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
+                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1) *  &
+                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
+
+                  !!case of only one layer in the ice (surface & ice equations are altered)
+                  IF ( nlay_i == 1 ) THEN
+                     ztrid(ji,numeqmin(ji),1)  =  0.0
+                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
+                     ztrid(ji,numeqmin(ji),3)  =  0.0
+                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
+                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
+                  ENDIF
+
+               ENDIF
+            ENDIF
+
+         END DO
+         !
+         !------------------------------------------------------------------------------|
+         ! 10) tridiagonal system solving                                               |
+         !------------------------------------------------------------------------------|
+         !
+         ! Solve the tridiagonal system with Gauss elimination method.
+         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, 
+         ! McGraw-Hill 1984.	
+
+         maxnumeqmax = 0
+         minnumeqmin = nlay_i+5
+
+         DO ji = 1, nidx
+            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
+            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
+            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
+            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
+         END DO
+
+         DO jk = minnumeqmin+1, maxnumeqmax
+            DO ji = 1, nidx
+               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
+               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
+               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
+            END DO
+         END DO
+
+         DO ji = 1, nidx
+            ! ice temperatures
+            t_i_1d(ji,nlay_i) =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
+         END DO
+
+         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
+            DO ji = 1, nidx
+               jk    =  numeq - nlay_s - 1
+               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
+            END DO
+         END DO
+
+         DO ji = 1, nidx
+            ! snow temperatures      
+            IF( ht_s_1d(ji) > 0._wp ) THEN
+               t_s_1d(ji,nlay_s) =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
+                  &                 / zdiagbis(ji,nlay_s+1)
+            ENDIF
+            ! surface temperature
+            ztsub(ji) = t_su_1d(ji)
+            IF( t_su_1d(ji) < rt0 ) THEN
+               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
+                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
+            ENDIF
+         END DO
+         !
+         !--------------------------------------------------------------------------
+         !  11) Has the scheme converged ?, end of the iterative procedure         |
+         !--------------------------------------------------------------------------
+         !
+         ! check that nowhere it has started to melt
+         ! zdti_max is a measure of error, it has to be under zdti_bnd
+         zdti_max = 0._wp
+         DO ji = 1, nidx
+            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp )
+            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )     
+         END DO
+
+         DO jk  =  1, nlay_s
+            DO ji = 1, nidx
+               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp )
+               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) )
+            END DO
+         END DO
+
+         DO jk  =  1, nlay_i
+            DO ji = 1, nidx
+               ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
+               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp )
+               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) )
+            END DO
+         END DO
+
+         ! Compute spatial maximum over all errors
+         ! note that this could be optimized substantially by iterating only the non-converging points
+         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
+
+      END DO  ! End of the do while iterative procedure
+
+      IF( ln_icectl .AND. lwp ) THEN
+         WRITE(numout,*) ' zdti_max : ', zdti_max
+         WRITE(numout,*) ' iconv    : ', iconv
+      ENDIF
+
+      !
+      !-------------------------------------------------------------------------!
+      !   12) Fluxes at the interfaces                                          !
+      !-------------------------------------------------------------------------!
+      DO ji = 1, nidx
+         !                                ! surface ice conduction flux
+         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
+            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
+         !                                ! bottom ice conduction flux
+         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
+      END DO
+
+      ! --- computes sea ice energy of melting compulsory for icethd_dh --- !
+      CALL ice_thd_enmelt
+
+      ! --- diagnose the change in non-solar flux due to surface temperature change --- !
+      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) 
+         END DO
+      END IF
+
+      ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
+      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
+      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
+      DO ji = 1, nidx
+         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
+            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
+
+         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
+            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
+         ELSE                          ! case T_su = 0degC
+            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
+         ENDIF
+         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
+
+         ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
+         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
+
+      END DO   
+
+      ! --- Diagnostics SIMIP --- !
+      DO ji = 1, nidx
+         !--- Conduction fluxes (positive downwards)
+         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
+         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
+
+         !--- Snow-ice interfacial temperature (diagnostic SIMIP)
+         zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
+         IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
+            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
+            t_si_1d(ji) = t_su_1d(ji)
+         ENDIF
+      END DO
+      !
+   END SUBROUTINE ice_thd_zdf
+
+
+   SUBROUTINE ice_thd_enmelt
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_thd_enmelt *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts  ! local scalar 
+      !!-------------------------------------------------------------------
+      !
+      DO jk = 1, nlay_i             ! Sea ice energy of melting
+         DO ji = 1, nidx
+            ztmelts      = - tmut  * s_i_1d(ji,jk)
+            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point
+                                                                !   (sometimes dif scheme produces abnormally high temperatures)   
+            e_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
+               &                    + lfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
+               &                    - rcp  *   ztmelts )
+         END DO
+      END DO
+      DO jk = 1, nlay_s             ! Snow energy of melting
+         DO ji = 1, nidx
+            e_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus )
+         END DO
+      END DO
+      !
+   END SUBROUTINE ice_thd_enmelt
+
+
+   SUBROUTINE ice_thd_zdf_init
+      !!-----------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_thd_zdf_init *** 
+      !!                 
+      !! ** Purpose :   Physical constants and parameters associated with
+      !!                ice thermodynamics
+      !!
+      !! ** Method  :   Read the namthd_zdf namelist and check the parameters
+      !!                called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namthd_zdf
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ios   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namthd_zdf/ ln_zdf_Beer, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i
+      !!-------------------------------------------------------------------
+      !
+      REWIND( numnam_ice_ref )              ! Namelist namthd_zdf in reference namelist : Ice thermodynamics
+      READ  ( numnam_ice_ref, namthd_zdf, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namthd_zdf in configuration namelist : Ice thermodynamics
+      READ  ( numnam_ice_cfg, namthd_zdf, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_zdf in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd_zdf )
+      !
+      !
+      IF(lwp) THEN                          ! control print
+         WRITE(numout,*) 'ice_thd_zdf_init: Ice vertical heat diffusion'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namthd_zdf:'
+         WRITE(numout,*) '      Diffusion follows a Beer Law                            ln_zdf_Beer  = ', ln_zdf_Beer
+         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,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s
+         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
+         WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i
+     ENDIF
+      !
+      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN
+         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conductivity (ln_cndi_U64 or ln_cndi_P07)' )
+      ENDIF
+      !
+   END SUBROUTINE ice_thd_zdf_init
+
+#else
+   !!----------------------------------------------------------------------
+   !!                   Dummy Module                 No ESIM sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!======================================================================
+END MODULE icethd_zdf
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8530)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8531)
@@ -79,21 +79,21 @@
       !!      over sea ice
       !!
-      !! ** Method  :  Read the namicemp  namelist and check the melt pond  
+      !! ** Method  :  Read the nammp  namelist and check the melt pond  
       !!       parameter values called at the first timestep (nit000)
       !!
-      !! ** input   :   Namelist namicemp  
+      !! ** input   :   Namelist nammp  
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicemp/  ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd
+      NAMELIST/nammp/  ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd
       !!-------------------------------------------------------------------
 
-      REWIND( numnam_ice_ref )              ! Namelist namicemp  in reference namelist : Melt Ponds  
-      READ  ( numnam_ice_ref, namicemp, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicemp  in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namicemp  in configuration namelist : Melt Ponds
-      READ  ( numnam_ice_cfg, namicemp, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicemp in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namicemp )
+      REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds  
+      READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds
+      READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, nammp )
       
       IF(lwp) THEN                        ! control print
