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 8283)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8285)
@@ -423,23 +423,6 @@
    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_dyn  !: Change in snow mass due to ice dynamics (kg/m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_ice !: X-component of ice mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_ice !: Y-component of ice mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_snw !: X-component of snow mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_snw !: Y-component of snow mass transport (kg/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xatrp    !: X-component of area transport (m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_yatrp    !: Y-component of area transport (m2/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)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_utau_oi  !: X-direction ocean-ice stress
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vtau_oi  !: Y-direction ocean-ice stress  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dx  !: X-direction sea-surface tilt term (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dy  !: X-direction sea-surface tilt term (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstrx  !: X-direction coriolis stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstry  !: Y-direction coriolis stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstrx  !: X-direction internal stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstry  !: Y-direction internal stress (N/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig1     !: Average normal stress in sea ice   
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig2     !: Maximum shear stress in sea ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_shear    !: Maximum shear of sea-ice velocity field
 
    !
@@ -544,13 +527,5 @@
       ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
                 diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    &
-                diag_xmtrp_ice(jpi,jpj), diag_ymtrp_ice(jpi,jpj),  &
-                diag_xmtrp_snw(jpi,jpj), diag_ymtrp_snw(jpi,jpj),  &
-                diag_xatrp(jpi,jpj)    , diag_yatrp(jpi,jpj)    ,  &
                 diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
-                diag_utau_oi(jpi,jpj) , diag_vtau_oi(jpi,jpj) ,    &
-                diag_dssh_dx(jpi,jpj) , diag_dssh_dy(jpi,jpj) ,    &
-                diag_corstrx(jpi,jpj) , diag_corstry(jpi,jpj) ,    &
-                diag_intstrx(jpi,jpj) , diag_intstry(jpi,jpj) ,    &
-                diag_sig1(jpi,jpj)    , diag_sig2(jpi,jpj)    ,    &
                 STAT = ierr(ii) )
 
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 8283)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8285)
@@ -38,4 +38,5 @@
    USE in_out_manager ! I/O manager
    USE prtctl         ! Print control
+   USE iom
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 #if defined key_agrif && defined key_lim2
@@ -129,5 +130,5 @@
       REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
       REAL(wp) ::   zintb, zintn                                             ! dummy argument
-      REAL(wp) ::   zswi, zfac_x, zfac_y
+      REAL(wp) ::   zfac_x, zfac_y
       
       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength
@@ -156,7 +157,28 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice
 
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_xmtrp_ice                 ! X-component of ice mass transport  (kg/s)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_ymtrp_ice                 ! Y-component of ice mass transport  (kg/s)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_xmtrp_snw                 ! X-component of snow mass transport (kg/s)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_ymtrp_snw                 ! Y-component of snow mass transport (kg/s)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_xatrp                     ! X-component of area transport (m2/s)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_yatrp                     ! Y-component of area transport (m2/s)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_utau_oi                   ! X-direction ocean-ice stress
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_vtau_oi                   ! Y-direction ocean-ice stress  
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_dssh_dx                   ! X-direction sea-surface tilt term (N/m2)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_dssh_dy                   ! X-direction sea-surface tilt term (N/m2)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_corstrx                   ! X-direction coriolis stress (N/m2)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_corstry                   ! Y-direction coriolis stress (N/m2)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_intstrx                   ! X-direction internal stress (N/m2)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_intstry                   ! Y-direction internal stress (N/m2)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_sig1                      ! Average normal stress in sea ice   
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiag_sig2                      ! Maximum shear stress in sea ice
+
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zswi, zmiss                     ! Switch & missing value array
+
       REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
       REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
       REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
+      REAL(wp), PARAMETER               ::   zmiss_val = 1.0e+20             ! missing value for outputs
+
       !!-------------------------------------------------------------------
 
@@ -168,4 +190,14 @@
       CALL wrk_alloc( jpi,jpj, zCorx, zCory)
       CALL wrk_alloc( jpi,jpj, ztaux_oi, ztauy_oi)
+
+      CALL wrk_alloc( jpi,jpj, zdiag_xmtrp_ice, zdiag_ymtrp_ice )
+      CALL wrk_alloc( jpi,jpj, zdiag_xmtrp_snw, zdiag_ymtrp_snw )
+      CALL wrk_alloc( jpi,jpj, zdiag_xatrp    , zdiag_yatrp     )
+      CALL wrk_alloc( jpi,jpj, zdiag_utau_oi  , zdiag_vtau_oi   )
+      CALL wrk_alloc( jpi,jpj, zdiag_dssh_dx  , zdiag_dssh_dy   )
+      CALL wrk_alloc( jpi,jpj, zdiag_corstrx  , zdiag_corstry   )
+      CALL wrk_alloc( jpi,jpj, zdiag_intstrx  , zdiag_intstry   )
+      CALL wrk_alloc( jpi,jpj, zdiag_sig1     , zdiag_sig2      )
+      CALL wrk_alloc( jpi,jpj, zswi           , zmiss           )
 
 #if  defined key_lim2 && ! defined key_lim2_vp
@@ -650,48 +682,86 @@
       ! 5) SIMIP diagnostics
       !------------------------------------------------------------------------------!
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+         END DO
+      END DO
+
+      zmiss(:,:)         = zmiss_val * ( 1. - zswi(:,:) )
                            
       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
+             zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * zswi(ji,jj)                                 ! normal stress
+             zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * zswi(ji,jj)   ! 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
-
-             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
+             zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * zswi(ji,jj)     ! sea surface slope stress term
+             zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * zswi(ji,jj)
+
+             zdiag_corstrx(ji,jj) = zCorx(ji,jj) * zswi(ji,jj)     ! Coriolis stress term
+             zdiag_corstry(ji,jj) = zCory(ji,jj) * zswi(ji,jj)
+
+             zdiag_intstrx(ji,jj) = zfU(ji,jj)   * zswi(ji,jj)     ! internal stress term
+             zdiag_intstry(ji,jj) = zfV(ji,jj)   * zswi(ji,jj)
            
-             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * zswi  ! oceanic stress
-             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * zswi
+             zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * zswi(ji,jj)  ! oceanic stress
+             zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * zswi(ji,jj)
 
              ! 2D ice mass, snow mass, area transport arrays (X, Y)
-             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zswi
-             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zswi
-
-             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
-             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
-
-             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
-             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
-
-             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
-             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
+             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zswi(ji,jj)
+             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zswi(ji,jj)
+
+             zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component (kg/s)
+             zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+             zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+             zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+             zdiag_xatrp(ji,jj)     = zfac_x         * ( at_i(ji+1,jj) + at_i(ji,jj) ) ! area transport,      X-component (m2/s)
+             zdiag_yatrp(ji,jj)     = zfac_y         * ( at_i(ji,jj+1) + at_i(ji,jj) ) !        ''            Y-   ''
 
          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_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
+      CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   &
+                 &          zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   &
+                 &          zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   & 
+                 &          zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    )
+
+      CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    )
+
+      CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1., &
+                 &          zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
+                 &          zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.  )
+
+      IF  ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  zdiag_xmtrp_ice(:,:)      )                            ! X-component of sea-ice mass transport (kg/s)
+      IF  ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice"     ,  zdiag_ymtrp_ice(:,:)      )                            ! Y-component of sea-ice mass transport 
+
+      IF  ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw"     ,  zdiag_xmtrp_snw(:,:)      )                            ! X-component of snow mass transport (kg/s)
+      IF  ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw"     ,  zdiag_ymtrp_snw(:,:)      )                            ! Y-component of snow mass transport
+
+      IF  ( iom_use( "xatrp"    ) ) CALL iom_put( "xatrp"        ,  zdiag_xatrp(:,:)          )                            ! X-component of ice area transport
+      IF  ( iom_use( "yatrp"    ) ) CALL iom_put( "yatrp"        ,  zdiag_yatrp(:,:)          )                            ! Y-component of ice area transport
+
+      IF  ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice"     ,  utau_ice(:,:)      * zswi(:,:) + zmiss(:,:) )          ! Wind stress term in force balance (x)
+      IF  ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice"     ,  vtau_ice(:,:)      * zswi(:,:) + zmiss(:,:) )          ! Wind stress term in force balance (y)
+
+      IF  ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   zdiag_utau_oi(:,:) * zswi(:,:) + zmiss(:,:) )          ! Ocean stress term in force balance (x)
+      IF  ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   zdiag_vtau_oi(:,:) * zswi(:,:) + zmiss(:,:) )          ! Ocean stress term in force balance (y)
+
+      IF  ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   zdiag_dssh_dx(:,:) * zswi(:,:) + zmiss(:,:) )          ! Sea-surface tilt term in force balance (x)
+      IF  ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   zdiag_dssh_dy(:,:) * zswi(:,:) + zmiss(:,:) )          ! Sea-surface tilt term in force balance (y)
+
+      IF  ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   zdiag_corstrx(:,:) * zswi(:,:) + zmiss(:,:) )          ! Coriolis force term in force balance (x)
+      IF  ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   zdiag_corstry(:,:) * zswi(:,:) + zmiss(:,:) )          ! Coriolis force term in force balance (y)
+
+      IF  ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   zdiag_intstrx(:,:) * zswi(:,:) + zmiss(:,:) )          ! Internal force term in force balance (x)
+      IF  ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   zdiag_intstry(:,:) * zswi(:,:) + zmiss(:,:) )          ! Internal force term in force balance (y)
+
+      IF  ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   zdiag_sig1(:,:)    * zswi(:,:) + zmiss(:,:) )          ! Normal stress
+      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   zdiag_sig2(:,:)    * zswi(:,:) + zmiss(:,:) )          ! Shear stress
 
       !
