Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ice_ref	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/cfgs/SHARED/namelist_ice_ref	(revision 11501)
@@ -233,4 +233,6 @@
 !------------------------------------------------------------------------------
    ln_icediachk     = .false.         !  check online the heat, mass & salt budgets (T) or not (F)
+      rn_icechk     =  1.             !     1. <=> 1mm of ice spuriously gained/lost during 1   year  (at any gridcell)
+                                      !                           --                        100 years (over the entire ice cover)
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/ice.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/ice.F90	(revision 11501)
@@ -196,4 +196,6 @@
    !                                     !!** ice-diagnostics namelist (namdia) **
    LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)
+   REAL(wp), PUBLIC ::   rn_icechk        !: 1. <=> 1mm of ice spuriously gained/lost during 1   year  (at any gridcell)
+                                          !                           --                     100 years (globally)
    LOGICAL , PUBLIC ::   ln_icediahsb     !: flag for ice diag (T) or not (F)
    LOGICAL , PUBLIC ::   ln_icectl        !: flag for sea-ice points output (T) or not (F)
@@ -351,7 +353,4 @@
    !! * Ice diagnostics
    !!----------------------------------------------------------------------
-   ! thd refers to changes induced by thermodynamics
-   ! trp   ''         ''     ''       advection (transport of ice)
-   !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume
@@ -365,4 +364,13 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s] 
 
+   !!----------------------------------------------------------------------
+   !! * Ice conservation
+   !!----------------------------------------------------------------------
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_v        !: conservation of ice volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_s        !: conservation of ice salt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_t        !: conservation of ice heat
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fv       !: conservation of ice volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fs       !: conservation of ice salt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_ft       !: conservation of ice heat
    !
    !!----------------------------------------------------------------------
@@ -389,5 +397,5 @@
       INTEGER :: ice_alloc
       !
-      INTEGER :: ierr(15), ii
+      INTEGER :: ierr(16), ii
       !!-----------------------------------------------------------------
       ierr(:) = 0
@@ -461,4 +469,9 @@
          &      diag_sice  (jpi,jpj) , diag_vice   (jpi,jpj) , diag_vsnw  (jpi,jpj), STAT=ierr(ii) )
 
+      ! * Ice conservation
+      ii = ii + 1
+      ALLOCATE( diag_v (jpi,jpj) , diag_s (jpi,jpj) , diag_t (jpi,jpj),   & 
+         &      diag_fv(jpi,jpj) , diag_fs(jpi,jpj) , diag_ft(jpi,jpj), STAT=ierr(ii) )
+      
       ! * SIMIP diagnostics
       ii = ii + 1
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icecor.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icecor.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icecor.F90	(revision 11501)
@@ -60,4 +60,5 @@
       IF( ln_timing    )   CALL timing_start('icecor')                                                             ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icecor',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp .AND. kn == 2 ) THEN
@@ -164,9 +165,10 @@
       !
       ! controls
-      IF( ln_icediachk   )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_ctl         )   CALL ice_prt3D   ('icecor')                                                             ! prints
+      IF( ln_ctl       )   CALL ice_prt3D   ('icecor')                                                             ! prints
       IF( ln_icectl .AND. kn == 2 ) &
-         &                   CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints
-      IF( ln_timing      )   CALL timing_stop ('icecor')                                                             ! timing
+         &                 CALL ice_prt     ( kt, iiceprt, jiceprt, 2, ' - Final state - ' )                       ! prints
+      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icecor', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icecor',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
+      IF( ln_timing    )   CALL timing_stop ('icecor')                                                             ! timing
       !
    END SUBROUTINE ice_cor
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icectl.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icectl.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icectl.F90	(revision 11501)
@@ -12,6 +12,7 @@
    !!   'key_si3'                                       SI3 sea-ice model
    !!----------------------------------------------------------------------
-   !!    ice_cons_hsm     : conservation tests on heat, salt and mass
-   !!    ice_cons_final   : conservation tests on heat, salt and mass at end of time step
+   !!    ice_cons_hsm     : conservation tests on heat, salt and mass during a  time step (global) 
+   !!    ice_cons_final   : conservation tests on heat, salt and mass at end of time step (global)
+   !!    ice_cons2D       : conservation tests on heat, salt and mass at each gridcell
    !!    ice_ctl          : control prints in case of crash
    !!    ice_prt          : control prints at a given grid point
@@ -27,4 +28,5 @@
    !
    USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
