Index: /NEMO/trunk/src/ICE/ice.F90
===================================================================
--- /NEMO/trunk/src/ICE/ice.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/ice.F90	(revision 9604)
@@ -4,6 +4,6 @@
    !!   sea-ice:  ice variables defined in memory
    !!======================================================================
-   !! History :  3.0  ! 2008-03  (M. Vancoppenolle) original code LIM-3
-   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
+   !! History :  3.0  !  2008-03  (M. Vancoppenolle) Original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -19,29 +19,4 @@
    PUBLIC   ice_alloc   ! called by icestp.F90
 
-   !!======================================================================
-   !! LIM3 by the use of sweat, agile fingers and sometimes brain juice, 
-   !!  was developed in Louvain-la-Neuve by : 
-   !!    * Martin Vancoppenolle (UCL-ASTR, Belgium)
-   !!    * Sylvain Bouillon (UCL-ASTR, Belgium)
-   !!    * Miguel Angel Morales Maqueda (NOC-L, UK)
-   !! 
-   !! Based on extremely valuable earlier work by
-   !!    * Thierry Fichefet
-   !!    * Hugues Goosse
-   !!
-   !! The following persons also contributed to the code in various ways
-   !!    * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France)
-   !!    * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
-   !!    * Bill Lipscomb (LANL), Cecilia Bitz (UWa) 
-   !!      and Elisabeth Hunke (LANL), USA.
-   !! 
-   !! For more info, the interested user is kindly invited to consult the following references
-   !!    For model description and validation :
-   !!    * Vancoppenolle et al., Ocean Modelling, 2008a.
-   !!    * Vancoppenolle et al., Ocean Modelling, 2008b.
-   !!    For a specific description of EVP :
-   !!    * Bouillon et al., Ocean Modelling 2009.
-   !!
-   !!    Or the reference manual, that should be available by 2011
    !!======================================================================
    !!                                                                     |
@@ -62,23 +37,4 @@
    !!   These are the variables that are either physically more           |
    !!   meaningful and/or used in ice thermodynamics                      |
-   !!                                                                     |
-   !! Routines in icevar.F90 perform conversions                          |
-   !!  - ice_var_glo2eqv  : from global to equivalent variables           |
-   !!  - ice_var_eqv2glo  : from equivalent to global variables           |
-   !!                                                                     |
-   !! For various purposes, the sea ice state variables have sometimes    |
-   !! to be aggregated over all ice thickness categories. This operation  |
-   !! is done in :                                                        |
-   !!  - ice_var_agg                                                      |
-   !!                                                                     |
-   !! in icestp.F90, the routines that compute the changes in the ice     |
-   !! state variables are called                                          |
-   !! - ice_rhg : ice dynamics                                            |
-   !! - ice_adv : ice advection                                           |
-   !! - ice_rdgrft : ice ridging and rafting                              |
-   !! - ice_thd : ice halo-thermodynamics and creation of new ice         |
-   !! - ice_itd : thermodynamic changes in ice thickness distribution     |
-   !!                                                                     |
-   !! See the associated routines for more information                    |
    !!                                                                     |
    !! List of ice state variables :                                       |
