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 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8518)
@@ -170,7 +170,4 @@
    
    !                                     !!** ice-dynamics namelist (namice_dyn) **
-   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
@@ -180,22 +177,9 @@
    REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction (landfast ice) 
    !
-   !                                     !!** ice-rdige/raft namelist (namice_rdgrft) **
-   LOGICAL , PUBLIC ::   ln_str_H79       !: ice strength parameterization (Hibler79)
-   REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength, Hibler JPO79
-   REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength
-   LOGICAL , PUBLIC ::   ln_str_R75       !: ice strength parameterization (Rothrock75)
-   REAL(wp), PUBLIC ::   rn_perdg         !: ridging work divided by pot. energy change in ridging
-   !
    !                                     !!** ice-rheology namelist (namice_rhg) **
-   LOGICAL , PUBLIC ::   ln_rhg_EVP       !: EVP rheology
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
    REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
-   !
-   !                                     !!** ice-advection namelist (namice_adv) **
-   LOGICAL , PUBLIC ::   ln_adv_Pra       !: Prather        advection scheme
-   LOGICAL , PUBLIC ::   ln_adv_UMx       !: Ultimate-Macho advection scheme
-   INTEGER , PUBLIC ::   nn_UMx           !: order of the UMx advection scheme   
    !
    !                                     !!** ice-thermodynamics namelist (namice_thd) **
@@ -230,5 +214,4 @@
 
    !                                     !!** ice-salinity namelist (namice_sal) **
-   LOGICAL , PUBLIC ::   ln_icedS         !: activate gravity drainage and flushing (T) or not (F)
    INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model
    !                                      ! 1 - constant salinity in both space and time
@@ -236,8 +219,4 @@
    !                                      ! 3 - salinity profile, constant in time
    REAL(wp), PUBLIC ::   rn_icesal        !: bulk salinity (ppt) in case of constant salinity
-   REAL(wp), PUBLIC ::   rn_sal_gd        !: restoring salinity for gravity drainage [PSU]
-   REAL(wp), PUBLIC ::   rn_time_gd       !: restoring time constant for gravity drainage (= 20 days) [s]
-   REAL(wp), PUBLIC ::   rn_sal_fl        !: restoring salinity for flushing [PSU]
-   REAL(wp), PUBLIC ::   rn_time_fl       !: restoring time constant for gravity drainage (= 10 days) [s]
    REAL(wp), PUBLIC ::   rn_simax         !: maximum ice salinity [PSU]
    REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU]
@@ -329,5 +308,4 @@
    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     !: heat flux error after heat diffusion             [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]
@@ -494,5 +472,5 @@
          &      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_err(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) ,     &
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8518)
@@ -47,5 +47,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d
@@ -178,5 +177,5 @@
          &      rn_amax_1d(jpij) ,                                         &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
-         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,   &
+         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) ,                      &
          &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , hfx_out_1d(jpij), STAT=ierr(ii) )
       !
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 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8518)
@@ -41,5 +41,10 @@
    INTEGER, PARAMETER ::   np_advPRA = 1   ! Prather scheme
    INTEGER, PARAMETER ::   np_advUMx = 2   ! Ultimate-Macho scheme
-
+   !
+   ! ** namelist (namice_adv) **
+   LOGICAL ::   ln_adv_Pra   ! Prather        advection scheme
+   LOGICAL ::   ln_adv_UMx   ! Ultimate-Macho advection scheme
+   INTEGER ::   nn_UMx       ! order of the UMx advection scheme   
+   !
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
@@ -85,5 +90,5 @@
       CASE( np_advUMx )                ! ULTIMATE-MACHO scheme !
          !                             !-----------------------!
-         CALL ice_adv_umx( kt, u_ice, v_ice,  &
+         CALL ice_adv_umx( nn_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 )
       !                                !-----------------------!
@@ -107,4 +112,8 @@
       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)
+
+      IF( lrst_ice ) THEN                       !* write Prather fields in the restart file
+         IF( ln_adv_Pra )   CALL adv_pra_rst( 'WRITE', kt )
+      ENDIF
                
       ! controls
@@ -158,4 +167,6 @@
       IF( ioptio /= 1 )   CALL ctl_stop( 'ice_adv_init: choose one and only one ice advection scheme (ln_adv_Pra or ln_adv_UMx)' )
       !
+      IF( ln_adv_Pra )   CALL adv_pra_rst( 'READ' )  !* read or initialize all required files
+      !
    END SUBROUTINE ice_adv_init
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_prather.F90	(revision 8518)
@@ -20,4 +20,5 @@
    USE lbclnk         ! lateral boundary condition - MPP exchanges
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O library
    USE prtctl         ! Print control
    USE lib_mpp        ! MPP library
@@ -28,4 +29,5 @@
 
    PUBLIC   ice_adv_prather   ! called by iceadv
+   PUBLIC   adv_pra_rst       ! called by iceadv
 
    !! * Substitutions
@@ -592,4 +594,405 @@
    END SUBROUTINE adv_y
 