@@ -37,8 +39,18 @@
    PUBLIC   ice_cons_hsm
    PUBLIC   ice_cons_final
+   PUBLIC   ice_cons2D
    PUBLIC   ice_ctl
    PUBLIC   ice_prt
    PUBLIC   ice_prt3D
 
+   ! thresold values for conservation
+   !    these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk
+   REAL(wp), PARAMETER ::   zchk_m   = 1.e-5   ! kg/m2/s <=> 1mm of ice per 1   year  spuriously gained/lost
+   REAL(wp), PARAMETER ::   zchk2D_m = 1.e-7   !                   --       100 years          --
+   REAL(wp), PARAMETER ::   zchk_s   = 1.e-4   ! g/m2/s  <=> 1mm of ice per 1   year  spuriously gained/lost (considering s=10g/kg)
+   REAL(wp), PARAMETER ::   zchk2D_s = 1.e-6   !                   --       100 years          --
+   REAL(wp), PARAMETER ::   zchk_t   = 3.      ! W/m2    <=> 1mm of ice per 1   year  spuriously gained/lost (considering Lf=3e5J/kg)
+   REAL(wp), PARAMETER ::   zchk2D_t = 0.03    !                   --       100 years          --
+   
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -59,5 +71,5 @@
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to
+      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine violations are set to
       !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
@@ -68,109 +80,86 @@
       REAL(wp)        , INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
       !!
-      REAL(wp) ::   zv, zs, zt, zfs, zfv, zft
-      REAL(wp) ::   zvmin, zamin, zamax, zeimin, zesmin, zsmin
+      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat, &
+         &          zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin
       REAL(wp) ::   zvtrp, zetrp
-      REAL(wp) ::   zarea, zv_sill, zs_sill, zt_sill
-      REAL(wp), PARAMETER ::   zconv = 1.e-9 ! convert W to GW and kg to Mt
+      REAL(wp) ::   zarea
       !!-------------------------------------------------------------------
       !
       IF( icount == 0 ) THEN
-         !                          ! water flux
-         pdiag_fv = glob_sum( 'icectl',                                                                       &
-            &                 -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  &
-            &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  &
-            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  &
-            &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  &
-            &                  ) * e1e2t(:,:) ) * zconv
+
+         pdiag_v = glob_sum( 'icectl',   SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t )
+         pdiag_s = glob_sum( 'icectl',   SUM( sv_i * rhoi            , dim=3 ) * e1e2t )
+         pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t )
+
+         ! mass flux
+         pdiag_fv = glob_sum( 'icectl',  &
+            &                         ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
+            &                           wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t )
+         ! salt flux
+         pdiag_fs = glob_sum( 'icectl',  &
+            &                         ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + &
+            &                           sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t )
+         ! heat flux
+         pdiag_ft = glob_sum( 'icectl',  &
+            &                         (   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  &
+            &                           - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )
+
+      ELSEIF( icount == 1 ) THEN
+
+         ! -- mass diag -- !
+         zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_rdtice       &
+            &         + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn +       &
+            &                                 wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + &
+            &                                 wfx_ice_sub + wfx_spr ) * e1e2t )                                           &
+            &         - pdiag_fv
          !
-         !                          ! salt flux
-         pdiag_fs = glob_sum( 'icectl',                                                                     &
-            &                  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
-            &                  ) *  e1e2t(:,:) ) * zconv 
+         ! -- salt diag -- !
+         zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice  &
+            &         + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni +           &
+            &                                 sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) &
+            &         - pdiag_fs
          !
