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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8150)
@@ -256,4 +256,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow/ice sublimation       [kg.m-2.s-1]
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: dynamical component of wfx_snw    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1]
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange                   [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice      [kg.m-2.s-1]
@@ -417,8 +420,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_si         !: mean temperature at the snow-ice interface (K) 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmi_dyn  !: Change in ice mass due to ice dynamics (kg/m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dms_mel  !: Change in snow mass due to melting (kg/m2/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dms_dyn  !: Change in snow mass due to ice dynamics (kg/m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmtx_dyn !: X-component of total snow + ice mass transport (kg/m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmty_dyn !: X-component of total snow + ice mass transport (kg/m2/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmtx_dyn !: X-component of total snow + ice mass transport (kg/s)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmty_dyn !: X-component of total snow + ice mass transport (kg/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2)
@@ -465,7 +467,8 @@
 
       ii = ii + 1
-      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                        &
+      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                                           &
          &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        &
-         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                        &
+         &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj), wfx_snw_sum(jpi,jpj)                   &
+         &      wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                                           &
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
@@ -534,6 +537,5 @@
       ii = ii + 1
       ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
-                diag_dmi_dyn(jpi,jpj)                         ,    &
-                diag_dms_mel(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    &
+                diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    &
                 diag_dmtx_dyn(jpi,jpj), diag_dmty_dyn(jpi,jpj),    &
                 diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
Index: /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 8150)
@@ -190,5 +190,6 @@
          ! water flux
          zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
-            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
+            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) +                                & 
+            &                  wfx_sub(:,:) + wfx_spr(:,:)                                                 &
             &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
 
@@ -215,5 +216,6 @@
          ! water flux
          zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
-            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
+            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) +                             +  &
+            &                wfx_sub(:,:) + wfx_spr(:,:)                                                 &
             &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfw_b
 
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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8150)
@@ -535,4 +535,5 @@
       INTEGER ::   icells            ! number of cells with a_i > puny
       REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
+      REAL(wp) ::   zwfx_snw         ! snow mass flux increment
 
       INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
@@ -648,6 +649,6 @@
             srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
             
-            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
-            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
+            sfx_dyn(ji,jj)     = sfx_dyn(ji,jj)     - smsw(ij) * rhoic * r1_rdtice
+            wfx_dyn(ji,jj)     = wfx_dyn(ji,jj)     - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
             
              ! virtual salt flux to keep salinity constant
@@ -666,11 +667,13 @@
             !  During the next time step, thermo_rates will determine whether
             !  the ocean cools or new ice grows.
-            wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
-               &                              + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
+            zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
+               &                 + 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) = diag_dms_dyn(ji,jj) - ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
-               &                                      - rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice
-            diag_dmi_dyn(ji,jj) = diag_dmi_dyn(ji,jj) + vsw (ij) * rhoic * r1_rdtice 
+            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)
 
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8150)
@@ -149,5 +149,5 @@
                                                                              !   ocean surface (ssh_m) if ice is not embedded
                                                                              !   ice top surface if ice is embedded   
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zCor                            ! Coriolis stress array (SIMIP)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zCorx, zCory                    ! Coriolis stress array (SIMIP)
       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays
       REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence
@@ -164,5 +164,5 @@
       CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice )
       CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf )
-      CALL wrk_alloc( jpi,jpj, zCor)
+      CALL wrk_alloc( jpi,jpj, zCorx, zCory)
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -449,10 +449,10 @@
 
                   ! Coriolis at V-points (energy conserving formulation)
-                  zCor(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor(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) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -465,9 +465,4 @@
             CALL lbc_lnk( v_ice, 'V', -1. )
 
-            ! SIMIP diag
-            IF ( jter .EQ. nn_nevp ) THEN
-               diag_corstry(:,:) = zCor(:,:) 
-            ENDIF
-            
 #if defined key_agrif && defined key_lim2
             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
@@ -485,10 +480,10 @@
 
                   ! Coriolis at U-points (energy conserving formulation)
-                  zCor(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                  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) + zCor(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) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -500,7 +495,4 @@
             END DO
             CALL lbc_lnk( u_ice, 'U', -1. )
-            IF ( jter .EQ. nn_nevp ) THEN
-               diag_corstrx(:,:) = zCor(:,:) 
-            ENDIF
             
 #if defined key_agrif && defined key_lim2
@@ -521,10 +513,10 @@
 
                   ! Coriolis at U-points (energy conserving formulation)
-                  zCor(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                  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) + zCor(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) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
 
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -552,10 +544,10 @@
 
                   ! Coriolis at V-points (energy conserving formulation)
-                  zCor(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
                      &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
                      &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
 
                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor(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) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
                   
                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
@@ -640,26 +632,39 @@
       stress12_i(:,:) = zs12(:,:)
 
-      ! SIMIP diagnostic internal stress terms (N/m2)
-      diag_dssh_dx(:,:) = zspgU(:,:)
-      diag_dssh_dy(:,:) = zspgV(:,:)
+      ! SIMIP diagnostics: sea surface sloop stress, coriolis and internal stress terms (N/m2) 
+      !                    stress tensor invariants (normal and shear stress N/m)
+      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
+
+             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
+
+             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * zswi    ! sea surface slope stress term
+             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * zswi
+
+             diag_corstrx(ji,jj) = zCorx(ji,jj) * zswi    ! Coriolis stress term
+             diag_corstry(ji,jj) = zCory(ji,jj) * zswi
+
+             diag_intstrx(ji,jj) = zfU(ji,jj)   * zswi    ! internal stress term
+             diag_intstry(ji,jj) = zfV(ji,jj)   * zswi
+
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi( diag_sig1   , 'T',  1., diag_sig2   , 'T',  1., &
+                 &        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. )
 
-      diag_intstrx(:,:) = zfU(:,:)
-      diag_intstry(:,:) = zfV(:,:)
       CALL lbc_lnk( diag_intstrx, 'U', -1. )
       CALL lbc_lnk( diag_intstry, 'V', -1. )
 
-      ! SIMIP diagnostic stress tensor invariants (normal and shear stress N/m)
-      DO jj = k_j1+1, k_jpj-1
-         DO ji = 2, jpim1
-             zswi  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 1.0e-6 ) ) ! 1 if ice, 0 if no ice
-             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * zswi
-             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * zswi
-         END DO
-      END DO
-      CALL lbc_lnk( diag_sig1, 'T', 1. )
-      CALL lbc_lnk( diag_sig2, 'T', 1. )
-
+      CALL lbc_lnk( diag_corstrx, 'U', -1. )
+      CALL lbc_lnk( diag_corstry, 'V', -1. )
+            
       !
       !------------------------------------------------------------------------------!