+   SUBROUTINE adv_pra_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE adv_pra_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER ::   jk, jl   ! dummy loop indices
+      INTEGER ::   iter     ! local integer
+      INTEGER ::   id1      ! local integer
+      CHARACTER(len=25) ::   znam
+      CHARACTER(len=2)  ::   zchar, zchar1
+      REAL(wp), DIMENSION(jpi,jpj) :: z2d
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )
+            !
+            IF( id1 > 0 ) THEN      ! fields exist
+               DO jl = 1, jpl 
+                  WRITE(zchar,'(I2.2)') jl
+                  znam = 'sxice'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxice(:,:,jl) = z2d(:,:)
+                  znam = 'syice'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syice(:,:,jl) = z2d(:,:)
+                  znam = 'sxxice'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxxice(:,:,jl) = z2d(:,:)
+                  znam = 'syyice'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syyice(:,:,jl) = z2d(:,:)
+                  znam = 'sxyice'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxyice(:,:,jl) = z2d(:,:)
+                  znam = 'sxsn'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxsn(:,:,jl) = z2d(:,:)
+                  znam = 'sysn'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sysn(:,:,jl) = z2d(:,:)
+                  znam = 'sxxsn'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxxsn(:,:,jl) = z2d(:,:)
+                  znam = 'syysn'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syysn(:,:,jl) = z2d(:,:)
+                  znam = 'sxysn'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxysn(:,:,jl) = z2d(:,:)
+                  znam = 'sxa'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxa(:,:,jl) = z2d(:,:)
+                  znam = 'sya'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sya(:,:,jl) = z2d(:,:)
+                  znam = 'sxxa'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxxa(:,:,jl) = z2d(:,:)
+                  znam = 'syya'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syya(:,:,jl) = z2d(:,:)
+                  znam = 'sxya'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxya(:,:,jl) = z2d(:,:)
+                  znam = 'sxc0'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxc0(:,:,jl) = z2d(:,:)
+                  znam = 'syc0'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syc0(:,:,jl) = z2d(:,:)
+                  znam = 'sxxc0'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxxc0(:,:,jl) = z2d(:,:)
+                  znam = 'syyc0'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syyc0(:,:,jl) = z2d(:,:)
+                  znam = 'sxyc0'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxyc0(:,:,jl) = z2d(:,:)
+                  znam = 'sxsal'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxsal(:,:,jl) = z2d(:,:)
+                  znam = 'sysal'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sysal(:,:,jl) = z2d(:,:)
+                  znam = 'sxxsal'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxxsal(:,:,jl) = z2d(:,:)
+                  znam = 'syysal'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syysal(:,:,jl) = z2d(:,:)
+                  znam = 'sxysal'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxysal(:,:,jl) = z2d(:,:)
+                  znam = 'sxage'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxage(:,:,jl) = z2d(:,:)
+                  znam = 'syage'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syage(:,:,jl) = z2d(:,:)
+                  znam = 'sxxage'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxxage(:,:,jl) = z2d(:,:)
+                  znam = 'syyage'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  syyage(:,:,jl) = z2d(:,:)
+                  znam = 'sxyage'//'_htc'//zchar
+                  CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                  sxyage(:,:,jl)= z2d(:,:)
+               END DO
+               ! MV MP 2016
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  DO jl = 1, jpl 
+                     WRITE(zchar,'(I2.2)') jl
+                     znam = 'sxap'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxap(:,:,jl) = z2d(:,:)
+                     znam = 'syap'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     syap(:,:,jl) = z2d(:,:)
+                     znam = 'sxxap'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxxap(:,:,jl) = z2d(:,:)
+                     znam = 'syyap'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     syyap(:,:,jl) = z2d(:,:)
+                     znam = 'sxyap'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxyap(:,:,jl) = z2d(:,:)
+
+                     znam = 'sxvp'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxvp(:,:,jl) = z2d(:,:)
+                     znam = 'syvp'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     syvp(:,:,jl) = z2d(:,:)
+                     znam = 'sxxvp'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxxvp(:,:,jl) = z2d(:,:)
+                     znam = 'syyvp'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     syyvp(:,:,jl) = z2d(:,:)
+                     znam = 'sxyvp'//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxyvp(:,:,jl) = z2d(:,:)
+                  END DO
+               ENDIF
+               ! END MV MP 2016
+
+               CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  )
+               CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw )
+               CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw )
+
+               DO jl = 1, jpl 
+                  WRITE(zchar,'(I2.2)') jl
+                  DO jk = 1, nlay_i 
+                     WRITE(zchar1,'(I2.2)') jk
+                     znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxe(:,:,jk,jl) = z2d(:,:)
+                     znam = 'sye'//'_il'//zchar1//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sye(:,:,jk,jl) = z2d(:,:)
+                     znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxxe(:,:,jk,jl) = z2d(:,:)
+                     znam = 'syye'//'_il'//zchar1//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     syye(:,:,jk,jl) = z2d(:,:)
+                     znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
+                     CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+                     sxye(:,:,jk,jl) = z2d(:,:)
+                  END DO
+               END DO
+               !
+            ELSE                                     ! start rheology from rest
+               IF(lwp) WRITE(numout,*) '   ==>>   previous run without Prather, set moments to 0'
+               sxopw (:,:) = 0._wp   ;   sxice (:,:,:) = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:) = 0._wp
+               syopw (:,:) = 0._wp   ;   syice (:,:,:) = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:) = 0._wp
+               sxxopw(:,:) = 0._wp   ;   sxxice(:,:,:) = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:) = 0._wp
+               syyopw(:,:) = 0._wp   ;   syyice(:,:,:) = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:) = 0._wp
+               sxyopw(:,:) = 0._wp   ;   sxyice(:,:,:) = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:) = 0._wp
+               !
+               sxc0  (:,:,:) = 0._wp   ;   sxe  (:,:,:,:) = 0._wp   ;   sxsal  (:,:,:) = 0._wp   ;   sxage  (:,:,:) = 0._wp
+               syc0  (:,:,:) = 0._wp   ;   sye  (:,:,:,:) = 0._wp   ;   sysal  (:,:,:) = 0._wp   ;   syage  (:,:,:) = 0._wp
+               sxxc0 (:,:,:) = 0._wp   ;   sxxe (:,:,:,:) = 0._wp   ;   sxxsal (:,:,:) = 0._wp   ;   sxxage (:,:,:) = 0._wp
+               syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
+               sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
+               ! MV MP 2016
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
+                  syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
+                  sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
+                  syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
+                  sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
+               ENDIF
+               ! END MV MP 2016
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set moments to 0'
+            sxopw (:,:) = 0._wp   ;   sxice (:,:,:) = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:) = 0._wp
+            syopw (:,:) = 0._wp   ;   syice (:,:,:) = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:) = 0._wp
+            sxxopw(:,:) = 0._wp   ;   sxxice(:,:,:) = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:) = 0._wp
+            syyopw(:,:) = 0._wp   ;   syyice(:,:,:) = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:) = 0._wp
+            sxyopw(:,:) = 0._wp   ;   sxyice(:,:,:) = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:) = 0._wp
+            !
+            sxc0  (:,:,:) = 0._wp   ;   sxe  (:,:,:,:) = 0._wp   ;   sxsal  (:,:,:) = 0._wp   ;   sxage  (:,:,:) = 0._wp
+            syc0  (:,:,:) = 0._wp   ;   sye  (:,:,:,:) = 0._wp   ;   sysal  (:,:,:) = 0._wp   ;   syage  (:,:,:) = 0._wp
+            sxxc0 (:,:,:) = 0._wp   ;   sxxe (:,:,:,:) = 0._wp   ;   sxxsal (:,:,:) = 0._wp   ;   sxxage (:,:,:) = 0._wp
+            syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
+            sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
+            ! MV MP 2016
+            IF ( nn_pnd_scheme > 0 ) THEN
+               sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
+               syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
+               sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
+               syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
+               sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
+            ENDIF
+            ! END MV MP 2016
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- adv-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         DO jl = 1, jpl 
+            WRITE(zchar,'(I2.2)') jl
+            znam = 'sxice'//'_htc'//zchar
+            z2d(:,:) = sxice(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syice'//'_htc'//zchar
+            z2d(:,:) = syice(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxice'//'_htc'//zchar
+            z2d(:,:) = sxxice(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syyice'//'_htc'//zchar
+            z2d(:,:) = syyice(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxyice'//'_htc'//zchar
+            z2d(:,:) = sxyice(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxsn'//'_htc'//zchar
+            z2d(:,:) = sxsn(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sysn'//'_htc'//zchar
+            z2d(:,:) = sysn(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxsn'//'_htc'//zchar
+            z2d(:,:) = sxxsn(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syysn'//'_htc'//zchar
+            z2d(:,:) = syysn(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxysn'//'_htc'//zchar
+            z2d(:,:) = sxysn(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxa'//'_htc'//zchar
+            z2d(:,:) = sxa(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sya'//'_htc'//zchar
+            z2d(:,:) = sya(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxa'//'_htc'//zchar
+            z2d(:,:) = sxxa(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syya'//'_htc'//zchar
+            z2d(:,:) = syya(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxya'//'_htc'//zchar
+            z2d(:,:) = sxya(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxc0'//'_htc'//zchar
+            z2d(:,:) = sxc0(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syc0'//'_htc'//zchar
+            z2d(:,:) = syc0(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxc0'//'_htc'//zchar
+            z2d(:,:) = sxxc0(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syyc0'//'_htc'//zchar
+            z2d(:,:) = syyc0(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxyc0'//'_htc'//zchar
+            z2d(:,:) = sxyc0(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxsal'//'_htc'//zchar
+            z2d(:,:) = sxsal(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sysal'//'_htc'//zchar
+            z2d(:,:) = sysal(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxsal'//'_htc'//zchar
+            z2d(:,:) = sxxsal(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syysal'//'_htc'//zchar
+            z2d(:,:) = syysal(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxysal'//'_htc'//zchar
+            z2d(:,:) = sxysal(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxage'//'_htc'//zchar
+            z2d(:,:) = sxage(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syage'//'_htc'//zchar
+            z2d(:,:) = syage(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxxage'//'_htc'//zchar
+            z2d(:,:) = sxxage(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'syyage'//'_htc'//zchar
+            z2d(:,:) = syyage(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            znam = 'sxyage'//'_htc'//zchar
+            z2d(:,:) = sxyage(:,:,jl)
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         END DO
+
+         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  )
+         CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw )
+         CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw )
+         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw )
+         
+         DO jl = 1, jpl 
+            WRITE(zchar,'(I2.2)') jl
+            DO jk = 1, nlay_i 
+               WRITE(zchar1,'(I2.2)') jk
+               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
+               z2d(:,:) = sxe(:,:,jk,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sye'//'_il'//zchar1//'_htc'//zchar
+               z2d(:,:) = sye(:,:,jk,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
+               z2d(:,:) = sxxe(:,:,jk,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'syye'//'_il'//zchar1//'_htc'//zchar
+               z2d(:,:) = syye(:,:,jk,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
+               z2d(:,:) = sxye(:,:,jk,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            END DO
+         END DO
+         ! MV MP 2016
+         IF ( nn_pnd_scheme > 0 ) THEN
+            DO jl = 1, jpl 
+               WRITE(zchar,'(I2.2)') jl
+               znam = 'sxap'//'_htc'//zchar
+               z2d(:,:) = sxap(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'syap'//'_htc'//zchar
+               z2d(:,:) = syap(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sxxap'//'_htc'//zchar
+               z2d(:,:) = sxxap(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'syyap'//'_htc'//zchar
+               z2d(:,:) = syyap(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sxyap'//'_htc'//zchar
+               z2d(:,:) = sxyap(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+   
+               znam = 'sxvp'//'_htc'//zchar
+               z2d(:,:) = sxvp(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'syvp'//'_htc'//zchar
+               z2d(:,:) = syvp(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sxxvp'//'_htc'//zchar
+               z2d(:,:) = sxxvp(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'syyvp'//'_htc'//zchar
+               z2d(:,:) = syyvp(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+               znam = 'sxyvp'//'_htc'//zchar
+               z2d(:,:) = sxyvp(:,:,jl)
+               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            END DO
+         ENDIF
+         !
+      ENDIF
+      !
+   END SUBROUTINE adv_pra_rst
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv_umx.F90	(revision 8518)
@@ -29,5 +29,5 @@
    PRIVATE
 
-   PUBLIC   ice_adv_umx    ! routine called by iceadv.F90
+   PUBLIC   ice_adv_umx    ! called by iceadv.F90
       
    REAL(wp) ::   z1_6   = 1._wp /   6._wp   ! =1/6
@@ -43,5 +43,5 @@
 CONTAINS
 
-   SUBROUTINE ice_adv_umx( kt, pu_ice, pv_ice,  &
+   SUBROUTINE ice_adv_umx( k_order, kt, pu_ice, pv_ice,  &
       &                    pato_i, pv_i, pv_s, psmv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
@@ -54,4 +54,5 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the scheme (1-5 or 20)
       INTEGER                     , INTENT(in   ) ::   kt         ! time step
       REAL(wp), DIMENSION(:,:)    , INTENT(in   ) ::   pu_ice     ! ice i-velocity
@@ -112,18 +113,18 @@
       !---------------!
       DO jt = 1, initad
-         CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area 
+         CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )             ! Open water area 
          DO jl = 1, jpl
-            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
-            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
-            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content
-            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content
+            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) )         ! Ice area
+            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) )         ! Ice  volume
+            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psmv_i(:,:,jl) )       ! Salt content
+            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i (:,:,jl) )       ! Age content
             DO jk = 1, nlay_i
-               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
-            END DO
-            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
-            CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
+               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) )   ! Ice  heat content
+            END DO
+            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
+            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
             IF ( nn_pnd_scheme > 0 ) THEN
-               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
-               CALL adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
+               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
+               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
             ENDIF
          END DO
@@ -134,5 +135,5 @@
    END SUBROUTINE ice_adv_umx
    
-   SUBROUTINE adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc )
+   SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE adv_umx  ***
@@ -147,4 +148,5 @@
       !! ** Action : - pt  the after advective tracer
       !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   k_order        ! order of the ULTIMATE scheme
       INTEGER                     , INTENT(in   ) ::   kt             ! number of iteration
       REAL(wp)                    , INTENT(in   ) ::   pdt            ! tracer time-step
@@ -189,5 +191,5 @@
       ! High order (_ho) fluxes 
       ! -----------------------
-      SELECT CASE( nn_UMx )
+      SELECT CASE( k_order )
       CASE ( 20 )                          ! centered second order
          DO jj = 2, jpjm1
@@ -199,5 +201,5 @@
          !
       CASE ( 1:5 )                      ! 1st to 5th order ULTIMATE-MACHO scheme
-         CALL macho( kt, nn_UMx, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
+         CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v )
          !
          DO jj = 2, jpjm1
@@ -240,5 +242,5 @@
 
 
-   SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
+   SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE ultimate_x  ***
@@ -251,6 +253,6 @@
       !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 
       !!----------------------------------------------------------------------
+      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
       INTEGER                     , INTENT(in   ) ::   kt         ! number of iteration
-      INTEGER                     , INTENT(in   ) ::   k_order    ! order of the ULTIMATE scheme
       REAL(wp)                    , INTENT(in   ) ::   pdt        ! tracer time-step
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   ptc        ! tracer fields
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8518)
@@ -45,4 +45,15 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
+   INTEGER FUNCTION ice_dia_alloc()
+      !!---------------------------------------------------------------------!
+      !!                ***  ROUTINE ice_rdgrft_alloc ***
+      !!---------------------------------------------------------------------!
+      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ice_dia_alloc )
+
+      IF( lk_mpp             )   CALL mpp_sum ( ice_dia_alloc )
+      IF( ice_dia_alloc /= 0 )   CALL ctl_warn( 'ice_dia_alloc: failed to allocate arrays' )
+      !
+   END FUNCTION ice_dia_alloc
 
    SUBROUTINE ice_dia( kt )
@@ -188,11 +199,6 @@
       !      
       IF( ln_icediahsb ) THEN
-         ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror )
-         IF( ierror > 0 )  THEN
-            CALL ctl_stop( 'ice_dia: unable to allocate vol_loc_ini' )
-            RETURN
-         ENDIF
-         !
-         CALL ice_dia_rst( 'READ' )  !* read or initialize all required files
+         IF( ice_dia_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_dia_init : unable to allocate arrays' )   ! allocate tke arrays
+         CALL ice_dia_rst( 'READ' )   ! read or initialize all required files
       ENDIF
       !
@@ -220,6 +226,6 @@
             CALL iom_get( numrir, 'kt_ice' , ziter )
             IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) ' ice_dia_rst read at time step = ', ziter
-            IF(lwp) WRITE(numout,*) '~~~~~~~'
+            IF(lwp) WRITE(numout,*) 'ice_dia_rst read at time step = ', ziter
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
             CALL iom_get( numrir, 'frc_voltop' , frc_voltop  )
             CALL iom_get( numrir, 'frc_volbot' , frc_volbot  )
@@ -252,6 +258,6 @@
          IF( iter == nitrst ) THEN
             IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) ' ice_dia_rst write at time step = ', kt
-            IF(lwp) WRITE(numout,*) '~~~~~~~'
+            IF(lwp) WRITE(numout,*) 'ice_dia_rst write at time step = ', kt
+            IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
          ENDIF
          !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8518)
@@ -42,6 +42,9 @@
    !
    ! ** namelist (namdyn) **
-   REAL(wp) ::   rn_uice          ! prescribed u-vel (case np_dynADV)
-   REAL(wp) ::   rn_vice          ! prescribed v-vel (case np_dynADV)
+   LOGICAL  ::   ln_dynFULL       ! full ice dynamics               (rheology + advection + ridging/rafting + correction)
+   LOGICAL  ::   ln_dynRHGADV     ! no ridge/raft & no corrections  (rheology + advection)
+   LOGICAL  ::   ln_dynADV        ! only advection w prescribed vel.(rn_uvice + advection)
+   REAL(wp) ::   rn_uice          !    prescribed u-vel (case np_dynADV)
+   REAL(wp) ::   rn_vice          !    prescribed v-vel (case np_dynADV)
    
    !! * Substitutions
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8518)
@@ -17,12 +17,13 @@
    !!   ice_istate_init  :  initialization of ice state and namelist read
    !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
    USE phycst         ! physical constant
    USE oce            ! dynamics and tracers variables
    USE dom_oce        ! ocean domain
-   USE sbc_oce , ONLY : sst_m, sss_m 
-   USE sbc_ice , ONLY : tn_ice
+   USE sbc_oce , ONLY : sst_m, sss_m, ln_ice_embd 
+   USE sbc_ice , ONLY : tn_ice, snwice_mass, snwice_mass_b
    USE eosbn2         ! equation of state
+   USE domvvl         ! Variable volume
    USE ice            ! sea-ice variables
-   USE par_oce        ! ocean parameters
    USE icevar         ! ice_var_salprof
    !
@@ -91,13 +92,14 @@
       !!              where there is no ice (clem: I do not know why, is it mandatory?) 
       !!--------------------------------------------------------------------
-      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
-      REAL(wp)   :: ztmelts, zdh
-      INTEGER    :: i_hemis, i_fill, jl0
-      REAL(wp)   :: zarg, zV, zconv, zdv
-      REAL(wp), DIMENSION(jpi,jpj)     :: zswitch    ! ice indicator
-      REAL(wp), DIMENSION(jpi,jpj)     :: zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj)     :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
-      REAL(wp), DIMENSION(jpi,jpj,jpl) :: zh_i_ini, za_i_ini                         !data by cattegories to fill
-      INTEGER , DIMENSION(4)           :: itest
+      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      REAL(wp) ::   ztmelts, zdh
+      INTEGER  ::   i_hemis, i_fill, jl0
+      REAL(wp) ::   zarg, zV, zconv, zdv
+      INTEGER , DIMENSION(4)           ::   itest
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zswitch    ! ice indicator
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, zvt_i_ini            !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zh_i_ini, za_i_ini                         !data by cattegories to fill
       !--------------------------------------------------------------------
 
@@ -119,5 +121,4 @@
       CALL eos_fzp( sss_m(:,:), t_bo(:,:) )
       t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 
-
 
       !--------------------------------------------------------------------
@@ -377,21 +378,15 @@
          ! Melt pond volume and fraction
          IF ( ln_pnd ) THEN
-
             DO jl = 1, jpl
-
                a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
                h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
                a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
                v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
-
-            END DO
-
+            END DO
          ELSE
-
             a_ip(:,:,:)      = 0._wp
             v_ip(:,:,:)      = 0._wp
             a_ip_frac(:,:,:) = 0._wp
             h_ip     (:,:,:) = 0._wp
-
          ENDIF
          ! END MV MP 2016
@@ -437,42 +432,57 @@
       u_ice (:,:)     = 0._wp
       v_ice (:,:)     = 0._wp
-      stress1_i(:,:)  = 0._wp
-      stress2_i(:,:)  = 0._wp
-      stress12_i(:,:) = 0._wp
-
-      !--------------------------------------------------------------------
-      ! 5) Moments for advection
-      !--------------------------------------------------------------------
-
-      sxopw (:,:) = 0._wp 
-      syopw (:,:) = 0._wp 
-      sxxopw(:,:) = 0._wp 
-      syyopw(:,:) = 0._wp 
-      sxyopw(:,:) = 0._wp
-
-      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
-      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
-      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
-      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
-      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
-
-      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
-      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
-      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
-      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
-      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
-
-      sxsal  (:,:,:)  = 0._wp
-      sysal  (:,:,:)  = 0._wp
-      sxxsal (:,:,:)  = 0._wp
-      syysal (:,:,:)  = 0._wp
-      sxysal (:,:,:)  = 0._wp
-
-      sxage  (:,:,:)  = 0._wp
-      syage  (:,:,:)  = 0._wp
-      sxxage (:,:,:)  = 0._wp
-      syyage (:,:,:)  = 0._wp
-      sxyage (:,:,:)  = 0._wp
-
+      !
+      !--------------------------------------------------------------------
+      ! Snow-ice mass (case ice is fully embedded)                    | 
+      !--------------------------------------------------------------------
+      snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass
+      snwice_mass_b(:,:) = snwice_mass(:,:)
+      !
+      IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
+
+         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+
+         IF( .NOT.ln_linssh ) THEN
+            
+            WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
+            ELSEWHERE                ;   z2d(:,:) = 1._wp   ;   END WHERE
+         
+            DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                
+               e3t_n(:,:,jk) = e3t_0(:,:,jk) * z2d(:,:)
+               e3t_b(:,:,jk) = e3t_n(:,:,jk)
+               e3t_a(:,:,jk) = e3t_n(:,:,jk)
+            END DO
+            
+            ! Reconstruction of all vertical scale factors at now and before time-steps
+            ! =========================================================================
+            ! Horizontal scale factor interpolations
+            ! --------------------------------------
+            CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
+            CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
+            CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
+            CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
+            CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
+            ! Vertical scale factor interpolations
+            ! ------------------------------------
+            CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
+            CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
+            CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
+            CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
+            CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
+            ! t- and w- points depth
+            ! ----------------------
+            !!gm not sure of that....
+            gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
+            gdepw_n(:,:,1) = 0.0_wp
+            gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
+            DO jk = 2, jpk
+               gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
+               gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
+               gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
+            END DO
+         ENDIF
+      ENDIF
+      
       !------------------------------------
       ! 6) store fields at before time-step
@@ -488,14 +498,4 @@
       u_ice_b(:,:)     = u_ice(:,:)
       v_ice_b(:,:)     = v_ice(:,:)
-
-      ! MV MP 2016
-      IF ( nn_pnd_scheme > 0 ) THEN
-         sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
-         syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
-         sxxap (:,:,:) = 0._wp    ; sxxvp (:,:,:) = 0._wp 
-         syyap (:,:,:) = 0._wp    ; syyvp (:,:,:) = 0._wp 
-         sxyap (:,:,:) = 0._wp    ; sxyvp (:,:,:) = 0._wp
-      ENDIF
-      ! END MV MP 2016
 
 !!!clem
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 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8518)
@@ -34,16 +34,15 @@
    PUBLIC   ice_rdgrft_strength      ! called by icerhg_evp
    PUBLIC   ice_rdgrft_init          ! called by ice_stp
-   PUBLIC   ice_rdgrft_alloc         ! called by ice_init 
 
    ! Variables shared among ridging subroutines
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
    !
    REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
@@ -52,4 +51,9 @@
    !
    ! ** namelist (namice_rdgrft) **
+   LOGICAL  ::   ln_str_H79       ! ice strength parameterization (Hibler79)
+   REAL(wp) ::   rn_pstar         ! determines ice strength, Hibler JPO79
+   REAL(wp) ::   rn_crhg          ! determines changes in ice strength
+   LOGICAL  ::   ln_str_R75       ! ice strength parameterization (Rothrock75)
+   REAL(wp) ::   rn_perdg         ! ridging work divided by pot. energy change in ridging
    REAL(wp) ::   rn_csrdg         ! fraction of shearing energy contributing to ridging            
    LOGICAL  ::   ln_partf_lin     ! participation function linear (Thorndike et al. (1975))
@@ -82,5 +86,6 @@
          &      hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl)   , aridge(jpi,jpj,jpl) ,     &
          &      hrmax(jpi,jpj,jpl) , krdg  (jpi,jpj,jpl)   , araft (jpi,jpj,jpl) , STAT=ice_rdgrft_alloc )
-         !
+
+      IF( lk_mpp                )   CALL mpp_sum ( ice_rdgrft_alloc )
       IF( ice_rdgrft_alloc /= 0 )   CALL ctl_warn( 'ice_rdgrft_alloc: failed to allocate arrays' )
       !
@@ -909,4 +914,6 @@
          CALL ctl_stop( 'ice_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
       ENDIF
+      !                              ! allocate tke arrays
+      IF( ice_rdgrft_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_rdgrft_init : unable to allocate arrays' )
       !
   END SUBROUTINE ice_rdgrft_init
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8518)
@@ -37,5 +37,8 @@
    INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
 !! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
-   
+
+   ! ** namelist (namrhg) **
+   LOGICAL ::   ln_rhg_EVP       ! EVP rheology
+   !
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -77,6 +80,10 @@
          !                             !------------------------!
          CALL ice_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
-
+         !         
       END SELECT
+      !
+      IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
+         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
+      ENDIF
       !
       ! controls
@@ -86,5 +93,4 @@
       !
    END SUBROUTINE ice_rhg
-
 
    SUBROUTINE ice_rhg_init
@@ -131,7 +137,6 @@
 !!    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' )
+      !
+      IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
       !
    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 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8518)
@@ -41,5 +41,6 @@
    PRIVATE
 
-   PUBLIC   ice_rhg_evp        ! routine called by icerhg.F90
+   PUBLIC   ice_rhg_evp   ! called by icerhg.F90
+   PUBLIC   rhg_evp_rst   ! called by icerhg.F90
 
    !! * Substitutions
@@ -108,6 +109,6 @@
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt      ! time step
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pu_ice, pv_ice, pshear_i, pdivu_i, pdelta_i 
       !!
       INTEGER ::   ji, jj       ! dummy loop indices
@@ -825,4 +826,57 @@
    END SUBROUTINE ice_rhg_evp
 
+   SUBROUTINE rhg_evp_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rhg_evp_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER  ::   iter            ! local integer
+      INTEGER  ::   id1, id2, id3   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
+            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
+            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
+               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+            ELSE                                     ! start rheology from rest
+               IF(lwp) WRITE(numout,*) '   ==>>   previous run without rheology, set stresses to 0'
+               stress1_i (:,:) = 0._wp
+               stress2_i (:,:) = 0._wp
+               stress12_i(:,:) = 0._wp
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set stresses to 0'
+            stress1_i (:,:) = 0._wp
+            stress2_i (:,:) = 0._wp
+            stress12_i(:,:) = 0._wp
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+         !
+      ENDIF
+      !
+   END SUBROUTINE rhg_evp_rst
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8518)
@@ -17,5 +17,4 @@
    !!----------------------------------------------------------------------
    USE ice            ! sea-ice variables
-   USE sbc_ice , ONLY :  snwice_mass, snwice_mass_b
    USE dom_oce        ! ocean domain
    USE sbc_oce , ONLY : nn_fsbc
@@ -30,10 +29,7 @@
    PRIVATE
 
-   PUBLIC   ice_rst_opn    ! routine called by icestep.F90
-   PUBLIC   ice_rst_write  ! routine called by icestep.F90
-   PUBLIC   ice_rst_read   ! routine called by ice_init
-
-   LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write 
-   INTEGER, PUBLIC ::   numrir, numriw   !: logical unit for ice restart (read and write)
+   PUBLIC   ice_rst_opn     ! called by icestp
+   PUBLIC   ice_rst_write   ! called by icestp
+   PUBLIC   ice_rst_read    ! called by ice_init
 
    !!----------------------------------------------------------------------
@@ -79,7 +75,8 @@
                   WRITE(numout,*) '             open ice restart NetCDF file: ',TRIM(clpath)//clname
                END SELECT
-               IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
-                  WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
-               ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
+               IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN
+                  WRITE(numout,*) '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
+               ELSE
+                  WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
                ENDIF
             ENDIF
@@ -102,5 +99,5 @@
       INTEGER, INTENT(in) ::   kt     ! number of iteration
       !!
-      INTEGER ::   ji, jj, jk ,jl   ! dummy loop indices
+      INTEGER ::   jk ,jl   ! dummy loop indices
       INTEGER ::   iter
       CHARACTER(len=25) ::   znam
@@ -127,5 +124,4 @@
 !!gm   "just" ask Sebatien
 
-
       ! Prognostic variables 
       DO jl = 1, jpl 
@@ -133,20 +129,20 @@
          znam = 'v_i'//'_htc'//zchar
          z2d(:,:) = v_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_i
          znam = 'v_s'//'_htc'//zchar
          z2d(:,:) = v_s(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_s
          znam = 'smv_i'//'_htc'//zchar
          z2d(:,:) = smv_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! smv_i
          znam = 'oa_i'//'_htc'//zchar
          z2d(:,:) = oa_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! oa_i
          znam = 'a_i'//'_htc'//zchar
          z2d(:,:) = a_i(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_i
          znam = 't_su'//'_htc'//zchar
          z2d(:,:) = t_su(:,:,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! t_su
       END DO
 
@@ -157,8 +153,8 @@
             znam = 'a_ip'//'_htc'//zchar
             z2d(:,:) = a_ip(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip
             znam = 'v_ip'//'_htc'//zchar
             z2d(:,:) = v_ip(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip
          END DO
       ENDIF
@@ -169,185 +165,15 @@
          znam = 'tempt_sl1'//'_htc'//zchar
          z2d(:,:) = e_s(:,:,1,jl)
-         CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
+         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_s
          DO jk = 1, nlay_i 
             WRITE(zchar1,'(I2.2)') jk
             znam = 'tempt'//'_il'//zchar1//'_htc'//zchar
             z2d(:,:) = e_i(:,:,jk,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-         END DO
-      END DO
-
-      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'        , u_ice      )
-      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'        , v_ice      )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'    , stress1_i  )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'    , stress2_i  )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i'   , stress12_i )
-      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass )
-      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
-
-      ! In case Prather scheme is used for advection, write second order moments
-      ! ------------------------------------------------------------------------
-      IF( ln_adv_Pra ) THEN
-         
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'sxice'//'_htc'//zchar
-            z2d(:,:) = sxice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syice'//'_htc'//zchar
-            z2d(:,:) = syice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxice'//'_htc'//zchar
-            z2d(:,:) = sxxice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syyice'//'_htc'//zchar
-            z2d(:,:) = syyice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxyice'//'_htc'//zchar
-            z2d(:,:) = sxyice(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxsn'//'_htc'//zchar
-            z2d(:,:) = sxsn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sysn'//'_htc'//zchar
-            z2d(:,:) = sysn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxsn'//'_htc'//zchar
-            z2d(:,:) = sxxsn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syysn'//'_htc'//zchar
-            z2d(:,:) = syysn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxysn'//'_htc'//zchar
-            z2d(:,:) = sxysn(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxa'//'_htc'//zchar
-            z2d(:,:) = sxa(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sya'//'_htc'//zchar
-            z2d(:,:) = sya(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxa'//'_htc'//zchar
-            z2d(:,:) = sxxa(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syya'//'_htc'//zchar
-            z2d(:,:) = syya(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxya'//'_htc'//zchar
-            z2d(:,:) = sxya(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxc0'//'_htc'//zchar
-            z2d(:,:) = sxc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syc0'//'_htc'//zchar
-            z2d(:,:) = syc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxc0'//'_htc'//zchar
-            z2d(:,:) = sxxc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syyc0'//'_htc'//zchar
-            z2d(:,:) = syyc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxyc0'//'_htc'//zchar
-            z2d(:,:) = sxyc0(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxsal'//'_htc'//zchar
-            z2d(:,:) = sxsal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sysal'//'_htc'//zchar
-            z2d(:,:) = sysal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxsal'//'_htc'//zchar
-            z2d(:,:) = sxxsal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syysal'//'_htc'//zchar
-            z2d(:,:) = syysal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxysal'//'_htc'//zchar
-            z2d(:,:) = sxysal(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxage'//'_htc'//zchar
-            z2d(:,:) = sxage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syage'//'_htc'//zchar
-            z2d(:,:) = syage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxxage'//'_htc'//zchar
-            z2d(:,:) = sxxage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'syyage'//'_htc'//zchar
-            z2d(:,:) = syyage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            znam = 'sxyage'//'_htc'//zchar
-            z2d(:,:) = sxyage(:,:,jl)
-            CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-         END DO
-
-         CALL iom_rstput( iter, nitrst, numriw, 'sxopw ' ,  sxopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syopw ' ,  syopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw' ,  sxxopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyopw' ,  syyopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw' ,  sxyopw )
-         
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            DO jk = 1, nlay_i 
-               WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sxe(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sye'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sye(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sxxe(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syye'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = syye(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
-               z2d(:,:) = sxye(:,:,jk,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            END DO
-         END DO
-         ! MV MP 2016
-         IF ( nn_pnd_scheme > 0 ) THEN
-            DO jl = 1, jpl 
-               WRITE(zchar,'(I2.2)') jl
-               znam = 'sxap'//'_htc'//zchar
-               z2d(:,:) = sxap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syap'//'_htc'//zchar
-               z2d(:,:) = syap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxxap'//'_htc'//zchar
-               z2d(:,:) = sxxap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syyap'//'_htc'//zchar
-               z2d(:,:) = syyap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxyap'//'_htc'//zchar
-               z2d(:,:) = sxyap(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-   
-               znam = 'sxvp'//'_htc'//zchar
-               z2d(:,:) = sxvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syvp'//'_htc'//zchar
-               z2d(:,:) = syvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxxvp'//'_htc'//zchar
-               z2d(:,:) = sxxvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'syyvp'//'_htc'//zchar
-               z2d(:,:) = syyvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-               znam = 'sxyvp'//'_htc'//zchar
-               z2d(:,:) = sxyvp(:,:,jl)
-               CALL iom_rstput( iter, nitrst, numriw, znam , z2d )
-            END DO
-         ENDIF
-
-      ENDIF
+            CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! e_i
+         END DO
+      END DO
+
+      CALL iom_rstput( iter, nitrst, numriw, 'u_ice', u_ice ) ! u_ice
+      CALL iom_rstput( iter, nitrst, numriw, 'v_ice', v_ice ) ! v_ice
       
       ! close restart file
@@ -368,5 +194,5 @@
       !! ** purpose  :   read of sea-ice variable restart in a netcdf file
       !!----------------------------------------------------------------------
-      INTEGER  :: ji, jj, jk, jl
+      INTEGER  :: jk, jl
       REAL(wp) :: zfice, ziter
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d
@@ -379,6 +205,6 @@
       IF(lwp) THEN
          WRITE(numout,*)
-         WRITE(numout,*) 'ice_rst_read : read ice NetCDF restart file'
-         WRITE(numout,*) '~~~~~~~~~~~~~'
+         WRITE(numout,*) 'ice_rst_read: read ice NetCDF restart file'
+         WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
 
@@ -390,6 +216,5 @@
       IF(lwp) WRITE(numout,*) '   in any case we force it to nit000 - 1 : ', nit000 - 1
 
-      !Control of date
-
+      ! Control of date
       IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   &
          &     CALL ctl_stop( 'ice_rst_read ===>>>> : problem with nit000 in ice restart',  &
@@ -451,191 +276,7 @@
       END DO
 
-      CALL iom_get( numrir, jpdom_autoglo, 'u_ice'     , u_ice      )
-      CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      )
-      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
-      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
-      CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
-      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass )
-      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
-
-      ! In case Prather scheme is used for advection, read second order moments
-      ! ------------------------------------------------------------------------
-      IF( ln_adv_Pra ) THEN
-
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            znam = 'sxice'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxice(:,:,jl) = z2d(:,:)
-            znam = 'syice'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syice(:,:,jl) = z2d(:,:)
-            znam = 'sxxice'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxxice(:,:,jl) = z2d(:,:)
-            znam = 'syyice'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syyice(:,:,jl) = z2d(:,:)
-            znam = 'sxyice'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxyice(:,:,jl) = z2d(:,:)
-            znam = 'sxsn'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxsn(:,:,jl) = z2d(:,:)
-            znam = 'sysn'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sysn(:,:,jl) = z2d(:,:)
-            znam = 'sxxsn'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxxsn(:,:,jl) = z2d(:,:)
-            znam = 'syysn'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syysn(:,:,jl) = z2d(:,:)
-            znam = 'sxysn'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxysn(:,:,jl) = z2d(:,:)
-            znam = 'sxa'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxa(:,:,jl) = z2d(:,:)
-            znam = 'sya'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sya(:,:,jl) = z2d(:,:)
-            znam = 'sxxa'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxxa(:,:,jl) = z2d(:,:)
-            znam = 'syya'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syya(:,:,jl) = z2d(:,:)
-            znam = 'sxya'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxya(:,:,jl) = z2d(:,:)
-            znam = 'sxc0'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxc0(:,:,jl) = z2d(:,:)
-            znam = 'syc0'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syc0(:,:,jl) = z2d(:,:)
-            znam = 'sxxc0'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxxc0(:,:,jl) = z2d(:,:)
-            znam = 'syyc0'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syyc0(:,:,jl) = z2d(:,:)
-            znam = 'sxyc0'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxyc0(:,:,jl) = z2d(:,:)
-            znam = 'sxsal'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxsal(:,:,jl) = z2d(:,:)
-            znam = 'sysal'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sysal(:,:,jl) = z2d(:,:)
-            znam = 'sxxsal'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxxsal(:,:,jl) = z2d(:,:)
-            znam = 'syysal'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syysal(:,:,jl) = z2d(:,:)
-            znam = 'sxysal'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxysal(:,:,jl) = z2d(:,:)
-            znam = 'sxage'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxage(:,:,jl) = z2d(:,:)
-            znam = 'syage'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syage(:,:,jl) = z2d(:,:)
-            znam = 'sxxage'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxxage(:,:,jl) = z2d(:,:)
-            znam = 'syyage'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            syyage(:,:,jl) = z2d(:,:)
-            znam = 'sxyage'//'_htc'//zchar
-            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-            sxyage(:,:,jl)= z2d(:,:)
-         END DO
-         ! MV MP 2016
-         IF ( nn_pnd_scheme > 0 ) THEN
-            DO jl = 1, jpl 
-               WRITE(zchar,'(I2.2)') jl
-               znam = 'sxap'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxap(:,:,jl) = z2d(:,:)
-               znam = 'syap'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               syap(:,:,jl) = z2d(:,:)
-               znam = 'sxxap'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxxap(:,:,jl) = z2d(:,:)
-               znam = 'syyap'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               syyap(:,:,jl) = z2d(:,:)
-               znam = 'sxyap'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxyap(:,:,jl) = z2d(:,:)
-   
-               znam = 'sxvp'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxvp(:,:,jl) = z2d(:,:)
-               znam = 'syvp'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               syvp(:,:,jl) = z2d(:,:)
-               znam = 'sxxvp'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxxvp(:,:,jl) = z2d(:,:)
-               znam = 'syyvp'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               syyvp(:,:,jl) = z2d(:,:)
-               znam = 'sxyvp'//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxyvp(:,:,jl) = z2d(:,:)
-            END DO
-         ENDIF
-         ! END MV MP 2016
-
-         CALL iom_get( numrir, jpdom_autoglo, 'sxopw ' ,  sxopw  )
-         CALL iom_get( numrir, jpdom_autoglo, 'syopw ' ,  syopw  )
-         CALL iom_get( numrir, jpdom_autoglo, 'sxxopw' ,  sxxopw )
-         CALL iom_get( numrir, jpdom_autoglo, 'syyopw' ,  syyopw )
-         CALL iom_get( numrir, jpdom_autoglo, 'sxyopw' ,  sxyopw )
-
-         DO jl = 1, jpl 
-            WRITE(zchar,'(I2.2)') jl
-            DO jk = 1, nlay_i 
-               WRITE(zchar1,'(I2.2)') jk
-               znam = 'sxe'//'_il'//zchar1//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxe(:,:,jk,jl) = z2d(:,:)
-               znam = 'sye'//'_il'//zchar1//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sye(:,:,jk,jl) = z2d(:,:)
-               znam = 'sxxe'//'_il'//zchar1//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxxe(:,:,jk,jl) = z2d(:,:)
-               znam = 'syye'//'_il'//zchar1//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               syye(:,:,jk,jl) = z2d(:,:)
-               znam = 'sxye'//'_il'//zchar1//'_htc'//zchar
-               CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-               sxye(:,:,jk,jl) = z2d(:,:)
-            END DO
-         END DO
-         !
-      END IF
-      
-      ! clem: I do not understand why the following IF is needed
-      !       I suspect something inconsistent in the main code with option nn_icesal=1
-      IF( nn_icesal == 1 ) THEN
-         DO jl = 1, jpl 
-            sm_i(:,:,jl) = rn_icesal
-            DO jk = 1, nlay_i 
-               s_i(:,:,jk,jl) = rn_icesal
-            END DO
-         END DO
-      ENDIF
-      !
-      !CALL iom_close( numrir ) !clem: closed in icestp.F90
-      !
+      CALL iom_get( numrir, jpdom_autoglo, 'u_ice', u_ice )
+      CALL iom_get( numrir, jpdom_autoglo, 'v_ice', v_ice )
+
    END SUBROUTINE ice_rst_read
 
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 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8518)
@@ -206,6 +206,4 @@
                                     CALL ice_wri( 1 )         ! -- Ice outputs 
          !
-         IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please
-            &                       CALL iom_close( numrir )  ! close input ice restart file
          !
          IF( lrst_ice )             CALL ice_rst_write( kt )  ! -- Ice restart file 
@@ -251,5 +249,4 @@
       ierr = ierr + sbc_ice_alloc    ()      ! surface forcing 
       ierr = ierr + ice1D_alloc      ()      ! thermodynamics
-      ierr = ierr + ice_rdgrft_alloc ()      ! ridging/rafting
       !
       IF( lk_mpp    )   CALL mpp_sum( ierr )
@@ -258,18 +255,6 @@
       CALL ice_itd_init                ! ice thickness distribution initialization
       !
-      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
+      CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here)
       ! END MV MP 2016
       !                                ! Initial sea-ice state
@@ -283,4 +268,16 @@
       CALL ice_var_glo2eqv
       !
+      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   
+      !
       CALL ice_update_init             ! ice surface boundary condition
       !
@@ -296,4 +293,6 @@
       ELSEWHERE                     ;   rn_amax_2d(:,:) = rn_amax_s  ! SH
       END WHERE
+
+      IF( ln_rstart )   CALL iom_close( numrir )  ! close input ice restart file
       !
    END SUBROUTINE ice_init
@@ -429,5 +428,5 @@
       hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp
       hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp
-      hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp
+      hfx_err_rem(:,:) = 0._wp
       hfx_err_dif(:,:) = 0._wp
       wfx_err_sub(:,:) = 0._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 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8518)
@@ -413,5 +413,4 @@
          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw          )
          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub          )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_1d (1:nidx), hfx_err          )
          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res          )
          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif   )
@@ -499,5 +498,4 @@
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_snw_1d (1:nidx), hfx_snw        )
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_sub_1d (1:nidx), hfx_sub        )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_1d (1:nidx), hfx_err        )
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_res_1d (1:nidx), hfx_res        )
          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_err_dif_1d(1:nidx), hfx_err_dif )
@@ -585,5 +583,5 @@
       ENDIF
       !
-      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'STOP', 'ice_thd_init : rn_hinew should be >= rn_himin' )
+      IF ( rn_hinew < rn_himin )   CALL ctl_stop( 'ice_thd_init : rn_hinew should be >= rn_himin' )
       !
       IF(lwp) WRITE(numout,*)
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8518)
@@ -94,11 +94,13 @@
       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered at 0C 
       REAL(wp) ::   ztmelt_i                  ! ice melting temperature
-      REAL(wp) ::   zhsu
+      REAL(wp) ::   z1_hsu
       REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
       REAL(wp) ::   zdti_bnd = 1.e-4_wp       ! maximal authorized error on temperature 
+      REAL(wp) ::   zcpi                      ! Ice specific heat
+      REAL(wp) ::   zhi                       ! 
+      REAL(wp) ::   zhfx_err, zdq             ! diag errors on heat
       
       REAL(wp), DIMENSION(jpij)     ::   isnow       ! switch for presence (1) or absence (0) of snow
-      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! old surface temperature (before the iterative procedure )
-      REAL(wp), DIMENSION(jpij)     ::   ztsubit     ! surface temperature at previous iteration
+      REAL(wp), DIMENSION(jpij)     ::   ztsub       ! surface temperature at previous iteration
       REAL(wp), DIMENSION(jpij)     ::   zh_i        ! ice layer thickness
       REAL(wp), DIMENSION(jpij)     ::   zh_s        ! snow layer thickness
@@ -106,31 +108,28 @@
       REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
-      REAL(wp), DIMENSION(jpij)     ::   dzf         ! derivative of the surface flux function
+      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
       REAL(wp), DIMENSION(jpij)     ::   zdti        ! current error on temperature
-      REAL(wp), DIMENSION(jpij)     ::   zdifcase    ! case of the equation resolution (1->4)
       REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
-      REAL(wp), DIMENSION(jpij)     ::   zihic
       
+      REAL(wp), DIMENSION(jpij,nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice
+      REAL(wp), DIMENSION(jpij,nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztib        ! Old temperature in the ice
+      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztsb        ! Temporary temperature in the snow
+      REAL(wp), DIMENSION(jpij,nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
+      REAL(wp), DIMENSION(jpij,nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztcond_i    ! Ice thermal conductivity
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradtr_i    ! Radiation transmitted through the ice
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zradab_i    ! Radiation absorbed in the ice
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztib        ! Old temperature in the ice
       REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   zspeche_i   ! Ice specific heat
-      REAL(wp), DIMENSION(jpij,0:nlay_i)   ::   z_i         ! Vertical cotes of the layers in the ice
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradtr_s    ! Radiation transmited through the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zradab_s    ! Radiation absorbed in the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zkappa_s    ! Kappa factor in the snow
       REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   zeta_s      ! Eta factor in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   ztsb        ! Temporary temperature in the snow
-      REAL(wp), DIMENSION(jpij,0:nlay_s)   ::   z_s         ! Vertical cotes of the layers in the snow
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindterm    ! 'Ind'ependent term
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zindtbis    ! Temporary 'ind'ependent term
       REAL(wp), DIMENSION(jpij,nlay_i+3)   ::   zdiagbis    ! Temporary 'dia'gonal term
       REAL(wp), DIMENSION(jpij,nlay_i+3,3) ::   ztrid       ! Tridiagonal system terms
-      REAL(wp), DIMENSION(jpij)            ::   zdq, zq_ini, zhfx_err ! diag errors on heat
+      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat
       REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat
       
@@ -148,5 +147,4 @@
 
       ! --- diag error on heat diffusion - PART 1 --- !
-      zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
       DO ji = 1, nidx
          zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
@@ -167,17 +165,15 @@
       ! Ice / snow layers
       !--------------------
-
-      z_s(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st snow layer
-      z_i(:,0) = 0._wp   ! vert. coord. of the up. lim. of the 1st ice layer
-
-      DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
-         DO ji = 1 , nidx
-            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
-         END DO
-      END DO
-
-      DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
-         DO ji = 1 , nidx
-            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
+      z_s(1:nidx,1) = zh_s(1:nidx)
+      z_i(1:nidx,1) = zh_i(1:nidx)
+      DO jk = 2, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
+         DO ji = 1 , nidx
+            z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji)
+         END DO
+      END DO
+
+      DO jk = 2, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
+         DO ji = 1 , nidx
+            z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji)
          END DO
       END DO
@@ -187,34 +183,28 @@
       !------------------------------------------------------------------------------|
       !
-      !-------------------
-      ! Computation of i0
-      !-------------------
-      ! i0 describes the fraction of solar radiation which does not contribute
-      ! to the surface energy budget but rather penetrates inside the ice.
-      ! We assume that no radiation is transmitted through the snow
-      ! If there is no no snow
-      ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
-      ! zftrice = io.qsr_ice       is below the surface 
-      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
-      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
-      zhsu = 0.1_wp ! threshold for the computation of i0
+      z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
       DO ji = 1 , nidx
-         ! switches
-         isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 
-         ! hs > 0, isnow = 1
-         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) )     
-
-         i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) )
-      END DO
-
-      !-------------------------------------------------------
-      ! Solar radiation absorbed / transmitted at the surface
-      ! Derivative of the non solar flux
-      !-------------------------------------------------------
-      DO ji = 1 , nidx
-         zfsw   (ji)    =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
-         zftrice(ji)    =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
-         dzf    (ji)    = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
-         zqns_ice_b(ji) = qns_ice_1d(ji)                     ! store previous qns_ice_1d value
+         !-------------------
+         ! Computation of i0
+         !-------------------
+         ! i0 describes the fraction of solar radiation which does not contribute
+         ! to the surface energy budget but rather penetrates inside the ice.
+         ! We assume that no radiation is transmitted through the snow
+         ! If there is no no snow
+         ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface 
+         ! zftrice = io.qsr_ice       is below the surface 
+         ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
+         ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
+         zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) )     
+         i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) )
+
+         !-------------------------------------------------------
+         ! Solar radiation absorbed / transmitted at the surface
+         ! Derivative of the non solar flux
+         !-------------------------------------------------------
+         zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface
+         zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer
+         zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux 
+         zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value
       END DO
 
@@ -222,26 +212,19 @@
       ! Transmission - absorption of solar radiation in the ice
       !---------------------------------------------------------
-
-      DO ji = 1, nidx           ! snow initialization
-         zradtr_s(ji,0) = zftrice(ji)     ! radiation penetrating through snow
-      END DO
-
-      DO jk = 1, nlay_s          ! Radiation through snow
+      zradtr_s(1:nidx,0) = zftrice(1:nidx)
+      DO jk = 1, nlay_s
          DO ji = 1, nidx
             !                             ! radiation transmitted below the layer-th snow layer
-            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) )
+            zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) )
             !                             ! radiation absorbed by the layer-th snow layer
             zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk)
          END DO
       END DO
-
-      DO ji = 1, nidx           ! ice initialization
-         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) )
-      END DO
-
-      DO jk = 1, nlay_i          ! Radiation through ice
+         
+      zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) )
+      DO jk = 1, nlay_i 
          DO ji = 1, nidx
             !                             ! radiation transmitted below the layer-th ice layer
-            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) )
+            zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) )
             !                             ! radiation absorbed by the layer-th ice layer
             zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk)
