Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8152)
@@ -534,5 +534,5 @@
          &      diag_smvi  (jpi,jpj), diag_vice   (jpi,jpj), diag_vsnw  (jpi,jpj), STAT=ierr(ii) )
 
-      ! * SIMIP data request diagnostics
+      ! * SIMIP diagnostics
       ii = ii + 1
       ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8152)
@@ -670,10 +670,6 @@
                &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
 
-            wfx_snw_dyn(ji,jj) = wfx_snw(ji,jj) + zwfx_snw
-            wfx_snw(ji,jj)     = wfx_snw(ji,jj) + zwfx_snw
-
-            ! SIMIP diagnostic
-            diag_dms_dyn(ji,jj) = - wfx_dyn(ji,jj)     + diag_trp_vi(ji,jj)
-            diag_dmi_dyn(ji,jj) = - wfx_snw_dyn(ji,jj) + diag_trp_vs(ji,jj)
+            wfx_snw_dyn(ji,jj)  =   wfx_snw_dyn(ji,jj) + zwfx_snw
+            wfx_snw(ji,jj)      =   wfx_snw(ji,jj)     + zwfx_snw
 
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
@@ -774,4 +770,8 @@
       END DO ! jl1 (deforming categories)
 
+      ! SIMIP diagnostics
+      diag_dmi_dyn(:,:) = - wfx_dyn(:,:)     + rhoic * diag_trp_vi(:,:)
+      diag_dms_dyn(:,:) = - wfx_snw_dyn(:,:) + rhosn * diag_trp_vs(:,:)
+      
       !
       CALL wrk_dealloc( jpij,        indxi, indxj )
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8152)
@@ -149,5 +149,7 @@
                                                                              !   ocean surface (ssh_m) if ice is not embedded
                                                                              !   ice top surface if ice is embedded   
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zCorx, zCory                    ! Coriolis stress array (SIMIP)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zCorx, zCory                    ! Coriolis stress array
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
+
       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
       REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
@@ -165,4 +167,5 @@
       CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
       CALL wrk_alloc( jpi,jpj, zCorx, zCory)
+      CALL wrk_alloc( jpi,jpj, ztaux_oi, ztauy_oi)
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -448,4 +451,7 @@
                      &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
                   ! Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
@@ -454,5 +460,5 @@
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -479,11 +485,14 @@
                      &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
                   ! Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
-                  
+
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -512,11 +521,14 @@
                      &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  
                   ! Coriolis at U-points (energy conserving formulation)
                   zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
                      &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
                      &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
-                  
+
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -543,4 +555,7 @@
                      &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
 
+                  ! Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
                   ! Coriolis at V-points (energy conserving formulation)
                   zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
@@ -549,5 +564,5 @@
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -632,13 +647,17 @@
       stress12_i(:,:) = zs12(:,:)
 
-      ! SIMIP diagnostics: sea surface sloop stress, coriolis and internal stress terms (N/m2) 
-      !                    stress tensor invariants (normal and shear stress N/m)
+      !------------------------------------------------------------------------------!
+      ! 5) SIMIP diagnostics
+      !------------------------------------------------------------------------------!
+                           
       DO jj = k_j1+1, k_jpj-1
          DO ji = 2, jpim1
              zswi  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
 
+             ! Stress tensor invariants (normal and shear stress N/m)
              diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * zswi                                 ! normal stress
              diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * zswi   ! shear stress
 
+             ! Stress terms of the momentum equation (N/m2)
              diag_dssh_dx(ji,jj) = zspgU(ji,jj) * zswi    ! sea surface slope stress term
              diag_dssh_dy(ji,jj) = zspgV(ji,jj) * zswi
@@ -649,4 +668,7 @@
              diag_intstrx(ji,jj) = zfU(ji,jj)   * zswi    ! internal stress term
              diag_intstry(ji,jj) = zfV(ji,jj)   * zswi
+           
+             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * zswi  ! oceanic stress
+             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * zswi
 
          END DO
@@ -656,18 +678,10 @@
                  &        diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1., &
                  &        diag_corstrx, 'U', -1., diag_corstry, 'V', -1., & 
