Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg	(revision 8517)
@@ -1,48 +1,64 @@
 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-!! LIM3 configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref
-!!              1 - Generic parameters                 (namicerun)
-!!              2 - Diagnostics                        (namicediag)
-!!              3 - Ice initialization                 (namiceini)
-!!              4 - Ice discretization                 (namiceitd)
-!!              5 - Ice dynamics and transport         (namicedyn)
-!!              6 - Ice thermodynamics                 (namicethd)
-!!              7 - Ice salinity                       (namicesal)
-!!              8 - Ice mechanical redistribution      (namiceitdme)
-!!              9 - Ice/snow albedos                   (namicealb)
-!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!! 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)
+!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+!
 !------------------------------------------------------------------------------
-&namicerun     !   Generic parameters
+&namice_run     !   Generic parameters
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicediag    !   Diagnostics
+&namice_itd     !   Ice discretization
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namiceini     !   Ice initialization
+&namice_dyn     !   Ice dynamics
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namiceitd     !   Ice discretization
+&namice_rdgrft  !   Ice ridging/rafting
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicedyn     !   Ice dynamics and transport
+&namice_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicethd     !   Ice thermodynamics
+&namice_adv     !   Ice advection
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namicesal     !   Ice salinity
+&namice_thd     !   Ice thermodynamics
 !------------------------------------------------------------------------------
 /
 !------------------------------------------------------------------------------
-&namiceitdme   !   Ice mechanical redistribution (ridging and rafting)
+&namice_sal     !   Ice salinity
 !------------------------------------------------------------------------------
 /
-!-----------------------------------------------------------------------
-&namicealb     !   albedo parameters
-!-----------------------------------------------------------------------
+!------------------------------------------------------------------------------
+&namicemp      !   Melt ponds
+!------------------------------------------------------------------------------
 /
+!------------------------------------------------------------------------------
+&namice_ini     !   Ice initialization
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namice_alb     !   albedo parameters
+!------------------------------------------------------------------------------
+/
+!------------------------------------------------------------------------------
+&namice_dia     !   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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8517)
@@ -18,11 +18,13 @@
 &namice_run     !   Generic parameters
 !------------------------------------------------------------------------------
-   jpl              =    5            !  number of ice  categories
-   nlay_i           =    2            !  number of ice  layers
-   nlay_s           =    1            !  number of snow layers (only 1 is working)
-   nn_monocat       =    0            !  virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0)
+   jpl              =   5             !  number of ice  categories
+   nlay_i           =   2             !  number of ice  layers
+   nlay_s           =   1             !  number of snow layers (only 1 is working)
+   nn_monocat       =   0             !  virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0)
                                       !     2: simple piling instead of ridging    --- temporary option
                                       !     3: activate G(he) only                 --- temporary option
                                       !     4: activate extra lateral melting only --- temporary option
+   ln_icedyn        = .true.          !  ice dynamics (T) or not (F)
+   ln_icethd        = .true.          !  ice thermo   (T) or not (F)
    rn_amax_n        =   0.997         !  maximum tolerated ice concentration NH
    rn_amax_s        =   0.997         !  maximum tolerated ice concentration SH
@@ -41,11 +43,9 @@
 &namice_dyn     !   Ice dynamics
 !------------------------------------------------------------------------------
-   ln_icedyn        = .true.          !  ice dynamics (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-      nn_icedyn     =   2             !     switch for ice dynamics   
-                                      !        2: total
-                                      !        1: advection only (no diffusion, no ridging/rafting)
-                                      !        0: advection only (as 1 but with prescribed velocity, bypass rheology)
-         rn_uice    =   0.00001       !           prescribed ice u-velocity
-         rn_vice    =  -0.00001       !           prescribed ice v-velocity
+   ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection)
+   ln_dynADV        = .false.         !  dyn.: only advection w prescribed vel.(rn_uvice + advection)
+      rn_uice       =   0.00001       !        prescribed ice u-velocity
+      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 (-)
@@ -53,8 +53,8 @@
       rn_gamma      =   0.15          !     fraction of ocean depth that ice must reach to initiate landfast
                                       !        recommended range: [0.1 ; 0.25]
-      rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact (N/m2)                 
+      rn_icebfr     =  10.            !     maximum bottom stress per unit area of contact [N/m2]                 
                                       !        a very large value ensures ice velocity=0 even with a small contact area
                                       !        recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2)
-      rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction (s-1)                 
+      rn_lfrelax    =   1.e-5         !     relaxation time scale to reach static friction [s-1]
 /
 !------------------------------------------------------------------------------
@@ -63,5 +63,5 @@
           ! -- ice_rdgrft_strength -- !
    ln_str_H79       = .true.          !  ice strength param.: Hibler_79   => P = pstar*<h>*exp(-c_rhg*A)                      
-      rn_pstar      =   2.0e+04       !     ice strength thickness parameter (N/m2) 
+      rn_pstar      =   2.0e+04       !     ice strength thickness parameter [N/m2]
       rn_crhg       =   20.0          !     ice strength conc. parameter (-)
    ln_str_R75       = .false.         !  ice strength param.: Rothrock_75 => P = Cf*coeff*integral(wr.h^2)    
@@ -74,11 +74,11 @@
    ln_partf_exp     = .true.          !  Exponential ridging participation function (Lipscomb, 2007)
       rn_astar      =   0.03          !     exponential measure of ridging ice fraction [set to 0.05 if hstar=100]
-   ln_ridging       = .true.          !  ridging activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-      rn_hstar      =  25.0           !     determines the maximum thickness of ridged ice (m) (Hibler, 1980)
+   ln_ridging       = .true.          !  ridging activated (T) or not (F)
+      rn_hstar      =  25.0           !     determines the maximum thickness of ridged ice [m] (Hibler, 1980)
       rn_porordg    =   0.3           !     porosity of newly ridged ice (Lepparanta et al., 1995)
       rn_fsnwrdg    =   0.5           !     snow volume fraction that survives in ridging
       rn_fpndrdg    =   1.0           !     pond fraction that survives in ridging (small a priori)
-   ln_rafting       = .true.          !  rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
-      rn_hraft      =   0.75          !     threshold thickness for rafting (m)
+   ln_rafting       = .true.          !  rafting activated (T) or not (F)
+      rn_hraft      =   0.75          !     threshold thickness for rafting [m]
       rn_craft      =   5.0           !     squeezing coefficient used in the rafting function
       rn_fsnwrft    =   0.5           !     snow volume fraction that survives in rafting