@@ -705,5 +710,5 @@
       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, zCor )
+      CALL wrk_dealloc( jpi,jpj, zCorx, zCory)
 
    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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8150)
@@ -169,6 +169,6 @@
 
             ! mass flux at the ocean/ice interface
-            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model
-            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
+            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )             ! F/M mass flux save at least for biogeochemical model
+            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)  ! mass flux + F/M mass flux (always ice/ocean mass exchange)
          END DO
       END DO
@@ -278,6 +278,6 @@
       !
       ! SIMIP diags
-      diag_utau_oi(:,:) = 0._wp
-      diag_vtau_oi(:,:) = 0._wp
+      !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,10 +290,11 @@
             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
+            
+            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( 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/limthd.F90
===================================================================
--- /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8150)
@@ -502,5 +502,4 @@
          
          CALL tab_2d_1d( nbpb, t_su_1d     (1:nbpb), t_su(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, t_si_1d     (1:nbpb), t_si(:,:,jl)    , jpi, jpj, npb(1:nbpb) ) ! MV SIMIP 2016
          CALL tab_2d_1d( nbpb, sm_i_1d     (1:nbpb), sm_i(:,:,jl)    , jpi, jpj, npb(1:nbpb) )
          DO jk = 1, nlay_s
@@ -530,4 +529,5 @@
          
          CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_snw_sum_1d(1:nbpb), wfx_snw_sum  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
          
@@ -562,5 +562,4 @@
 
          ! SIMIP diagnostics
-         CALL tab_2d_1d( nbpb, diag_dms_mel_1d (1:nbpb), diag_dms_mel, jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, diag_fc_bo_1d   (1:nbpb), diag_fc_bo  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, diag_fc_su_1d   (1:nbpb), diag_fc_su  , jpi, jpj, npb(1:nbpb) )
@@ -573,5 +572,4 @@
          CALL tab_1d_2d( nbpb, a_i (:,:,jl)  , npb, a_i_1d     (1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, t_su(:,:,jl)  , npb, t_su_1d    (1:nbpb)   , jpi, jpj )
-         CALL tab_1d_2d( nbpb, t_si(:,:,jl)  , npb, t_si_1d    (1:nbpb)   , jpi, jpj ) ! MV 2016
          CALL tab_1d_2d( nbpb, sm_i(:,:,jl)  , npb, sm_i_1d    (1:nbpb)   , jpi, jpj )
          DO jk = 1, nlay_s
@@ -587,4 +585,5 @@
          
          CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_snw_sum   , npb, wfx_snw_sum_1d(1:nbpb),jpi, jpj )
          CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
          
@@ -622,5 +621,5 @@
 
          ! SIMIP diagnostics         
-         CALL tab_1d_2d( nbpb, diag_dms_mel   , npb, diag_dms_mel_1d(1:nbpb) , jpi, jpj )
+         CALL tab_1d_2d( nbpb, t_si(:,:,jl)   , npb, t_si_1d    (1:nbpb)     , jpi, jpj )
          CALL tab_1d_2d( nbpb, diag_fc_bo     , npb, diag_fc_bo_1d(1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, diag_fc_su     , npb, diag_fc_su_1d(1:nbpb)   , jpi, jpj )
Index: /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8150)
@@ -88,5 +88,5 @@
       REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg)
       REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean
-      REAL(wp) ::   zsstK        ! SST in Kelvin
+      REAL(wp) ::   zsstK        ! SST (K)
 
       REAL(wp), POINTER, DIMENSION(:) ::   zqprec      ! energy of fallen snow                       (J.m-3)
@@ -174,7 +174,5 @@
             hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! Contribution to mass flux
-            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
-            ! SIMIP snow melt diagnostic
-            diag_dms_mel_1d(ji) = diag_dms_mel_1d(ji) - rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
+            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! updates
             ht_s_1d(ji)   = 0._wp
@@ -240,7 +238,5 @@
          hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
          ! snow melting only = water into the ocean (then without snow precip), >0
-         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice    
-         ! SIMIP snow melt diagnostic
-         diag_dms_mel_1d(ji) = diag_dms_mel_1d(ji) + rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
+         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice    
          ! updates available heat + precipitations after melting
          zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
@@ -264,7 +260,5 @@
             hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
-            ! SIMIP snow melt diagnostic
-            diag_dms_mel_1d(ji) = diag_dms_mel_1d(ji) + rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
+            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
             ! updates available heat + thickness
             zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
@@ -618,7 +612,5 @@
          hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
          ! Contribution to mass flux
-         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
-         ! SIMIP snow melt diagnostic
-         diag_dms_mel_1d(ji) = diag_dms_mel_1d(ji) + rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
+         wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
          !    
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
@@ -628,5 +620,10 @@
          IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
       END DO
-      
+
+      ! add snow melt to mass snow mass flux to the ocean
+      DO ji = kideb, kiut
+         wfx_snw_1d(ji) = wfx_snw_1d(ji) + wfx_snw_sum_1d(ji)
+      END DO
+
       !
       !------------------------------------------------------------------------------|
@@ -708,5 +705,5 @@
       ! --- ensure that a_i = 0 where ht_i = 0 ---
       WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
-      
+
       CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
Index: /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8150)
@@ -743,5 +743,5 @@
       END DO
 
-      WHERE( zh_s .LT. epsi10) ; t_si_1d(:) = t_su_1d(:) ; END WHERE
+      WHERE( ( zh_s .LT. 1.0e-3 ) ) ; t_si_1d(:) = t_su_1d(:) ; END WHERE
       ! END MV SIMIP 2016
 
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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 8150)
@@ -102,4 +102,7 @@
 
       zsm(:,:) = e12t(:,:)