@@ -249,7 +232,5 @@
       END DO
 
-      DO ji = 1, nidx           ! Radiation transmitted below the ice
-         ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 
-      END DO
+      ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i)   ! record radiation transmitted below the ice
 
       !------------------------------------------------------------------------------|
@@ -257,22 +238,9 @@
       !------------------------------------------------------------------------------|
       !
-      DO ji = 1, nidx        ! Old surface temperature
-         ztsub  (ji) =  t_su_1d(ji)                              ! temperature at the beg of iter pr.
-         ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
-         t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary
-         zdti   (ji) =  1000._wp                                 ! initial value of error
-      END DO
-
-      DO jk = 1, nlay_s       ! Old snow temperature
-         DO ji = 1 , nidx
-            ztsb(ji,jk) =  t_s_1d(ji,jk)
-         END DO
-      END DO
-
-      DO jk = 1, nlay_i       ! Old ice temperature
-         DO ji = 1 , nidx
-            ztib(ji,jk) =  t_i_1d(ji,jk)
-         END DO
-      END DO
+      ztsub  (1:nidx)   =  t_su_1d(1:nidx)                          ! temperature at the previous iter
+      t_su_1d(1:nidx)   =  MIN( t_su_1d(1:nidx), rt0 - ztsu_err )   ! necessary
+      zdti   (1:nidx)   =  1000._wp                                 ! initial value of error
+      ztsb   (1:nidx,:) =  t_s_1d(1:nidx,:)                         ! Old snow temperature
+      ztib   (1:nidx,:) =  t_i_1d(1:nidx,:)                         ! Old ice temperature
 
       iconv    =  0           ! number of iterations