@@ -89,5 +89,5 @@
 !------------------------------------------------------------------------------
    ln_rhg_EVP       = .true.          !  EVP rheology
-      rn_creepl     =   1.0e-12       !     creep limit (s-1)
+      rn_creepl     =   1.0e-12       !     creep limit [1/s]
       rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve          
       nn_nevp       = 120             !     number of EVP subcycles                             
@@ -105,7 +105,6 @@
 &namice_thd     !   Ice thermodynamics
 !------------------------------------------------------------------------------
-   ln_icethd        = .true.            !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
                      ! -- icethd_dif -- !
-   rn_kappa_i       =   1.0             !  radiation attenuation coefficient in sea ice (m-1)
+   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)
@@ -114,9 +113,9 @@
                                         !     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) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+   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) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+   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)
@@ -127,5 +126,5 @@
                                         !         10 vs 8m = -20% melting
                      ! -- icethd_lac -- !
-   ln_icedO         = .true.            !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+   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)
@@ -148,5 +147,5 @@
 !------------------------------------------------------------------------------
                     ! -- icethd_sal -- !
-   ln_icedS       = .true.             !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
+   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)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8517)
@@ -429,17 +429,4 @@
 /
 !-----------------------------------------------------------------------
-&namsbc_alb    !   albedo parameters
-!-----------------------------------------------------------------------
-   nn_ice_alb   =    1   !  parameterization of ice/snow albedo
-                         !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
-                         !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004),
-                         !        giving cloud-sky albedo
-   rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky)
-   rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' )
-   rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' )
-   rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' )
-   rn_alb_dpnd  =  0.27  !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' )
-/
-!-----------------------------------------------------------------------
 &namsbc_wave   ! External fields from wave model                        (ln_wave=T)
 !-----------------------------------------------------------------------
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8517)
@@ -157,4 +157,6 @@
    INTEGER           , PUBLIC ::   nlay_i          !: number of ice  layers 
    INTEGER           , PUBLIC ::   nlay_s          !: number of snow layers 
+   LOGICAL           , PUBLIC ::   ln_icedyn       !: flag for ice dynamics (T) or not (F)
+   LOGICAL           , PUBLIC ::   ln_icethd       !: flag for ice thermo   (T) or not (F)
    REAL(wp)          , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere
    REAL(wp)          , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere
@@ -168,14 +170,13 @@
    
    !                                     !!** ice-dynamics namelist (namice_dyn) **
-   LOGICAL , PUBLIC ::   ln_icedyn        !: flag for ice dynamics (T) or not (F)
-   INTEGER , PUBLIC ::   nn_icedyn        !: flag for ice dynamics
-   REAL(wp), PUBLIC ::   rn_uice          !: prescribed u-vel (case nn_icedyn=0)
-   REAL(wp), PUBLIC ::   rn_vice          !: prescribed v-vel (case nn_icedyn=0)
+   LOGICAL,  PUBLIC ::   ln_dynFULL       !:  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   LOGICAL,  PUBLIC ::   ln_dynRHGADV     !:  dyn.: no ridge/raft & no corrections  (rheology + advection)
+   LOGICAL,  PUBLIC ::   ln_dynADV        !:  dyn.: only advection w prescribed vel.(rn_uvice + advection)
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
    LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F) 
-   REAL(wp), PUBLIC ::   rn_gamma         !: fraction of ocean depth that ice must reach to initiate landfast ice
-   REAL(wp), PUBLIC ::   rn_icebfr        !: maximum bottom stress per unit area of contact (landfast ice) 
-   REAL(wp), PUBLIC ::   rn_lfrelax       !: relaxation time scale (s-1) to reach static friction (landfast ice) 
+   REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice
+   REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (landfast ice) 
+   REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice) 
    !
    !                                     !!** ice-rdige/raft namelist (namice_rdgrft) **
@@ -199,5 +200,4 @@
    !
    !                                     !!** ice-thermodynamics namelist (namice_thd) **
-   LOGICAL , PUBLIC ::   ln_icethd        !: flag for ice thermo (T) or not (F)
                                           ! -- icethd_dif -- !
    REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
@@ -266,9 +266,9 @@
    REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s 
    REAL(wp), PUBLIC ::   rswitch          !: switch for the presence of ice (1) or not (0)
+   REAL(wp), PUBLIC ::   rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft   !: conservation diagnostics
    REAL(wp), PUBLIC, PARAMETER ::   epsi06 = 1.e-06_wp  !: small number 
    REAL(wp), PUBLIC, PARAMETER ::   epsi10 = 1.e-10_wp  !: small number 
    REAL(wp), PUBLIC, PARAMETER ::   epsi20 = 1.e-20_wp  !: small number 
-   !
-   LOGICAL , PUBLIC ::   l_piling         !: =T simple conservative piling, comparable with LIM2
+
 
    !                                     !!** define arrays
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8517)
@@ -37,4 +37,9 @@
    PUBLIC   ice_adv_init   ! called by icestp
 
+   INTEGER ::              nice_adv   ! choice of the type of advection scheme
+   !                                        ! associated indices:
+   INTEGER, PARAMETER ::   np_advPRA = 1   ! Prather scheme
+   INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme
+
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
@@ -61,14 +66,4 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt   ! number of iteration
-      !
-      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
-      INTEGER  ::   initad                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zcfl , zusnit           !   -      -
-      CHARACTER(len=80) :: cltmp
-      !
-      REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
-      REAL(wp) ::    zdv
-      REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold 
-      REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('iceadv')
@@ -80,69 +75,31 @@
       ENDIF
       
-      CALL ice_var_agg( 1 ) ! integrated values + ato_i
-
       ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      
-      ! store old values for diag
-      zviold (:,:,:) = v_i(:,:,:)
-      zvsold (:,:,:) = v_s(:,:,:)
-      zsmvold(:,:)   = SUM( smv_i(:,:,:), dim=3 )
-      zeiold (:,:)   = et_i(:,:)
-      zesold (:,:)   = et_s(:,:)
-
-      ! Thickness correction init.
-      zatold(:,:) = at_i
-      WHERE( a_i(:,:,:) >= epsi20 )
-         ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
-         ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:)
-      ELSEWHERE
-         ht_i(:,:,:) = 0._wp
-         ht_s(:,:,:) = 0._wp         
-      END WHERE
-         
-      ! Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick
-      zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
-      DO jl = 1, jpl
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
-               zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) )
-            END DO
-         END DO
-      END DO
-      CALL lbc_lnk( zhimax(:,:,:), 'T', 1. )
-      
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+                     
       !----------
       ! Advection
       !----------
