Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml	(revision 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml	(revision 8562)
@@ -129,5 +129,5 @@
          <field id="albedo"       long_name="Mean albedo over sea ice and ocean"                                                                                              unit=""             />
 
-		 <field id="iceamp"       long_name="melt pond fraction"                                           standard_name="sea_ice_meltpond_fraction"                          unit="%"            /> 
+	 <field id="iceamp"       long_name="melt pond fraction"                                           standard_name="sea_ice_meltpond_fraction"                          unit="%"            /> 
          <field id="icevmp"       long_name="melt pond volume"                                             standard_name="sea_ice_meltpond_volume"                            unit="m"            /> 
 
@@ -154,36 +154,36 @@
          <field id="tau_icebfr"   long_name="ice friction on ocean bottom for landfast ice"                unit="N/2"      />
 
-         <field id="icetrp"       long_name="ice volume transport"                                         unit="m/day"          />
-         <field id="snwtrp"       long_name="snw volume transport"                                         unit="m/day"          />
-         <field id="saltrp"       long_name="salt content transport"                                       unit="1e-3*kg/m2/day" />
+         <field id="icetrp"       long_name="ice mass transport"                                           unit="kg/m2/s"          />
+         <field id="snwtrp"       long_name="snw mass transport"                                           unit="kg/m2/s"          />
+         <field id="saltrp"       long_name="salt     transport"                                           unit="1e-3*kg/m2/s" />
          <field id="deitrp"       long_name="advected ice enthalpy"                                        unit="W/m2"           />
          <field id="destrp"       long_name="advected snw enthalpy"                                        unit="W/m2"           />
 
-         <field id="sfxbri"       long_name="brine salt flux"                                              unit="1e-3*kg/m2/day" />
-         <field id="sfxdyn"       long_name="salt flux from ridging rafting"                               unit="1e-3*kg/m2/day" />
-         <field id="sfxres"       long_name="salt flux from lipupdate (resultant)"                         unit="1e-3*kg/m2/day" />
-         <field id="sfxbog"       long_name="salt flux from bot growth"                                    unit="1e-3*kg/m2/day" />
-         <field id="sfxbom"       long_name="salt flux from bot melt"                                      unit="1e-3*kg/m2/day" />
-         <field id="sfxsum"       long_name="salt flux from surf melt"                                     unit="1e-3*kg/m2/day" />
-         <field id="sfxlam"       long_name="salt flux from lateral melt"                                  unit="1e-3*kg/m2/day" />
-         <field id="sfxsni"       long_name="salt flux from snow-ice formation"                            unit="1e-3*kg/m2/day" />
-         <field id="sfxopw"       long_name="salt flux from open water ice formation"                      unit="1e-3*kg/m2/day" />
-         <field id="sfxsub"       long_name="salt flux from sublimation"                                   unit="1e-3*kg/m2/day" />
-         <field id="sfx"          long_name="Salt flux from sea ice"                                       unit="1e-3*kg/m2/day" />
-
-         <field id="vfxbog"       long_name="daily bottom thermo ice prod."                                unit="m/day"   />
-         <field id="vfxdyn"       long_name="daily  dynamic ice prod."                                     unit="m/day"   />
-         <field id="vfxopw"       long_name="daily lateral thermo ice prod."                               unit="m/day"   />
-         <field id="vfxsni"       long_name="daily snowice ice prod."                                      unit="m/day"   />
-         <field id="vfxsum"       long_name="surface melt"                                                 unit="m/day"   />
-         <field id="vfxlam"       long_name="lateral melt"                                                 unit="m/day"   />
-         <field id="vfxbom"       long_name="bottom melt"                                                  unit="m/day"   />
-         <field id="vfxres"       long_name="daily resultant ice prod./melting from limupdate"             unit="m/day"   />
-         <field id="vfxice"       long_name="ice melt/growth"                                              unit="m/day"   />
-         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   />
-         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   />
-         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   />
-         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   />
-         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(20cm) + open water"       unit="m/day"   />
+         <field id="sfxbri"       long_name="salt flux from brines"                                        unit="1e-3*kg/m2/s" />
+         <field id="sfxdyn"       long_name="salt flux from ridging rafting"                               unit="1e-3*kg/m2/s" />
+         <field id="sfxres"       long_name="salt flux from lipupdate (resultant)"                         unit="1e-3*kg/m2/s" />
+         <field id="sfxbog"       long_name="salt flux from bot growth"                                    unit="1e-3*kg/m2/s" />
+         <field id="sfxbom"       long_name="salt flux from bot melt"                                      unit="1e-3*kg/m2/s" />
+         <field id="sfxsum"       long_name="salt flux from surf melt"                                     unit="1e-3*kg/m2/s" />
+         <field id="sfxlam"       long_name="salt flux from lateral melt"                                  unit="1e-3*kg/m2/s" />
+         <field id="sfxsni"       long_name="salt flux from snow-ice formation"                            unit="1e-3*kg/m2/s" />
+         <field id="sfxopw"       long_name="salt flux from open water ice formation"                      unit="1e-3*kg/m2/s" />
+         <field id="sfxsub"       long_name="salt flux from sublimation"                                   unit="1e-3*kg/m2/s" />
+         <field id="sfx"          long_name="Salt flux from sea ice"                                       unit="1e-3*kg/m2/s" />
+
+         <field id="vfxbog"       long_name="bottom thermo ice prod."                                      unit="kg/m2/s"   />
+         <field id="vfxdyn"       long_name="dynamic ice prod."                                            unit="kg/m2/s"   />
+         <field id="vfxopw"       long_name="lateral thermo ice prod."                                     unit="kg/m2/s"   />
+         <field id="vfxsni"       long_name="snowice ice prod."                                            unit="kg/m2/s"   />
+         <field id="vfxsum"       long_name="surface melt"                                                 unit="kg/m2/s"   />
+         <field id="vfxlam"       long_name="lateral melt"                                                 unit="kg/m2/s"   />
+         <field id="vfxbom"       long_name="bottom melt"                                                  unit="kg/m2/s"   />
+         <field id="vfxres"       long_name="resultant ice prod./melting"                                  unit="kg/m2/s"   />
+         <field id="vfxice"       long_name="ice melt/growth"                                              unit="kg/m2/s"   />
+         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="kg/m2/s"   />
+         <field id="vfxsub"       long_name="snw sublimation"                                              unit="kg/m2/s"   />
+         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="kg/m2/s"   />
+         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="kg/m2/s"   />
+         <field id="vfxthin"      long_name="thermo ice prod. for thin ice(20cm) + open water"             unit="kg/m2/s"   />
 
          <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   />
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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8562)
@@ -172,5 +172,4 @@
    !                                     !!** ice-dynamics namelist (namdyn) **
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
-   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
    LOGICAL , PUBLIC ::   ln_landfast      !: landfast ice parameterization (T or F) 
    REAL(wp), PUBLIC ::   rn_gamma         !:    fraction of ocean depth that ice must reach to initiate landfast ice