@@ -740,4 +810,14 @@
       CALL wrk_dealloc( jpi,jpj, zCorx, zCory )
       CALL wrk_dealloc( jpi,jpj, ztaux_oi, ztauy_oi )
+      
+      CALL wrk_dealloc( jpi,jpj, zdiag_xmtrp_ice, zdiag_ymtrp_ice )
+      CALL wrk_dealloc( jpi,jpj, zdiag_xmtrp_snw, zdiag_ymtrp_snw )
+      CALL wrk_dealloc( jpi,jpj, zdiag_xatrp    , zdiag_yatrp     )
+      CALL wrk_dealloc( jpi,jpj, zdiag_utau_oi  , zdiag_vtau_oi   )
+      CALL wrk_dealloc( jpi,jpj, zdiag_dssh_dx  , zdiag_dssh_dy   )
+      CALL wrk_dealloc( jpi,jpj, zdiag_corstrx  , zdiag_corstry   )
+      CALL wrk_dealloc( jpi,jpj, zdiag_intstrx  , zdiag_intstry   )
+      CALL wrk_dealloc( jpi,jpj, zdiag_sig1     , zdiag_sig2      )
+      CALL wrk_dealloc( jpi,jpj, zswi           , zmiss           )
 
    END SUBROUTINE lim_rhg
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 8283)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8285)
@@ -300,33 +300,6 @@
       IF  ( iom_use( "uice_mv"  ) ) CALL iom_put( "uice_mv"     ,   u_ice(:,:)               * zswi(:,:) + zmiss(:,:) )   ! ice velocity u component
       IF  ( iom_use( "vice_mv"  ) ) CALL iom_put( "vice_mv"     ,   v_ice(:,:)               * zswi(:,:) + zmiss(:,:) )   ! ice velocity v component