-                 &        diag_intstrx, 'U', -1., diag_intstry, 'V', -1.  )
-
-!     CALL lbc_lnk( diag_dssh_dx, 'U', -1. )
-!     CALL lbc_lnk( diag_dssh_dy, 'V', -1. )
-
-!     CALL lbc_lnk( diag_intstrx, 'U', -1. )
-!     CALL lbc_lnk( diag_intstry, 'V', -1. )
-
-!     CALL lbc_lnk( diag_corstrx, 'U', -1. )
-!     CALL lbc_lnk( diag_corstry, 'V', -1. )
-            
-      !
-      !------------------------------------------------------------------------------!
-      ! 5) Control prints of residual and charge ellipse
+                 &        diag_intstrx, 'U', -1., diag_intstry, 'V', -1.  &
+                 &        diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
+
+      !
+      !------------------------------------------------------------------------------!
+      ! 6) Control prints of residual and charge ellipse
       !------------------------------------------------------------------------------!
       !
@@ -710,5 +724,6 @@
       CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
-      CALL wrk_dealloc( jpi,jpj, zCorx, zCory)
+      CALL wrk_dealloc( jpi,jpj, zCorx, zCory )
+      CALL wrk_dealloc( jpi,jpj, ztaux_oi, ztauy_oi )
 
    END SUBROUTINE lim_rhg
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8152)
@@ -277,7 +277,4 @@
       !                                      !==  every ocean time-step  ==!
       !
-      ! SIMIP diags
-      !diag_utau_oi(:,:) = 0._wp
-      !diag_vtau_oi(:,:) = 0._wp
       DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
          DO ji = fs_2, fs_jpim1   ! Vect. Opt.
@@ -290,11 +287,7 @@
             utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
             vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
-            
-            diag_utau_oi(ji,jj) = - zutau_ice                   ! diagnostic oceanic stress on sea ice
-            diag_vtau_oi(ji,jj) = - zvtau_ice
          END DO
       END DO
       CALL lbc_lnk( utau, 'U', -1. )           ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
-      CALL lbc_lnk( diag_utau_oi, 'U', -1. )   ;   CALL lbc_lnk( diag_vtau_oi, 'V', -1. )   ! lateral boundary condition
       !
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 8152)
@@ -103,7 +103,4 @@
       zsm(:,:) = e12t(:,:)
 
-      ! SIMIP diags
-      diag_dmtx_dyn(:,:) = 0._wp ; diag_dmty_dyn(:,:) = 0._wp
-      
       !                             !-------------------------------------!
       IF( ln_limdyn ) THEN          !   Advection of sea ice properties   !
@@ -119,4 +116,7 @@
          zeiold(:,:)    = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
          zesold(:,:)    = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
+
+         ! SIMIP diags init
+         diag_dmtx_dyn(:,:) = 0._wp ; diag_dmty_dyn(:,:) = 0._wp
 
          !--- Thickness correction init. -------------------------------
@@ -192,19 +192,24 @@
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
+
                   ! SIMIP mass transport diags
                   diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
+
                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
+
                   diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
 
                   diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
+
                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
                   CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
+
                   diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
-                  ! END SIMIP mass transport diags
+
                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
@@ -240,4 +245,5 @@
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
+
                   ! SIMIP mass transport diags
                   diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
@@ -249,10 +255,9 @@
 
                   diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
+                  CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
                   diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
-                  ! END SIMIP mass transport diags
 
                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
@@ -386,8 +391,4 @@
                diag_trp_vs (ji,jj) = SUM(   v_s(ji,jj,:) -  zvsold(ji,jj,:) ) * r1_rdtice
                diag_trp_smv(ji,jj) = SUM( smv_i(ji,jj,:) - zsmvold(ji,jj,:) ) * r1_rdtice
-
-               ! SIMIP diagnostics
-               diag_dms_dyn(ji,jj) = diag_dms_dyn(ji,jj) + rhosn * diag_trp_vs(ji,jj)
-               diag_dmi_dyn(ji,jj) = diag_dmi_dyn(ji,jj) + rhoic * diag_trp_vi(ji,jj)
             END DO
          END DO
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8152)
@@ -116,5 +116,5 @@
          tm_i(:,:)  = 0._wp
          tm_su(:,:) = 0._wp