-      IF( ln_adv_UMx ) THEN                    !-- ULTIMATE-MACHO scheme
+      SELECT CASE( nice_adv )
+      !                                !-----------------------!
+      CASE( np_advUMx )                ! ULTIMATE-MACHO scheme !
+         !                             !-----------------------!
          CALL ice_adv_umx( kt, u_ice, v_ice,  &
             &              ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
-         
-      ELSEIF( ln_adv_Pra ) THEN                !-- PRATHER scheme
+      !                                !-----------------------!
+      CASE( np_advPRA )                ! PRATHER scheme        !
+         !                             !-----------------------!
          CALL ice_adv_prather( kt, u_ice, v_ice,  &
             &                  ato_i, v_i, v_s, smv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
-         
-      ENDIF
+      END SELECT
 
-      ! total ice fraction
-      at_i(:,:) = a_i(:,:,1)
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-      
       !------------
       ! diagnostics
       !------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            diag_trp_ei (ji,jj) = ( SUM( e_i  (ji,jj,1:nlay_i,:) ) -  zeiold(ji,jj) ) * r1_rdtice
-            diag_trp_es (ji,jj) = ( SUM( e_s  (ji,jj,1:nlay_s,:) ) -  zesold(ji,jj) ) * r1_rdtice
-            diag_trp_smv(ji,jj) = ( SUM( smv_i(ji,jj,:)          ) - zsmvold(ji,jj) ) * r1_rdtice
-            diag_trp_vi (ji,jj) =   SUM(   v_i(ji,jj,:)            -  zviold(ji,jj,:) ) * r1_rdtice
-            diag_trp_vs (ji,jj) =   SUM(   v_s(ji,jj,:)            -  zvsold(ji,jj,:) ) * r1_rdtice
-         END DO
-      END DO
+      diag_trp_ei (:,:) = SUM(SUM( e_i  (:,:,1:nlay_i,:) - e_i_b  (:,:,1:nlay_i,:), dim=4 ), dim=3 ) * r1_rdtice
+      diag_trp_es (:,:) = SUM(SUM( e_s  (:,:,1:nlay_s,:) - e_s_b  (:,:,1:nlay_s,:), dim=4 ), dim=3 ) * r1_rdtice
+      diag_trp_smv(:,:) = SUM(     smv_i(:,:,:)          - smv_i_b(:,:,:)                  , dim=3 ) * r1_rdtice
+      diag_trp_vi (:,:) = SUM(     v_i  (:,:,:)          - v_i_b  (:,:,:)                  , dim=3 ) * r1_rdtice
+      diag_trp_vs (:,:) = SUM(     v_s  (:,:,:)          - v_s_b  (:,:,:)                  , dim=3 ) * r1_rdtice
       IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi * rday  )         ! ice volume transport
       IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs * rday  )         ! snw volume transport
@@ -150,79 +107,9 @@
       IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei         )         ! advected ice enthalpy (W/m2)
       IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2)
-      
-      !--------------------------------------
-      ! Thickness correction in case too high
-      !--------------------------------------
-      IF( nn_icedyn == 2 ) THEN
-         !
-         CALL ice_var_zapsmall                       !-- zap small areas
-         !
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to zhimax
-                     !
-                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
-                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl)
-                     zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)  
-                     !
-                     IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
-                        & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
-                        a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl)
-                        ht_i(ji,jj,jl) =   v_i(ji,jj,jl) / a_i(ji,jj,jl)
-                     ENDIF
-                     !
-                  ENDIF
-               END DO
-            END DO
-         END DO
-                  
-         WHERE( ht_i(:,:,jpl) > hi_max(jpl) )        !-- bound ht_i to hi_max (99 m)
-            ht_i(:,:,jpl) = hi_max(jpl)
-            a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl)
-         END WHERE
-
-         IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
-            WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-         ENDIF
-         !
-      ENDIF
-         
-      !------------------------------------------------------------
-      ! Impose a_i < amax if no ridging/rafting or in mono-category
-      !------------------------------------------------------------
-      IF( l_piling ) THEN                            !-- simple conservative piling, comparable with 1-cat models
-         at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
-         DO jl = 1, jpl
-            WHERE( at_i(:,:) > epsi20 )
-               a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
-            END WHERE
-         END DO
-      ENDIF
-      
-      ! agglomerate variables 
-      vt_i(:,:) = SUM( v_i(:,:,:), dim=3 )
-      vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
-      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
-
-      ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i)
-      IF ( nn_pnd_scheme > 0 ) THEN
-          at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 )
-          vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 )
-      ENDIF
-      ! END MP 2016
-      
-      ! open water = 1 if at_i=0 
-      WHERE( at_i == 0._wp ) ato_i = 1._wp 
-      
-      ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceadv', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-        
-      ! --------------
-      ! control prints
-      ! --------------
-      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ' )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('iceadv')
+               
+      ! controls
+      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'iceadv', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icectl      )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn & trp - ')                       ! prints
+      IF( nn_timing == 1 )   CALL timing_stop ('iceadv')                                                             ! timing
       !
    END SUBROUTINE ice_adv
@@ -241,5 +128,5 @@
       !! ** input   :   Namelist namice_adv
       !!-------------------------------------------------------------------
-      INTEGER ::   ios   ! Local integer output status for namelist read
+      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
       NAMELIST/namice_adv/ ln_adv_Pra, ln_adv_UMx, nn_UMx
@@ -260,13 +147,14 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_adv'
-         WRITE(numout,*) '      advection scheme for ice transport (limtrp)'
-         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,*) '      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
       !