-         !                          ! heat flux
-         pdiag_ft = glob_sum( 'icectl',                                                                    &
-            &                  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &                  - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                  ) *  e1e2t(:,:) ) * zconv
-
-         pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t * zconv )
-
-         pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi            , dim=3 ) * e1e2t * zconv )
-
-         pdiag_t = glob_sum( 'icectl', (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
-            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv
-
-      ELSEIF( icount == 1 ) THEN
-
-         ! water flux
-         zfv = glob_sum( 'icectl',                                                                        &
-            &             -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  &
-            &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  &
-            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  &
-            &                wfx_ice_sub(:,:) + wfx_spr(:,:)  &
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv
-
-         ! salt flux
-         zfs = glob_sum( 'icectl',                                                                       &
-            &              ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_fs
-
-         ! heat flux
-         zft = glob_sum( 'icectl',                                                                      &
-            &              ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_ft
- 
-         ! outputs
-         zv = ( ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) * zconv  &
-            &     - pdiag_v ) * r1_rdtice - zfv ) * rday
-
-         zs = ( ( glob_sum( 'icectl', SUM( sv_i * rhoi            , dim=3 ) * e1e2t ) * zconv  &
-            &     - pdiag_s ) * r1_rdtice + zfs ) * rday
-
-         zt = ( glob_sum( 'icectl',                                                                &
-            &             (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )                       &
-            &              + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv   &
-            &   - pdiag_t ) * r1_rdtice + zft
-
-         ! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
-         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t  ) * zconv * rday 
-         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t  ) * zconv
-
-         zamax  = glob_max( 'icectl', SUM( a_i, dim=3 ) )
-         zvmin  = glob_min( 'icectl', v_i )
-         zamin  = glob_min( 'icectl', a_i )
-         zsmin  = glob_min( 'icectl', sv_i )
-         zeimin = glob_min( 'icectl', SUM( e_i, dim=3 ) )
-         zesmin = glob_min( 'icectl', SUM( e_s, dim=3 ) )
-
-         ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
-         zarea   = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2
-         zv_sill = zarea * 2.5e-5
-         zs_sill = zarea * 25.e-5
-         zt_sill = zarea * 10.e-5
-
-         IF(lwp) THEN
+         ! -- heat diag -- !
+         zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t &
+            &         ) * r1_rdtice                                                                                           &
+            &         + glob_sum( 'icectl', (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                      &
+            &                                - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )                    &
+            &         - pdiag_ft
+
+         ! -- min/max diag -- !
+         zdiag_amax  = glob_max( 'icectl', SUM( a_i, dim=3 ) )
+         zdiag_vmin  = glob_min( 'icectl', v_i )
+         zdiag_amin  = glob_min( 'icectl', a_i )
+         zdiag_smin  = glob_min( 'icectl', sv_i )
+         zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) )
+         zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) )
+
+         ! -- advection scheme is conservative? -- !
+         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0
+         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0
+
+         ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
+         zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t )
+
+         IF( lwp ) THEN
             ! check conservation issues