-         tm_si(:,:) = 0._wp ! MV SIMIP 2016
+         tm_si(:,:) = 0._wp
          om_i (:,:) = 0._wp
          DO jl = 1, jpl
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8152)
@@ -159,7 +159,7 @@
       CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
 
-      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
-      CALL iom_put( "idive"       , divu_i              )    ! divergence
-      CALL iom_put( "ishear"      , shear_i             )    ! shear
+      CALL iom_put( "icestr"      , strength * zswi )            ! ice strength
+      CALL iom_put( "idive"       , divu_i              )        ! divergence
+      CALL iom_put( "ishear"      , shear_i             )        ! shear
       CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume
       
@@ -248,15 +248,14 @@
       zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn / rau0
 
-      IF  ( iom_use( "icethic"  ) ) CALL iom_put( "icethic"     , htm_i * zswi               )          ! ice thickness 
-      IF  ( iom_use( "icepres"  ) ) CALL iom_put( "icepres"     , zswi                       )          ! ice presence (1 or 0) 
-      IF  ( iom_use( "snowthic" ) ) CALL iom_put( "snowthic"    , htm_s * zswi               )          ! snow thickness       
-      IF  ( iom_use( "icemass"  ) ) CALL iom_put( "icemass"     , rhoic * vt_i(:,:) * zswi   )          ! ice mass per cell area 
-      IF  ( iom_use( "snomass"  ) ) CALL iom_put( "snomass"     , rhosn * vt_s(:,:) * zswi   )          ! snow mass per cell area
-      IF  ( iom_use( "icesnt"   ) ) CALL iom_put( "icesnt"      , ( tm_si - rt0 ) * zswi     )          ! snow-ice interface temperature
-      IF  ( iom_use( "icebot"   ) ) CALL iom_put( "icebot"      , ( t_bo  - rt0 ) * zswi     )          ! ice bottom temperature
-      IF  ( iom_use( "icesmass" ) ) CALL iom_put( "icesmass"    , SUM( smv_i, DIM = 3 ) * rhoic / 1000. * zswi )         ! mass of salt in sea ice per cell area
-      IF  ( iom_use( "icefb"    ) ) CALL iom_put( "icefb"       , ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) ) * zswi )   ! ice freeboard
-
-      IF  ( iom_use( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum                  )          ! Freshwater flux from sea-ice surface
+      IF  ( iom_use( "icethic"  ) ) CALL iom_put( "icethic"     , htm_i * zswi               )          ! Ice thickness 
+      IF  ( iom_use( "icepres"  ) ) CALL iom_put( "icepres"     , zswi                       )          ! Ice presence (1 or 0) 
+      IF  ( iom_use( "snowthic" ) ) CALL iom_put( "snowthic"    , htm_s * zswi               )          ! Snow thickness       
+      IF  ( iom_use( "icemass"  ) ) CALL iom_put( "icemass"     , rhoic * vt_i(:,:) * zswi   )          ! Ice mass per cell area 
+      IF  ( iom_use( "snomass"  ) ) CALL iom_put( "snomass"     , rhosn * vt_s(:,:) * zswi   )          ! Snow mass per cell area
+      IF  ( iom_use( "icesnt"   ) ) CALL iom_put( "icesnt"      , ( tm_si - rt0 ) * zswi     )          ! Snow-ice interface temperature
+      IF  ( iom_use( "icebot"   ) ) CALL iom_put( "icebot"      , ( t_bo  - rt0 ) * zswi     )          ! Ice bottom temperature
+      IF  ( iom_use( "icesmass" ) ) CALL iom_put( "icesmass"    , SUM( smv_i, DIM = 3 ) * rhoic / 1000. * zswi )         ! Mass of salt in sea ice per cell area
+      IF  ( iom_use( "icefb"    ) ) CALL iom_put( "icefb"       , ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) ) * zswi )   ! Ice freeboard
+
       IF  ( iom_use( "dmithd"   ) ) CALL iom_put( "dmithd"      , - wfx_bog - wfx_bom - wfx_sum   &     ! Sea-ice mass change from thermodynamics
               &                     - wfx_sni - wfx_opw - wfx_res )