@@ -249,6 +205,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   delta_i     !: ice rheology elta factor (Flato & Hibler 95) [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   divu_i      !: Divergence of the velocity field [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   divu_i      !: Divergence of the velocity field             [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field                  [s-1]
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
@@ -257,14 +213,14 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: snow-ocean mass exchange   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: snow-ocean mass exchange             [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sni !: snow ice growth component of wfx_snw [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: melt pond-ocean mass exchange   [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice  [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: sublimation of snow/ice    [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sub !: snow sublimation           [kg.m-2.s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice_sub !: ice sublimation            [kg.m-2.s-1]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: dynamical component of wfx_snw    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: melt pond-ocean mass exchange        [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: snow precipitation on ice            [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: sublimation of snow/ice              [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sub !: snow sublimation                     [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice_sub !: ice sublimation                      [kg.m-2.s-1]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: dynamical component of wfx_snw       [kg.m-2.s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice     !: ice-ocean mass exchange                   [kg.m-2.s-1]
@@ -278,5 +234,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res     !: residual component of wfx_ice             [kg.m-2.s-1]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)        [s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -292,18 +248,18 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth        [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt          [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt         [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice   [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                          [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth           [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt             [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt            [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation    [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice      [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                             [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping             [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations   [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping                [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations      [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  [W.m-2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation                   [kg.m-2.s-1]
    
    ! heat flux associated with ice-atmosphere mass exchange
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation            [W.m-2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  [W.m-2]
 
@@ -421,15 +377,12 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(18), ii
+      INTEGER :: ierr(15), ii
       !!-----------------------------------------------------------------
-
       ierr(:) = 0
 
-      ! What could be one huge allocate statement is broken-up to try to
-      ! stay within Fortran's max-line length limit.
       ii = 1
-      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new (jpi,jpj) ,                        &
-         &      strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  &
-         &      delta_i (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) , STAT=ierr(ii) )
+      ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new  (jpi,jpj) , strength(jpi,jpj) ,  &
+         &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      &
+         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) )
 
       ii = ii + 1
@@ -437,12 +390,12 @@
          &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) ,  &
          &      wfx_ice(jpi,jpj) , wfx_sub    (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam    (jpi,jpj) ,  &
-         &      wfx_pnd(jpi,jpj) ,                                                              &
-         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
-         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
-         &      afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj),                                         &
-         &      fhtur  (jpi,jpj) , qlead  (jpi,jpj) ,                                           &
+         &      wfx_pnd(jpi,jpj) ,                                                                       &
+         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,              &
+         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,              &
+         &      afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj),                                                  &
+         &      fhtur  (jpi,jpj) , qlead  (jpi,jpj) ,                                                    &
          &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  &
-         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) ,     & 
+         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) ,                        & 
          &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld   (jpi,jpj) ,                        &
          &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     &
@@ -452,17 +405,20 @@
       ! * Ice global state variables
       ii = ii + 1
-      ALLOCATE( ftr_ice(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , &
-         &      h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) ,     &
-         &      v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , t_su(jpi,jpj,jpl) ,     &
-         &      s_i (jpi,jpj,jpl) , sv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) ,     &
-         &      oa_i(jpi,jpj,jpl) , bv_i(jpi,jpj,jpl) ,  STAT=ierr(ii) )
-      ii = ii + 1
-      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                      &
-         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,     &
-         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s (jpi,jpj) , bvm_i(jpi,jpj) ,  &
-         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s (jpi,jpj) ,     &
-         &      om_i (jpi,jpj) , tau_icebfr(jpi,jpj)                             , STAT=ierr(ii) )
+      ALLOCATE( ftr_ice(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  &
+         &      h_i    (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  &
+         &      v_s    (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  &
+         &      s_i    (jpi,jpj,jpl) , sv_i   (jpi,jpj,jpl) , o_i   (jpi,jpj,jpl) ,  &
+         &      oa_i   (jpi,jpj,jpl) , bv_i   (jpi,jpj,jpl) , STAT=ierr(ii) )
+
+      ii = ii + 1
+      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,                                   &
+         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i(jpi,jpj) , ato_i(jpi,jpj) ,  &
+         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s (jpi,jpj) ,  &
+         &      sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s (jpi,jpj) ,  &
+         &      om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj)            , STAT=ierr(ii) )
+
       ii = ii + 1
       ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) )
+
       ii = ii + 1
       ALLOCATE( t_i(jpi,jpj,nlay_i,jpl) , e_i(jpi,jpj,nlay_i,jpl) , sz_i(jpi,jpj,nlay_i,jpl) , STAT=ierr(ii) )
@@ -470,4 +426,5 @@
       ii = ii + 1
       ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
+
       ii = ii + 1
       ALLOCATE( at_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) )
@@ -475,7 +432,8 @@
       ! * Old values of global variables
       ii = ii + 1
-      ALLOCATE( v_s_b  (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   &
-         &      a_i_b  (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   &
-         &      oa_i_b (jpi,jpj,jpl)                                                   , STAT=ierr(ii) )
+      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl)        ,   &
+         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,   &
+         &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) )
+
       ii = ii + 1
       ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) )
Index: /NEMO/trunk/src/ICE/ice1d.F90
===================================================================
--- /NEMO/trunk/src/ICE/ice1d.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/ice1d.F90	(revision 9604)
@@ -4,5 +4,6 @@
    !! sea-ice :   Ice thermodynamics variables in 1D
    !!=====================================================================
-   !! History :  3.0  !  2002-11  (C. Ethe)  F90: Free form and module
+   !! History :  3.0  !  2002-11  (C. Ethe)          original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -179,10 +180,9 @@
       ii = 1
       ALLOCATE( nptidx    (jpij) ,   &
-         &      qlead_1d  (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d (jpij) ,   &
-         &      qns_ice_1d(jpij) ,   &
-         &      qml_ice_1d(jpij) , qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , &
-         &      cnd_ice_1d(jpij) , t1_ice_1d(jpij)  , t_bo_1d   (jpij) ,      &
-         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) , hfx_bog_1d(jpij) ,    & 
-         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) , hfx_dyn_1d(jpij) ,    &
+         &      qlead_1d  (jpij) , ftr_ice_1d(jpij) , qsr_ice_1d(jpij) ,   &
+         &      qns_ice_1d(jpij) , qml_ice_1d(jpij) , qcn_ice_1d(jpij) , qsr_ice_tr_1d(jpij) , &
+         &      cnd_ice_1d(jpij) , t1_ice_1d (jpij) , t_bo_1d   (jpij) ,   &
+         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) , hfx_bog_1d(jpij) ,   & 
+         &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) , hfx_dyn_1d(jpij) ,   &
          &      rn_amax_1d(jpij) ,                                         &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
@@ -191,27 +191,27 @@
       !
       ii = ii + 1
-      ALLOCATE( sprecip_1d (jpij) , at_i_1d    (jpij) , ato_i_1d(jpij) ,                    &
-         &      fhtur_1d   (jpij) , wfx_snw_sni_1d (jpij) , wfx_spr_1d (jpij) , wfx_snw_sum_1d(jpij) , &
-         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  &
-         &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
-         &      wfx_snw_sub_1d(jpij), wfx_snw_dyn_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,   &
-         &      wfx_lam_1d(jpij)  , wfx_dyn_1d(jpij), wfx_pnd_1d(jpij),dqns_ice_1d(jpij) , evap_ice_1d (jpij),   &
-         &      qprec_ice_1d(jpij),                                                             &  
-         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
-         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
-         &      sfx_lam_1d (jpij) , sfx_dyn_1d(jpij)  , STAT=ierr(ii) )
-      !
-      ii = ii + 1
-      ALLOCATE( t_su_1d  (jpij) , t_si_1d   (jpij) , a_i_1d  (jpij) , a_ib_1d(jpij) ,                  &
-         &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &    
-         &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &    
-         &      dh_i_bott(jpij) , dh_s_mlt(jpij), dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , &
-         &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , &
-         &      h_ip_1d  (jpij) , a_ip_frac_1d(jpij) , &
-         &      sv_i_1d (jpij) , oa_i_1d   (jpij) , STAT=ierr(ii) )
+      ALLOCATE( sprecip_1d    (jpij) , at_i_1d       (jpij) , ato_i_1d      (jpij) ,                         &
+         &      fhtur_1d      (jpij) , wfx_snw_sni_1d(jpij) , wfx_spr_1d    (jpij) , wfx_snw_sum_1d(jpij) ,  &
+         &      fhld_1d       (jpij) , wfx_sub_1d    (jpij) , wfx_bog_1d    (jpij) , wfx_bom_1d    (jpij) ,  &
+         &      wfx_sum_1d    (jpij) , wfx_sni_1d    (jpij) , wfx_opw_1d    (jpij) , wfx_res_1d    (jpij) ,  &
+         &      wfx_snw_sub_1d(jpij) , wfx_snw_dyn_1d(jpij) , wfx_ice_sub_1d(jpij) , wfx_err_sub_1d(jpij) ,  &
+         &      wfx_lam_1d    (jpij) , wfx_dyn_1d    (jpij) , wfx_pnd_1d    (jpij) , dqns_ice_1d   (jpij) , evap_ice_1d (jpij) , &
+         &      qprec_ice_1d  (jpij) ,                                                             &  
+         &      sfx_bri_1d    (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
+         &      sfx_sni_1d    (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
+         &      sfx_lam_1d    (jpij) , sfx_dyn_1d(jpij)  , STAT=ierr(ii) )
+      !
+      ii = ii + 1
+      ALLOCATE( t_su_1d  (jpij) , t_si_1d  (jpij) , a_i_1d    (jpij) , a_ib_1d(jpij) ,                  &
+         &      h_i_1d   (jpij) , h_ib_1d  (jpij) , h_s_1d    (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &    
+         &      dh_s_tot (jpij) , dh_i_surf(jpij) , dh_i_sub  (jpij) ,                                  &    
+         &      dh_i_bott(jpij) , dh_s_mlt (jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) ,  &
+         &      a_ip_1d  (jpij) , v_ip_1d  (jpij) , v_i_1d    (jpij) , v_s_1d (jpij) ,                  &
+         &      h_ip_1d  (jpij) , a_ip_frac_1d(jpij) ,                                                  &
+         &      sv_i_1d  (jpij) , oa_i_1d  (jpij) , STAT=ierr(ii) )
       !
       ii = ii + 1
       ALLOCATE( t_s_1d  (jpij,nlay_s)     , t_i_1d (jpij,nlay_i)     , sz_i_1d(jpij,nlay_i) ,  &            
-         &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                        &
+         &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                         &
          &      eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) )
       !
@@ -223,7 +223,7 @@
       !
       ii = ii + 1
-      ALLOCATE( a_i_2d(jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d(jpij,jpl) , h_ib_2d(jpij,jpl) , &
-         &      v_i_2d(jpij,jpl) ,v_s_2d(jpij,jpl) ,oa_i_2d(jpij,jpl) ,sv_i_2d(jpij,jpl) ,  &
-         &      a_ip_2d(jpij,jpl) ,v_ip_2d(jpij,jpl) ,t_su_2d(jpij,jpl) ,  &
+      ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  &
+         &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  &
+         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      &
          &      STAT=ierr(ii) )
 
Index: /NEMO/trunk/src/ICE/icealb.F90
===================================================================
--- /NEMO/trunk/src/ICE/icealb.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icealb.F90	(revision 9604)
@@ -4,5 +4,6 @@
    !! Atmospheric forcing:  Albedo over sea ice
    !!=====================================================================
-   !! History :  4.0  ! 2017-07  (C. Rousset) Split ice and ocean albedos
+   !! History :  4.0  !  2017-07  (C. Rousset)       Split ice and ocean albedos
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -28,10 +29,4 @@
 
    REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
-
-   REAL(wp) , PARAMETER ::   ppc1    = 0.05    ! snow thickness (only for nn_ice_alb=0)
-   REAL(wp) , PARAMETER ::   ppc2    = 0.10    !  "        "
-   REAL(wp) , PARAMETER ::   ppcloud = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0)
-   REAL(wp) , PARAMETER ::   pp1_c1 = 1. / ppc1
-   REAL(wp) , PARAMETER ::   pp1_c2 = 1. / ppc2
    !
    !                             !!* albedo namelist (namalb)
Index: /NEMO/trunk/src/ICE/icecor.F90
===================================================================
--- /NEMO/trunk/src/ICE/icecor.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icecor.F90	(revision 9604)
@@ -2,8 +2,9 @@
    !!======================================================================
    !!                     ***  MODULE  icecor  ***
-   !!   SI3 : Corrections on sea-ice variables at the end of the time step
+   !!   sea-ice: Corrections on sea-ice variables at the end of the time step
    !!======================================================================
    !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
    !!            3.5  !  2014-06  (C. Rousset)       Complete rewriting/cleaning
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icectl.F90
===================================================================
--- /NEMO/trunk/src/ICE/icectl.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icectl.F90	(revision 9604)
@@ -6,4 +6,5 @@
    !! History :  3.5  !  2015-01  (M. Vancoppenolle) Original code
    !!            3.7  !  2016-10  (C. Rousset)       Add routine ice_prt3D
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedia.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedia.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedia.F90	(revision 9604)
@@ -2,8 +2,8 @@
    !!======================================================================
    !!                       ***  MODULE icedia  ***
-   !!  Sea-Ice model :   global budgets 
+   !!  Sea-Ice:   global budgets 
    !!======================================================================
-   !! History :  3.4  ! 2012-10  (C. Rousset)  original code
-   !!            4.0  ! 2017-08  (C. Rousset)  fits nemo4.0 standards
+   !! History :  3.4  !  2012-10  (C. Rousset)       original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedyn.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn.F90	(revision 9604)
@@ -4,5 +4,5 @@
    !!   Sea-Ice dynamics : master routine for sea ice dynamics 
    !!======================================================================
-   !! history :  4.0  ! 2017-09  (C. Rousset)  original code 
+   !! history :  4.0  ! 2018  (C. Rousset)  original code SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -63,6 +63,9 @@
       !! ** Purpose :   this routine manages sea ice dynamics
       !!
-      !! ** Action : - Initialisation of some variables
-      !!             - call ice_rhg
+      !! ** Action : - calculation of friction in case of landfast ice
+      !!             - call ice_dyn_rhg    = rheology
+      !!             - call ice_dyn_adv    = advection
+      !!             - call ice_dyn_rdgrft = ridging/rafting
+      !!             - call ice_cor        = corrections if fields are out of bounds
       !!--------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ice time step
Index: /NEMO/trunk/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn_adv.F90	(revision 9604)
@@ -2,9 +2,7 @@
    !!======================================================================
    !!                       ***  MODULE icedyn_adv   ***
-   !!   SI3 : sea-ice advection
+   !!   sea-ice: advection
    !!======================================================================
-   !! History : LIM-2 ! 2000-01 (M.A. Morales Maqueda, H. Goosse, and T. Fichefet)  Original code
-   !!            3.0  ! 2005-11 (M. Vancoppenolle)   Multi-layer sea ice, salinity variations
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !! History :  4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn_adv_pra.F90	(revision 9604)
@@ -4,7 +4,6 @@
    !!   sea-ice : advection => Prather scheme
    !!======================================================================
-   !! History :  LIM  ! 2008-03 (M. Vancoppenolle)  LIM-3 from LIM-2 code
-   !!            3.2  ! 2009-06 (F. Dupont)  correct a error in the North fold b.c.
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !! History :       !  2008-03  (M. Vancoppenolle) LIM-3 from LIM-2 code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!--------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn_adv_umx.F90	(revision 9604)
@@ -5,4 +5,5 @@
    !!==============================================================================
    !! History :  3.6  !  2014-11  (C. Rousset, G. Madec)  Original code
+   !!            4.0  !  2018     (many people)           SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn_rdgrft.F90	(revision 9604)
@@ -4,7 +4,6 @@
    !!    sea-ice : Mechanical impact on ice thickness distribution      
    !!======================================================================
-   !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code 
-   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn
-   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation
+   !! History :       !  2006-02  (M. Vancoppenolle) Original code 
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rhg.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn_rhg.F90	(revision 9604)
@@ -4,7 +4,5 @@
    !!   Sea-Ice dynamics : master routine for rheology 
    !!======================================================================
-   !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code 
-   !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
-   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
+   !! history :  4.0  !  2018     (C. Rousset)      Original code
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icedyn_rhg_evp.F90	(revision 9604)
@@ -7,8 +7,10 @@
    !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
    !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
-   !!            3.3  !  2009-05  (G.Garric) addition of the evp cas
-   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
-   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
-   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
+   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
+   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation 
+   !!            3.5  !  2012-08  (R. Benshila) AGRIF
+   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
+   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
+   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/iceforcing.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceforcing.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/iceforcing.F90	(revision 9604)
@@ -4,5 +4,6 @@
    !! Sea-Ice :   air-ice forcing fields
    !!=====================================================================
-   !! History :  4.0  ! 2017-08  (C. Rousset) Original code
+   !! History :  4.0  !  2017-08  (C. Rousset)       Original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/iceistate.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceistate.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/iceistate.F90	(revision 9604)
@@ -4,9 +4,7 @@
    !!   sea-ice : Initialization of ice variables
    !!======================================================================
-   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code
-   !!            3.0  ! 2007    (M. Vancoppenolle)  Rewrite for ice cats
-   !!            3.0  ! 2009-11 (M. Vancoppenolle)  Enhanced version for ice cats
-   !!            3.0  ! 2011-02 (G. Madec)  dynamical allocation
-   !!             -   ! 2014    (C. Rousset)  add N/S initializations
+   !! History :  2.0  !  2004-01  (C. Ethe, G. Madec) Original code
+   !!            3.0  !  2007     (M. Vancoppenolle)  Rewrite for ice cats
+   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/iceitd.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceitd.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/iceitd.F90	(revision 9604)
@@ -4,8 +4,6 @@
    !!   sea-ice : ice thickness distribution
    !!======================================================================
-   !! History :   -   !          (W. H. Lipscomb and E.C. Hunke) CICE (c) original code
-   !!            3.0  ! 2005-12  (M. Vancoppenolle) adaptation to LIM-3
-   !!             -   ! 2006-06  (M. Vancoppenolle) adaptation to include salt, age
-   !!             -   ! 2007-04  (M. Vancoppenolle) Mass conservation checked
+   !! History :  3.0  !  2005-12  (M. Vancoppenolle) original code (based on CICE)
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icerst.F90
===================================================================
--- /NEMO/trunk/src/ICE/icerst.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icerst.F90	(revision 9604)
@@ -4,8 +4,5 @@
    !!   sea-ice :  write/read the ice restart file
    !!======================================================================
-   !! History:   3.0  ! 2005-04 (M. Vancoppenolle) Original code
-   !!             -   ! 2008-03 (C. Ethe) restart files in using IOM interface
-   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
-   !!            4.0  ! 2017-11 (M. Andrejczuk) Extend IOM interface to write 3D ice fields
+   !! History:   4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icestp.F90
===================================================================
--- /NEMO/trunk/src/ICE/icestp.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icestp.F90	(revision 9604)
@@ -2,15 +2,27 @@
    !!======================================================================
    !!                       ***  MODULE  icestp  ***
-   !! Surface module :  update the ocean surface boundary condition over ice
-   !!                   covered area using SI3 sea-ice model
+   !! sea ice : Master routine for all the sea ice model
    !!=====================================================================
-   !! History :  2.0  ! 2006-12  (M. Vancoppenolle) Original code
-   !!            3.0  ! 2008-02  (C. Talandier)  Surface module from icestp.F90
-   !!             -   ! 2008-04  (G. Madec)  sltyle and ice_ctl routine
-   !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
-   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
-   !!             -   ! 2012-10  (C. Rousset)  add ice_dia
-   !!            3.6  ! 2014-07  (M. Vancoppenolle, G. Madec, O. Marti) revise coupled interface
-   !!            4.0  ! 2016-06  (L. Brodeau) new unified bulk routine (based on AeroBulk)
+   !!
+   !! The sea ice model SI3 (Sea Ice modelling Integrated Initiative),
+   !!                        aka Sea Ice cube for its nickname
+   !!
+   !!    is originally based on LIM3, developed in Louvain-la-Neuve by: 
+   !!       * Martin Vancoppenolle (UCL-ASTR, Belgium)
+   !!       * Sylvain Bouillon (UCL-ASTR, Belgium)
+   !!       * Miguel Angel Morales Maqueda (NOC-L, UK)
+   !!      thanks to valuable earlier work by
+   !!       * Thierry Fichefet
+   !!       * Hugues Goosse
+   !!      thanks also to the following persons who contributed
+   !!       * Gurvan Madec, Claude Talandier, Christian Ethe (LOCEAN, France)
+   !!       * Xavier Fettweis (UCL-ASTR), Ralph Timmermann (AWI, Germany)
+   !!       * Bill Lipscomb (LANL), Cecilia Bitz (UWa) and Elisabeth Hunke (LANL), USA.
+   !!
+   !! SI3 has been made possible by a handful of persons who met as working group
+   !!     in France and UK:
+   !!    * ...
+   !!======================================================================
+   !! History :  4.0  !  2018     (C. Rousset)      Original code SI3
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -315,7 +327,7 @@
       ENDIF
       !
-	  IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
-	     CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
-	  ENDIF
+      IF( ln_cpl .AND. nn_cats_cpl /= 1 .AND. nn_cats_cpl /= jpl ) THEN
+         CALL ctl_stop( 'STOP', 'par_init: in coupled mode, nn_cats_cpl should be either 1 or jpl' )
+      ENDIF
       !
       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY')
Index: /NEMO/trunk/src/ICE/icetab.F90
===================================================================
--- /NEMO/trunk/src/ICE/icetab.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icetab.F90	(revision 9604)
@@ -4,4 +4,6 @@
    !!   sea-ice : transform 1D (2D) array to a 2D (1D) table
    !!======================================================================
+   !! History :  4.0  !  2018     (C. Rousset)       Original code SI3
+   !!----------------------------------------------------------------------
 #if defined key_si3
    !!----------------------------------------------------------------------
Index: /NEMO/trunk/src/ICE/icethd.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd.F90	(revision 9604)
@@ -4,12 +4,8 @@
    !!   sea-ice : master routine for thermodynamics
    !!======================================================================
-   !! History :  LIM  ! 2000-01 (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1
-   !!            2.0  ! 2002-07 (C. Ethe, G. Madec)  LIM-2 (F90 rewriting)
-   !!            3.0  ! 2005-11 (M. Vancoppenolle)  LIM-3 : Multi-layer thermodynamics + salinity variations
-   !!             -   ! 2007-04 (M. Vancoppenolle) add ice_thd_glohec, ice_thd_con_dh and ice_thd_con_dif
-   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw
-   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas)
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
-   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
+   !! History :  1.0  !  2000-01  (M.A. Morales Maqueda, H. Goosse, T. Fichefet) LIM-1
+   !!            2.0  !  2002-07  (C. Ethe, G. Madec) LIM-2
+   !!            3.0  !  2005-11  (M. Vancoppenolle)  LIM-3
+   !!            4.0  !  2018     (many people)       SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -73,18 +69,19 @@
       !! ** Purpose : This routine manages ice thermodynamics
       !!         
-      !! ** Action : - Initialisation of some variables
-      !!             - Some preliminary computation (oceanic heat flux
-      !!               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_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 active virtual thickness distribution
-      !!             - call ice_thd_da   for lateral ice melt
+      !! ** Action : - computation of oceanic sensible heat flux at the ice base
+      !!                              energy budget in the leads
+      !!                              net fluxes on top of ice and of ocean
+      !!             - selection of grid cells with ice
+      !!                - call ice_thd_zdf  for vertical heat diffusion
+      !!                - call ice_thd_dh   for vertical ice growth and melt
+      !!                - call ice_thd_pnd  for melt ponds
+      !!                - 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 active virtual thickness distribution
+      !!                - call ice_thd_da   for lateral ice melt
       !!             - back to the geographic grid
-      !!             - call ice_thd_rem  for remapping thickness distribution
-      !!             - call ice_thd_do   for ice growth in leads
+      !!                - call ice_thd_rem  for remapping thickness distribution
+      !!                - call ice_thd_do   for ice growth in leads
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt    ! number of iteration
Index: /NEMO/trunk/src/ICE/icethd_da.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_da.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_da.F90	(revision 9604)
@@ -4,5 +4,6 @@
    !!   sea-ice : lateral melting
    !!======================================================================
-   !! History :  4.0  ! 2016-03 (C. Rousset)  original code
+   !! History :  3.7  !  2016-03  (C. Rousset)       Original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!---------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_dh.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_dh.F90	(revision 9604)
@@ -4,9 +4,7 @@
    !!   seaice : thermodynamic growth and melt 
    !!======================================================================
-   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
-   !!                 ! 2005-06 (M. Vancoppenolle) 3D version 
-   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
-   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
-   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes 
+   !! History :       !  2003-05  (M. Vancoppenolle) Original code in 1D
+   !!                 !  2005-06  (M. Vancoppenolle) 3D version 
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -125,8 +123,7 @@
       END DO
       !
-      ! ---------------------------------------------- !
-      ! Available heat for surface and bottom ablation !
-      ! ---------------------------------------------- !
-      !
+      !                       ! ============================================== !
+      !                       ! Available heat for surface and bottom ablation !
+      !                       ! ============================================== !
       SELECT CASE( nice_jules )
       !
Index: /NEMO/trunk/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_do.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_do.F90	(revision 9604)
@@ -4,8 +4,6 @@
    !!   sea-ice: sea ice growth in the leads (open water)  
    !!======================================================================
-   !! 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
+   !! History :       !  2005-12  (M. Vancoppenolle) Original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -58,16 +56,15 @@
       !!  
       !! ** Purpose : Computation of the evolution of the ice thickness and 
-      !!      concentration as a function of the heat balance in the leads.
+      !!              concentration as a function of the heat balance in the leads
       !!       
-      !! ** 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:
+      !! ** Method  : Ice is formed in the open water when ocean looses heat
+      !!              (heat budget of open water is negative) following
+      !!
+      !!       (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.)
+      !!                - exld is the exponent closure rate (=2 default val.)
       !! 
       !! ** Action : - Adjustment of snow and ice thicknesses and heat
@@ -79,6 +76,6 @@
       !!------------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  ::   iter     !   -       -
-      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                          ! local scalars
+      INTEGER  ::   iter             !   -       -
+      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                               ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
@@ -105,12 +102,12 @@
       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) ::   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,nlay_i,jpl) ::   ze_i_2d !: 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)
+      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)
       !!-----------------------------------------------------------------------!
 
@@ -121,8 +118,8 @@
 
       !------------------------------------------------------------------------------!
-      ! 3) Collection thickness of ice formed in leads and polynyas
+      ! 1) Collection thickness of ice formed in leads and polynyas
       !------------------------------------------------------------------------------!    
       ! ht_i_new is the thickness of new ice formed in open water
-      ! ht_i_new can be either prescribed (frazswi = 0) or computed (frazswi = 1)
+      ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)
       ! Frazil ice forms in open water, is transported by wind
       ! accumulates at the edge of the consolidated ice edge
@@ -130,7 +127,4 @@
       ! collection thickness.
 
-      ! Note : the following algorithm currently breaks vectorization
-      ! 
-
       zvrel(:,:) = 0._wp
 
@@ -142,12 +136,10 @@
       IF( ln_frazil ) THEN
          !
-         !--------------------
-         ! Physical constants
-         !--------------------
          ht_i_new(:,:) = 0._wp
          !
-         zhicrit = 0.04 ! frazil ice thickness
+         ! Physical constants
+         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)
+         zsqcd   = 1.0 / SQRT( 1.3 * zcai )                      ! 1/SQRT(airdensity*drag)
          zgamafr = 0.03
          !
@@ -155,8 +147,5 @@
             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
+                  ! -- Wind stress -- !
                   ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
                      &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
@@ -166,22 +155,14 @@
                   ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
 
-                  !---------------------
-                  ! Frazil ice velocity
-                  !---------------------
+                  ! -- 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
+                  ! -- Pack 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
+                  ! -- Relative frazil/pack ice velocity -- !
                   rswitch      = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
                   zvrel2       = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
@@ -189,7 +170,5 @@
                   zvrel(ji,jj) = SQRT( zvrel2 )
 
-                  !---------------------
-                  ! Iterative procedure
-                  !---------------------
+                  ! -- new ice thickness (iterative loop) -- !
                   ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
                      &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
@@ -205,5 +184,5 @@
                   END DO
                   !
-               ENDIF ! end of selection of pixels where ice forms
+               ENDIF
                !
             END DO 
@@ -212,14 +191,12 @@
          CALL lbc_lnk_multi( zvrel, 'T', 1., ht_i_new, 'T', 1.  )
 
-      ENDIF ! End of computation of frazil ice collection thickness
+      ENDIF
 
       !------------------------------------------------------------------------------!
-      ! 4) Identify grid points where new ice forms
+      ! 2) Compute thickness, salinity, enthalpy, age, area and volume of new ice
       !------------------------------------------------------------------------------!
-
-      !-------------------------------------
-      ! Select points for new ice formation
-      !-------------------------------------
       ! This occurs if open water energy budget is negative (cooling) and there is no landfast ice
+
+      ! Identify grid points where new ice forms
       npti = 0   ;   nptidx(:) = 0
       DO jj = 1, jpj
@@ -232,9 +209,5 @@
       END DO
 
-      !------------------------------
       ! Move from 2-D to 1-D vectors
-      !------------------------------
-      ! If ocean gains heat do nothing. Otherwise compute new ice formation
-
       IF ( npti > 0 ) THEN
 
@@ -248,19 +221,17 @@
             END DO
          END DO
-         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead       )
-         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo        )
-         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new    )
-         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel       )
-
-         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d  )
-         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m       )
-
-         !------------------------------------------------------------------------------|
-         ! 2) Convert units for ice internal energy
-         !------------------------------------------------------------------------------|
+         CALL tab_2d_1d( npti, nptidx(1:npti), qlead_1d  (1:npti) , qlead      )
+         CALL tab_2d_1d( npti, nptidx(1:npti), t_bo_1d   (1:npti) , t_bo       )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_opw_1d(1:npti) , sfx_opw    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_opw_1d(1:npti) , wfx_opw    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zh_newice (1:npti) , ht_i_new   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zvrel_1d  (1:npti) , zvrel      )
+
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d(1:npti) , hfx_thd    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_opw_1d(1:npti) , hfx_opw    )
+         CALL tab_2d_1d( npti, nptidx(1:npti), rn_amax_1d(1:npti) , rn_amax_2d )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d    (1:npti) , sss_m      )
+
+         ! Convert units for ice internal energy
          DO jl = 1, jpl
             DO jk = 1, nlay_i               