-            IF ( ABS( zv ) > zv_sill )   WRITE(numout,*) 'violation volume [Mt/day]     (',cd_routine,') = ',zv
-            IF ( ABS( zs ) > zs_sill )   WRITE(numout,*) 'violation saline [Mkg/day]    (',cd_routine,') = ',zs
-            IF ( ABS( zt ) > zt_sill )   WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',zt
+            IF( ABS(zdiag_mass) > zchk_m * rn_icechk * zarea ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice
+            IF( ABS(zdiag_salt) > zchk_s * rn_icechk * zarea ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice
+            IF( ABS(zdiag_heat) > zchk_t * rn_icechk * zarea ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice
+            ! check negative values
+            IF( zdiag_vmin  < 0. ) WRITE(numout,*)   cd_routine,' : violation v_i < 0         = ',zdiag_vmin
+            IF( zdiag_amin  < 0. ) WRITE(numout,*)   cd_routine,' : violation a_i < 0         = ',zdiag_amin
+            IF( zdiag_smin  < 0. ) WRITE(numout,*)   cd_routine,' : violation s_i < 0         = ',zdiag_smin
+            IF( zdiag_eimin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_i < 0         = ',zdiag_eimin
+            IF( zdiag_esmin < 0. ) WRITE(numout,*)   cd_routine,' : violation e_s < 0         = ',zdiag_esmin
             ! check maximum ice concentration
-            IF ( zamax > MAX( rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' )  &
-               &                         WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
-            ! check negative values
-            IF ( zvmin  < 0. )           WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
-            IF ( zamin  < 0. )           WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
-            IF ( zsmin  < 0. )           WRITE(numout,*) 'violation s_i<0               (',cd_routine,') = ',zsmin
-            IF ( zeimin < 0. )           WRITE(numout,*) 'violation e_i<0               (',cd_routine,') = ',zeimin
-            IF ( zesmin < 0. )           WRITE(numout,*) 'violation e_s<0               (',cd_routine,') = ',zesmin
-!clem: the following check fails (I think...)
-!            IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
-!                                           WRITE(numout,*) 'violation vtrp [Mt/day]       (',cd_routine,') = ',zvtrp
-!                                           WRITE(numout,*) 'violation etrp [GW]           (',cd_routine,') = ',zetrp
-!            ENDIF
+            IF( zdiag_amax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax
+            ! check if advection scheme is conservative
+            IF( ABS(zvtrp) > zchk_m*rn_icechk*zarea .AND. cd_routine == 'icedyn_adv' ) &
+               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
          ENDIF
          !
@@ -179,5 +168,4 @@
    END SUBROUTINE ice_cons_hsm
 
-
    SUBROUTINE ice_cons_final( cd_routine )
       !!-------------------------------------------------------------------
@@ -188,40 +176,175 @@
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to
+      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine the violation are set to
       !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
       !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
       !!-------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in)    :: cd_routine    ! name of the routine
-      REAL(wp)                        :: zqmass, zhfx, zsfx, zvfx
-      REAL(wp)                        :: zarea, zv_sill, zs_sill, zt_sill
-      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt
+      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine
+      REAL(wp) ::   zdiag_mass, zdiag_salt, zdiag_heat
+      REAL(wp) ::   zarea
       !!-------------------------------------------------------------------
 
       ! water flux
-      zvfx  = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday
-
-      ! salt flux
-      zsfx  = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t ) * zconv * rday
-
-      ! heat flux
+      ! -- mass diag -- !
+      zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t )
+
+      ! -- salt diag -- !
+      zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice ) * e1e2t )
+
+      ! -- heat diag -- !
       ! clem: not the good formulation
-      !!zhfx  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  &
-      !!   &                        ) * e1e2t ) * zconv
-
-      ! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice) 
-      zarea   = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2
-      zv_sill = zarea * 2.5e-5
-      zs_sill = zarea * 25.e-5
-      zt_sill = zarea * 10.e-5
-
-      IF(lwp) THEN
-         IF( ABS( zvfx ) > zv_sill )   WRITE(numout,*) 'violation vfx  [Mt/day]       (',cd_routine,') = ',zvfx
-         IF( ABS( zsfx ) > zs_sill )   WRITE(numout,*) 'violation sfx  [Mkg/day]      (',cd_routine,') = ',zsfx
-         !!IF( ABS( zhfx ) > zt_sill )   WRITE(numout,*) 'violation hfx  [GW]           (',cd_routine,') = ',zhfx
+      !!zdiag_heat  = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat + hfx_thd + hfx_dyn + hfx_res + hfx_sub + hfx_spr  &
+      !!   &                              ) * e1e2t )
+
+      ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
+      zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t )
+
+      IF( lwp ) THEN
+         IF( ABS(zdiag_mass) > zchk_m * rn_icechk * zarea ) &
+            &                   WRITE(numout,*) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass * rdt_ice
+         IF( ABS(zdiag_salt) > zchk_s * rn_icechk * zarea ) &
+            &                   WRITE(numout,*) cd_routine,' : violation salt cons. [g]  = ',zdiag_salt * rdt_ice
+         !!IF( ABS(zdiag_heat) > zchk_t * rn_icechk * zarea ) WRITE(numout,*) cd_routine,' : violation heat cons. [J]  = ',zdiag_heat * rdt_ice
       ENDIF
       !
    END SUBROUTINE ice_cons_final
 