@@ -289,12 +257,10 @@
          IF( ln_cndi_U64 ) THEN         !-- Untersteiner (1964) formula: k = k0 + beta.S/T
             DO ji = 1 , nidx
-               ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )
-               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
+               ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1 , nidx
-                  ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
-                     MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0)
-                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
+                  ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) /  &
+                     &                           MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) )
                END DO
             END DO
@@ -302,21 +268,18 @@
          ELSEIF( ln_cndi_P07 ) THEN     !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T
             DO ji = 1 , nidx
-               ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
-                  &                   - 0.011_wp * ( t_i_1d(ji,1) - rt0 )  
-               ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin )
+               ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )   &
+                  &                               - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) )
             END DO
             DO jk = 1, nlay_i-1
                DO ji = 1 , nidx
-                  ztcond_i(ji,jk) = rcdic +                                                                       & 
+                  ztcond_i(ji,jk) = MAX( zkimin, rcdic +                                                           & 
                      &                 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) )                              &
                      &                 / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 )   &
-                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  
-                  ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin )
+                     &               - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 )  )
                END DO
             END DO
             DO ji = 1 , nidx
-               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
-                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 )  
-               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin )
+               ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )   &
+                  &                        - 0.011_wp * ( t_bo_1d(ji) - rt0 ) )  
             END DO
          ENDIF