-      IF ( ( ln_adv_Pra .AND. ln_adv_UMx ) .OR. ( .NOT.ln_adv_Pra .AND. .NOT.ln_adv_UMx ) ) THEN
-         CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' )
-      ENDIF
+      !                             !== set the choice of ice advection ==!
+      ioptio = 0 
+      IF( ln_adv_Pra ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advPRA    ;   ENDIF
+      IF( ln_adv_UMx ) THEN   ;   ioptio = ioptio + 1   ;   nice_adv = np_advUMx    ;   ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' )
       !
    END SUBROUTINE ice_adv_init
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8517)
@@ -321,5 +321,5 @@
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namice_alb in reference namelist', lwp )
 
-      REWIND( numnam_ice_cfg )              ! Namelist namsbc_alb in configuration namelist : Albedo parameters
+      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 )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90	(revision 8517)
@@ -53,8 +53,10 @@
       !
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc
+      REAL(wp) ::   zsal, zzc
       REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag
       !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )   CALL timing_start('icecor')
+      ! controls
+      IF( nn_timing == 1 )   CALL timing_start('icecor')                                                             ! timing
+      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
@@ -63,7 +65,4 @@
          WRITE(numout,*) '~~~~~~~'
       ENDIF
-      !                             !--- conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      !
       !
       IF( kn == 2 ) THEN
@@ -175,12 +174,9 @@
       END SELECT
       !
-      !                                !--- conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      !
-      !                                !--- control prints
-      IF( ln_ctl )                    CALL ice_prt3D( 'icecor' )
-      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('icecor')
+      ! controls
+      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints
+      IF( ln_icectl .AND. kn == 2 )   CALL ice_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                   ! prints
+      IF( nn_timing == 1 )   CALL timing_stop ('icecor')                                                             ! timing
       !
    END SUBROUTINE ice_cor
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8517)
@@ -47,5 +47,5 @@
 CONTAINS
 
-   SUBROUTINE ice_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
+   SUBROUTINE ice_cons_hsm( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft )
       !!----------------------------------------------------------------------
       !!                       ***  ROUTINE ice_cons_hsm ***
@@ -56,5 +56,5 @@
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to
+      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to
       !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
@@ -63,65 +63,65 @@
       INTEGER         , INTENT(in)    ::   icount        ! called at: =0 the begining of the routine, =1  the end
       CHARACTER(len=*), INTENT(in)    ::   cd_routine    ! name of the routine
-      REAL(wp)        , INTENT(inout) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b   ! ????
-      !!
-      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
-      REAL(wp)                        :: zvmin, zamin, zamax 
-      REAL(wp)                        :: zvtrp, zetrp
-      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
-      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
+      REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
+      !!
+      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft
+      REAL(wp) ::   zvmin, zamin, zamax 
+      REAL(wp) ::   zvtrp, zetrp
+      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill
+      REAL(wp), PARAMETER ::   zconv = 1.e-9 ! convert W to GW and kg to Mt
       !!----------------------------------------------------------------------
       !
       IF( icount == 0 ) THEN
+         !                          ! water flux
+         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
+            &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
+            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    &
+            &                  ) * e1e2t(:,:) ) * zconv
+         !
          !                          ! salt flux
-         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
-            &                ) *  e1e2t(:,:) ) * zconv 
-         !
-         !                          ! water flux
-         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                  wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                  wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    &
-            &                ) * e1e2t(:,:) ) * zconv
+         pdiag_fs = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
+            &                  ) *  e1e2t(:,:) ) * zconv 
          !
          !                          ! heat flux
-         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                ) *  e1e2t(:,:) ) * zconv
-
-         zvi_b  = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )
-
-         zsmv_b = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv )
-
-         zei_b  = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
-            &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv
+         pdiag_ft = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &                  ) *  e1e2t(:,:) ) * zconv
+
+         pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )
+
+         pdiag_s = glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t * zconv )
+
+         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
+            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv
 
       ELSEIF( icount == 1 ) THEN
+
+         ! water flux
+         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
+            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
+            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
+            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv
 
          ! salt flux
          zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
             &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
-            &              ) * e1e2t(:,:) ) * zconv - zfs_b
-
-         ! water flux
-         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
-            &              ) * e1e2t(:,:) ) * zconv - zfw_b
+            &              ) * e1e2t(:,:) ) * zconv - pdiag_fs
 
          ! heat flux
          zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
             &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &              ) * e1e2t(:,:) ) * zconv - zft_b
+            &              ) * e1e2t(:,:) ) * zconv - pdiag_ft
  
          ! outputs
-         zvi  = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  &
-            &       - zvi_b ) * r1_rdtice - zfw ) * rday
-
-         zsmv = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  &
-            &       - zsmv_b ) * r1_rdtice + zfs ) * rday
-
-         zei  = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   &
+         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  &
+            &     - pdiag_v ) * r1_rdtice - zfv ) * rday
+
+         zs = ( ( glob_sum( SUM( smv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  &
+            &     - pdiag_s ) * r1_rdtice + zfs ) * rday
+
+         zt = ( glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )   &
             &                + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   &
-            &   - zei_b ) * r1_rdtice + zft
+            &   - pdiag_t ) * r1_rdtice + zft
 
          ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
@@ -137,21 +137,18 @@
          zv_sill = zarea * 2.5e-5
          zs_sill = zarea * 25.e-5
-         zh_sill = zarea * 10.e-5
+         zt_sill = zarea * 10.e-5
 
          IF(lwp) THEN