@@ -272,17 +243,10 @@
             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:npti,:) = v_i_2d(1:npti,:) 
          za_b(1:npti,:) = a_i_2d(1:npti,:)
 
-         !----------------------
-         ! Salinity of new ice 
-         !----------------------
+         ! --- Salinity of new ice --- ! 
          SELECT CASE ( nn_icesal )
          CASE ( 1 )                    ! Sice = constant 
@@ -296,7 +260,5 @@
          END SELECT
 
-         !-------------------------
-         ! Heat content of new ice
-         !-------------------------
+         ! --- Heat content of new ice --- !
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = 1, npti
@@ -307,12 +269,8 @@
          END DO
 
-         !----------------
-         ! Age of new ice
-         !----------------
+         ! --- Age of new ice --- !
          zo_newice(1:npti) = 0._wp
 
-         !-------------------
-         ! Volume of new ice
-         !-------------------
+         ! --- Volume of new ice --- !
          DO ji = 1, npti
 
@@ -351,7 +309,5 @@
          END IF
          
-         !-----------------
-         ! Area of new ice
-         !-----------------
+         ! --- Area of new ice --- !
          DO ji = 1, npti
             za_newice(ji) = zv_newice(ji) / zh_newice(ji)
@@ -359,10 +315,8 @@
 
          !------------------------------------------------------------------------------!