@@ -184,8 +183,9 @@
    REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 
    !
-   !                                     !!** ice-thermodynamics namelist (namthd) **
+   !                                     !!** ice-surface forcing namelist (namforcing) **
                                           ! -- icethd_dh -- !
    REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice
                                           ! -- icethd -- !
+   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
    INTEGER , PUBLIC ::   nn_iceflx        !: Redistribute heat flux over ice categories
    !                                      !   =-1  Do nothing (needs N(cat) fluxes)
@@ -235,5 +235,5 @@
    !                                     !!** define arrays
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce,v_oce !: surface ocean velocity used in ice dynamics
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hicol       !: ice collection thickness accreted in leads
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_i_new    !: ice collection thickness accreted in leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strength    !: ice strength
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element
@@ -437,5 +437,5 @@
       ! stay within Fortran's max-line length limit.
       ii = 1
-      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) , hicol    (jpi,jpj) ,                        &
+      ALLOCATE( u_oce   (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new (jpi,jpj) ,                        &
          &      strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,  &
          &      delta_i (jpi,jpj) , divu_i   (jpi,jpj) , shear_i  (jpi,jpj) , STAT=ierr(ii) )
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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8562)
@@ -208,5 +208,5 @@
       ii = ii + 1
       ALLOCATE( t_s_1d  (jpij,nlay_s)     , t_i_1d (jpij,nlay_i)     , s_i_1d(jpij,nlay_i) ,  &            
-         &      e_i_1d  (jpij,nlay_i+1)   , e_s_1d (jpij,nlay_s)     ,                        &
+         &      e_i_1d  (jpij,nlay_i)     , e_s_1d (jpij,nlay_s)     ,                        &
          &      eh_i_old(jpij,0:nlay_i+1) , h_i_old(jpij,0:nlay_i+1) , STAT=ierr(ii) )
       !
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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedia.F90	(revision 8562)
@@ -76,5 +76,5 @@
       IF( kt == nit000 .AND. lwp ) THEN
          WRITE(numout,*)
-         WRITE(numout,*)'icedia : outpout ice diagnostics (integrated over the domain)'
+         WRITE(numout,*)'icedia: output ice diagnostics (integrated over the domain)'
          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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8562)
@@ -87,4 +87,5 @@
             WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )   tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
          END DO
+         IF( iom_use('tau_icebfr') )   CALL iom_put( 'tau_icebfr', tau_icebfr )  
       ENDIF
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv.F90	(revision 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv.F90	(revision 8562)
@@ -105,9 +105,9 @@
       diag_trp_vi (:,:) = SUM(     v_i  (:,:,:)          - v_i_b  (:,:,:)                  , dim=3 ) * r1_rdtice
       diag_trp_vs (:,:) = SUM(     v_s  (:,:,:)          - v_s_b  (:,:,:)                  , dim=3 ) * r1_rdtice
-      IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi * rday  )         ! ice volume transport
-      IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs * rday  )         ! snw volume transport
-      IF( iom_use('saltrp') )   CALL iom_put( "saltrp" , diag_trp_smv * rday * rhoic ) ! salt content transport
-      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( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi )          ! ice volume transport
+      IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs )          ! snw volume transport
+      IF( iom_use('saltrp') )   CALL iom_put( "saltrp" , diag_trp_smv * rhoic ) ! salt content transport
+      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
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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8562)
@@ -128,6 +128,5 @@
          !    utau_ice, vtau_ice = surface ice stress [N/m2]
          !------------------------------------------------!
-         CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )
-                                      
+                                        CALL ice_forcing_tau( kt, ksbc, utau_ice, vtau_ice )          
          !-------------------------------------!
          ! --- ice dynamics and advection  --- !
@@ -136,15 +135,15 @@
          CALL ice_rst_opn( kt )     ! Open Ice restart file (if necessary) 
          !
-         IF( ln_icedyn .AND. .NOT.lk_c1d )   CALL ice_dyn( kt )            ! -- Ice dynamics
+         IF( ln_icedyn .AND. .NOT.lk_c1d )   &
+            &                           CALL ice_dyn( kt )            ! -- Ice dynamics
          !
          !                          !==  lateral boundary conditions  ==!
 #if defined key_agrif
-         IF( .NOT. Agrif_Root()     )        CALL agrif_interp_lim3('T')   ! -- AGRIF 
+         IF( .NOT. Agrif_Root()     )   CALL agrif_interp_lim3('T')   ! -- AGRIF 
 #endif
-         IF( ln_icethd .AND. ln_bdy )        CALL bdy_ice( kt )            ! -- bdy ice thermo
+         IF( ln_icethd .AND. ln_bdy )   CALL bdy_ice( kt )            ! -- bdy ice thermo
          !
          !
          !                          !==  previous lead fraction and ice volume for flux calculations