+
+      ! SIMIP diags
+      diag_dmtx_dyn(:,:) = 0._wp ; diag_dmty_dyn(:,:) = 0._wp
       
       !                             !-------------------------------------!
@@ -190,27 +193,18 @@
                DO jl = 1, jpl
                   ! SIMIP mass transport diags
-                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
+                  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) ) * r1_e12t(:,:)
-
-                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
+                  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) ) * r1_e12t(:,:)
+                  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, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-!                    &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-!                 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
-!                    &                                       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)  )
-!                 CALL lim_adv_y( zusnit, v_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
-!                    &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
                   CALL lim_adv_x( zusnit, u_ice, 1._wp, zsm, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
                      &                                       sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
@@ -246,26 +240,19 @@
                   &                                       sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
                DO jl = 1, jpl
-                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
+                  ! SIMIP mass transport diags
+                  diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) )
                   CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
                      &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
                   CALL lim_adv_y( zusnit, v_ice, 1._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) ) * r1_e12t(:,:)
-
-                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) - ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
-                  CALL lim_adv_y( zusnit, v_ice, 1._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) )
+
+                  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),   &
+                     &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
                   CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
                      &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-                  diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) + ( rhoic * z0ice(:,:,jl) + rhosn * z0snw(:,:,jl) ) * r1_e12t(:,:)
-
-!                 CALL lim_adv_y( zusnit, v_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, 0._wp, zsm, z0ice (:,:,jl), sxice(:,:,jl),   &
-!                    &                                       sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-!                 CALL lim_adv_y( zusnit, v_ice, 1._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-!                    &                                       sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-!                 CALL lim_adv_x( zusnit, u_ice, 0._wp, zsm, z0snw (:,:,jl), sxsn (:,:,jl),   &
-!                    &                                       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 ---
@@ -298,6 +285,6 @@
 
          ! SIMIP diags
-         diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) / ( rdt_ice * zusnit )
-         diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) / ( rdt_ice * zusnit )
+         diag_dmtx_dyn(:,:) = diag_dmtx_dyn(:,:) / rdt_ice 
+         diag_dmty_dyn(:,:) = diag_dmty_dyn(:,:) / rdt_ice
 
          !-------------------------------------------
@@ -401,6 +388,6 @@
 
                ! SIMIP diagnostics
-               diag_dms_dyn(ji,jj) = rhosn * diag_trp_vs(ji,jj)
-               diag_dmi_dyn(ji,jj) = rhoic * diag_trp_vi(ji,jj)
+               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
@@ -447,7 +434,4 @@
                         hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
                         hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
-
-                        ! SIMIP diagnostic
-                        diag_dms_dyn(ji,jj) = diag_dms_dyn(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
 
                      ENDIF
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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8150)
@@ -600,7 +600,4 @@
                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
 
-               ! SIMIP snow melt diagnostic
-               diag_dms_mel_1d(ji) = diag_dms_mel(ji,jj) + ( v_s(ji,jj,jl) - zvs  ) * rhosn * r1_rdtice
-
             END DO
          END DO
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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8150)
@@ -80,5 +80,5 @@
       INTEGER  :: Nu, Nv                                                       ! passage size
 
-      INTEGER, PARAMETER :: i_grid = 2                                         ! grid type (eORCA1 = 1, ORCA2 = 2)
+      INTEGER, DIMENSION(4)  :: ji0, ji1, jj0, jj1
       
       !!-------------------------------------------------------------------
@@ -142,6 +142,4 @@
       !
       IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age
-      IF ( iom_use( "icethic_cea" ) )   CALL iom_put( "icethic_cea" , htm_i * zswi           )  ! ice thickness mean
-      IF ( iom_use( "snowthic_cea" ) )  CALL iom_put( "snowthic_cea", htm_s * zswi           )  ! snow thickness mean
       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
@@ -228,337 +226,4 @@
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
 