-            IF ( ABS( zvi  ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zvi
-            IF ( ABS( zsmv ) > zs_sill ) WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zsmv
-            IF ( ABS( zei  ) > zh_sill ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zei
+            IF ( ABS( zv   ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv
+            IF ( ABS( zs   ) > zs_sill )   WRITE(numout,*) 'violation saline [psu*Mt/day] (',cd_routine,') = ',zs
+            IF ( ABS( zt   ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt
             IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'iceadv' ) THEN
-                                         WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
-                                         WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
-            ENDIF
-            IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
-            IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. &
-               &                         cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' ) THEN
-                                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
-            IF (     zamax   > 1._wp   ) WRITE(numout,*) 'violation a_i>1               (',cd_routine,') = ',zamax
-            ENDIF
-            IF (      zamin  < -epsi10 ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
+                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
+            ENDIF
+            IF ( zvmin < -epsi10 )         WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
+            IF ( zamax > MAX(rn_amax_n,rn_amax_s) + epsi10 .AND. cd_routine /= 'iceadv' .AND. cd_routine /= 'icerdgrft' )  &
+               &                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
+            IF ( zamin < -epsi10 )         WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
          ENDIF
          !
@@ -169,5 +166,5 @@
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to
+      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to
       !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
@@ -176,7 +173,13 @@
       CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
       REAL(wp)                        :: zhfx, zsfx, zvfx
-      REAL(wp)                        :: zarea, zv_sill, zs_sill, zh_sill
+      REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill
       REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
       !!----------------------------------------------------------------------
+
+      ! water flux
+      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday
+
+      ! salt flux
+      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday
 
       ! heat flux
@@ -184,8 +187,4 @@
       !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check)
          &              ) * e1e2t ) * zconv
-      ! salt flux
-      zsfx  = glob_sum( ( sfx + diag_smvi ) * e1e2t ) * zconv * rday
-      ! water flux
-      zvfx  = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday
 
       ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
@@ -193,10 +192,10 @@
       zv_sill = zarea * 2.5e-5
       zs_sill = zarea * 25.e-5
-      zh_sill = zarea * 10.e-5
+      zt_sill = zarea * 10.e-5
 
       IF(lwp) THEN
-         IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx    [Mt/day]       (',cd_routine,')  = ',(zvfx)
-         IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx    [psu*Mt/day]   (',cd_routine,')  = ',(zsfx)
-         IF( ABS( zhfx ) > zh_sill )   WRITE(numout,*) 'violation hfx    [GW]           (',cd_routine,')  = ',(zhfx)
+         IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx  [Mt/day]       (',cd_routine,') = ',zvfx
+         IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx  [psu*Mt/day]   (',cd_routine,') = ',zsfx
+         IF( ABS( zhfx ) > zt_sill )   WRITE(numout,*) 'violation hfx  [GW]           (',cd_routine,') = ',zhfx
       ENDIF
       !
@@ -215,6 +214,6 @@
       INTEGER  ::   ialert_id         ! number of the current alert
       REAL(wp) ::   ztmelts           ! ice layer melting point
-      CHARACTER (len=30), DIMENSION(20)      ::   cl_alname   ! name of alert
-      INTEGER           , DIMENSION(20)      ::   inb_alp     ! number of alerts positive
+      CHARACTER (len=30), DIMENSION(20) ::   cl_alname   ! name of alert
+      INTEGER           , DIMENSION(20) ::   inb_alp     ! number of alerts positive
       !!-------------------------------------------------------------------
 
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8517)
@@ -34,11 +34,15 @@
    PUBLIC   ice_dyn_init   ! called by icestp.F90
    
-   INTEGER ::              nice_dyn   ! choice of the type of advection scheme
+   INTEGER ::              nice_dyn   ! choice of the type of dynamics
    !                                        ! associated indices:
-   INTEGER, PARAMETER ::   np_dynNO   = 0   ! no ice dynamics and ice advection
-   INTEGER, PARAMETER ::   np_dynFULL = 1   ! full ice dynamics  (rheology + advection + ridging/rafting + correction)
-   INTEGER, PARAMETER ::   np_dyn     = 2   ! no ridging/rafting (rheology + advection                   + correction)
-   INTEGER, PARAMETER ::   np_dynPURE = 3   ! pure dynamics      (rheology + advection) 
-
+   INTEGER, PARAMETER ::   np_dynFULL    = 1   ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   INTEGER, PARAMETER ::   np_dynRHGADV1 = 2   ! no ridging/rafting              (rheology + advection                   + correction)
+   INTEGER, PARAMETER ::   np_dynRHGADV2 = 3   ! pure dynamics                   (rheology + advection) 
+   INTEGER, PARAMETER ::   np_dynADV     = 4   ! only advection w prescribed vel.(rn_uvice + advection)
+   !
+   ! ** namelist (namdyn) **
+   REAL(wp) ::   rn_uice          ! prescribed u-vel (case np_dynADV)
+   REAL(wp) ::   rn_vice          ! prescribed v-vel (case np_dynADV)
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -61,5 +65,6 @@
       INTEGER, INTENT(in) ::   kt     ! ice time step
       !!
-      INTEGER  ::   jl   ! dummy loop indices
+      INTEGER ::   ji, jj, jl         ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhmax
       !!--------------------------------------------------------------------
       !
@@ -72,7 +77,6 @@
       ENDIF
 
-      CALL ice_var_agg(1)           ! -- aggregate ice categories
       !                      
-      IF( ln_landfast ) THEN        ! -- Landfast ice parameterization: define max bottom friction
+      IF( ln_landfast ) THEN            !-- Landfast ice parameterization: define max bottom friction
          tau_icebfr(:,:) = 0._wp
          DO jl = 1, jpl
@@ -81,26 +85,41 @@
       ENDIF
 
-      SELECT CASE( nice_dyn )       ! -- Set which dynamics is running
+      zhmax(:,:,:) = ht_i_b(:,:,:)      !-- Record max of the surrounding 9-pts ice thick. (for CALL Hbig)
+      DO jl = 1, jpl
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+!!gm use of MAXVAL here is very probably less efficient than expending the 9 values
+               zhmax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i_b(ji-1:ji+1,jj-1:jj+1,jl) ) )
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( zhmax(:,:,:), 'T', 1. )
+      !
+      !
+      SELECT CASE( nice_dyn )           !-- Set which dynamics is running
 
       CASE ( np_dynFULL )          !==  all dynamical processes  ==!
-         CALL ice_rhg   ( kt )          ! -- rheology  
-         CALL ice_adv   ( kt )          ! -- advection of ice
-         CALL ice_rdgrft( kt )          ! -- ridging/rafting 
-         CALL ice_cor   ( kt , 1 )      ! -- Corrections
-
-      CASE ( np_dyn )              !==  pure dynamics only ==!   (no ridging/rafting)   (nono cat. case 2)
-         CALL ice_rhg   ( kt )          ! -- rheology  
-         CALL ice_adv   ( kt )          ! -- advection of ice
-         CALL ice_cor   ( kt , 1 )      ! -- Corrections
-
-      CASE ( np_dynPURE )          !==  pure dynamics only ==!   (nn_icedyn= 1 )
-         CALL ice_rhg   ( kt )          ! -- rheology  
-         CALL ice_adv   ( kt )          ! -- advection of ice
-
-      CASE ( np_dynNO )            !==  prescribed ice velocities ==!   (nn_icedyn= 0 )
+         CALL ice_rhg   ( kt )                            ! -- rheology  
+         CALL ice_adv   ( kt )   ;   CALL Hbig( zhmax )   ! -- advection of ice + correction on ice thickness
+         CALL ice_rdgrft( kt )                            ! -- ridging/rafting 
+         CALL ice_cor   ( kt , 1 )                        ! -- Corrections
+
+      CASE ( np_dynRHGADV1 )       !==  no ridge/raft ==!   (mono cat. case 2)
+         CALL ice_rhg   ( kt )                            ! -- rheology  
+         CALL ice_adv   ( kt )                            ! -- advection of ice
+         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting)
+         CALL ice_cor   ( kt , 1 )                        ! -- Corrections
+
+      CASE ( np_dynRHGADV2 )       !==  no ridge/raft & no corrections ==!
+         CALL ice_rhg   ( kt )                            ! -- rheology  
+         CALL ice_adv   ( kt )                            ! -- advection of ice
+         CALL Hpiling                                     ! -- simple pile-up (replaces ridging/rafting)
+
+      CASE ( np_dynADV )           !==  pure advection ==!   (prescribed velocities)
          u_ice(:,:) = rn_uice * umask(:,:,1)
          v_ice(:,:) = rn_vice * vmask(:,:,1)
          !!CALL RANDOM_NUMBER(u_ice(:,:))
          !!CALL RANDOM_NUMBER(v_ice(:,:))