-         !
          CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation
          CALL ice_var_agg(1)             ! at_i for coupling 
@@ -165,40 +164,35 @@
          !    qprec_ice, qevap_ice
          !------------------------------------------------------!
-                                    CALL ice_forcing_flx( kt, ksbc )
-
+                                        CALL ice_forcing_flx( kt, ksbc )
          !----------------------------!
          ! --- ice thermodynamics --- !
          !----------------------------!
-         IF( ln_icethd )            CALL ice_thd( kt )          ! -- Ice thermodynamics      
-
-         ! MV MP 2016
-         IF ( ln_pnd )              CALL lim_mp( kt )           ! -- Melt ponds
-         ! END MV MP 2016
-
-         IF( ln_icethd )            CALL ice_cor( kt , 2 )      ! -- Corrections
-         ! ---
+         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
+         !
+         IF ( ln_pnd )                  CALL lim_mp( kt )             ! -- Melt ponds
+         !
+         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
 # if defined key_agrif
-         IF( .NOT. Agrif_Root() )   CALL agrif_update_lim3( kt )
+         IF( .NOT. Agrif_Root() )       CALL agrif_update_lim3( kt )
 # endif
-                                    CALL ice_var_glo2eqv        ! necessary calls (at least for coupling)
-                                    CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling)
-                                    !
+         CALL ice_var_glo2eqv        ! necessary calls (at least for coupling)
+         CALL ice_var_agg( 2 )       ! necessary calls (at least for coupling)
+                                     !
 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
 !!       moreover it should only be called at the update frequency (as in agrif_lim3_update.F90)
 !!# if defined key_agrif
-!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ChildGrid_To_ParentGrid()
+!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ChildGrid_To_ParentGrid()
 !!# endif
-                                    CALL ice_update_flx( kt )   ! -- Update ocean surface mass, heat and salt fluxes
+                                        CALL ice_update_flx( kt )     ! -- Update ocean surface mass, heat and salt fluxes
 !!# if defined key_agrif
-!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid()
+!!         IF( .NOT. Agrif_Root() )     CALL Agrif_ParentGrid_To_ChildGrid()
 !!# endif
-         IF( ln_icediahsb )         CALL ice_dia( kt )          ! -- Diagnostics and outputs 
-         !
-                                    CALL ice_wri( 1 )           ! -- Ice outputs 
-         !
-         !
-         IF( lrst_ice )             CALL ice_rst_write( kt )    ! -- Ice restart file 
-         !
-         IF( ln_icectl )            CALL ice_ctl( kt )          ! -- alerts in case of model crash
+         IF( ln_icediahsb )             CALL ice_dia( kt )            ! -- Diagnostics outputs 
+         !
+                                        CALL ice_wri( kt )            ! -- Ice outputs 
+         !
+         IF( lrst_ice )                 CALL ice_rst_write( kt )      ! -- Ice restart file 
+         !
+         IF( ln_icectl )                CALL ice_ctl( kt )            ! -- alerts in case of model crash
          !
       ENDIF   ! End sea-ice time step only
@@ -207,10 +201,8 @@
       ! --- Ocean time step --- !
       !-------------------------!
-      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere
-      !    using before instantaneous surf. currents
-      IF( ln_icedyn )               CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )
+      IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
-      IF( nn_timing == 1 ) CALL timing_stop('ice_stp')
+      IF( nn_timing == 1 )   CALL timing_stop('ice_stp')
       !
    END SUBROUTINE ice_stp
@@ -246,7 +238,6 @@
       CALL ice_itd_init                ! ice thickness distribution initialization
       !
-      ! MV MP 2016
       CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here)
-      ! END MV MP 2016
+      !
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8562)
@@ -116,10 +116,4 @@
       END SELECT
 
-      DO ji = 1, nidx
-         icount (ji,:) = 0
-         zdh_s_mel(ji) = 0._wp
-         e_i_1d(ji,nlay_i+1) = 0._wp ! Initialize enthalpy at nlay_i+1
-      END DO
-
       ! initialize layer thicknesses and enthalpies
       h_i_old (1:nidx,0:nlay_i+1) = 0._wp
@@ -227,4 +221,5 @@
       ! If heat still available (zq_su > 0), then melt more snow
       zdeltah(1:nidx,:) = 0._wp
+      zdh_s_mel(1:nidx) = 0._wp
       DO jk = 1, nlay_s
          DO ji = 1, nidx
@@ -443,6 +438,4 @@
 
                dh_i_bott(ji)      = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
-
-               e_i_1d(ji,nlay_i+1) = -zEi * rhoic                        ! New ice energy of melting (J/m3, >0)
                
             END DO
@@ -475,5 +468,5 @@
 
             ! update heat content (J.m-2) and layer thickness
-            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * e_i_1d(ji,nlay_i+1)
+            eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * (-zEi * rhoic)
             h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
 
@@ -654,7 +647,5 @@
 
       ! --- ensure that a_i = 0 where ht_i = 0 ---
-      DO ji = 1, nidx
-         IF( ht_i_1d(ji) == 0._wp )   a_i_1d(ji) = 0._wp
-      END DO         
+      WHERE( ht_i_1d(1:nidx) == 0._wp )   a_i_1d(1:nidx) = 0._wp
       !
    END SUBROUTINE ice_thd_dh
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8562)
@@ -123,6 +123,6 @@
       ! 3) Collection thickness of ice formed in leads and polynyas
       !------------------------------------------------------------------------------!    
-      ! hicol is the thickness of new ice formed in open water
-      ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1)
+      ! ht_i_new is the thickness of new ice formed in open water
+      ! ht_i_new can be either prescribed (frazswi = 0) or computed (frazswi = 1)
       ! Frazil ice forms in open water, is transported by wind
       ! accumulates at the edge of the consolidated ice edge