@@ -268,14 +267,15 @@
       IF  ( iom_use( "dmibom"   ) ) CALL iom_put( "dmibom"      , - wfx_bom                  )          ! Sea-ice mass change through bottom melting
       IF  ( iom_use( "dmtsub"   ) ) CALL iom_put( "dmtsub"      , - wfx_sub                  )          ! Sea-ice mass change through evaporation and sublimation
-      IF  ( iom_use( "dmsspr"   ) ) CALL iom_put( "dmsspr"      , - wfx_spr                  )          ! snow mass change through snow fall
-      IF  ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn/rhoic      )          ! snow mass change through snow-to-ice conversion
-
-      IF  ( iom_use( "dmsmel"   ) ) CALL iom_put( "dmsmel"      , - wfx_snw_sum              )          ! snow mass change through melt
-      IF  ( iom_use( "dmsdyn"   ) ) CALL iom_put( "dmsdyn"      ,   diag_dms_dyn             )          ! snow mass change through dynamics
-
-      IF  ( iom_use( "hfxconbo" ) ) CALL iom_put( "hfxconbo"    ,   diag_fc_bo               )          ! bottom conduction flux
-      IF  ( iom_use( "hfxconsu" ) ) CALL iom_put( "hfxconsu"    ,   diag_fc_su               )          ! surface conduction flux
-
-      IF  ( iom_use( "wfxtot"   ) ) CALL iom_put( "wfxtot"      ,   -wfx_ice                 )          ! total freshwater flux from sea ice
+      IF  ( iom_use( "dmsspr"   ) ) CALL iom_put( "dmsspr"      , - wfx_spr                  )          ! Snow mass change through snow fall
+      IF  ( iom_use( "dmsssi"   ) ) CALL iom_put( "dmsssi"      ,   wfx_sni*rhosn/rhoic      )          ! Snow mass change through snow-to-ice conversion
+
+      IF  ( iom_use( "dmsmel"   ) ) CALL iom_put( "dmsmel"      , - wfx_snw_sum              )          ! Snow mass change through melt
+      IF  ( iom_use( "dmsdyn"   ) ) CALL iom_put( "dmsdyn"      ,   diag_dms_dyn             )          ! Snow mass change through dynamics
+
+      IF  ( iom_use( "hfxconbo" ) ) CALL iom_put( "hfxconbo"    ,   diag_fc_bo               )          ! Bottom conduction flux
+      IF  ( iom_use( "hfxconsu" ) ) CALL iom_put( "hfxconsu"    ,   diag_fc_su               )          ! Surface conduction flux
+
+      IF  ( iom_use( "wfxtot"   ) ) CALL iom_put( "wfxtot"      ,   wfx_ice                  )          ! Total freshwater flux from sea ice
+      IF  ( iom_use( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum                  )          ! Freshwater flux from sea-ice surface
 
       IF  ( iom_use( "dmtxdyn"  ) ) CALL iom_put( "dmtxdyn"     ,   diag_dmtx_dyn            )          ! X-component of sea-ice mass transport
Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8151)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8152)
@@ -105,5 +105,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_si_1d       !: <==> the 2D  t_si MV SIMIP 2016
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_si_1d       !: <==> the 2D  t_si
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ht_i_1d       !: <==> the 2D  ht_s
@@ -128,9 +128,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old  !: ice thickness layer (m)
 
-   INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point
-
-   ! Extra SIMIP data request diagnostics
+   ! Conduction flux diagnostics (SIMIP)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   diag_fc_bo_1d      !: <==> the 2D  diag_fc_bo
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   diag_fc_su_1d      !: <==> the 2D  diag_fc_su
+
+   INTEGER , PUBLIC ::   jiindex_1d   ! 1D index of debugging point
 
    !!----------------------------------------------------------------------
@@ -178,5 +178,5 @@
          &      q_i_1d(jpij,nlay_i+1) , q_s_1d(jpij,nlay_s) ,                        &
          &      qh_i_old(jpij,0:nlay_i+1), h_i_old(jpij,0:nlay_i+1) , STAT=ierr(3) )
-      ! Extra SIMIP diagnostics
+
       ALLOCATE( diag_fc_bo_1d(jpij)      , diag_fc_su_1d(jpij)      , STAT=ierr(4) )
 