+         CALL ice_adv   ( kt )                            ! -- advection of ice
 
       END SELECT
@@ -110,4 +129,73 @@
    END SUBROUTINE ice_dyn
 
+   SUBROUTINE Hbig( phmax )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hbig  ***
+      !!
+      !! ** Purpose : Thickness correction in case advection scheme creates
+      !!              abnormally tick ice
+      !!
+      !! ** Method  : 1- check whether ice thickness resulting from advection is
+      !!                 larger than the surrounding 9-points before advection
+      !!                 and reduce it if a) divergence or b) convergence & at_i>0.8
+      !!              2- bound ice thickness with hi_max (99m)
+      !!
+      !! ** input   : Max thickness of the surrounding 9-points
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(in) ::   phmax   ! max ice thick from surrounding 9-pts
+      !
+      INTEGER  ::   ji, jj, jl         ! dummy loop indices
+      REAL(wp) ::   zh, zdv
+      !!-------------------------------------------------------------------
+      !
+      CALL ice_var_zapsmall                       !-- zap small areas
+      !
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( v_i(ji,jj,jl) > 0._wp ) THEN  !-- bound to hmax
+                  !
+                  zh  = v_i (ji,jj,jl) / a_i(ji,jj,jl)
+                  zdv = v_i(ji,jj,jl) - v_i_b(ji,jj,jl)  
+                  !
+                  IF ( ( zdv >  0.0 .AND. zh > phmax(ji,jj,jl) .AND. at_i_b(ji,jj) < 0.80 ) .OR. &
+                     & ( zdv <= 0.0 .AND. zh > phmax(ji,jj,jl) ) ) THEN
+                     a_i (ji,jj,jl) = v_i(ji,jj,jl) / MIN( phmax(ji,jj,jl), hi_max(jpl) )   !-- bound ht_i to hi_max (99 m)
+                  ENDIF
+                  !
+               ENDIF
+            END DO
+         END DO
+      END DO
+            
+      IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
+         WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
+      ENDIF
+      !
+   END SUBROUTINE Hbig
+
+   SUBROUTINE Hpiling
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hpiling  ***
+      !!
+      !! ** Purpose : Simple conservative piling comparable with 1-cat models
+      !!
+      !! ** Method  : pile-up ice when no ridging/rafting
+      !!
+      !! ** input   : a_i
+      !!-------------------------------------------------------------------
+      INTEGER ::   jl         ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+      CALL ice_var_zapsmall                       !-- zap small areas
+      !
+      at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+      DO jl = 1, jpl
+         WHERE( at_i(:,:) > epsi20 )
+            a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
+         END WHERE
+      END DO
+      !
+   END SUBROUTINE Hpiling
 
    SUBROUTINE ice_dyn_init
@@ -123,7 +211,7 @@
       !! ** input   :   Namelist namice_dyn
       !!-------------------------------------------------------------------
-      INTEGER ::   ios   ! Local integer output status for namelist read
-      !!
-      NAMELIST/namice_dyn/ ln_icedyn  , nn_icedyn, rn_uice  , rn_vice ,    &
+      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
@@ -144,9 +232,8 @@
          WRITE(numout,*) '~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_dyn'
-         WRITE(numout,*) '      Ice dynamics       (T) or not (F)                         ln_icedyn  = ', ln_icedyn
-         WRITE(numout,*) '         associated switch                                      nn_icedyn  = ', nn_icedyn
-         WRITE(numout,*) '            =2 all processes (default option)'
-         WRITE(numout,*) '            =1 advection only (no ridging/rafting)'
-         WRITE(numout,*) '            =0 advection only with prescribed velocity given by '
+         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
+         WRITE(numout,*) '      Advection only         (rn_uvice + adv)                 ln_dynADV    = ', ln_dynADV
+         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
@@ -158,16 +245,15 @@
       ENDIF
       !                             !== set the choice of ice dynamics ==!