@@ -136,6 +136,6 @@
 
       ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp )   ;   hicol(:,:) = rn_hinew
-      ELSEWHERE                     ;   hicol(:,:) = 0._wp
+      WHERE( qlead(:,:) < 0._wp )   ;   ht_i_new(:,:) = rn_hinew
+      ELSEWHERE                     ;   ht_i_new(:,:) = 0._wp
       END WHERE
 
@@ -145,5 +145,5 @@
          ! Physical constants
          !--------------------
-         hicol(:,:) = 0._wp
+         ht_i_new(:,:) = 0._wp
 
          zhicrit = 0.04 ! frazil ice thickness
@@ -192,14 +192,14 @@
                   ! Iterative procedure
                   !---------------------
-                  hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
+                  ht_i_new(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
                      &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
 
                   iter = 1
                   DO WHILE ( iter < 20 ) 
-                     zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   &
-                        &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2
-                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
-
-                     hicol(ji,jj) = hicol(ji,jj) - zf / MAX( zfp, epsi20 )
+                     zf  = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) -   &
+                        &    ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2
+                     zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0 * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
+
+                     ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )
                      iter = iter + 1
                   END DO
@@ -210,5 +210,5 @@
          END DO 
          ! 
-         CALL lbc_lnk_multi( zvrel, 'T', 1., hicol, 'T', 1.  )
+         CALL lbc_lnk_multi( zvrel, 'T', 1., ht_i_new, 'T', 1.  )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -252,5 +252,5 @@
          CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw     )
          CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , hicol       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , ht_i_new    )
          CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx) , zvrel       )
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8562)
@@ -78,6 +78,6 @@
       !  1) Cumulative integral of old enthalpy * thickness and layers interfaces
       !--------------------------------------------------------------------------
-      zeh_cum0(:,0:nlay_i+2) = 0._wp 
-      zh_cum0 (:,0:nlay_i+2) = 0._wp
+      zeh_cum0(1:nidx,0) = 0._wp 
+      zh_cum0 (1:nidx,0) = 0._wp
       DO jk0 = 1, nlay_i+2
          DO ji = 1, nidx
@@ -96,5 +96,5 @@
 
       ! new layers interfaces
-      zh_cum1(:,0:nlay_i) = 0._wp
+      zh_cum1(1:nidx,0) = 0._wp
       DO jk1 = 1, nlay_i
          DO ji = 1, nidx
@@ -103,5 +103,5 @@
       END DO
 
-      zeh_cum1(:,0:nlay_i) = 0._wp 
+      zeh_cum1(1:nidx,0) = 0._wp 
       ! new cumulative q*h => linear interpolation
       DO jk0 = 1, nlay_i+1
@@ -117,5 +117,5 @@
       END DO
       ! to ensure that total heat content is strictly conserved, set:
-      zeh_cum1(:,nlay_i) = zeh_cum0(:,nlay_i+2) 
+      zeh_cum1(1:nidx,nlay_i) = zeh_cum0(1:nidx,nlay_i+2) 
 
       ! new enthalpies
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90	(revision 8562)
@@ -84,11 +84,11 @@
       !!-------------------------------------------------------------------
       INTEGER ::   ji, jk         ! spatial loop index
-      INTEGER ::   numeq          ! current reference number of equation
-      INTEGER ::   minnumeqmin, maxnumeqmax
+      INTEGER ::   jm             ! current reference number of equation
+      INTEGER ::   jm_mint, jm_maxt
       INTEGER ::   iconv          ! number of iterations in iterative procedure
       INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
       
-      INTEGER, DIMENSION(jpij) ::   numeqmin   ! reference number of top equation
-      INTEGER, DIMENSION(jpij) ::   numeqmax   ! reference number of bottom equation
+      INTEGER, DIMENSION(jpij) ::   jm_min    ! reference number of top equation
+      INTEGER, DIMENSION(jpij) ::   jm_max    ! reference number of bottom equation
       
       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
@@ -113,5 +113,5 @@
       REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface
       REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface
-      REAL(wp), DIMENSION(jpij)     ::   zf          ! surface flux function
+      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function
       REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function
       REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice
@@ -176,5 +176,4 @@
       ! 2) Radiation
       !-------------
-      !
       z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0
       DO ji = 1, nidx
@@ -224,8 +223,7 @@
       iconv    =  0          ! number of iterations
       zdti_max =  1000._wp   ! maximal value of error on all points
-      !                                                          !----------------------------!
+      !                                                          !============================!
       DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )   ! Iterative procedure begins !
-         !                                                       !----------------------------!
-         !
+         !                                                       !============================!
          iconv = iconv + 1
          !
@@ -272,7 +270,6 @@
          ! Used in mono-category case only to simulate an ITD implicitly
          ! Fichefet and Morales Maqueda, JGR 1997
-
          zghe(1:nidx) = 1._wp
-         
+         !
          SELECT CASE ( nn_monocat )
 
@@ -281,16 +278,9 @@
             zepsilon = 0.1_wp
             DO ji = 1, nidx
-   
-               ! Mean sea ice thermal conductivity
-               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )
-
-               ! Effective thickness he (zhe)
-               zhe  = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )
-
-               ! G(he)
+               zcnd_i = SUM( ztcond_i(ji,:) ) / REAL( nlay_i+1, wp )                             ! Mean sea ice thermal conductivity
+               zhe = ( rn_cnd_s * ht_i_1d(ji) + zcnd_i * ht_s_1d(ji) ) / ( rn_cnd_s + zcnd_i )   ! Effective thickness he (zhe)
                IF( zhe >=  zepsilon * 0.5_wp * EXP(1._wp) ) THEN
-                  zghe(ji) =  MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )
+                  zghe(ji) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) )    ! G(he)
                ENDIF