+   SUBROUTINE ice_cons2D( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft )
+      !!-------------------------------------------------------------------
+      !!                       ***  ROUTINE ice_cons2D ***
+      !!
+      !! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
+      !!                     + test if ice concentration and volume are > 0
+      !!
+      !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
+      !!              It stops the code if there is a violation of conservation at any gridcell
+      !!-------------------------------------------------------------------
+      INTEGER         , INTENT(in) ::   icount        ! called at: =0 the begining of the routine, =1  the end
+      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine
+      REAL(wp)        , DIMENSION(jpi,jpj), INTENT(inout) ::   pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
+      !!
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdiag_mass, zdiag_salt, zdiag_heat, &
+         &                              zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax  
+      INTEGER ::   jl, jk
+      LOGICAL ::   ll_stop_m = .FALSE.
+      LOGICAL ::   ll_stop_s = .FALSE.
+      LOGICAL ::   ll_stop_t = .FALSE.
+      CHARACTER(len=120) ::   clnam   ! filename for the output
+      !!-------------------------------------------------------------------
+      !
+      IF( icount == 0 ) THEN
+
+         pdiag_v = SUM( v_i  * rhoi + v_s * rhos, dim=3 )
+         pdiag_s = SUM( sv_i * rhoi             , dim=3 )
+         pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 )
+
+         ! mass flux
+         pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd  +  &
+            &       wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr
+         ! salt flux
+         pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam 
+         ! heat flux
+         pdiag_ft =   hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw  & 
+            &       - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr
+
+      ELSEIF( icount == 1 ) THEN
+
+         ! -- mass diag -- !
+         zdiag_mass =   ( SUM( v_i * rhoi + v_s * rhos, dim=3 ) - pdiag_v ) * r1_rdtice                             &
+            &         + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
+            &             wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr )           &
+            &         - pdiag_fv
+         IF( MAXVAL( ABS(zdiag_mass) ) > zchk2D_m * rn_icechk )   ll_stop_m = .TRUE.
+         !
+         ! -- salt diag -- !
+         zdiag_salt =   ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice                                                  &
+            &         + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) &
+            &         - pdiag_fs
+         IF( MAXVAL( ABS(zdiag_salt) ) > zchk2D_s * rn_icechk )   ll_stop_s = .TRUE.
+         !
+         ! -- heat diag -- !
+         zdiag_heat =   ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_rdtice &
+            &         + (  hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw                                & 
+            &            - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr )                                        &
+            &         - pdiag_ft
+         IF( MAXVAL( ABS(zdiag_heat) ) > zchk2D_t * rn_icechk )   ll_stop_t = .TRUE.
+         !
+         ! -- other diags -- !
+         ! a_i < 0
+         zdiag_amin(:,:) = 0._wp
+         DO jl = 1, jpl
+            WHERE( a_i(:,:,jl) < 0._wp )   zdiag_amin(:,:) = 1._wp
+         ENDDO
+         ! v_i < 0
+         zdiag_vmin(:,:) = 0._wp
+         DO jl = 1, jpl
+            WHERE( v_i(:,:,jl) < 0._wp )   zdiag_vmin(:,:) = 1._wp
+         ENDDO
+         ! s_i < 0
+         zdiag_smin(:,:) = 0._wp
+         DO jl = 1, jpl
+            WHERE( s_i(:,:,jl) < 0._wp )   zdiag_smin(:,:) = 1._wp
+         ENDDO
+         ! e_i < 0
+         zdiag_emin(:,:) = 0._wp
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               WHERE( e_i(:,:,jk,jl) < 0._wp )   zdiag_emin(:,:) = 1._wp
+            ENDDO
+         ENDDO
+         ! a_i > amax
+         !WHERE( SUM( a_i, dim=3 ) > ( MAX( rn_amax_n, rn_amax_s ) + epsi10 )   ;   zdiag_amax(:,:) = 1._wp
+         !ELSEWHERE                                                             ;   zdiag_amax(:,:) = 0._wp
+         !END WHERE
+
+         IF( ll_stop_m .OR. ll_stop_s .OR. ll_stop_t ) THEN
+            clnam = 'diag_ice_conservation_'//cd_routine
+            CALL ice_cons_wri( clnam, zdiag_mass, zdiag_salt, zdiag_heat, zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin )
+         ENDIF
+
+         IF( ll_stop_m )   CALL ctl_stop( 'STOP', cd_routine//': ice mass conservation issue' )
+         IF( ll_stop_s )   CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' )
+         IF( ll_stop_t )   CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' )
+         
+      ENDIF
+
+   END SUBROUTINE ice_cons2D
+
+   SUBROUTINE ice_cons_wri( cdfile_name, pdiag_mass, pdiag_salt, pdiag_heat, pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE ice_cons_wri  ***
+      !!        
+      !! ** Purpose :   create a NetCDF file named cdfile_name which contains 
+      !!                the instantaneous fields when conservation issue occurs
+      !!
+      !! ** Method  :   NetCDF files using ioipsl
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*), INTENT( in ) ::   cdfile_name      ! name of the file created
+      REAL(wp), DIMENSION(:,:), INTENT( in ) ::   pdiag_mass, pdiag_salt, pdiag_heat, &
+         &                                        pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax  
+      !!
+      INTEGER ::   inum
+      !!----------------------------------------------------------------------
+      ! 
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~  named :', cdfile_name, '...nc'
+      IF(lwp) WRITE(numout,*)                
+
+      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
+      
+      CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 )    ! ice mass spurious lost/gain
+      CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 )    ! ice salt spurious lost/gain
+      CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 )    ! ice heat spurious lost/gain
+      ! other diags
+      CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 )    ! 
+      CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 )    ! 
+      CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 )    ! 
+      CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 )    ! 
+      
+      CALL iom_close( inum )
+
+   END SUBROUTINE ice_cons_wri
    
    SUBROUTINE ice_ctl( kt )