-      SELECT CASE( nn_icedyn )
-      CASE( 2 )                    
-         IF( nn_monocat /= 2 ) THEN          !--- full dynamics (rheology + advection + ridging/rafting + correction)
-            nice_dyn = np_dynFULL
-         ELSE
-            nice_dyn = np_dyn                !--- dynamics without ridging/rafting
-         ENDIF
-      CASE( 1 )                              !--- dynamics without ridging/rafting and correction 
-         nice_dyn = np_dynPURE
-      CASE( 0 )                              !--- prescribed ice velocities (from namelist) 
-         nice_dyn = np_dynNO
-      END SELECT
+      ioptio = 0 
+      !      !--- full dynamics                               (rheology + advection + ridging/rafting + correction)
+      IF( ln_dynFULL .AND. nn_monocat /= 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynFULL      ;   ENDIF
+      !      !--- dynamics without ridging/rafting            (rheology + advection                   + correction)
+      IF( ln_dynFULL .AND. nn_monocat == 2 ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV1   ;   ENDIF
+      !      !--- dynamics without ridging/rafting and corr   (rheology + advection)
+      IF( ln_dynRHGADV                     ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynRHGADV2   ;   ENDIF
+      !      !--- advection only with prescribed ice velocities (from namelist)
+      IF( ln_dynADV                        ) THEN   ;   ioptio = ioptio + 1   ;   nice_dyn = np_dynADV       ;   ENDIF
+      !
+      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_dyn_init: one and only one ice dynamics option has to be defined ' )
       !
       !                                      !--- Lateral boundary conditions
@@ -177,11 +263,7 @@
       ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ===>>>   ice lateral  strong-slip'
       ENDIF
-      !
       !                                      !--- NO Landfast ice : set to zero once for all
       IF( .NOT. ln_landfast )   tau_icebfr(:,:) = 0._wp 
       !
-      !                                      !--- simple conservative piling, comparable with LIM2
-      l_piling = nn_icedyn == 1 .OR. ( nn_monocat == 2  .AND.  jpl == 1 )
-      !
    END SUBROUTINE ice_dyn_init
 
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8517)
@@ -67,5 +67,4 @@
       REAL(wp) ::   zx3        
       REAL(wp) ::   zslope          ! used to compute local thermodynamic "speeds"
-      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
 
       INTEGER , DIMENSION(jpij)       ::   idxice2         ! compute remapping or not
@@ -81,5 +80,5 @@
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' 
 
-      IF( ln_icediachk ) CALL ice_cons_hsm(0, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
 
       !-----------------------------------------------------------------------------------------------
@@ -311,5 +310,5 @@
       ENDIF
       !
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
       !
    END SUBROUTINE ice_itd_rem
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8517)
@@ -127,8 +127,8 @@
       !
       INTEGER, PARAMETER ::   nitermax = 20    
-      !
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!-----------------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('icerdgrft')
+      ! controls
+      IF( nn_timing == 1 )   CALL timing_start('icerdgrft')                                                             ! timing
+      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
 
       IF( kt == nit000 ) THEN
@@ -140,6 +140,4 @@
          !
       ENDIF      
-      !                    ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
       !-----------------------------------------------------------------------------!
@@ -284,14 +282,8 @@
       CALL ice_var_agg( 1 ) 
 
-      !-----------------------------------------------------------------------------!
-      ! control prints
-      !-----------------------------------------------------------------------------!
-      !                    ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      !                    ! control prints
-      IF( ln_ctl       )   CALL ice_prt3D( 'icerdgrft' )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('icerdgrft')
+      ! controls
+      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_ctl         )   CALL ice_prt3D   ('icerdgrft')                                                             ! prints
+      IF( nn_timing == 1 )   CALL timing_stop ('icerdgrft')                                                             ! timing
       !
    END SUBROUTINE ice_rdgrft
@@ -795,5 +787,5 @@
       ELSEIF( ln_str_H79 ) THEN      ! Ice strength => Hibler (1979) method             !
       !                              !--------------------------------------------------!
-         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1)
+         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) )
          !
          ismooth = 1
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8517)
@@ -32,4 +32,9 @@
    PUBLIC   ice_rhg        ! called by icestp.F90
    PUBLIC   ice_rhg_init   ! called by icestp.F90
+
+   INTEGER ::              nice_rhg   ! choice of the type of rheology
+   !                                        ! associated indices:
+   INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
+!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
    
    !! * Substitutions
@@ -53,9 +58,8 @@
       !!--------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ice time step
-      !!
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
       !!--------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )  CALL timing_start('icerhg')
+      ! controls
+      IF( nn_timing == 1 )   CALL timing_start('icerhg')                                                             ! timing
+      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icerhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -64,16 +68,20 @@
          WRITE(numout,*)'~~~~~~~'
       ENDIF
-      !                             ! -- conservation test
-      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      ! -----------------------
-      ! Rheology (ice dynamics)
-      ! -----------------------   
-      CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+      ! --------
+      ! Rheology
+      ! --------   
+      SELECT CASE( nice_rhg )
+      !                                !------------------------!
+      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic !
+         !                             !------------------------!
+         CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+
+      END SELECT
       !
-      !                             !- conservation test
-      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      IF( ln_ctl         )   CALL ice_prt3D  ('icerhg')   !- Control prints
-      IF( nn_timing == 1 )   CALL timing_stop('icerhg')   !- timing
+      ! controls
+      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icerhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_ctl         )   CALL ice_prt3D   ('icerhg')                                                             ! prints
+      IF( nn_timing == 1 )   CALL timing_stop ('icerhg')                                                             ! timing
       !
    END SUBROUTINE ice_rhg
@@ -92,5 +100,5 @@
       !! ** input   :   Namelist namice_rhg
       !!-------------------------------------------------------------------
-      INTEGER ::   ios   ! Local integer output status for namelist read
+      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
       NAMELIST/namice_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
@@ -118,4 +126,13 @@
       ENDIF
       !
+      !                             !== set the choice of ice advection ==!
+      ioptio = 0 
+      IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF
+!!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
+      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_rhg_init: choose one and only one ice rheology' )
+
+     !                              ! allocate tke arrays
+!!clem example      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
+      !
    END SUBROUTINE ice_rhg_init
 
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8517)
@@ -23,6 +23,6 @@
    USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
    USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
-   USE ice            ! ice variables
-   USE icerdgrft      ! ice strength
+   USE ice            ! sea-ice: ice variables
+   USE icerdgrft      ! sea-ice: ice strength
    !
    USE lbclnk         ! Lateral Boundary Condition / MPP link
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8517)
@@ -258,16 +258,16 @@
       CALL ice_itd_init                ! ice thickness distribution initialization
       !
-      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
-      !
-      CALL ice_thd_init                ! set ice thermodynics parameters
-      !
-      CALL ice_thd_sal_init            ! set ice salinity 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
+   
       ! MV MP 2016
       CALL lim_mp_init                 ! set melt ponds parameters
@@ -283,13 +283,14 @@
       CALL ice_var_glo2eqv
       !