-   
             END DO
 
@@ -303,5 +293,5 @@
          DO jk = 0, nlay_s-1
             DO ji = 1, nidx
-               zkappa_s(ji,jk)  = zghe(ji) * rn_cnd_s * z1_h_s(ji)
+               zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji)
             END DO
          END DO
@@ -322,5 +312,5 @@
          END DO
          DO ji = 1, nidx   ! Snow-ice interface
-            zkappa_i(ji,0)     = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
+            zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) )
          END DO
          !
@@ -337,5 +327,5 @@
          DO jk = 1, nlay_s
             DO ji = 1, nidx
-               zeta_s(ji,jk)  = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
+               zeta_s(ji,jk) = rdt_ice * r1_rhosn * r1_cpic * z1_h_s(ji)
             END DO
          END DO
@@ -352,5 +342,5 @@
 
          DO ji = 1, nidx
-            zf(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
+            zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH)
          END DO
          !
@@ -359,39 +349,32 @@
          !----------------------------
          !!layer denotes the number of the layer in the snow or in the ice
-         !!numeq denotes the reference number of the equation in the tridiagonal
+         !!jm denotes the reference number of the equation in the tridiagonal
          !!system, terms of tridiagonal system are indexed as following :
          !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one
 
          !!ice interior terms (top equation has the same form as the others)
-
-         DO numeq=1,nlay_i+3
-            DO ji = 1, nidx
-               ztrid(ji,numeq,1) = 0.
-               ztrid(ji,numeq,2) = 0.
-               ztrid(ji,numeq,3) = 0.
-               zindterm(ji,numeq)= 0.
-               zindtbis(ji,numeq)= 0.
-               zdiagbis(ji,numeq)= 0.
-            ENDDO
+         ztrid   (1:nidx,:,:) = 0._wp
+         zindterm(1:nidx,:)   = 0._wp
+         zindtbis(1:nidx,:)   = 0._wp
+         zdiagbis(1:nidx,:)   = 0._wp
+
+         DO jm = nlay_s + 2, nlay_s + nlay_i 
+            DO ji = 1, nidx
+               jk = jm - nlay_s - 1
+               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
+               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
+               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk)
+               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
+            END DO
          ENDDO
 
-         DO numeq = nlay_s + 2, nlay_s + nlay_i 
-            DO ji = 1, nidx
-               jk                 = numeq - nlay_s - 1
-               ztrid(ji,numeq,1)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk-1)
-               ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) )
-               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk) * zkappa_i(ji,jk)
-               zindterm(ji,numeq) =  ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk)
-            END DO
-         ENDDO
-
-         numeq =  nlay_s + nlay_i + 1
+         jm =  nlay_s + nlay_i + 1
          DO ji = 1, nidx
             !!ice bottom term
-            ztrid(ji,numeq,1)  =  - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
-            ztrid(ji,numeq,2)  =  1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
-            ztrid(ji,numeq,3)  =  0.0
-            zindterm(ji,numeq) =  ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
-               &                  ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
+            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)   
+            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) )
+            ztrid(ji,jm,3)  = 0.0
+            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  &
+               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) ) 
          ENDDO
 
@@ -401,48 +384,47 @@
             IF ( ht_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells !
                !                            !---------------------!
-               !
-               !!snow interior terms (bottom equation has the same form as the others)
-               DO numeq = 3, nlay_s + 1
-                  jk                  =  numeq - 1
-                  ztrid(ji,numeq,1)   =  - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
-                  ztrid(ji,numeq,2)   =  1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
-                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
-                  zindterm(ji,numeq)  =  ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
+               ! snow interior terms (bottom equation has the same form as the others)
+               DO jm = 3, nlay_s + 1
+                  jk = jm - 1
+                  ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1)
+                  ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) )
+                  ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk)
+                  zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk)
                END DO
 
-               !!case of only one layer in the ice (ice equation is altered)
+               ! case of only one layer in the ice (ice equation is altered)
                IF ( nlay_i == 1 ) THEN
-                  ztrid(ji,nlay_s+2,3)    =  0.0
-                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
+                  ztrid(ji,nlay_s+2,3)  = 0.0
+                  zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 
                ENDIF
 
                IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
 
-                  numeqmin(ji)    =  1
-                  numeqmax(ji)    =  nlay_i + nlay_s + 1
-
-                  !!surface equation
+                  jm_min(ji) = 1
+                  jm_max(ji) = nlay_i + nlay_s + 1
+
+                  ! surface equation
                   ztrid(ji,1,1)  = 0.0
                   ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0)
                   ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0)
-                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji)
-
-                  !!first layer of snow equation
-                  ztrid(ji,2,1)  =  - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
-                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zindterm(ji,2) =  ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
+                  zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji)
+
+                  ! first layer of snow equation
+                  ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1)
+                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+                  ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1)
+                  zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1)
 
                ELSE                            !--  case 2 : surface is melting
                   !
-                  numeqmin(ji)    =  2
-                  numeqmax(ji)    =  nlay_i + nlay_s + 1
-
-                  !!first layer of snow equation
-                  ztrid(ji,2,1)  =  0.0
-                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
-                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
+                  jm_min(ji) = 2
+                  jm_max(ji) = nlay_i + nlay_s + 1
+
+                  ! first layer of snow equation
+                  ztrid(ji,2,1)  = 0.0
+                  ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s )
+                  ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1) 
                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  &
-                     &             ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
+                     &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
                ENDIF
             !                               !---------------------!
@@ -452,57 +434,58 @@
                IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting
                   !