@@ -335,5 +298,5 @@
          SELECT CASE ( nn_monocat )
 
-         CASE (1,3) ! LIM3
+         CASE (1,3)
 
             zepsilon = 0.1_wp
@@ -404,7 +367,7 @@
          DO jk = 1, nlay_i
             DO ji = 1 , nidx
-               ztitemp(ji,jk)   = t_i_1d(ji,jk)
-               zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
-               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 )
+               ztitemp(ji,jk) = t_i_1d(ji,jk)
+               zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )
+               zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 )
             END DO
          END DO
@@ -425,5 +388,5 @@
             DO ji = 1 , nidx
                ! update of the non solar flux according to the update in T_su
-               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) )
+               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) )
             END DO
          ENDIF
@@ -507,5 +470,4 @@
                   !  case 1 : no surface melting - snow present                                  |
                   !------------------------------------------------------------------------------|
-                  zdifcase(ji)    =  1.0
                   numeqmin(ji)    =  1
                   numeqmax(ji)    =  nlay_i + nlay_s + 1
@@ -513,7 +475,7 @@
                   !!surface equation
                   ztrid(ji,1,1)  = 0.0
-                  ztrid(ji,1,2)  = dzf(ji) - zg1s * zkappa_s(ji,0)
+                  ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
                   ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