-      CALL ice_update_init                 ! ice surface boundary condition
-      !
-      CALL ice_alb_init                    ! ice surface albedo
-      !
-      CALL ice_dia_init                    ! initialization for diags
-      !
-      fr_i  (:,:)   = at_i(:,:)         ! initialisation of sea-ice fraction
-      tn_ice(:,:,:) = t_su(:,:,:)       ! initialisation of surface temp for coupled simu
-      !
+      CALL ice_update_init             ! ice surface boundary condition
+      !
+      CALL ice_alb_init                ! ice surface albedo
+      !
+      CALL ice_dia_init                ! initialization for diags
+      !
+      fr_i  (:,:)   = at_i(:,:)        ! initialisation of sea-ice fraction
+      tn_ice(:,:,:) = t_su(:,:,:)      ! initialisation of surface temp for coupled simu
+      !
+      !                                ! set max concentration in both hemispheres
       WHERE( gphit(:,:) > 0._wp )   ;   rn_amax_2d(:,:) = rn_amax_n  ! NH
       ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH
@@ -306,11 +307,11 @@
       !!
       !! ** Method  :   Read the namice_run namelist and check the parameter
-      !!              values called at the first timestep (nit000)
+      !!                values called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist namice_run
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namice_run/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
-         &                cn_icerst_out, cn_icerst_outdir
+      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
       !!-------------------------------------------------------------------
       !
@@ -333,9 +334,11 @@
          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,*) '      maximum ice concentration for NH                              = ', rn_amax_n 
          WRITE(numout,*) '      maximum ice concentration for SH                              = ', rn_amax_s
       ENDIF
       !
-      !                                         !--- check consistency
+      !                                        !--- check consistency
       IF ( jpl > 1 .AND. nn_monocat == 1 ) THEN
          nn_monocat = 0
@@ -349,9 +352,9 @@
       IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('ice_run_init: online conservation check does not work with BDY')
       !
-      rdt_ice   = REAL(nn_fsbc) * rdt           !--- sea-ice timestep and inverse
+      rdt_ice   = REAL(nn_fsbc) * rdt          !--- sea-ice timestep and inverse
       r1_rdtice = 1._wp / rdt_ice
-      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice  = ', rdt_ice
-      !
-      r1_nlay_i = 1._wp / REAL( nlay_i, wp )    !--- inverse of nlay_i and nlay_s
+      IF( lwp ) WRITE(numout,*) '   ice timestep rdt_ice = ', rdt_ice
+      !
+      r1_nlay_i = 1._wp / REAL( nlay_i, wp )   !--- inverse of nlay_i and nlay_s
       r1_nlay_s = 1._wp / REAL( nlay_s, wp )
       !
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 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8517)
@@ -81,5 +81,4 @@
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
       REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
       REAL(wp), PARAMETER :: zfric_umin = 0._wp           ! lower bound for the friction velocity (cice value=5.e-04)
       REAL(wp), PARAMETER :: zch        = 0.0057_wp       ! heat transfer coefficient
@@ -87,6 +86,7 @@
       !
       !!-------------------------------------------------------------------
-
-      IF( nn_timing == 1 )   CALL timing_start('icethd')
+      ! controls
+      IF( nn_timing == 1 )   CALL timing_start('icethd')                                                             ! timing
+      IF( ln_icediachk   )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
 
       IF( kt == nit000 .AND. lwp ) THEN
@@ -96,7 +96,4 @@
       ENDIF
       
-      ! conservation test
-      IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
-
       CALL ice_var_glo2eqv
 
@@ -257,5 +254,5 @@
       oa_i(:,:,:) = o_i(:,:,:) * a_i(:,:,:)
 
-      IF( ln_icediachk )   CALL ice_cons_hsm( 1, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b )
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
       !
                            CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
@@ -266,8 +263,8 @@
       IF( ln_icedO )       CALL ice_thd_lac                ! --- frazil ice growing in leads --- !
       !
-      IF( ln_icectl )      CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ' )   ! control print
-      IF( ln_ctl )         CALL ice_prt3D( 'icethd' )      ! Control print
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('icethd')
+      ! controls
+      IF( ln_icectl      )   CALL ice_prt    (kt, iiceprt, jiceprt, 1, ' - ice thermodyn. - ') ! prints
+      IF( ln_ctl         )   CALL ice_prt3D  ('icethd')                                        ! prints
+      IF( nn_timing == 1 )   CALL timing_stop('icethd')                                        ! timing
       !
    END SUBROUTINE ice_thd 
@@ -539,5 +536,5 @@
       INTEGER  ::   ios   ! Local integer output status for namelist read
       !!
-      NAMELIST/namice_thd/ ln_icethd, rn_kappa_i, ln_cndi_U64, ln_cndi_P07, ln_dqns_i, rn_cnd_s,   &
+      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,                                             &
@@ -560,5 +557,4 @@
          WRITE(numout,*) '~~~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namice_thd'
-         WRITE(numout,*) '      Ice thermodynamics (T) or not (F)                       ln_icethd    = ', ln_icethd
          WRITE(numout,*) '   -- icethd_dif --'
          WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8517)
@@ -83,5 +83,4 @@
      
       REAL(wp) ::   zv_newfra
-      REAL(wp) ::   zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b ! conservation check
   
       INTEGER , DIMENSION(jpij) ::   jcat        ! indexes of categories where new ice grows
@@ -113,5 +112,5 @@
       !!-----------------------------------------------------------------------!
 
-      IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
 
       CALL ice_var_agg(1)
@@ -479,5 +478,5 @@
       ENDIF ! nidx > 0
       !
-      IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_lac', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
       !
    END SUBROUTINE ice_thd_lac
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8516)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8517)
@@ -104,12 +104,5 @@
       ! END MP 2016
 
-      DO jj = 1, jpj                         ! open water fraction
-         DO ji = 1, jpi
-            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   
-         END DO
-      END DO
-!!gm I think this should do the work :
-!      ato_i(:,:) = MAX( 1._wp - at_i(:,:), 0._wp )  
-!!gm end
+      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction  
 
       IF( kn > 1 ) THEN
@@ -543,18 +536,9 @@
 
       ! to be sure that at_i is the sum of a_i(jl)
-      at_i (:,:) = a_i(:,:,1)
-      vt_i (:,:) = v_i(:,:,1)
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
-      END DO
-
-      ! open water = 1 if at_i=0 (no re-calculation of ato_i here)
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
-            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
-         END DO
-      END DO
+      at_i (:,:) = SUM( a_i(:,:,:), dim=3 )
+      vt_i (:,:) = SUM( v_i(:,:,:), dim=3 )
+
+      ! open water = 1 if at_i=0
+      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
       !
    END SUBROUTINE ice_var_zapsmall