@@ -246,14 +369,9 @@
       ialert_id = 2 ! reference number of this alert
       cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert
-
       DO jl = 1, jpl
          DO jj = 1, jpj
             DO ji = 1, jpi
                IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN
-                  !WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
-                  !WRITE(numout,*) ' at_i     ', at_i(ji,jj)
-                  !WRITE(numout,*) ' Point - category', ji, jj, jl
-                  !WRITE(numout,*) ' a_i *** a_i_b   ', a_i      (ji,jj,jl), a_i_b  (ji,jj,jl)
-                  !WRITE(numout,*) ' v_i *** v_i_b   ', v_i      (ji,jj,jl), v_i_b  (ji,jj,jl)
+                  WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                ENDIF
@@ -269,4 +387,5 @@
          DO ji = 1, jpi
             IF(   h_i(ji,jj,jl)  >  50._wp   ) THEN
+               WRITE(numout,*) ' ALERTE 3 :   Very thick ice'
                !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -280,15 +399,21 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5  .AND.  &
+            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2.  .AND.  &
                &  at_i(ji,jj) > 0._wp   ) THEN
+               WRITE(numout,*) ' ALERTE 4 :   Very fast ice'
                !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
-               !WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
-               !WRITE(numout,*) ' oceanic stress utau      : ', utau(ji,jj) 
-               !WRITE(numout,*) ' oceanic stress vtau      : ', vtau(ji,jj)
-               !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj) 
-               !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj)
-               !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj)
-               !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj)
-               !WRITE(numout,*) 
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+            ENDIF
+         END DO
+      END DO
+
+      ! Alert on salt flux
+      ialert_id = 5 ! reference number of this alert
+      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth
+               WRITE(numout,*) ' ALERTE 5 :   High salt flux'
+               !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
             ENDIF
@@ -302,14 +427,6 @@
          DO ji = 1, jpi
             IF(   tmask(ji,jj,1) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN 
+               WRITE(numout,*) ' ALERTE 6 :   Ice on continents'
                !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
-               !WRITE(numout,*) ' masks s, u, v        : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1) 
-               !WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
-               !WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
-               !WRITE(numout,*) ' at_i(ji,jj)          : ', at_i(ji,jj)
-               !WRITE(numout,*) ' v_ice(ji,jj)         : ', v_ice(ji,jj)
-               !WRITE(numout,*) ' v_ice(ji,jj-1)       : ', v_ice(ji,jj-1)
-               !WRITE(numout,*) ' u_ice(ji-1,jj)       : ', u_ice(ji-1,jj)
-               !WRITE(numout,*) ' u_ice(ji,jj)         : ', v_ice(ji,jj)
-               !
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
             ENDIF
@@ -325,8 +442,6 @@
             DO ji = 1, jpi
                IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN
+                  WRITE(numout,*) ' ALERTE 7 :   Very fresh ice'
 !                 CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
-!                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
-!                 WRITE(numout,*) ' sss                  : ', sss_m(ji,jj)
-!                 WRITE(numout,*) 
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                ENDIF
@@ -335,4 +450,19 @@
       END DO
 !
+      ! Alert if qns very big
+      ialert_id = 8 ! reference number of this alert
+      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN
+               !
+               WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
+               !CALL ice_prt( kt, ji, jj, 2, '   ')
+               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+               !
+            ENDIF
+         END DO
+      END DO
+      !+++++
 
 !     ! Alert if too old ice
@@ -345,4 +475,5 @@
                       ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. &
                              ( a_i(ji,jj,jl) > 0._wp ) ) THEN
+                  WRITE(numout,*) ' ALERTE 9 :   Wrong ice age'
                   !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -351,44 +482,5 @@
          END DO
       END DO