-         ! 6) Redistribute new ice area and volume into ice categories                  !
+         ! 3) Redistribute new ice area and volume into ice categories                  !
          !------------------------------------------------------------------------------!
 
-         !------------------------
-         ! 6.1) lateral ice growth
-         !------------------------
+         ! --- 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
@@ -407,7 +361,5 @@
          END DO
 
-         !------------------------------------------------
-         ! 6.2) bottom ice growth + ice enthalpy remapping
-         !------------------------------------------------
+         ! --- bottom ice growth + ice enthalpy remapping --- !
          DO jl = 1, jpl
 
@@ -436,7 +388,5 @@
          END DO
 
-         !-----------------
-         ! Update salinity
-         !-----------------
+         ! --- Update salinity --- !
          DO jl = 1, jpl
             DO ji = 1, npti
@@ -445,7 +395,5 @@
          END DO
 
-         !------------------------------------------------------------------------------!
-         ! 8) Change units for e_i
-         !------------------------------------------------------------------------------!    
+         ! Change units for e_i
          DO jl = 1, jpl
             DO jk = 1, nlay_i
@@ -453,7 +401,6 @@
             END DO
          END DO
-         !------------------------------------------------------------------------------!
-         ! 7) Change 2D vectors to 1D vectors 
-         !------------------------------------------------------------------------------!
+
+         ! Move 2D vectors to 1D vectors 
          CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) )