-      !--------------------------------
-      ! Add-ons for SIMIP
-      !--------------------------------
-      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"    , smt_i * vt_i * 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 ) ! mass of salt in sea ice per cell area
-
-      IF  ( iom_use( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum                  )          ! Freshwater flux from sea-ice surface
-      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 )
-      IF  ( iom_use( "dmidyn"   ) ) CALL iom_put( "dmidyn"      ,   diag_dmi_dyn             )          ! Sea-ice mass change from dynamics
-      IF  ( iom_use( "dmiopw"   ) ) CALL iom_put( "dmiopw"      , - wfx_opw                  )          ! Sea-ice mass change through growth in open water
-      IF  ( iom_use( "dmibog"   ) ) CALL iom_put( "dmibog"      , - wfx_bog                  )          ! Sea-ice mass change through basal growth
-      IF  ( iom_use( "dmisni"   ) ) CALL iom_put( "dmisni"      , - wfx_sni                  )          ! Sea-ice mass change through snow-to-ice conversion
-      IF  ( iom_use( "dmisum"   ) ) CALL iom_put( "dmisum"      , - wfx_sum                  )          ! Sea-ice mass change through surface melting
-      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"      ,   diag_dms_mel             )          ! 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( "dmtxdyn"  ) ) CALL iom_put( "dmtxdyn"     ,   diag_dmtx_dyn            )          ! X-component of sea-ice mass transport
-      IF  ( iom_use( "dmtydyn"  ) ) CALL iom_put( "dmtydyn"     ,   diag_dmty_dyn            )          ! Y-component of sea-ice mass transport
-
-      IF  ( iom_use( "utau_oce" ) ) CALL iom_put( "utau_oce"    ,   diag_utau_oi*zswi        )          ! X-component of ocean stress on sea ice
-      IF  ( iom_use( "vtau_oce" ) ) CALL iom_put( "vtau_oce"    ,   diag_vtau_oi*zswi        )          ! Y-component of ocean stress on sea ice
-
-      IF  ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx*zswi        )          ! Sea-surface tilt term in force balance (x-component)
-      IF  ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy*zswi        )          ! Sea-surface tilt term in force balance (y-component)
-
-      IF  ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx*zswi        )          ! Coriolis force term in force balance (x-component)
-      IF  ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry*zswi        )          ! Coriolis force term in force balance (y-component)
-
-      IF  ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx*zswi        )          ! Internal force term in force balance (x-component)
-      IF  ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry*zswi        )          ! Internal force term in force balance (y-component)
-
-      IF  ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1   *zswi        )          ! Normal stress
-      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2   *zswi        )          ! Shear stress
-
-      !--------------------------------
-      ! Global ice diagnostics (SIMIP)
-      !--------------------------------
-      !
-      ! What follows could be a separate routine
-      !
-
-       IF ( iom_use( "NH_icearea" ) .OR. iom_use( "NH_icevolu" ) ) THEN    ! NH ice area 
- 
-         WHERE( fcor > 0 ); zswi(:,:) = 1.0e-12; 
-         ELSEWHERE        ; zswi(:,:) = 0.
-         END WHERE 
-
-           IF ( iom_use( "NH_icearea" ) ) THEN
-              zdiag_area_nh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) )
-              CALL iom_put( "NH_icearea", zdiag_area_nh )
-           ENDIF
-           IF ( iom_use( "NH_icevolu" ) ) THEN
-              zdiag_volu_nh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
-              CALL iom_put( "NH_icevolu", zdiag_volu_nh )
-           ENDIF
-
-      ENDIF 
-
-      IF ( iom_use( "NH_iceextt" ) ) THEN                                   ! NH ice extt 
-
-         WHERE( fcor > 0 .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12; 
-         ELSEWHERE                          ; zswi(:,:) = 0.
-         END WHERE 
-
-         zdiag_extt_nh = glob_sum( zswi(:,:) * e12t(:,:) )
-         CALL iom_put( "NH_iceextt", zdiag_extt_nh )
-
-      ENDIF 
-
-      IF ( iom_use( "SH_icearea" ) .OR. iom_use( "SH_icevolu" ) ) THEN    ! SH ice area  / volume
-
-         WHERE( fcor < 0 ); zswi(:,:) = 1.0e-12; 
-         ELSEWHERE        ; zswi(:,:) = 0.
-         END WHERE 
-
-         IF ( iom_use( "SH_icearea" ) ) THEN
-            zdiag_area_sh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) ) 
-            CALL iom_put( "SH_icearea", zdiag_area_sh )
-         ENDIF
-         IF ( iom_use( "SH_icevolu" ) ) THEN
-            zdiag_volu_sh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
-            CALL iom_put( "SH_icevolu", zdiag_volu_sh )
-         ENDIF
-
-      ENDIF 
-
-      IF ( iom_use( "SH_iceextt" ) ) THEN                                   ! SH ice extt 
-
-         WHERE( fcor < 0 .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12; 
-         ELSEWHERE                          ; zswi(:,:) = 0.
-         END WHERE 
-
-         zdiag_extt_sh = glob_sum( zswi(:,:) * e12t(:,:) )
-
-         CALL iom_put( "SH_iceextt", zdiag_extt_sh )
-
-      ENDIF 
-
-      !------------------------
-      ! Fluxes through straits
-      !------------------------
-      !
-      ! This piece of code is quite awful and should probably be a separate routine
-      !
-      ! See Notz et al 2016 for definitions
-      ! 4 Arctic passages are considered (Fram, CAA, Barents, Bering)
-      !
-
-      IF ( iom_use("strait_arfl") .OR. iom_use("strait_mifl") .OR. iom_use("strait_msfl") ) THEN
-
-         zdiag_area_strait(:) = 0._wp
-         zdiag_mice_strait(:) = 0._wp
-         zdiag_msno_strait(:) = 0._wp
-   
-         ! === Fram Strait ===
-         IF ( i_grid == 2 ) THEN     ! ORCA2
-            Nv  = 4
-            zvi(1:Nv) = (/ 133, 134, 135, 136 /)
-            zvj(1:Nv) = (/ 136, 136, 136, 136 /)
-         ENDIF
-   
-         IF ( i_grid == 1 ) THEN    ! eORCA1
-            Nv  = 10
-            zvi(1:Nv) = (/ 268,269,270,271,272,273,274,275,276,277 /)
-            zvj(1:Nv) = (/ 311,311,311,311,311,311,311,311,311,311 /)
-         ENDIF
-   
-         DO ii = 1, Nv
-            ji = zvi(ii)    
-            jj = zvj(ii)
-            zdiag_area_strait(1) = zdiag_area_strait(1)                               & ! --- ice area flux ---
-                &                + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         & ! northwards (positive) flow
-                &                + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 )           ! southwards (negative) flow
-   
-            zdiag_mice_strait(1) = zdiag_mice_strait(1) + rhoic *                     & ! --- ice mass flux  ---
-                &                ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         & !
-                &                + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 ) )         !
-   
-            zdiag_msno_strait(1) = zdiag_msno_strait(1) + rhosn *                     & ! --- snow mass flux ---
-                &                ( vt_s(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         & 
-                &                + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 ) )  
-         END DO
-   
-         ! === Bering Strait ===
-         IF ( i_grid == 1 ) THEN ! eORCA1
-            Nv  = 3
-            zvi(1:Nv) = (/ 113,114,115 /)
-            zvj(1:Nv) = (/ 285,285,285 /)
-         ENDIF
-   
-         IF ( i_grid == 2 ) THEN ! ORCA2
-            Nv  = 2
-            zvi(1:Nv) = (/ 55 , 56 /)
-            zvj(1:Nv) = (/ 122,122 /)
-         ENDIF
-   
-         DO ii = 1, Nv
-            ji = zvi(ii)     
-            jj = zvj(ii)
-            zdiag_area_strait(4) = zdiag_area_strait(4)                               & ! --- ice area flux ---
-                &                + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         & ! northwards (positive) flow
-                &                + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 )           ! southwards (negative) flow
-   
-            zdiag_mice_strait(4) = zdiag_mice_strait(4) + rhoic *                     & ! --- ice mass flux  ---
-                &                ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )            & !
-                &                + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 ) )            !
-   
-            zdiag_msno_strait(4) = zdiag_msno_strait(4) + rhosn *                     & ! --- snow mass flux ---
-                &                ( vt_s(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )            & 
-                &                + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 ) )  
-         END DO
-   
-         ! === Barents throughflow (eORCA1)
-
-         ! U-flow
-         IF ( i_grid == 1 ) THEN ! 'eORCA1'
-            Nu = 11
-            zui(1:Nu) = (/ 282,283,284,285,286,286,287,288,289,290,292/)
-            zuj(1:Nu) = (/ 308,307,306,305,304,303,302,301,300,299,298/)
-         ENDIF
-         IF ( i_grid == 2 ) THEN ! 'ORCA2'
-            Nu = 5
-            zui(1:Nu) = (/ 141,142,142,143,144 /)
-            zuj(1:Nu) = (/ 134,133,132,131,130 /)
-         ENDIF
-         zfarea_u = 0._wp
-         zfmice_u = 0._wp
-         zfmsno_u = 0._wp
-   
-         DO ii = 1, Nu
-            ji = zui(ii)     
-            jj = zuj(ii)
-            zfarea_u          = zfarea_u                                           & ! --- ice area zonal flux ---
-                &             + at_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         & ! --- eastward
-                &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 )           ! --- westward
-            zfmice_u          = zfmice_u + rhoic *                                 & ! --- ice mass zonal flux --- 
-                &             ( vt_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
-                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
-            zfmsno_u          = zfmsno_u + rhosn *                                 & ! --- snow mass zonal flux --- 
-                &             ( vt_s(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
-                &             + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
-         END DO
-   
-         ! V-flow
-         IF ( i_grid == 1 ) THEN ! 'eORCA1'
-            Nv = 9
-            zvi(1:Nv) = (/ 282,283,284,285,286,287,288,289,290/)
-            zvj(1:Nv) = (/ 308,307,306,305,303,302,301,300,299/)
-         ENDIF
-         IF ( i_grid == 2 ) THEN ! 'ORCA2'
-            Nv = 4
-            zvi(1:Nv) = (/ 140,141,142,143 /)
-            zvj(1:Nv) = (/ 135,134,132,131 /)
-         ENDIF
-         zfarea_v = 0._wp
-         zfmice_v = 0._wp
-         zfmsno_v = 0._wp
-   
-         DO ii  = 1, Nv  
-            ji = zvi(ii)     
-            jj = zvj(ii)
-            zfarea_v          = zfarea_v                                           & ! --- ice area meridian flux ---
-                &             + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & ! --- eastward
-                &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 )           ! --- westward
-            zfmice_v          = zfmice_v + rhoic *                                 & ! --- ice mass meridian flux ---
-                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
-                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
-            zfmsno_v          = zfmsno_v + rhosn *                                 & ! --- snow mass meridian flux ---
-                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
-                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
-         END DO
-   
-         ! Sum U/V contributions
-         zdiag_area_strait(3) = zfarea_u + zfarea_v 
-         zdiag_mice_strait(3) = zfmice_u + zfmice_v
-         zdiag_msno_strait(3) = zfmsno_u + zfmsno_v
-   
-         ! === CAA throughflow ===
-         ! U-flow through Queen Elisabeth Islands and McClure straits
-         IF ( i_grid == 1 ) THEN  ! eORCA1
-            Nu = 8 
-            zui(1:Nu) = (/ 231,231,231,  132,132,132,132,132  /)
-            zuj(1:Nu) = (/ 328,329,330,  318,319,320,321,322  /)
-            zfarea_u = 0._wp
-            zfmice_u = 0._wp
-            zfmsno_u = 0._wp
-      
-            DO ii = 1, Nu
-               ji = zui(ii)     
-               jj = zuj(ii)
-               zfarea_u          = zfarea_u                                           & ! --- ice area zonal flux ---
-                   &             + at_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         & ! --- eastward
-                   &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 )           ! --- westward
-               zfmice_u          = zfmice_u + rhoic *                                 & ! --- ice mass zonal flux --- 
-                   &             ( vt_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
-                   &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
-               zfmsno_u          = zfmsno_u + rhosn *                                 & ! --- snow mass zonal flux --- 
-                   &             ( vt_s(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
-                   &             + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
-            END DO
-         ENDIF
-
-         IF ( i_grid == 2 ) THEN   ! ORCA2
-            zfarea_u = 0._wp       ! QEI and McClure straits are not resolved in ORCA2
-            zfmice_u = 0._wp
-            zfmsno_u = 0._wp
-         ENDIF
-   
-         ! V-flow through Nares Strait
-         IF ( i_grid == 1 ) THEN   ! eORCA1
-            Nv = 4
-            zvi(1:Nv) = (/ 254,255,256,257 /)
-            zvj(1:Nv) = (/ 317,317,317,317 /)
-         ENDIF
-         IF ( i_grid == 2 ) THEN   ! ORCA2
-            Nv = 2
-            zvi(1:Nv) = (/ 117,118 /)
-            zvj(1:Nv) = (/ 145,145 /)
-         ENDIF
-         zfarea_v = 0._wp
-         zfmice_v = 0._wp
-         zfmsno_v = 0._wp
-   
-         DO ii = 1, Nv  
-            ji = zvi(ii)     
-            jj = zvj(ii)
-            zfarea_v          = zfarea_v                                           & ! --- ice area meridian flux ---
-                &             + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & ! --- eastward
-                &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 )           ! --- westward
-            zfmice_v          = zfmice_v + rhoic *                                 & ! --- ice mass meridian flux ---
-                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
-                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
-            zfmsno_v          = zfmsno_v + rhosn *                                 & ! --- snow mass meridian flux ---
-                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
-                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
-         END DO
-
-         ! Sum U/V contributions
-         zdiag_area_strait(2) = zfarea_u + zfarea_v 
-         zdiag_mice_strait(2) = zfmice_u + zfmice_v
-         zdiag_msno_strait(2) = zfmsno_u + zfmsno_v
-   
-         ! === Write in file
-         IF ( iom_use("strait_arfl") ) CALL iom_put( "strait_arfl", zdiag_area_strait )
-         IF ( iom_use("strait_mifl") ) CALL iom_put( "strait_mifl", zdiag_mice_strait )
-         IF ( iom_use("strait_msfl") ) CALL iom_put( "strait_msfl", zdiag_msno_strait ) 
-
-         WRITE(numout,*) " area flx ", zdiag_area_strait(:)
-         WRITE(numout,*) " mice flx ", zdiag_mice_strait(:)
-         WRITE(numout,*) " msno flx ", zdiag_msno_strait(:)
-
-      ENDIF
-
       !----------------------------------
       ! Output category-dependent fields
@@ -576,4 +241,326 @@
       ! brine volume
       IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )
+
+      !--------------------------------
+      ! Add-ons for SIMIP
+      !--------------------------------
+      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( "dmithd"   ) ) CALL iom_put( "dmithd"      , - wfx_bog - wfx_bom - wfx_sum   &     ! Sea-ice mass change from thermodynamics
+              &                     - wfx_sni - wfx_opw - wfx_res )
+      IF  ( iom_use( "dmidyn"   ) ) CALL iom_put( "dmidyn"      ,   diag_dmi_dyn             )          ! Sea-ice mass change from dynamics
+      IF  ( iom_use( "dmiopw"   ) ) CALL iom_put( "dmiopw"      , - wfx_opw                  )          ! Sea-ice mass change through growth in open water
+      IF  ( iom_use( "dmibog"   ) ) CALL iom_put( "dmibog"      , - wfx_bog                  )          ! Sea-ice mass change through basal growth
+      IF  ( iom_use( "dmisni"   ) ) CALL iom_put( "dmisni"      , - wfx_sni                  )          ! Sea-ice mass change through snow-to-ice conversion
+      IF  ( iom_use( "dmisum"   ) ) CALL iom_put( "dmisum"      , - wfx_sum                  )          ! Sea-ice mass change through surface melting
+      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( "dmtxdyn"  ) ) CALL iom_put( "dmtxdyn"     ,   diag_dmtx_dyn            )          ! X-component of sea-ice mass transport
+      IF  ( iom_use( "dmtydyn"  ) ) CALL iom_put( "dmtydyn"     ,   diag_dmty_dyn            )          ! Y-component of sea-ice mass transport
+
+      IF  ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   diag_utau_oi*zswi        )          ! X-component of ocean stress on sea ice
+      IF  ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   diag_vtau_oi*zswi        )          ! Y-component of ocean stress on sea ice
+
+      IF  ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx*zswi        )          ! Sea-surface tilt term in force balance (x-component)
+      IF  ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy*zswi        )          ! Sea-surface tilt term in force balance (y-component)
+
+      IF  ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx*zswi        )          ! Coriolis force term in force balance (x-component)
+      IF  ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry*zswi        )          ! Coriolis force term in force balance (y-component)
+
+      IF  ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx*zswi        )          ! Internal force term in force balance (x-component)
+      IF  ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry*zswi        )          ! Internal force term in force balance (y-component)
+
+      IF  ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1   *zswi        )          ! Normal stress
+      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2   *zswi        )          ! Shear stress
+
+      !--------------------------------
+      ! Global ice diagnostics (SIMIP)
+      !--------------------------------
+
+      IF ( iom_use( "NH_icearea" ) .OR. iom_use( "NH_icevolu" ) .OR. iom_use( "NH_iceextt" ) )   THEN   ! NH integrated diagnostics
+ 
+         WHERE( fcor > 0._wp ); zswi(:,:) = 1.0e-12
+         ELSEWHERE            ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_area_nh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) )
+         zdiag_volu_nh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
+
+         WHERE( fcor > 0._wp .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12
+         ELSEWHERE                              ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_extt_nh = glob_sum( zswi(:,:) * e12t(:,:) )
+
+         IF ( iom_use( "NH_icearea" ) ) CALL iom_put( "NH_icearea" ,  zdiag_area_nh  )
+         IF ( iom_use( "NH_icevolu" ) ) CALL iom_put( "NH_icevolu" ,  zdiag_volu_nh  )
+         IF ( iom_use( "NH_iceextt" ) ) CALL iom_put( "NH_iceextt" ,  zdiag_extt_nh  )
+
+      ENDIF
+
+      IF ( iom_use( "SH_icearea" ) .OR. iom_use( "SH_icevolu" ) .OR. iom_use( "SH_iceextt" ) )   THEN   ! SH integrated diagnostics
+
+         WHERE( fcor < 0._wp ); zswi(:,:) = 1.0e-12; 
+         ELSEWHERE            ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_area_sh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) ) 
+         zdiag_volu_sh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
+
+         WHERE( fcor < 0._wp .AND. at_i > 0.15 ); zswi(:,:) = 1.0e-12
+         ELSEWHERE                              ; zswi(:,:) = 0.
+         END WHERE 
+
+         zdiag_extt_sh = glob_sum( zswi(:,:) * e12t(:,:) )
+
+         IF ( iom_use( "SH_icearea" ) ) CALL iom_put( "SH_icearea", zdiag_area_sh )
+         IF ( iom_use( "SH_icevolu" ) ) CALL iom_put( "SH_icevolu", zdiag_volu_sh )
+         IF ( iom_use( "SH_iceextt" ) ) CALL iom_put( "SH_iceextt", zdiag_extt_sh )
+
+      ENDIF 
+
+      !--------------------------------
+      ! Fluxes through straits (SIMIP)
+      !--------------------------------
+      !
+      ! Valid only for ORCA-like grids
+      !
+      ! 4 Arctic passages are considered (Fram, CAA, Barents, Bering; see Notz et al (GMD 2016) for definitions)
+      !
+      ! Fram and Bering  straits are easy because they follow parallels
+      ! Barents and Canadian Arctic Archipelago are less easy because they do not, which is why they look so awful.
+      ! 
+
+      IF ( iom_use( "strait_arfl" ) .OR. iom_use( "strait_mifl" ) .OR. iom_use( "strait_msfl" ) .AND. cp_cfg == "orca" ) THEN
+
+         zdiag_area_strait(:) = 0._wp   ;   zdiag_mice_strait(:) = 0._wp   ;   zdiag_msno_strait(:) = 0._wp
+   
+         !------------------------------
+         ! === Fram & Bering Straits ===
+         !------------------------------
+          
+         SELECT CASE ( jp_cfg ) 
+          
+         CASE ( 2 )   ! --- ORCA2
+          
+            ! Fram Strait   (i_strait = 1)
+            ji0(1) = 133   ;   ji1(1) = 136
+            jj0(1) = 136 
+
+            ! Bering Strait (i_strait = 4)
+            ji0(4) = 55    ;   ji1(4) = 56
+            jj0(4) = 122
+          
+         CASE ( 1 )   ! --- eORCA1
+          
+            ! Fram Strait
+            ji0(1) = 268   ;   ji1(1) = 277
+            jj0(1) = 311
+
+            ! Bering Strait
+            ji0(4) = 113   ;   jj1(4) = 115
+            jj0(4) = 285
+           
+         END SELECT
+          
+         DO i_strait = 1, 4, 3
+
+            DO ji = mi0(ji0), mi1(ji1)
+               jj = mj0(jj0)
+   
+               zdiag_area_strait(i_strait) = zdiag_area_strait(i_strait)                                 &     ! --- ice area flux ---
+                   &                + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         &     ! northwards (positive) flow
+                   &                + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 )               ! southwards (negative) flow
+      
+               zdiag_mice_strait(i_strait) = zdiag_mice_strait(i_strait) + rhoic *                       &     ! --- ice mass flux  ---
+                   &                ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         & 
+                   &                + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 ) )         
+      
+               zdiag_msno_strait(i_strait) = zdiag_msno_strait(i_strait) + rhosn *                       &     ! --- snow mass flux ---
+                   &                ( vt_s(ji,jj-1) * e12t(ji,jj-1) * MAX( v_ice(ji,jj-1), 0.0 )         & 
+                   &                + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( v_ice(ji,jj-1), 0.0 ) )  
+   
+            END DO
+
+         END DO
+
+         !---------------------
+         ! === Barents opening
+         !---------------------
+   
+         SELECT CASE ( jp_cfg ) 
+
+            CASE ( 1 )   ! 'eORCA1'
+
+               Nu = 11   ! U-Flow
+               zui(1:Nu) = (/ 282,283,284,285,286,286,287,288,289,290,292/)
+               zuj(1:Nu) = (/ 308,307,306,305,304,303,302,301,300,299,298/)
+
+               Nv = 9    ! V-Flow
+               zvi(1:Nv) = (/ 282,283,284,285,286,287,288,289,290/)
+               zvj(1:Nv) = (/ 308,307,306,305,303,302,301,300,299/)
+
+            CASE ( 2 )   ! 'ORCA2'
+
+               Nu = 5    ! U-Flow
+               zui(1:Nu) = (/ 141,142,142,143,144 /)
+               zuj(1:Nu) = (/ 134,133,132,131,130 /)
+
+               Nv = 4    ! V-Flow
+               zvi(1:Nv) = (/ 140,141,142,143 /)
+               zvj(1:Nv) = (/ 135,134,132,131 /)
+
+         END SELECT
+
+         ! Barents U-flow
+         zfarea_u = 0._wp   ;   zfmice_u = 0._wp   ;   zfmsno_u = 0._wp
+   
+         DO ii = 1, Nu
+
+            ji = mi0(zui(ii))
+            jj = mj0(zuj(ii))
+
+            zfarea_u          = zfarea_u                                           & ! --- ice area zonal flux ---
+                &             + at_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         & ! --- northward
+                &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 )           ! --- southward
+            zfmice_u          = zfmice_u + rhoic *                                 & ! --- ice mass zonal flux --- 
+                &             ( vt_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
+                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
+            zfmsno_u          = zfmsno_u + rhosn *                                 & ! --- snow mass zonal flux --- 
+                &             ( vt_s(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
+                &             + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
+         END DO
+   
+         ! Barents V-flow
+         zfarea_v = 0._wp   ;   zfmice_v = 0._wp   ;   zfmsno_v = 0._wp
+
+         DO ii  = 1, Nv  
+
+            ji = mi0(zvi(ii))
+            jj = mj0(zvj(ii))
+
+            zfarea_v          = zfarea_v                                           & ! --- ice area meridian flux ---
+                &             + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & ! --- eastward
+                &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 )           ! --- westward
+            zfmice_v          = zfmice_v + rhoic *                                 & ! --- ice mass meridian flux ---
+                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
+                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
+            zfmsno_v          = zfmsno_v + rhosn *                                 & ! --- snow mass meridian flux ---
+                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
+                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
+         END DO
+   
+         ! Sum Barents U-/V- contributions
+         zdiag_area_strait(3) = zfarea_u + zfarea_v 
+         zdiag_mice_strait(3) = zfmice_u + zfmice_v
+         zdiag_msno_strait(3) = zfmsno_u + zfmsno_v
+
+         !---------------------
+         ! === CAA throughflow
+         !---------------------
+   
+         SELECT CASE ( jp_cfg ) 
+
+            CASE ( 1 )   ! eORCA1
+
+               ! V-flow through Nares Strait
+               Nv = 4
+               zvi(1:Nv) = (/ 254,255,256,257 /)
+               zvj(1:Nv) = (/ 317,317,317,317 /)
+
+               ! U-flow through Queen Elisabeth Islands and McClure straits
+               Nu = 8 
+               zui(1:Nu) = (/ 231,231,231,  132,132,132,132,132  /)
+               zuj(1:Nu) = (/ 328,329,330,  318,319,320,321,322  /)
+
+               zfarea_u = 0._wp   ;   zfmice_u = 0._wp   ;   zfmsno_u = 0._wp
+
+               DO ii = 1, Nu
+
+                  ji = mi0(zui(ii))
+                  jj = mj0(zuj(ii))
+
+                  zfarea_u          = zfarea_u                                                           & ! --- ice area zonal flux ---
+                      &             + at_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         & ! --- eastward
+                      &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 )           ! --- westward
+                  zfmice_u          = zfmice_u + rhoic *                                                 & ! --- ice mass zonal flux --- 
+                      &             ( vt_i(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
+                      &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
+                  zfmsno_u          = zfmsno_u + rhosn *                                                 & ! --- snow mass zonal flux --- 
+                      &             ( vt_s(ji-1,jj) * e12t(ji-1,jj) * MAX( u_ice(ji-1,jj), 0.0 )         &   
+                      &             + vt_s(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji-1,jj), 0.0 ) )          
+
+               END DO
+
+
+            CASE ( 2 )   ! ORCA2
+
+               ! V-flow through Nares Strait
+               Nv = 2
+               zvi(1:Nv) = (/ 117,118 /)
+               zvj(1:Nv) = (/ 145,145 /)
+
+               ! U-flow through Queen Elisabeth Islands and McClure straits (not resolved in ORCA2)
+               zfarea_u = 0._wp   ;   zfmice_u = 0._wp   ;   zfmsno_u = 0._wp
+
+            END SELECT
+   
+         ! V-flow through Nares Strait
+         zfarea_v = 0._wp   ;   zfmice_v = 0._wp   ;   zfmsno_v = 0._wp
+   
+         DO ii = 1, Nv  
+
+            ji = mi0(zvi(ii))
+            jj = mj0(zvj(ii))
+
+            zfarea_v          = zfarea_v                                           & ! --- ice area meridian flux ---
+                &             + at_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & ! --- eastward
+                &             + at_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 )           ! --- westward
+            zfmice_v          = zfmice_v + rhoic *                                 & ! --- ice mass meridian flux ---
+                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
+                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
+            zfmsno_v          = zfmsno_v + rhosn *                                 & ! --- snow mass meridian flux ---
+                &             ( vt_i(ji,jj-1) * e12t(ji,jj-1) * MAX( u_ice(ji,jj-1), 0.0 )         & !
+                &             + vt_i(ji,jj  ) * e12t(ji,jj)   * MIN( u_ice(ji,jj-1), 0.0 ) )         !
+
+         END DO
+
+         ! Sum U/V contributions
+         zdiag_area_strait(2) = zfarea_u + zfarea_v 
+         zdiag_mice_strait(2) = zfmice_u + zfmice_v
+         zdiag_msno_strait(2) = zfmsno_u + zfmsno_v
+   
+         ! === Ncdf output
+         IF ( iom_use("strait_arfl") ) CALL iom_put( "strait_arfl", zdiag_area_strait )
+         IF ( iom_use("strait_mifl") ) CALL iom_put( "strait_mifl", zdiag_mice_strait )
+         IF ( iom_use("strait_msfl") ) CALL iom_put( "strait_msfl", zdiag_msno_strait ) 
+
+         WRITE(numout,*) " area flx ", zdiag_area_strait(:)
+         WRITE(numout,*) " mice flx ", zdiag_mice_strait(:)
+         WRITE(numout,*) " msno flx ", zdiag_msno_strait(:)
+
+      ENDIF
 
       !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s
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 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8150)
@@ -66,4 +66,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sum_1d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d
 