- 
-      ! Alert on salt flux
-      ialert_id = 5 ! reference number of this alert
-      cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth
-               !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
-               !DO jl = 1, jpl
-                  !WRITE(numout,*) ' Category no: ', jl
-                  !WRITE(numout,*) ' a_i        : ', a_i      (ji,jj,jl) , ' a_i_b      : ', a_i_b  (ji,jj,jl)   
-                  !WRITE(numout,*) ' v_i        : ', v_i      (ji,jj,jl) , ' v_i_b      : ', v_i_b  (ji,jj,jl)   
-                  !WRITE(numout,*) ' '
-               !END DO
-               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
-            ENDIF
-         END DO
-      END DO
-
-      ! Alert if qns very big
-      ialert_id = 8 ! reference number of this alert
-      cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN
-               !
-               !WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
-               !WRITE(numout,*) ' ji, jj    : ', ji, jj
-               !WRITE(numout,*) ' qns       : ', qns(ji,jj)
-               !WRITE(numout,*) ' sst       : ', sst_m(ji,jj)
-               !WRITE(numout,*) ' sss       : ', sss_m(ji,jj)
-               !
-               !CALL ice_prt( kt, ji, jj, 2, '   ')
-               inb_alp(ialert_id) = inb_alp(ialert_id) + 1
-               !
-            ENDIF
-         END DO
-      END DO
-      !+++++
- 
+  
       ! Alert if very warm ice
       ialert_id = 10 ! reference number of this alert
@@ -400,13 +492,8 @@
                DO ji = 1, jpi
                   ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0
-                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   &
-                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
-                     !WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
-                     !WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
-                     !WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
-                     !WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
-                     !WRITE(numout,*) ' sz_i: ', sz_i(ji,jj,jk,jl)
-                     !WRITE(numout,*) ' ztmelts : ', ztmelts
-                     inb_alp(ialert_id) = inb_alp(ialert_id) + 1
+                  IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   &
+                     &                            .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
+                     WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
+                    inb_alp(ialert_id) = inb_alp(ialert_id) + 1
                   ENDIF
                END DO
@@ -435,5 +522,4 @@
    END SUBROUTINE ice_ctl
  
-   
    SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 )
       !!-------------------------------------------------------------------
@@ -443,5 +529,5 @@
       !!                in ocean.ouput 
       !!                3 possibilities exist 
-      !!                n = 1/-1 -> simple ice state (plus Mechanical Check if -1)
+      !!                n = 1/-1 -> simple ice state
       !!                n = 2    -> exhaustive state
       !!                n = 3    -> ice/ocean salt fluxes
@@ -482,6 +568,6 @@
                WRITE(numout,*) ' - Cell values '
                WRITE(numout,*) '   ~~~~~~~~~~~ '
-               WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj)
                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
+               WRITE(numout,*) ' ato_i         : ', ato_i(ji,jj)       
                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
                WRITE(numout,*) ' vt_s          : ', vt_s(ji,jj)       
@@ -503,13 +589,4 @@
                END DO
             ENDIF
-            IF( kn == -1 ) THEN
-               WRITE(numout,*) ' Mechanical Check ************** '
-               WRITE(numout,*) ' Check what means ice divergence '
-               WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj)
-               WRITE(numout,*) ' Total lead fraction     ', ato_i(ji,jj)
-               WRITE(numout,*) ' Sum of both             ', ato_i(ji,jj) + at_i(ji,jj)
-               WRITE(numout,*) ' Sum of both minus 1     ', ato_i(ji,jj) + at_i(ji,jj) - 1.00
-            ENDIF
-            
 
             !--------------------
@@ -525,5 +602,4 @@
                WRITE(numout,*) ' - Cell values '
                WRITE(numout,*) '   ~~~~~~~~~~~ '
-               WRITE(numout,*) ' cell area     : ', e1e2t(ji,jj)
                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
@@ -624,6 +700,6 @@
       !!
       !!-------------------------------------------------------------------
-      CHARACTER(len=*), INTENT(in)  :: cd_routine    ! name of the routine
-      INTEGER                       :: jk, jl        ! dummy loop indices
+      CHARACTER(len=*), INTENT(in) ::   cd_routine    ! name of the routine
+      INTEGER                      ::   jk, jl        ! dummy loop indices
       
       CALL prt_ctl_info(' ========== ')
@@ -684,5 +760,5 @@
       
    END SUBROUTINE ice_prt3D
-
+      
 #else
    !!----------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedia.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedia.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedia.F90	(revision 11501)
@@ -175,5 +175,5 @@
       INTEGER            ::   ios, ierror   ! local integer
       !!
-      NAMELIST/namdia/ ln_icediachk, ln_icediahsb, ln_icectl, iiceprt, jiceprt  
+      NAMELIST/namdia/ ln_icediachk, rn_icechk, ln_icediahsb, ln_icectl, iiceprt, jiceprt  
       !!----------------------------------------------------------------------
       !
