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 8257)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8258)
@@ -54,9 +54,10 @@
       !
       INTEGER  ::  ii, ji, jj, jk, jl  ! dummy loop indices
-      REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2
+      REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val
       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2
       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zfb          ! ice freeboard
       REAL(wp), POINTER, DIMENSION(:,:)   ::  zamask, zamask15 ! 15% concentration mask
+      REAL(wp), POINTER, DIMENSION(:,:)   ::  zmiss        ! missing value array
 
       ! Global ice diagnostics (SIMIP)
@@ -75,4 +76,11 @@
       CALL wrk_alloc( jpi, jpj     , z2d, zswi )
       CALL wrk_alloc( jpi, jpj     , zfb, zamask, zamask15 )
+      CALL wrk_alloc( jpi, jpj     , zmiss                 )
+
+      !-----------------------------
+      ! Missing value array 
+      !-----------------------------
+      zmiss_val  = -1.0e20
+      zmiss(:,:) = zmiss_val * ( 1. - zswi(:,:) )
 
       !-----------------------------
@@ -119,5 +127,8 @@
 
       ! velocity
-      IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 
+      IF ( iom_use('uice_ipa') ) CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
+      IF ( iom_use('vice_ipa') ) CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
+
+      IF ( ( iom_use( "icevel" ) .OR. ( iom_use( "icevel_mv" ) ) ) THEN 
          DO jj = 2 , jpjm1
             DO ji = 2 , jpim1
@@ -128,7 +139,6 @@
          END DO
          CALL lbc_lnk( z2d, 'T', 1. )
-         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component
-         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component
-         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module
+         IF ( iom_use( "icevel" )    CALL iom_put( "icevel"       , z2d        )                                            ! ice velocity module
+         IF ( iom_use( "icevel_mv" ) CALL iom_put( "icevel_mv"    , z2d(:,:) * zswi(:,:) + zmiss(:,:)                       ! ice velocity module (missing value)
       ENDIF
       !
@@ -239,65 +249,81 @@
       zrho1 = ( rau0 - rhoic ) / rau0; zrho2 = rhosn / rau0
 
-      IF  ( iom_use( "icethic"  ) ) CALL iom_put( "icethic"     , htm_i * zamask             )          ! 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 * zamask             )          ! 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 * 1.0e-3 * zswi )   ! Mass of salt in sea ice per cell area
+      IF  ( iom_use( "icepres"  ) ) CALL iom_put( "icepres"     , zswi(:,:)                     )                         ! Ice presence (1 or 0) 
+      IF  ( iom_use( "icemass"  ) ) CALL iom_put( "icemass"     , rhoic * vt_i(:,:) * zswi(:,:) )                         ! Ice mass per cell area 
+      IF  ( iom_use( "icethic"  ) ) CALL iom_put( "icethic"     , htm_i(:,:) * zamask(:,:)  * zswi(:,:) + zmiss(:,:) )    ! Ice thickness 
+      IF  ( iom_use( "snomass"  ) ) CALL iom_put( "snomass"     , rhosn * vt_s(:,:)         * zswi(:,:) + zmiss(:,:) )    ! Snow mass per cell area
+      IF  ( iom_use( "snowthic" ) ) CALL iom_put( "snowthic"    , htm_s(:,:) * zamask(:,:)  * zswi(:,:) + zmiss(:,:) )    ! Snow thickness       
+      IF  ( iom_use( "icesnt"   ) ) CALL iom_put( "icesnt"      , ( tm_si(:,:) - rt0 )      * zswi(:,:) + zmiss(:,:) )    ! Snow-ice interface temperature
+      IF  ( iom_use( "icebot"   ) ) CALL iom_put( "icebot"      , ( t_bo(:,:)  - rt0 )      * zswi(:,:) + zmiss(:,:) )    ! Ice bottom temperature
+      IF  ( iom_use( "iceage_mv") ) CALL iom_put( "iceage_mv"   , om_i(:,:) * zamask15(:,:) * zswi(:,:) + zmiss(:,:) )    ! Ice age
+      IF  ( iom_use( "icesmass" ) ) CALL iom_put( "icesmass"    , SUM( smv_i, DIM = 3 ) * rhoic * 1.0e-3 * zswi(:,:) )    ! Mass of salt in sea ice per cell area
+      IF  ( iom_use( "icesal_mv") ) CALL iom_put( "icesal_mv"   , smt_i(:,:)                * zswi(:,:) + zmiss(:,:) )    ! Ice salinity
+
       IF  ( iom_use( "icefb"    ) ) THEN
-         zfb(:,:) = ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) ) * zswi(:,:)                             
+         zfb(:,:) = ( zrho1 * htm_i(:,:) - zrho2 * htm_s(:,:) )                                         
          WHERE( zfb < 0._wp ) ;   zfb = 0._wp ;   END WHERE