Index: /NEMO/trunk/src/ICE/icethd_ent.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_ent.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_ent.F90	(revision 9604)
@@ -5,10 +5,8 @@
    !!                       after vertical growth/melt
    !!======================================================================
-   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D
-   !!                 ! 2005-07 (M. Vancoppenolle) 3D version 
-   !!                 ! 2006-11 (X. Fettweis) Vectorized 
-   !!            3.0  ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code
-   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
-   !!             -   ! 2014-05 (C. Rousset) complete rewriting
+   !! History :       !  2003-05  (M. Vancoppenolle) Original code in 1D
+   !!                 !  2005-07  (M. Vancoppenolle) 3D version 
+   !!            3.6  !  2014-05  (C. Rousset)       New version
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_pnd.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_pnd.F90	(revision 9604)
@@ -2,9 +2,9 @@
    !!======================================================================
    !!                     ***  MODULE  icethd_pnd   ***
-   !!   Melt ponds
+   !!   sea-ice: Melt ponds on top of sea ice
    !!======================================================================
-   !! history :       ! Original code by Daniela Flocco and Adrian Turner
-   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
-   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO4
+   !! history :       !  2012     (O. Lecomte)       Adaptation to LIM from Flocco and Turner
+   !!                 !  2017     (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -74,10 +74,10 @@
       !!
       !! ** Method  :   Melt pond fraction and thickness are prescribed 
-      !!              to non-zero values when t_su = 0C
+      !!                to non-zero values when t_su = 0C
       !!
       !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
       !!                
       !! ** Note   : Coupling with such melt ponds is only radiative
-      !!         Advection, ridging, rafting... are bypassed
+      !!             Advection, ridging, rafting... are bypassed
       !!
       !! ** References : Bush, G.W., and Trump, D.J. (2017)
Index: /NEMO/trunk/src/ICE/icethd_sal.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_sal.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_sal.F90	(revision 9604)
@@ -4,7 +4,7 @@
    !!   sea-ice : computation of salinity variations in the ice
    !!======================================================================
-   !! History :   -   ! 2003-05 (M. Vancoppenolle) UCL-ASTR first coding for LIM3-1D
-   !!            3.0  ! 2005-12 (M. Vancoppenolle) adapted to the 3-D version
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !! History :   -   !  2003-05  (M. Vancoppenolle) UCL-ASTR first coding for LIM3-1D
+   !!            3.0  !  2005-12  (M. Vancoppenolle) adapted to the 3-D version
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!---------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icethd_zdf.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_zdf.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_zdf.F90	(revision 9604)
@@ -2,12 +2,7 @@
    !!======================================================================
    !!                       ***  MODULE icethd_zdf ***
-   !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 
+   !!   sea-ice: master routine for vertical heat diffusion in sea ice 
    !!======================================================================
-   !! 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
+   !! History :  4.0  !  2018     (C. Rousset)      Original code SI3
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -16,6 +11,6 @@
    !!----------------------------------------------------------------------
    !!  ice_thd_zdf      : select the appropriate routine for vertical heat diffusion calculation
-   !!  ice_thd_zdf_BL99 : 
-   !!  ice_thd_zdf_init :
+   !!  ice_thd_zdf_BL99 : heat diffusion from Bitz and Lipscomb 1999 
+   !!  ice_thd_zdf_init : initialization
    !!----------------------------------------------------------------------
    USE dom_oce         ! ocean space and time domain
Index: /NEMO/trunk/src/ICE/icethd_zdf_bl99.F90
===================================================================
--- /NEMO/trunk/src/ICE/icethd_zdf_bl99.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icethd_zdf_bl99.F90	(revision 9604)
@@ -4,10 +4,7 @@
    !!   sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 
    !!======================================================================
-   !! 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
+   !! History :       !  2003-02  (M. Vancoppenolle) original 1D code LIM
+   !!                 !  2005-06  (M. Vancoppenolle) 3d version
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/trunk/src/ICE/iceupdate.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/iceupdate.F90	(revision 9604)
@@ -4,13 +4,5 @@
    !!  Sea-ice :   computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code
-   !!            3.0  ! 2008-03 (C. Tallandier)  surface module
-   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
-   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
-   !!                 !                  + simplification of the ice-ocean stress calculation
-   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
-   !!             -   ! 2012    (D. Iovino) salt flux change
-   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux
-   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
+   !! History :  4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
    !!----------------------------------------------------------------------
 #if defined key_si3
Index: /NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- /NEMO/trunk/src/ICE/icevar.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icevar.F90	(revision 9604)
@@ -2,6 +2,13 @@
    !!======================================================================
    !!                       ***  MODULE icevar ***
-   !!   sea-ice:     Different sets of ice model variables 
-   !!                   how to switch from one to another
+   !!   sea-ice:  series of functions to transform or compute ice variables
+   !!======================================================================
+   !! History :   -   !  2006-01  (M. Vancoppenolle) Original code
+   !!            4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
+   !!----------------------------------------------------------------------
+#if defined key_si3
+   !!----------------------------------------------------------------------
+   !!   'key_si3'                                       SI3 sea-ice model
+   !!----------------------------------------------------------------------
    !!
    !!                 There are three sets of variables
@@ -12,4 +19,5 @@
    !!                        - t_s (jpi,jpj,jpl)
    !!                        - e_i (jpi,jpj,nlay_i,jpl)
+   !!                        - e_s (jpi,jpj,nlay_s,jpl)
    !!                        - sv_i(jpi,jpj,jpl)
    !!                        - oa_i(jpi,jpj,jpl)
@@ -29,13 +37,4 @@
    !!                        - tm_i(jpi,jpj)  mean ice temperature
    !!                        - tm_s(jpi,jpj)  mean snw temperature
-   !!======================================================================
-   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
-   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
-   !!            3.5  ! 2012    (M. Vancoppenolle)  add ice_var_itd
-   !!            3.6  ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
-   !!----------------------------------------------------------------------
-#if defined key_si3
-   !!----------------------------------------------------------------------
-   !!   'key_si3'                                       SI3 sea-ice model
    !!----------------------------------------------------------------------
    !!   ice_var_agg       : integrate variables over layers and categories
Index: /NEMO/trunk/src/ICE/icewri.F90
===================================================================
--- /NEMO/trunk/src/ICE/icewri.F90	(revision 9603)
+++ /NEMO/trunk/src/ICE/icewri.F90	(revision 9604)
@@ -4,4 +4,6 @@
    !!   sea-ice : output ice variables
    !!======================================================================
+   !! History :  4.0  !  2018     (many people)      SI3 [aka Sea Ice cube]
+   !!----------------------------------------------------------------------
 #if defined key_si3
    !!----------------------------------------------------------------------
@@ -42,6 +44,5 @@
    SUBROUTINE ice_wri( kt )
       !!-------------------------------------------------------------------
-      !!  This routine computes the average of some variables and write it
-      !!  on the ouput files.
+      !!  This routine ouputs some (most?) of the sea ice fields
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! time-step
@@ -60,9 +61,6 @@
       IF( ln_timing )   CALL timing_start('icewri')
 
-      !----------------------------------------
-      ! Brine volume, switches, missing values
-      !----------------------------------------
-
-      CALL ice_var_bv      ! brine volume
+      ! brine volume
+      CALL ice_var_bv
 
       ! tresholds for outputs