-                  numeqmin(ji)      =  nlay_s + 1
-                  numeqmax(ji)      =  nlay_i + nlay_s + 1
-
-                  !!surface equation	
-                  ztrid(ji,numeqmin(ji),1)   =  0.0
-                  ztrid(ji,numeqmin(ji),2)   =  zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
-                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zindterm(ji,numeqmin(ji))  =  zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji)
-
-                  !!first layer of ice equation
-                  ztrid(ji,numeqmin(ji)+1,1) =  - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
-                  ztrid(ji,numeqmin(ji)+1,2) =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
-                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1) * zkappa_i(ji,1)  
-                  zindterm(ji,numeqmin(ji)+1)=  ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
-
-                  !!case of only one layer in the ice (surface & ice equations are altered)
-
+                  jm_min(ji) = nlay_s + 1
+                  jm_max(ji) = nlay_i + nlay_s + 1
+
+                  ! surface equation	
+                  ztrid(ji,jm_min(ji),1)  = 0.0
+                  ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1    
+                  ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1
+                  zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji)
+
+                  ! first layer of ice equation
+                  ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1)
+                  ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )
+                  ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)  
+                  zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)  
+
+                  ! case of only one layer in the ice (surface & ice equations are altered)
                   IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,numeqmin(ji),1)    =  0.0
-                     ztrid(ji,numeqmin(ji),2)    =  zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
-                     ztrid(ji,numeqmin(ji),3)    =  zkappa_i(ji,0) * 2.0
-                     ztrid(ji,numeqmin(ji)+1,1)  =  -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
-                     ztrid(ji,numeqmin(ji)+1,2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0
-
-                     zindterm(ji,numeqmin(ji)+1) =  ztiold(ji,1) + zeta_i(ji,1) *  &
-                        &                          ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
+                     ztrid(ji,jm_min(ji),1)    = 0.0
+                     ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0
+                     ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0
+                     ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1)
+                     ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
+                     ztrid(ji,jm_min(ji)+1,3)  = 0.0
+                     zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  &
+                        &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )
                   ENDIF
 
                ELSE                            !--  case 2 : surface is melting
 
-                  numeqmin(ji)    =  nlay_s + 2
-                  numeqmax(ji)    =  nlay_i + nlay_s + 1
-
-                  !!first layer of ice equation
-                  ztrid(ji,numeqmin(ji),1)      =  0.0
-                  ztrid(ji,numeqmin(ji),2)      =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
-                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zindterm(ji,numeqmin(ji))     =  ztiold(ji,1) + zeta_i(ji,1) *  &
-                     &                             ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
-
-                  !!case of only one layer in the ice (surface & ice equations are altered)
+                  jm_min(ji)    =  nlay_s + 2
+                  jm_max(ji)    =  nlay_i + nlay_s + 1
+
+                  ! first layer of ice equation
+                  ztrid(ji,jm_min(ji),1)  = 0.0
+                  ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )  
+                  ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1)
+                  zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  &
+                     &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
+
+                  ! case of only one layer in the ice (surface & ice equations are altered)
                   IF ( nlay_i == 1 ) THEN
-                     ztrid(ji,numeqmin(ji),1)  =  0.0
-                     ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
-                     ztrid(ji,numeqmin(ji),3)  =  0.0
-                     zindterm(ji,numeqmin(ji)) =  ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
-                        &                       + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
+                     ztrid(ji,jm_min(ji),1)  = 0.0
+                     ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
+                     ztrid(ji,jm_min(ji),3)  = 0.0
+                     zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  &
+                        &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0
                   ENDIF
 
                ENDIF
             ENDIF
-
+            !
+            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji))
+            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2)
+            !
          END DO
          !
@@ -511,21 +494,17 @@
          !------------------------------
          ! Solve the tridiagonal system with Gauss elimination method.
-         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984.	
-
-         maxnumeqmax = 0
-         minnumeqmin = nlay_i+5
-
-         DO ji = 1, nidx
-            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
-            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
-            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
-            maxnumeqmax               =  MAX(numeqmax(ji),maxnumeqmax)
-         END DO
-
-         DO jk = minnumeqmin+1, maxnumeqmax
-            DO ji = 1, nidx
-               numeq               =  min(max(numeqmin(ji)+1,jk),numeqmax(ji))
-               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2)  - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3)  / zdiagbis(ji,numeq-1)
-               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1)
+         ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984
+         jm_maxt = 0
+         jm_mint = nlay_i+5
+         DO ji = 1, nidx
+            jm_mint = MIN(jm_min(ji),jm_mint)
+            jm_maxt = MAX(jm_max(ji),jm_maxt)
+         END DO
+
+         DO jk = jm_mint+1, jm_maxt
+            DO ji = 1, nidx
+               jm = min(max(jm_min(ji)+1,jk),jm_max(ji))
+               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1)
+               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1)
             END DO
          END DO
@@ -533,11 +512,11 @@
          DO ji = 1, nidx
             ! ice temperatures
-            t_i_1d(ji,nlay_i) =  zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji))
-         END DO
-
-         DO numeq = nlay_i + nlay_s, nlay_s + 2, -1
-            DO ji = 1, nidx
-               jk    =  numeq - nlay_s - 1
-               t_i_1d(ji,jk)  =  ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq)
+            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji))
+         END DO
+
+         DO jm = nlay_i + nlay_s, nlay_s + 2, -1
+            DO ji = 1, nidx
+               jk    =  jm - nlay_s - 1
+               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm)
             END DO
          END DO
@@ -546,12 +525,12 @@
             ! snow temperatures      
             IF( ht_s_1d(ji) > 0._wp ) THEN
-               t_s_1d(ji,nlay_s) =  ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
-                  &                 / zdiagbis(ji,nlay_s+1)
+               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  &
+                  &                / zdiagbis(ji,nlay_s+1)
             ENDIF
             ! surface temperature
             ztsub(ji) = t_su_1d(ji)
             IF( t_su_1d(ji) < rt0 ) THEN
-               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) *  &
-                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))
+               t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  &
+                  &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))
             ENDIF
          END DO
@@ -599,8 +578,8 @@
       DO ji = 1, nidx
          !                                ! surface ice conduction flux