-                  zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)
+                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji)
 
                   !!first layer of snow equation
@@ -529,5 +491,4 @@
                   !------------------------------------------------------------------------------|
                   !
-                  zdifcase(ji)    =  2.0
                   numeqmin(ji)    =  2
                   numeqmax(ji)    =  nlay_i + nlay_s + 1
@@ -552,5 +513,4 @@
                   !------------------------------------------------------------------------------|
                   !
-                  zdifcase(ji)      =  3.0
                   numeqmin(ji)      =  nlay_s + 1
                   numeqmax(ji)      =  nlay_i + nlay_s + 1
@@ -558,7 +518,7 @@
                   !!surface equation	
                   ztrid(ji,numeqmin(ji),1)   =  0.0
-                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1    
+                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
                   ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
+                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji)
 
                   !!first layer of ice equation
@@ -572,5 +532,5 @@
                   IF ( nlay_i == 1 ) THEN
                      ztrid(ji,numeqmin(ji),1)    =  0.0
-                     ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
+                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
                      ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
                      ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
@@ -589,5 +549,4 @@
                   !------------------------------------------------------------------------------|
                   !
-                  zdifcase(ji)    =  4.0
                   numeqmin(ji)    =  nlay_s + 2
                   numeqmax(ji)    =  nlay_i + nlay_s + 1