@@ -130,5 +131,4 @@
 
    ! Extra SIMIP data request diagnostics
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   diag_dms_mel_1d    !: <==> the 2D  diag_dms_mel
    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
@@ -177,9 +177,7 @@
          &      t_s_1d(jpij,nlay_s) , t_i_1d(jpij,nlay_i) , s_i_1d(jpij,nlay_i) ,  &            
          &      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))
+         &      qh_i_old(jpij,0:nlay_i+1), h_i_old(jpij,0:nlay_i+1) , STAT=ierr(3) )
       ! Extra SIMIP diagnostics
-      ALLOCATE( diag_dms_mel_1d(jpij)                               ,                & 
-         &      diag_fc_bo_1d(jpij)      , diag_fc_su_1d(jpij)      ,                &
-         &      STAT=ierr(4) )
+      ALLOCATE( diag_fc_bo_1d(jpij)      , diag_fc_su_1d(jpij)      , STAT=ierr(4) )
 
       thd_ice_alloc = MAXVAL( ierr )
Index: /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 8150)
@@ -230,4 +230,6 @@
       ! automatic definitions of some of the xml attributs
       CALL set_xmlatt
+
+      CALL set_1point
 
       ! end file definition
@@ -1573,4 +1575,20 @@
 
 
+   SUBROUTINE set_1point 
+   !!---------------------------------------------------------------------- 
+   !!                     ***  ROUTINE set_1point  *** 
+   !! 
+   !! ** Purpose :   define zoom grid for scalar fields 
+   !! 
+   !!---------------------------------------------------------------------- 
+   REAL(wp), DIMENSION(1)   ::   zz = 1. 
+   INTEGER  :: ix, iy 
+   !!---------------------------------------------------------------------- 
+   CALL dom_ngb( 180., 90., ix, iy, 'T' ) !  Nearest point to north pole should be ocean 
+   CALL iom_set_domain_attr('1point', zoom_ibegin=ix, zoom_jbegin=iy) 
+
+   END SUBROUTINE set_1point 
+
+
    SUBROUTINE set_xmlatt
       !!----------------------------------------------------------------------
Index: /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8149)
+++ /branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8150)
@@ -595,4 +595,6 @@
       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
       wfx_spr(:,:) = 0._wp   ;   
+
+      wfx_snw_dyn(:,:) = 0._wp ; wfx_snw_sum(:,:) = 0._wp
       
       hfx_thd(:,:) = 0._wp   ;   
@@ -613,6 +615,6 @@
 
       ! SIMIP diagnostics
-      diag_dms_mel(:,:)  = 0._wp ; diag_dms_dyn(:,:)  = 0._wp
-      diag_dmtx_dyn(:,:) = 0._wp ; diag_dmty_dyn(:,:) = 0._wp
+      diag_dms_sum(:,:)  = 0._wp ; diag_dms_dyn(:,:)  = 0._wp
+      diag_dmi_dyn(:,:)  = 0._wp ;
       diag_fc_bo(:,:)    = 0._wp ; diag_fc_su(:,:)    = 0._wp
       