-         fc_su(ji)       =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
-            &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
+         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   &
+            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji))
          !                                ! bottom ice conduction flux
-         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
+         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) )
       END DO
 
@@ -637,11 +616,11 @@
       DO ji = 1, nidx
          !--- Conduction fluxes (positive downwards)
-         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
-         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
+         diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
+         diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
 
          !--- Snow-ice interfacial temperature (diagnostic SIMIP)
          zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji)
          IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN
-            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) + &
+            t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   &
                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac
          ELSE
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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8562)
@@ -29,8 +29,6 @@
 !!gm  It should be probably better to pass these variable in argument of the routine, 
 !!gm  rather than having this long list in USE. This will also highlight what is updated, and what is just use.
-   USE sbc_ice , ONLY : emp_oce, qns_oce, qsr_oce , qemp_oce ,                             &
-      &                 emp_ice, qsr_ice, qemp_ice, qevap_ice, alb_ice, tn_ice, cldf_ice,  &
-      &                 snwice_mass, snwice_mass_b, snwice_fmass
-   USE sbc_oce , ONLY : nn_fsbc, ln_ice_embd, sfx, fr_i, qsr_tot, qns, qsr, fmmflx, emp, taum, utau, vtau
+   USE sbc_ice
+   USE sbc_oce
 !!gm end
    USE sbccpl         ! Surface boundary condition: coupled interface
@@ -107,4 +105,5 @@
       REAL(wp) ::   zqmass           ! Heat flux associated with mass exchange ice->ocean (W.m-2)
       REAL(wp) ::   zqsr             ! New solar flux received by the ocean
+      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d                  ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_cs, zalb_os     ! 3D workspace
       !!---------------------------------------------------------------------
@@ -208,5 +207,78 @@
       ENDIF
       !
+      ! output all fluxes
+      !------------------
+      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface
+      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )      ! non-solar flux at ocean surface
+      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
+      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
+      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
+      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )
+      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
+         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
+!!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ??? 
+      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
+      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
+      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce (:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
+      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice (:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
+
+      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation [m/day]
+
+      CALL iom_put( "sfxbog"      , sfx_bog             )        ! salt flux from bottom growth
+      CALL iom_put( "sfxbom"      , sfx_bom             )        ! salt flux from bottom melting
+      CALL iom_put( "sfxsum"      , sfx_sum             )        ! salt flux from surface melting
+      CALL iom_put( "sfxlam"      , sfx_lam             )        ! salt flux from lateral melting
+      CALL iom_put( "sfxsni"      , sfx_sni             )        ! salt flux from snow ice formation
+      CALL iom_put( "sfxopw"      , sfx_opw             )        ! salt flux from open water formation
+      CALL iom_put( "sfxdyn"      , sfx_dyn             )        ! salt flux from ridging rafting
+      CALL iom_put( "sfxres"      , sfx_res             )        ! salt flux from corrections (resultant)
+      CALL iom_put( "sfxbri"      , sfx_bri             )        ! salt flux from brines
+      CALL iom_put( "sfxsub"      , sfx_sub             )        ! salt flux from sublimation
+      CALL iom_put( "sfx"         , sfx                 )        ! total salt flux
+
+      CALL iom_put( "vfxres"     , wfx_res              )        ! prod./melting due to corrections 
+      CALL iom_put( "vfxopw"     , wfx_opw              )        ! lateral thermodynamic ice production
+      CALL iom_put( "vfxsni"     , wfx_sni              )        ! snowice ice production
+      CALL iom_put( "vfxbog"     , wfx_bog              )        ! bottom thermodynamic ice production
+      CALL iom_put( "vfxdyn"     , wfx_dyn              )        ! dynamic ice production (rid/raft)
+      CALL iom_put( "vfxsum"     , wfx_sum              )        ! surface melt 
+      CALL iom_put( "vfxbom"     , wfx_bom              )        ! bottom melt 
+      CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt 
+      CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt 
+      IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
+
+      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
+         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
+         ELSEWHERE                                       ; z2d = 0._wp
+         END WHERE
+         CALL iom_put( "vfxthin", ( wfx_opw + z2d )     )
+      ENDIF
+
+      CALL iom_put( "vfxspr"     , wfx_spr              )        ! precip (snow)
+      CALL iom_put( "vfxsnw"     , wfx_snw              )        ! total snw growth/melt 
+      CALL iom_put( "vfxsub"     , wfx_sub              )        ! sublimation (snow/ice) 
+      CALL iom_put( "vfxsub_err" , wfx_err_sub          )        ! "excess" of sublimation sent to ocean      
+ 
+      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !  
+      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !  
+      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   !  
+      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   !  
+      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   !  
+      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !  
+      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !  
+      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   !  
+      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !  
+      
+      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   !  
+      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   !  
+      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   !  
+      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !  
+      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !  
+      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base 
+      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
+      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
+      !
       ! 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
@@ -282,5 +354,5 @@
       !                                      !==  every ocean time-step  ==!
       !
-      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
+      DO jj = 2, jpjm1                                !* update the stress WITHOUT an ice-ocean rotation angle
          DO ji = fs_2, fs_jpim1   ! Vect. Opt.
             zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
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 8559)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8562)
@@ -40,5 +40,5 @@
 CONTAINS
 
-   SUBROUTINE ice_wri( kindic )
+   SUBROUTINE ice_wri( kt )
       !!-------------------------------------------------------------------
       !!  This routine computes the average of some variables and write it
@@ -50,8 +50,8 @@
       !!  modif : 03/06/98
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kindic   ! if kindic < 0 there has been an error somewhere
+      INTEGER, INTENT(in) ::   kt   ! time-step
       !
       INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
-      REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val
+      REAL(wp) ::  z2da, z2db, zrho1, zrho2, zmiss_val
       REAL(wp), DIMENSION(jpi,jpj)     ::  z2d, zswi, zmiss !  2D workspace
       REAL(wp), DIMENSION(jpi,jpj)     ::  zfb              ! ice freeboard