@@ -662,5 +621,5 @@
             ! surface temperature
             isnow(ji)   = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) )
-            ztsubit(ji) = t_su_1d(ji)
+            ztsub(ji) = t_su_1d(ji)
             IF( t_su_1d(ji) < rt0 ) &
                t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
@@ -676,5 +635,5 @@
          DO ji = 1 , nidx
             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
-            zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
+            zdti   (ji) =  ABS( t_su_1d(ji) - ztsub(ji) )     
          END DO
 
@@ -754,31 +713,20 @@
 
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
+      !     hfx_dif = Heat flux used to warm/cool ice in W.m-2
+      !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
       DO ji = 1, nidx
-         zdq(ji)        = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
-            &                              SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
+         zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
+            &                   SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
+
          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-            zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
+            zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
          ELSE                          ! case T_su = 0degC
-            zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
+            zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)
          ENDIF
-         hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
+         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji)
 
          ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation)
-         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji)
-      END DO 
-
-      !-----------------------------------------
-      ! Heat flux used to warm/cool ice in W.m-2
-      !-----------------------------------------
-      DO ji = 1, nidx
-         IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC
-            hfx_dif_1d(ji) = hfx_dif_1d(ji)  +   &
-               &            ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
-         ELSE                          ! case T_su = 0degC
-            hfx_dif_1d(ji) = hfx_dif_1d(ji) +    &
-               &             ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji)
-         ENDIF
-         ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation
-         hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji)
+         hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err
+
       END DO   
       !
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8518)
@@ -29,4 +29,11 @@
    PUBLIC   ice_thd_sal        ! called by icethd module
    PUBLIC   ice_thd_sal_init   ! called by ice_init