-      
-      IF  ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  diag_xmtrp_ice(:,:)      )                            ! X-component of sea-ice mass transport (kg/s)
-      IF  ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice"     ,  diag_ymtrp_ice(:,:)      )                            ! Y-component of sea-ice mass transport 
-
-      IF  ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw"     ,  diag_xmtrp_snw(:,:)      )                            ! X-component of snow mass transport (kg/s)
-      IF  ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw"     ,  diag_ymtrp_snw(:,:)      )                            ! Y-component of snow mass transport
-
-      IF  ( iom_use( "xatrp"    ) ) CALL iom_put( "xatrp"        ,  diag_xatrp(:,:)              )                        ! X-component of ice area transport
-      IF  ( iom_use( "yatrp"    ) ) CALL iom_put( "yatrp"        ,  diag_yatrp(:,:)              )                        ! Y-component of ice area transport
-
-      IF  ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice"     ,  utau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (x)
-      IF  ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice"     ,  vtau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (y)
-
-      IF  ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   diag_utau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (x)
-      IF  ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   diag_vtau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (y)
 
       IF  ( iom_use( "icestr"   ) ) CALL iom_put( "icestr"      ,   strength(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Ice strength
-
-      IF  ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (x)
-      IF  ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (y)
-
-      IF  ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (x)
-      IF  ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (y)
-
-      IF  ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (x)
-      IF  ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (y)
-
-      IF  ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Normal stress
-      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Shear stress
 
       !--------------------------------