@@ -191,8 +191,9 @@
          WRITE(numout,*) ' ~~~~~~~~~~~'
          WRITE(numout,*) '   Namelist namdia:'
-         WRITE(numout,*) '      Diagnose online heat/mass/salt budget      ln_icediachk = ', ln_icediachk
-         WRITE(numout,*) '      Output          heat/mass/salt budget      ln_icediahsb = ', ln_icediahsb
-         WRITE(numout,*) '      control prints for a given grid point      ln_icectl    = ', ln_icectl
-         WRITE(numout,*) '         chosen grid point position         (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')'
+         WRITE(numout,*) '      Diagnose online heat/mass/salt conservation ln_icediachk = ', ln_icediachk
+         WRITE(numout,*) '         threshold for conservation               rn_icechk    = ', rn_icechk
+         WRITE(numout,*) '      Output          heat/mass/salt budget       ln_icediahsb = ', ln_icediahsb
+         WRITE(numout,*) '      control prints for a given grid point       ln_icectl    = ', ln_icectl
+         WRITE(numout,*) '         chosen grid point position          (iiceprt,jiceprt) = (', iiceprt,',', jiceprt,')'
       ENDIF
       !      
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rdgrft.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rdgrft.F90	(revision 11501)
@@ -145,4 +145,5 @@
       IF( ln_timing    )   CALL timing_start('icedyn_rdgrft')                                                             ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
 
       IF( kt == nit000 ) THEN
@@ -276,6 +277,7 @@
 
       ! controls
+      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       IF( ln_timing    )   CALL timing_stop ('icedyn_rdgrft')                                                             ! timing
       !
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icedyn_rhg.F90	(revision 11501)
@@ -64,4 +64,5 @@
       IF( ln_timing    )   CALL timing_start('icedyn_rhg')                                                             ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       !
       IF( kt == nit000 .AND. lwp ) THEN
@@ -87,6 +88,7 @@
       !
       ! controls
+      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rhg', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rhg')                                                             ! prints
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rhg',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       IF( ln_timing    )   CALL timing_stop ('icedyn_rhg')                                                             ! timing
       !
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/iceitd.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/iceitd.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/iceitd.F90	(revision 11501)
@@ -88,4 +88,5 @@
 
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
 
       !-----------------------------------------------------------------------------------------------
@@ -316,4 +317,5 @@
       !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_rem',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !
    END SUBROUTINE ice_itd_rem
@@ -586,4 +588,5 @@
       !
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !
       jdonor(:,:) = 0
@@ -664,4 +667,5 @@
       !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'iceitd_reb',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !
    END SUBROUTINE ice_itd_reb
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd.F90	(revision 11501)
@@ -95,4 +95,5 @@
       IF( ln_timing    )   CALL timing_start('icethd')                                                             ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
 
       IF( kt == nit000 .AND. lwp ) THEN
@@ -243,4 +244,5 @@
       !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !                   
       IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd_do.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd_do.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/ICE/icethd_do.F90	(revision 11501)
@@ -113,4 +113,5 @@
 
       IF( ln_icediachk )   CALL ice_cons_hsm( 0, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft )
+      IF( ln_icediachk )   CALL ice_cons2D  ( 0, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft )
 
       at_i(:,:) = SUM( a_i, dim=3 )
@@ -420,4 +421,5 @@
       !
       IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icethd_do', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft)
+      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd_do',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !
    END SUBROUTINE ice_thd_do
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyice.F90	(revision 11500)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyice.F90	(revision 11501)
@@ -63,4 +63,5 @@
       IF( ln_timing    )   CALL timing_start('bdy_ice_thd')                                                            ! timing
       IF( ln_icediachk )   CALL ice_cons_hsm(0,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
+      IF( ln_icediachk )   CALL ice_cons2D  (0,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       !
       CALL ice_var_glo2eqv
@@ -109,6 +110,7 @@
       !
       ! controls
+      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
       IF( ln_icediachk )   CALL ice_cons_hsm(1,'bdy_ice_thd', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation
-      IF( ln_icectl    )   CALL ice_prt     ( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )                        ! prints
+      IF( ln_icediachk )   CALL ice_cons2D  (1,'bdy_ice_thd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation
       IF( ln_timing    )   CALL timing_stop ('bdy_ice_thd')                                                            ! timing
       !