+   
+   ! ** namelist (namsal) **
+   LOGICAL  ::   ln_icedS         ! activate gravity drainage and flushing (T) or not (F)
+   REAL(wp) ::   rn_sal_gd        !    restoring salinity for gravity drainage [PSU]
+   REAL(wp) ::   rn_time_gd       !    restoring time constant for gravity drainage (= 20 days) [s]
+   REAL(wp) ::   rn_sal_fl        !    restoring salinity for flushing [PSU]
+   REAL(wp) ::   rn_time_fl       !    restoring time constant for gravity drainage (= 10 days) [s]
 
    !!----------------------------------------------------------------------
@@ -115,5 +122,5 @@
       !!
       NAMELIST/namice_sal/ ln_icedS , nn_icesal , rn_icesal, rn_sal_gd, rn_time_gd,   &
-         &                rn_sal_fl, rn_time_fl, rn_simax , rn_simin 
+         &                 rn_sal_fl, rn_time_fl, rn_simax , rn_simin 
       !!-------------------------------------------------------------------
       !
@@ -143,4 +150,9 @@
       ENDIF
       !
+      IF( ln_icedS .AND. nn_icesal == 1 ) THEN
+         ln_icedS = .FALSE.
+         CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)')
+      ENDIF
+      !
    END SUBROUTINE ice_thd_sal_init
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8518)
@@ -36,8 +36,7 @@
 !!gm end
    USE sbccpl         ! Surface boundary condition: coupled interface
-   USE icealb         ! albedo parameters
+   USE icealb         ! sea-ice: albedo parameters
    USE traqsr         ! add penetration of solar flux in the calculation of heat budget
-   USE domvvl         ! Variable volume
-   USE icectl         ! ???
+   USE icectl         ! sea-ice: control prints
    USE bdy_oce , ONLY : ln_bdy
    !
@@ -58,5 +57,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress     [N/m2]
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmod_io              ! modulus of the ice-ocean velocity   [m/s]
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   soce_0  , sice_0     ! cst SSS and ice salinity (levitating sea-ice) 
 
    !! * Substitutions
@@ -73,9 +71,9 @@
       !!             ***  ROUTINE ice_update_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( soce_0(jpi,jpj) , utau_oce(jpi,jpj) ,                       &
-         &      sice_0(jpi,jpj) , vtau_oce(jpi,jpj) , tmod_io(jpi,jpj), STAT=ice_update_alloc)
-         !
+      ALLOCATE( utau_oce(jpi,jpj), vtau_oce(jpi,jpj), tmod_io(jpi,jpj), STAT=ice_update_alloc )
+      !
       IF( lk_mpp                )   CALL mpp_sum( ice_update_alloc )
       IF( ice_update_alloc /= 0 )   CALL ctl_warn('ice_update_alloc: failed to allocate arrays')
+      !
    END FUNCTION ice_update_alloc
 
@@ -206,12 +204,14 @@
       !
       alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-
-      !                    ! conservation test
-      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final( 'iceupdate' )
-      !                    ! control prints
-      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 3, ' - Final state ice_update - ' )
-      IF( ln_ctl    )   CALL ice_prt3D( 'iceupdate' )
-      !
-      IF( nn_timing == 1 )   CALL timing_stop('ice_update_flx')
+      !
+      IF( lrst_ice ) THEN                       !* write snwice_mass fields in the restart file
+         CALL update_rst( 'WRITE', kt )
+      ENDIF
+      !
+      ! controls
+      IF( ln_icediachk .AND. .NOT. ln_bdy)   CALL ice_cons_final('iceupdate')                                       ! conservation
+      IF( ln_icectl                      )   CALL ice_prt       (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints
+      IF( ln_ctl                         )   CALL ice_prt3D     ('iceupdate')                                       ! prints
+      IF( nn_timing == 1                 )   CALL timing_stop   ('ice_update_flx')                                  ! timing
       !
    END SUBROUTINE ice_update_flx
@@ -320,81 +320,54 @@
       IF( ice_update_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'ice_update_init : unable to allocate standard arrays' )
       !
-      soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating case 0 (i.e. virtual salt flux)
-      sice_0(:,:) = sice
-      WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &   ! reduced values in the Baltic Sea area
-         &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
-         soce_0(:,:) = 4._wp
-         sice_0(:,:) = 2._wp
-      END WHERE
-      !
-      IF( .NOT.ln_rstart ) THEN              ! set 
-         !
-         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass
-         snwice_mass_b(:,:) = snwice_mass(:,:)
-         !
-         IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
-            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
-            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-
-!!gm I really don't like this stuff here...  Find a way to put that elsewhere or differently
-!!gm
-            IF( .NOT.ln_linssh ) THEN
-               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  )
-                  e3t_b(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshb(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1) )  )
-               END DO
-               e3t_a(:,:,:) = e3t_b(:,:,:)
-!!gm  we are in no-restart case, so sshn=sshb   ==>> faster calculation:
-!!    REAL(wp) ::   ze3t   ! local scalar
-!!    REAL(wp), DIMENSION(jpi,jpj) ::   z2d   ! workspace
-!!
-!!             WHERE( ht_0(:,:) > 0 )   ;   z2d(:,:) = 1._wp + sshn(:,:)*tmask(:,:,1) / ht_0(:,:)
-!!             ELSEWHERE                ;   z2d(:,:) = 1._wp
-!!             END WHERE
-!!             DO jk = 1,jpkm1                     ! adjust initial vertical scale factors                
-!!                ze3t = e3t_0(:,:,jk) * z2d(:,:)
-!!                e3t_n(:,:,jk) = ze3t
-!!                e3t_b(:,:,jk) = ze3t
-!!                e3t_a(:,:,jk) = ze3t
-!!             END DO
-!!gm  but since it is only done at the initialisation....  just the following can be acceptable:
-!               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors
-!                  e3t_n(:,:,jk) = e3t_0(:,:,jk) * (  1._wp + sshn(:,:)*tmask(:,:,1) / (ht_0(:,:) + 1._wp - tmask(:,:,1))  )
-!               END DO
-!               e3t_b(:,:,:) = e3t_n(:,:,:)
-!               e3t_a(:,:,:) = e3t_n(:,:,:)
-!!gm end               
-               ! Reconstruction of all vertical scale factors at now and before time-steps
-               ! =========================================================================
-               ! Horizontal scale factor interpolations
-               ! --------------------------------------
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' )
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' )
-               CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' )
-               ! Vertical scale factor interpolations
-                 ! ------------------------------------
-               CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W'  )
-               CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
-               CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
-               CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
-               CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
-               ! t- and w- points depth
-               ! ----------------------
-!!gm not sure of that....
-               gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
-               gdepw_n(:,:,1) = 0.0_wp
-               gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
-               DO jk = 2, jpk
-                  gdept_n(:,:,jk) = gdept_n(:,:,jk-1) + e3w_n(:,:,jk  )
-                  gdepw_n(:,:,jk) = gdepw_n(:,:,jk-1) + e3t_n(:,:,jk-1)
-                  gde3w_n(:,:,jk) = gdept_n(:,:,jk  ) - sshn (:,:)
-               END DO
+      CALL update_rst( 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE ice_update_init
+
+   SUBROUTINE update_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rhg_evp_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER  ::   iter   ! local integer
+      INTEGER  ::   id1    ! local integer
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'snwice_mass' , ldstop = .FALSE. )
+            !
+            IF( id1 > 0 ) THEN                       ! fields exist
+               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass   )
+               CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b )
+            ELSE                                     ! start from rest
+               IF(lwp) WRITE(numout,*) '   ==>>   previous run without snow-ice mass output then set it'
+               snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
+               snwice_mass_b(:,:) = snwice_mass(:,:)
             ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set the snow-ice mass'
+            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:) )
+            snwice_mass_b(:,:) = snwice_mass(:,:)
          ENDIF
-      ENDIF ! .NOT. ln_rstart
-      !
-   END SUBROUTINE ice_update_init
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- update-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass   )
+         CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b )
+         !
+      ENDIF
+      !
+   END SUBROUTINE update_rst
 
 #else
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8518)
@@ -192,5 +192,5 @@
       CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !  
       CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !  
-      CALL iom_put ('hfxerr'     , hfx_err(:,:)         )   !  
+      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   !  
       CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !  
       
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 8517)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/IOM/in_out_manager.F90	(revision 8518)
@@ -88,6 +88,9 @@
    INTEGER ::   nitrst                !: time step at which restart file should be written
    LOGICAL ::   lrst_oce              !: logical to control the oce restart write 
+   LOGICAL ::   lrst_ice              !: logical to control the ice restart write 
    INTEGER ::   numror = 0            !: logical unit for ocean restart (read). Init to 0 is needed for SAS (in daymod.F90)
+   INTEGER ::   numrir                !: logical unit for ice   restart (read)
    INTEGER ::   numrow                !: logical unit for ocean restart (write)
+   INTEGER ::   numriw                !: logical unit for ice   restart (write)
    INTEGER ::   nrst_lst              !: number of restart to output next
 