@@ -95,19 +95,4 @@
       ! Standard outputs
       !----------------------------------------
-      ! fluxes 
-      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface
-      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )      ! non-solar flux at ocean surface
-      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
-      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
-      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
-      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )
-      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
-         &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
-!!gm I don't understand the variable below.... why not multiplied by a_i_b or (1-a_i_b) ??? 
-      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
-      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
-      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce (:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
-      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice (:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
-
       ! velocity
       IF( iom_use('uice_ipa') )  CALL iom_put( "uice_ipa" , u_ice         )   ! ice velocity u component
@@ -126,11 +111,9 @@
          IF( iom_use('icevel_mv') )   CALL iom_put( "icevel_mv"    , z2d(:,:) * zswi(:,:) + zmiss(:,:) )   ! ice velocity module (missing value)
       ENDIF
-
-      IF( iom_use('tau_icebfr') )     CALL iom_put( "tau_icebfr"  , tau_icebfr             )  ! ice friction with ocean bottom (landfast ice)  
       !
       IF( iom_use('miceage')  )       CALL iom_put( "miceage"     , om_i * zswi * zamask15 )  ! mean ice age
       IF( iom_use('micet')    )       CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature
       IF( iom_use('icest')    )       CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature
-      IF( iom_use('icecolf')  )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness
+      IF( iom_use('icecolf')  )       CALL iom_put( "icecolf"     , ht_i_new               )  ! new ice thickness formed in the leads
       !
       CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature
@@ -142,74 +125,11 @@
       CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content
       CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume
-      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation [m/day]
       CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
       CALL iom_put( "snowvol"     , vt_s    * zswi      )        ! snow volume
       
-      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth
-      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting
-      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting
-      CALL iom_put( "sfxlam"      , sfx_lam * rday      )        ! salt flux from lateral melting
-      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation
-      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation
-      CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting
-      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from corrections (resultant)
-      CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines
-      CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation
-      CALL iom_put( "sfx"         , sfx     * rday      )        ! total salt flux
-
-      ztmp = rday / rhoic
-      CALL iom_put( "vfxres"     , wfx_res * ztmp       )        ! daily prod./melting due to corrections 
-      CALL iom_put( "vfxopw"     , wfx_opw * ztmp       )        ! daily lateral thermodynamic ice production
-      CALL iom_put( "vfxsni"     , wfx_sni * ztmp       )        ! daily snowice ice production
-      CALL iom_put( "vfxbog"     , wfx_bog * ztmp       )        ! daily bottom thermodynamic ice production
-      CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp       )        ! daily dynamic ice production (rid/raft)
-      CALL iom_put( "vfxsum"     , wfx_sum * ztmp       )        ! surface melt 
-      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt 
-      CALL iom_put( "vfxlam"     , wfx_lam * ztmp       )        ! lateral melt 
-      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt 
-
-      IF ( ln_pnd ) &
-         CALL iom_put( "vfxpnd"  , wfx_pnd * ztmp       )        ! melt pond water flux
-
-      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
-         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog
-         ELSEWHERE                                       ; z2d = 0._wp
-         END WHERE
-         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp )
-      ENDIF
-
-      ztmp = rday * r1_rhosn
-      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
-      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
-      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice) 
-      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean      
- 
-      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !  
-      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !  
-      CALL iom_put ('hfxres'     , hfx_res(:,:)         )   !  
-      CALL iom_put ('hfxout'     , hfx_out(:,:)         )   !  
-      CALL iom_put ('hfxin'      , hfx_in(:,:)          )   !  
-      CALL iom_put ('hfxsnw'     , hfx_snw(:,:)         )   !  
-      CALL iom_put ('hfxsub'     , hfx_sub(:,:)         )   !  
-      CALL iom_put ('hfxerr'     , hfx_err_dif(:,:)     )   !  
-      CALL iom_put ('hfxerr_rem' , hfx_err_rem(:,:)     )   !  
-      
-      CALL iom_put ('hfxsum'     , hfx_sum(:,:)         )   !  
-      CALL iom_put ('hfxbom'     , hfx_bom(:,:)         )   !  
-      CALL iom_put ('hfxbog'     , hfx_bog(:,:)         )   !  
-      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !  
-      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !  
-      CALL iom_put ('hfxtur'     , fhtur(:,:) * at_i_b(:,:) ) ! turbulent heat flux at ice base 
-      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
-      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
-
-!!gm ====>>>>>  THIS should be moved where at_ip, vt_ip are computed fro the last time in the time-step  (limmpd)
-      ! MV MP 2016
       IF ( ln_pnd ) THEN
          CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction
          CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area
       ENDIF
-      ! END MV MP 2016
-!!gm  <<<<<<======= end
 
       !----------------------------------
@@ -225,5 +145,4 @@
       IF ( iom_use('brinevol_cat') )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )          ! brine volume
 
-      ! MV MP 2016
       IF ( ln_pnd ) THEN
          IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories
@@ -232,5 +151,4 @@
          IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories
       ENDIF
-      ! END MV MP 2016
 
       !--------------------------------
@@ -360,5 +278,5 @@
       !! History :   4.0  !  2013-06  (C. Rousset)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index)
+      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index
       INTEGER, INTENT( in )   ::   kid , kh_i
       INTEGER                 ::   nz_i, jl
@@ -476,8 +394,6 @@
       CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )    
 
-      ! Close the file
-      ! -----------------
-!!gm I don't understand why the file is not closed !
-      !CALL histclo( kid )
+      !! The file is closed in dia_wri_state (ocean routine)
+      !! CALL histclo( kid )
       !
     END SUBROUTINE ice_wri_state