-                                    CALL iom_put( "icefb"       , zfb                        )          ! Ice freeboard
+                                    CALL iom_put( "icefb"       , zfb(:,:)                  * zswi(:,:) + zmiss(:,:) )    ! Ice freeboard
       ENDIF
-      IF  ( iom_use( "dmithd"   ) ) CALL iom_put( "dmithd"      , - wfx_bog - wfx_bom - wfx_sum   &     ! Sea-ice mass change from thermodynamics
+
+      IF  ( iom_use( "snhc_mv"  ) ) CALL iom_put( "snhc_mv"     , et_s(:,:)                 * zswi(:,:) + zmiss(:,:) )    ! Snow total heat content
+
+      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( "dmssub"   ) ) CALL iom_put( "dmssub"      , - wfx_snw_sub              )          ! Snow mass change through sublimation
-      IF  ( iom_use( "dmisub"   ) ) CALL iom_put( "dmisub"      , - wfx_ice_sub              )          ! Sea-ice mass change through 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( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum                  )          ! Freshwater flux from sea-ice surface
-
-      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
-
-      IF  ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  diag_xmtrp_ice           )          ! X-component of sea-ice mass transport
-      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
-      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( "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( "dmssub"   ) ) CALL iom_put( "dmssub"      , - wfx_snw_sub              )                            ! Snow mass change through sublimation
+      IF  ( iom_use( "dmisub"   ) ) CALL iom_put( "dmisub"      , - wfx_ice_sub              )                            ! Sea-ice mass change through 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( "hfxconsu" ) ) CALL iom_put( "hfxsenso"    ,   -fhtur(:,:)              * zswi(:,:) + zmiss(:,:) )   ! Sensible oceanic heat flux
+      IF  ( iom_use( "hfxconbo" ) ) CALL iom_put( "hfxconbo"    ,   diag_fc_bo               * zswi(:,:) + zmiss(:,:) )   ! Bottom conduction flux
+      IF  ( iom_use( "hfxconsu" ) ) CALL iom_put( "hfxconsu"    ,   diag_fc_su               * zswi(:,:) + zmiss(:,:) )   ! Surface conduction flux
+
+      IF  ( iom_use( "wfxtot"   ) ) CALL iom_put( "wfxtot"      ,   wfx_ice(:,:)             * zswi(:,:) + zmiss(:,:) )   ! Total freshwater flux from sea ice
+      IF  ( iom_use( "wfxsum"   ) ) CALL iom_put( "wfxsum"      ,   wfx_sum(:,:)             * zswi(:,:) + zmiss(:,:) )   ! Freshwater flux from sea-ice surface
+      IF  ( iom_use( "sfx_mv"   ) ) CALL iom_put( "sfx_mv"      ,   sfx(:,:)                 * zswi(:,:) + zmiss(:,:) )   ! Total salt flux
+
+      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_ai_mv") ) CALL iom_put( "utau_ai_mv"   ,  utau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (x)
+      IF  ( iom_use("vtau_ai_mv") ) CALL iom_put( "vtau_ai_mv"   ,  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_mv") ) CALL iom_put( "icestr_mv"   ,   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(:,:)   )                        ! Normal stress
+      IF  ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2(:,:) * zswi(:,:)   )                        ! Shear stress
 
       !--------------------------------
@@ -354,4 +380,5 @@
       CALL wrk_dealloc( jpi, jpj     , z2d, zswi )
       CALL wrk_dealloc( jpi, jpj     , zfb, zamask, zamask15 )
+      CALL wrk_dealloc( jpi, jpj     , zmiss                 )
 
       IF( nn_timing == 1 )  CALL timing_stop('limwri')
