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 7506)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7517)
@@ -430,4 +430,7 @@
    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
 
    !
@@ -538,4 +541,5 @@
                 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 7506)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 7517)
@@ -129,4 +129,5 @@
       REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
       REAL(wp) ::   zintb, zintn                                             ! dummy argument
+      REAL(wp) ::   zswi
       
       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength
@@ -639,5 +640,5 @@
       stress12_i(:,:) = zs12(:,:)
 
-      ! SIMIP diagnostic internal stress
+      ! SIMIP diagnostic internal stress terms (N/m2)
       diag_dssh_dx(:,:) = zspgU(:,:)
       diag_dssh_dy(:,:) = zspgV(:,:)
@@ -649,4 +650,15 @@
       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. )
 
       !
@@ -686,4 +698,5 @@
          ENDIF
       ENDIF
+
       !
       CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt )
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 7506)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 7517)
@@ -18,4 +18,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE ice
+   USE dom_ice
    USE limvar
    USE in_out_manager
@@ -52,12 +53,33 @@
       INTEGER, INTENT(in) ::   kindic   ! if kindic < 0 there has been an error somewhere
       !
-      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices
+      INTEGER  ::  ii, ji, jj, jk, jl  ! dummy loop indices
       REAL(wp) ::  z1_365
       REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2
-      REAL(wp) ::  zdiag_area_nh, zdiag_extt_nh, zdiag_area_sh, zdiag_extt_sh, zdiag_volu_nh, zdiag_volu_sh
-      
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2
       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace
 
+      ! Global ice diagnostics (SIMIP)
+      REAL(wp) ::  zdiag_area_nh, &   ! area, extent, volume
+         &         zdiag_extt_nh, &
+         &         zdiag_area_sh, & 
+         &         zdiag_extt_sh, & 
+         &         zdiag_volu_nh, & 
+         &         zdiag_volu_sh 
+
+      ! Strait / passage fluxes (SIMIP)
+      REAL(wp), DIMENSION(4) ::  &                                    ! Strait fluxes for output
+         &         zdiag_area_strait  ,   &                                    ! 1=Fram Strait, 2=CAA, 3= Barents, 4 = Bering
+         &         zdiag_mice_strait  ,   &
+         &         zdiag_msno_strait
+
+      REAL(wp) :: zfarea_u, zfmice_u, zfmsno_u, zfarea_v, zfmice_v, zfmsno_v   ! dummy fluxes
+
+      REAL(wp), DIMENSION(11) :: &
+         &         zui, zuj, zvi, zvj                                          ! strait addresses
+
+      INTEGER  :: Nu, Nv                                                       ! passage size
+
+      INTEGER, PARAMETER :: i_grid = 2                                         ! grid type (eORCA1 = 1, ORCA2 = 2)
+      
       !!-------------------------------------------------------------------
 
@@ -139,6 +161,6 @@
 
       CALL iom_put( "icestr"      , strength * zswi )    ! ice strength
-      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence
-      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear
+      CALL iom_put( "idive"       , divu_i              )    ! divergence
+      CALL iom_put( "ishear"      , shear_i             )    ! shear
       CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume
       
@@ -257,20 +279,26 @@
       CALL iom_put( "intstry"     ,   diag_intstry*zswi        )          ! Internal force term in force balance (y-component)
 
+      CALL iom_put( "normstr"     ,   diag_sig1   *zswi        )          ! Normal stress
+      CALL iom_put( "sheastr"     ,   diag_sig2   *zswi        )          ! Shear stress
+
       !--------------------------------
-      ! Global ice diagnostics
+      ! Global ice diagnostics (SIMIP)
       !--------------------------------
-
-      IF ( iom_use ( "NH_icearea" ) .OR. ( iom_use( "NH_icevolu" ) ) THEN             ! NH ice area 
-
-           WHERE( fcor > 0 ); zswi(:,:) = 1.; 
-           ELSEWHERE        ; zswi(:,:) = 0.
-           END WHERE 
-
-           IF ( iom_use ( "NH_icearea" ) THEN
-              zdiag_area_nh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) ) / 1.0e12
+      !
+      ! 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(:,:) ) / 1.0e12
+           IF ( iom_use( "NH_icevolu" ) ) THEN
+              zdiag_volu_nh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) )
               CALL iom_put( "NH_icevolu", zdiag_volu_nh )
            ENDIF
@@ -278,47 +306,262 @@
       ENDIF 
 
-      IF ( iom_use ( "NH_iceextt" ) ) THEN                                            ! NH ice extt 
-
-           WHERE( fcor > 0 .AND. at_i > 0.15 ); zswi(:,:) = 1.; 
-           ELSEWHERE                          ; zswi(:,:) = 0.
-           END WHERE 
-
-           zdiag_extt_nh = glob_sum( zswi(:,:) * e12t(:,:) ) / 1.0e12
-           CALL iom_put( "NH_iceextt", zdiag_extt_nh )
+      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.; 
-           ELSEWHERE        ; zswi(:,:) = 0.
-           END WHERE 
-
-           IF ( iom_use ( "SH_icearea" ) THEN
-              zdiag_area_sh = glob_sum( at_i(:,:) * zswi(:,:) * e12t(:,:) ) / 1.0e12
-              CALL iom_put( "SH_icearea", zdiag_area_sh )
-           ENDIF
-           IF ( iom_use ( "SH_icevolu" ) THEN
-              zdiag_volu_sh = glob_sum( vt_i(:,:) * zswi(:,:) * e12t(:,:) ) / 1.0e12
-              CALL iom_put( "SH_icevolu", zdiag_volu_sh )
-           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.; 
-           ELSEWHERE                          ; zswi(:,:) = 0.
-           END WHERE 
-
-           zdiag_extt_sh = glob_sum( zswi(:,:) * e12t(:,:) )/ 1.0e12
-
-           CALL iom_put( "SH_iceextt", zdiag_extt_sh )
+      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 
 
-      !--------------------------------
-      ! Output values for each category
-      !--------------------------------
+      !------------------------
+      ! 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
+      !----------------------------------
       IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories
       IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories
